{
let (x, setX) = React.useState(() => 0.);
let (state, dispatch) =
React.useReducer(DistPlusPlotReducer.reducer, DistPlusPlotReducer.init);
+
{state.distributions
|> E.L.fmapi((index, config) =>
diff --git a/src/components/charts/DistributionPlot/distPlotD3.js b/src/components/charts/DistributionPlot/distPlotD3.js
index 03cc671c..6c58419b 100644
--- a/src/components/charts/DistributionPlot/distPlotD3.js
+++ b/src/components/charts/DistributionPlot/distPlotD3.js
@@ -427,7 +427,7 @@ export class DistPlotD3 {
addLollipopsChart(common) {
const data = this.getDataPoints('discrete');
- const yMin = 0.; //d3.min(this.attrs.data.discrete.ys);
+ const yMin = 0.;
const yMax = d3.max(this.attrs.data.discrete.ys);
// X axis.
diff --git a/src/distPlus/distribution/AlgebraicShapeCombination.re b/src/distPlus/distribution/AlgebraicShapeCombination.re
index a4239200..fd6f7453 100644
--- a/src/distPlus/distribution/AlgebraicShapeCombination.re
+++ b/src/distPlus/distribution/AlgebraicShapeCombination.re
@@ -17,25 +17,25 @@ let toDiscretePointMassesFromTriangulars =
let n = s |> XYShape.T.length;
// first, double up the leftmost and rightmost points:
let {xs, ys}: XYShape.T.t = s;
- let _ = Js.Array.unshift(xs[0], xs);
- let _ = Js.Array.unshift(ys[0], ys);
- let _ = Js.Array.push(xs[n - 1], xs);
- let _ = Js.Array.push(ys[n - 1], ys);
+ Js.Array.unshift(xs[0], xs) |> ignore;
+ Js.Array.unshift(ys[0], ys) |> ignore;
+ Js.Array.push(xs[n - 1], xs) |> ignore;
+ Js.Array.push(ys[n - 1], ys) |> ignore;
let n = E.A.length(xs);
// squares and neighbourly products of the xs
let xsSq: array(float) = Belt.Array.makeUninitializedUnsafe(n);
let xsProdN1: array(float) = Belt.Array.makeUninitializedUnsafe(n - 1);
let xsProdN2: array(float) = Belt.Array.makeUninitializedUnsafe(n - 2);
for (i in 0 to n - 1) {
- let _ = Belt.Array.set(xsSq, i, xs[i] *. xs[i]);
+ Belt.Array.set(xsSq, i, xs[i] *. xs[i]) |> ignore;
();
};
for (i in 0 to n - 2) {
- let _ = Belt.Array.set(xsProdN1, i, xs[i] *. xs[i + 1]);
+ Belt.Array.set(xsProdN1, i, xs[i] *. xs[i + 1]) |> ignore;
();
};
for (i in 0 to n - 3) {
- let _ = Belt.Array.set(xsProdN2, i, xs[i] *. xs[i + 2]);
+ Belt.Array.set(xsProdN2, i, xs[i] *. xs[i + 2]) |> ignore;
();
};
// means and variances
@@ -45,12 +45,11 @@ let toDiscretePointMassesFromTriangulars =
if (inverse) {
for (i in 1 to n - 2) {
- let _ =
- Belt.Array.set(
- masses,
- i - 1,
- (xs[i + 1] -. xs[i - 1]) *. ys[i] /. 2.,
- );
+ Belt.Array.set(
+ masses,
+ i - 1,
+ (xs[i + 1] -. xs[i - 1]) *. ys[i] /. 2.,
+ ) |> ignore;
// this only works when the whole triange is either on the left or on the right of zero
let a = xs[i - 1];
@@ -71,43 +70,39 @@ let toDiscretePointMassesFromTriangulars =
-. inverseMean
** 2.;
- let _ = Belt.Array.set(means, i - 1, inverseMean);
+ Belt.Array.set(means, i - 1, inverseMean) |> ignore;
- let _ = Belt.Array.set(variances, i - 1, inverseVar);
+ Belt.Array.set(variances, i - 1, inverseVar) |> ignore;
();
};
{n: n - 2, masses, means, variances};
} else {
for (i in 1 to n - 2) {
-
// area of triangle = width * height / 2
- let _ =
- Belt.Array.set(
- masses,
- i - 1,
- (xs[i + 1] -. xs[i - 1]) *. ys[i] /. 2.,
- );
+ Belt.Array.set(
+ masses,
+ i - 1,
+ (xs[i + 1] -. xs[i - 1]) *. ys[i] /. 2.,
+ ) |> ignore;
// means of triangle = (a + b + c) / 3
- let _ =
- Belt.Array.set(means, i - 1, (xs[i - 1] +. xs[i] +. xs[i + 1]) /. 3.);
+ Belt.Array.set(means, i - 1, (xs[i - 1] +. xs[i] +. xs[i + 1]) /. 3.) |> ignore;
// variance of triangle = (a^2 + b^2 + c^2 - ab - ac - bc) / 18
- let _ =
- Belt.Array.set(
- variances,
- i - 1,
- (
- xsSq[i - 1]
- +. xsSq[i]
- +. xsSq[i + 1]
- -. xsProdN1[i - 1]
- -. xsProdN1[i]
- -. xsProdN2[i - 1]
- )
- /. 18.,
- );
+ Belt.Array.set(
+ variances,
+ i - 1,
+ (
+ xsSq[i - 1]
+ +. xsSq[i]
+ +. xsSq[i + 1]
+ -. xsProdN1[i - 1]
+ -. xsProdN1[i]
+ -. xsProdN2[i - 1]
+ )
+ /. 18.,
+ ) |> ignore;
();
};
{n: n - 2, masses, means, variances};
@@ -115,7 +110,11 @@ let toDiscretePointMassesFromTriangulars =
};
let combineShapesContinuousContinuous =
- (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape)
+ (
+ op: ExpressionTypes.algebraicOperation,
+ s1: DistTypes.xyShape,
+ s2: DistTypes.xyShape,
+ )
: DistTypes.xyShape => {
let t1n = s1 |> XYShape.T.length;
let t2n = s2 |> XYShape.T.length;
@@ -123,10 +122,11 @@ let combineShapesContinuousContinuous =
// if we add the two distributions, we should probably use normal filters.
// if we multiply the two distributions, we should probably use lognormal filters.
let t1m = toDiscretePointMassesFromTriangulars(s1);
- let t2m = switch (op) {
+ let t2m =
+ switch (op) {
| `Divide => toDiscretePointMassesFromTriangulars(~inverse=true, s2)
| _ => toDiscretePointMassesFromTriangulars(~inverse=false, s2)
- };
+ };
let combineMeansFn =
switch (op) {
@@ -163,7 +163,7 @@ let combineShapesContinuousContinuous =
for (i in 0 to t1m.n - 1) {
for (j in 0 to t2m.n - 1) {
let k = i * t2m.n + j;
- let _ = Belt.Array.set(masses, k, t1m.masses[i] *. t2m.masses[j]);
+ Belt.Array.set(masses, k, t1m.masses[i] *. t2m.masses[j]) |> ignore;
let mean = combineMeansFn(t1m.means[i], t2m.means[j]);
let variance =
@@ -173,8 +173,8 @@ let combineShapesContinuousContinuous =
t1m.means[i],
t2m.means[j],
);
- let _ = Belt.Array.set(means, k, mean);
- let _ = Belt.Array.set(variances, k, variance);
+ Belt.Array.set(means, k, mean) |> ignore;
+ Belt.Array.set(variances, k, variance) |> ignore;
// update bounds
let minX = mean -. 2. *. sqrt(variance) *. 1.644854;
let maxX = mean +. 2. *. sqrt(variance) *. 1.644854;
@@ -190,66 +190,46 @@ let combineShapesContinuousContinuous =
// we now want to create a set of target points. For now, let's just evenly distribute 200 points between
// between the outputMinX and outputMaxX
let nOut = 300;
- let outputXs: array(float) = E.A.Floats.range(outputMinX^, outputMaxX^, nOut);
+ let outputXs: array(float) =
+ E.A.Floats.range(outputMinX^, outputMaxX^, nOut);
let outputYs: array(float) = Belt.Array.make(nOut, 0.0);
// now, for each of the outputYs, accumulate from a Gaussian kernel over each input point.
- for (j in 0 to E.A.length(masses) - 1) { // go through all of the result points
- let _ = if (variances[j] > 0. && masses[j] > 0.) {
- for (i in 0 to E.A.length(outputXs) - 1) { // go through all of the target points
+ for (j in 0 to E.A.length(masses) - 1) {
+ // go through all of the result points
+ if (variances[j] > 0. && masses[j] > 0.) {
+ for (i in 0 to E.A.length(outputXs) - 1) {
+ // go through all of the target points
let dx = outputXs[i] -. means[j];
- let contribution = masses[j] *. exp(-. (dx ** 2.) /. (2. *. variances[j])) /. (sqrt(2. *. 3.14159276 *. variances[j]));
- let _ = Belt.Array.set(outputYs, i, outputYs[i] +. contribution);
- ();
+ let contribution =
+ masses[j]
+ *. exp(-. (dx ** 2.) /. (2. *. variances[j]))
+ /. sqrt(2. *. 3.14159276 *. variances[j]);
+ Belt.Array.set(outputYs, i, outputYs[i] +. contribution) |> ignore;
};
- ();
};
- ();
};
{xs: outputXs, ys: outputYs};
};
-let toDiscretePointMassesFromDiscrete = (s: DistTypes.xyShape): pointMassesWithMoments => {
- let n = s |> XYShape.T.length;
+let toDiscretePointMassesFromDiscrete =
+ (s: DistTypes.xyShape): pointMassesWithMoments => {
let {xs, ys}: XYShape.T.t = s;
let n = E.A.length(xs);
- let masses: array(float) = Belt.Array.makeUninitializedUnsafe(n); // doesn't include the fake first and last points
- let means: array(float) = Belt.Array.makeUninitializedUnsafe(n);
- let variances: array(float) = Belt.Array.makeUninitializedUnsafe(n);
-
- for (i in 0 to n - 1) {
- let _ =
- Belt.Array.set(
- masses,
- i,
- ys[i]
- );
-
- let _ =
- Belt.Array.set(
- means,
- i,
- xs[i]
- );
-
- let _ =
- Belt.Array.set(
- variances,
- i,
- 0.0
- );
- ();
- };
+ let masses: array(float) = Belt.Array.makeBy(n, i => ys[i]);
+ let means: array(float) = Belt.Array.makeBy(n, i => xs[i]);
+ let variances: array(float) = Belt.Array.makeBy(n, i => 0.0);
{n, masses, means, variances};
};
-let combineShapesContinuousDiscreteAdd =
- (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape)
- : DistTypes.xyShape => {
- let t1n = s1 |> XYShape.T.length;
- let t2n = s2 |> XYShape.T.length;
+let combineShapesContinuousDiscrete =
+ (op: ExpressionTypes.algebraicOperation, continuousShape: DistTypes.xyShape, discreteShape: DistTypes.xyShape)
+ : DistTypes.xyShape => {
+
+ let t1n = continuousShape |> XYShape.T.length;
+ let t2n = discreteShape |> XYShape.T.length;
// each x pair is added/subtracted
let fn = Operation.Algebraic.toFn(op);
@@ -257,126 +237,48 @@ let combineShapesContinuousDiscreteAdd =
let outXYShapes: array(array((float, float))) =
Belt.Array.makeUninitializedUnsafe(t2n);
- for (j in 0 to t2n - 1) {
- // for each one of the discrete points
- // create a new distribution, as long as the original continuous one
-
- let dxyShape: array((float, float)) =
- Belt.Array.makeUninitializedUnsafe(t1n);
- for (i in 0 to t1n - 1) {
- let _ =
+ switch (op) {
+ | `Add
+ | `Subtract =>
+ for (j in 0 to t2n - 1) {
+ // creates a new continuous shape for each one of the discrete points, and collects them in outXYShapes.
+ let dxyShape: array((float, float)) =
+ Belt.Array.makeUninitializedUnsafe(t1n);
+ for (i in 0 to t1n - 1) {
Belt.Array.set(
dxyShape,
i,
- (fn(s1.xs[i], s2.xs[j]), s1.ys[i] *. s2.ys[j]),
- );
+ (fn(continuousShape.xs[i], discreteShape.xs[j]),
+ continuousShape.ys[i] *. discreteShape.ys[j]),
+ ) |> ignore;
+ ();
+ };
+ Belt.Array.set(outXYShapes, j, dxyShape) |> ignore;
();
- };
-
- let _ = Belt.Array.set(outXYShapes, j, dxyShape);
- ();
+ }
+ | `Multiply
+ | `Divide =>
+ for (j in 0 to t2n - 1) {
+ // creates a new continuous shape for each one of the discrete points, and collects them in outXYShapes.
+ let dxyShape: array((float, float)) =
+ Belt.Array.makeUninitializedUnsafe(t1n);
+ for (i in 0 to t1n - 1) {
+ Belt.Array.set(
+ dxyShape,
+ i,
+ (fn(continuousShape.xs[i], discreteShape.xs[j]), continuousShape.ys[i] *. discreteShape.ys[j] /. discreteShape.xs[j]),
+ ) |> ignore;
+ ();
+ };
+ Belt.Array.set(outXYShapes, j, dxyShape) |> ignore;
+ ();
+ }
};
outXYShapes
- |> E.A.fold_left(XYShape.PointwiseCombination.combineLinear((+.)), XYShape.T.empty);
-};
-
-let combineShapesContinuousDiscreteMultiply =
- (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape)
- : DistTypes.xyShape => {
- let t1n = s1 |> XYShape.T.length;
- let t2n = s2 |> XYShape.T.length;
-
- let t1m = toDiscretePointMassesFromTriangulars(s1);
- let t2m = toDiscretePointMassesFromDiscrete(s2);
-
- let combineMeansFn =
- switch (op) {
- | `Add => ((m1, m2) => m1 +. m2)
- | `Subtract => ((m1, m2) => m1 -. m2)
- | `Multiply => ((m1, m2) => m1 *. m2)
- | `Divide => ((m1, m2) => m1 /. m2)
- };
-
- let combineVariancesFn =
- switch (op) {
- | `Add
- | `Subtract => ((v1, v2, _, _) => v1 +. v2)
- | `Multiply
- | `Divide => (
- (v1, v2, m1, m2) => v1 *. m2 ** 2.
- )
- };
-
- let outputMinX: ref(float) = ref(infinity);
- let outputMaxX: ref(float) = ref(neg_infinity);
- let masses: array(float) =
- Belt.Array.makeUninitializedUnsafe(t1m.n * t2m.n);
- let means: array(float) =
- Belt.Array.makeUninitializedUnsafe(t1m.n * t2m.n);
- let variances: array(float) =
- Belt.Array.makeUninitializedUnsafe(t1m.n * t2m.n);
- // then convolve the two sets of pointMassesWithMoments
- for (i in 0 to t1m.n - 1) {
- for (j in 0 to t2m.n - 1) {
- let k = i * t2m.n + j;
- let _ = Belt.Array.set(masses, k, t1m.masses[i] *. t2m.masses[j]);
-
- let mean = combineMeansFn(t1m.means[i], t2m.means[j]);
- let variance =
- combineVariancesFn(
- t1m.variances[i],
- t2m.variances[j],
- t1m.means[i],
- t2m.means[j],
- );
- let _ = Belt.Array.set(means, k, mean);
- let _ = Belt.Array.set(variances, k, variance);
-
- // update bounds
- let minX = mean -. 2. *. sqrt(variance) *. 1.644854;
- let maxX = mean +. 2. *. sqrt(variance) *. 1.644854;
- if (minX < outputMinX^) {
- outputMinX := minX;
- };
- if (maxX > outputMaxX^) {
- outputMaxX := maxX;
- };
- };
- };
-
-
- // we now want to create a set of target points. For now, let's just evenly distribute 200 points between
- // between the outputMinX and outputMaxX
- let nOut = 300;
- let outputXs: array(float) = E.A.Floats.range(outputMinX^, outputMaxX^, nOut);
- let outputYs: array(float) = Belt.Array.make(nOut, 0.0);
- // now, for each of the outputYs, accumulate from a Gaussian kernel over each input point.
- for (j in 0 to E.A.length(masses) - 1) { // go through all of the result points
- let _ = if (variances[j] > 0. && masses[j] > 0.) {
- for (i in 0 to E.A.length(outputXs) - 1) { // go through all of the target points
- let dx = outputXs[i] -. means[j];
- let contribution = masses[j] *. exp(-. (dx ** 2.) /. (2. *. variances[j])) /. (sqrt(2. *. 3.14159276 *. variances[j]));
- let _ = Belt.Array.set(outputYs, i, outputYs[i] +. contribution);
- ();
- };
- ();
- };
- ();
- };
-
- {xs: outputXs, ys: outputYs};
-};
-
-let combineShapesContinuousDiscrete =
- (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape)
- : DistTypes.xyShape => {
-
- switch (op) {
- | `Add
- | `Subtract => combineShapesContinuousDiscreteAdd(op, s1, s2);
- | `Multiply
- | `Divide => combineShapesContinuousDiscreteMultiply(op, s1, s2);
- };
-
+ |> E.A.fmap(XYShape.T.fromZippedArray)
+ |> E.A.fold_left(
+ XYShape.PointwiseCombination.combine((+.),
+ XYShape.XtoY.continuousInterpolator(`Linear, `UseZero)),
+ XYShape.T.empty);
};
diff --git a/src/distPlus/distribution/Continuous.re b/src/distPlus/distribution/Continuous.re
index 173c65bc..1adb60c0 100644
--- a/src/distPlus/distribution/Continuous.re
+++ b/src/distPlus/distribution/Continuous.re
@@ -3,30 +3,45 @@ open Distributions;
type t = DistTypes.continuousShape;
let getShape = (t: t) => t.xyShape;
let interpolation = (t: t) => t.interpolation;
-let make = (interpolation, xyShape, knownIntegralSum): t => {
+let make = (~interpolation=`Linear, ~integralSumCache=None, ~integralCache=None, xyShape): t => {
xyShape,
interpolation,
- knownIntegralSum,
+ integralSumCache,
+ integralCache,
};
-let shapeMap = (fn, {xyShape, interpolation, knownIntegralSum}: t): t => {
+let shapeMap = (fn, {xyShape, interpolation, integralSumCache, integralCache}: t): t => {
xyShape: fn(xyShape),
interpolation,
- knownIntegralSum,
+ integralSumCache,
+ integralCache,
};
let lastY = (t: t) => t |> getShape |> XYShape.T.lastY;
let oShapeMap =
- (fn, {xyShape, interpolation, knownIntegralSum}: t)
+ (fn, {xyShape, interpolation, integralSumCache, integralCache}: t)
: option(DistTypes.continuousShape) =>
- fn(xyShape) |> E.O.fmap(make(interpolation, _, knownIntegralSum));
+ fn(xyShape) |> E.O.fmap(make(~interpolation, ~integralSumCache, ~integralCache));
+let emptyIntegral: DistTypes.continuousShape = {
+ xyShape: {xs: [|neg_infinity|], ys: [|0.0|]},
+ interpolation: `Linear,
+ integralSumCache: Some(0.0),
+ integralCache: None,
+};
let empty: DistTypes.continuousShape = {
xyShape: XYShape.T.empty,
interpolation: `Linear,
- knownIntegralSum: Some(0.0),
+ integralSumCache: Some(0.0),
+ integralCache: Some(emptyIntegral),
};
+
+let stepwiseToLinear = (t: t): t =>
+ make(~integralSumCache=t.integralSumCache, ~integralCache=t.integralCache, XYShape.Range.stepwiseToLinear(t.xyShape));
+
let combinePointwise =
(
- ~knownIntegralSumsFn,
+ ~integralSumCachesFn=(_, _) => None,
+ ~integralCachesFn: (t, t) => option(t) =(_, _) => None,
+ ~distributionType: DistTypes.distributionType = `PDF,
fn: (float, float) => float,
t1: DistTypes.continuousShape,
t2: DistTypes.continuousShape,
@@ -36,59 +51,91 @@ let combinePointwise =
// can just sum them up. Otherwise, all bets are off.
let combinedIntegralSum =
Common.combineIntegralSums(
- knownIntegralSumsFn,
- t1.knownIntegralSum,
- t2.knownIntegralSum,
+ integralSumCachesFn,
+ t1.integralSumCache,
+ t2.integralSumCache,
);
+ // TODO: does it ever make sense to pointwise combine the integrals here?
+ // It could be done for pointwise additions, but is that ever needed?
+
+ // If combining stepwise and linear, we must convert the stepwise to linear first,
+ // i.e. add a point at the bottom of each step
+ let (t1, t2) = switch (t1.interpolation, t2.interpolation) {
+ | (`Linear, `Linear) => (t1, t2);
+ | (`Stepwise, `Stepwise) => (t1, t2);
+ | (`Linear, `Stepwise) => (t1, stepwiseToLinear(t2));
+ | (`Stepwise, `Linear) => (stepwiseToLinear(t1), t2);
+ };
+
+ let extrapolation = switch (distributionType) {
+ | `PDF => `UseZero
+ | `CDF => `UseOutermostPoints
+ };
+
+ let interpolator = XYShape.XtoY.continuousInterpolator(t1.interpolation, extrapolation);
+
make(
- `Linear,
- XYShape.PointwiseCombination.combineLinear(
- ~fn=(+.),
+ ~integralSumCache=combinedIntegralSum,
+ XYShape.PointwiseCombination.combine(
+ (+.),
+ interpolator,
t1.xyShape,
t2.xyShape,
),
- combinedIntegralSum,
);
};
let toLinear = (t: t): option(t) => {
switch (t) {
- | {interpolation: `Stepwise, xyShape, knownIntegralSum} =>
+ | {interpolation: `Stepwise, xyShape, integralSumCache, integralCache} =>
xyShape
|> XYShape.Range.stepsToContinuous
- |> E.O.fmap(make(`Linear, _, knownIntegralSum))
+ |> E.O.fmap(make(~integralSumCache, ~integralCache))
| {interpolation: `Linear} => Some(t)
};
};
let shapeFn = (fn, t: t) => t |> getShape |> fn;
-let updateKnownIntegralSum = (knownIntegralSum, t: t): t => {
+
+let updateIntegralSumCache = (integralSumCache, t: t): t => {
...t,
- knownIntegralSum,
+ integralSumCache,
+};
+
+let updateIntegralCache = (integralCache, t: t): t => {
+ ...t,
+ integralCache,
};
let reduce =
(
- ~knownIntegralSumsFn: (float, float) => option(float)=(_, _) => None,
+ ~integralSumCachesFn: (float, float) => option(float)=(_, _) => None,
+ ~integralCachesFn: (t, t) => option(t)=(_, _) => None,
fn,
continuousShapes,
) =>
continuousShapes
- |> E.A.fold_left(combinePointwise(~knownIntegralSumsFn, fn), empty);
+ |> E.A.fold_left(combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn), empty);
-let mapY = (~knownIntegralSumFn=_ => None, fn, t: t) => {
- let u = E.O.bind(_, knownIntegralSumFn);
- let yMapFn = shapeMap(XYShape.T.mapY(fn));
-
- t |> yMapFn |> updateKnownIntegralSum(u(t.knownIntegralSum));
+let mapY = (~integralSumCacheFn=_ => None,
+ ~integralCacheFn=_ => None,
+ ~fn, t: t) => {
+ make(
+ ~interpolation=t.interpolation,
+ ~integralSumCache=t.integralSumCache |> E.O.bind(_, integralSumCacheFn),
+ ~integralCache=t.integralCache |> E.O.bind(_, integralCacheFn),
+ t |> getShape |> XYShape.T.mapY(fn),
+ );
};
-let scaleBy = (~scale=1.0, t: t): t => {
+let rec scaleBy = (~scale=1.0, t: t): t => {
+ let scaledIntegralSumCache = E.O.bind(t.integralSumCache, v => Some(scale *. v));
+ let scaledIntegralCache = E.O.bind(t.integralCache, v => Some(scaleBy(~scale, v)));
+
t
- |> mapY((r: float) => r *. scale)
- |> updateKnownIntegralSum(
- E.O.bind(t.knownIntegralSum, v => Some(scale *. v)),
- );
+ |> mapY(~fn=(r: float) => r *. scale)
+ |> updateIntegralSumCache(scaledIntegralSumCache)
+ |> updateIntegralCache(scaledIntegralCache)
};
module T =
@@ -98,6 +145,7 @@ module T =
let minX = shapeFn(XYShape.T.minX);
let maxX = shapeFn(XYShape.T.maxX);
let mapY = mapY;
+ let updateIntegralCache = updateIntegralCache;
let toDiscreteProbabilityMassFraction = _ => 0.0;
let toShape = (t: t): DistTypes.shape => Continuous(t);
let xToY = (f, {interpolation, xyShape}: t) => {
@@ -121,63 +169,56 @@ module T =
|> XYShape.T.zip
|> XYShape.Zipped.filterByX(x => x >= lc && x <= rc);
- let eps = (t |> getShape |> XYShape.T.xTotalRange) *. 0.0001;
-
let leftNewPoint =
- leftCutoff |> E.O.dimap(lc => [|(lc -. eps, 0.)|], _ => [||]);
+ leftCutoff |> E.O.dimap(lc => [|(lc -. epsilon_float, 0.)|], _ => [||]);
let rightNewPoint =
- rightCutoff |> E.O.dimap(rc => [|(rc +. eps, 0.)|], _ => [||]);
+ rightCutoff |> E.O.dimap(rc => [|(rc +. epsilon_float, 0.)|], _ => [||]);
let truncatedZippedPairsWithNewPoints =
E.A.concatMany([|leftNewPoint, truncatedZippedPairs, rightNewPoint|]);
let truncatedShape =
XYShape.T.fromZippedArray(truncatedZippedPairsWithNewPoints);
- make(`Linear, truncatedShape, None);
+ make(truncatedShape)
};
// TODO: This should work with stepwise plots.
- let integral = (~cache, t) =>
- if (t |> getShape |> XYShape.T.length > 0) {
- switch (cache) {
- | Some(cache) => cache
- | None =>
- t
- |> getShape
- |> XYShape.Range.integrateWithTriangles
- |> E.O.toExt("This should not have happened")
- |> make(`Linear, _, None)
- };
- } else {
- make(`Linear, {xs: [|neg_infinity|], ys: [|0.0|]}, None);
+ let integral = (t) =>
+ switch (getShape(t) |> XYShape.T.isEmpty, t.integralCache) {
+ | (true, _) => emptyIntegral
+ | (false, Some(cache)) => cache
+ | (false, None) =>
+ t
+ |> getShape
+ |> XYShape.Range.integrateWithTriangles
+ |> E.O.toExt("This should not have happened")
+ |> make
};
- let downsample = (~cache=None, length, t): t =>
+ let downsample = (length, t): t =>
t
|> shapeMap(
XYShape.XsConversion.proportionByProbabilityMass(
length,
- integral(~cache, t).xyShape,
+ integral(t).xyShape,
),
);
- let integralEndY = (~cache, t: t) =>
- t.knownIntegralSum |> E.O.default(t |> integral(~cache) |> lastY);
- let integralXtoY = (~cache, f, t: t) =>
- t |> integral(~cache) |> shapeFn(XYShape.XtoY.linear(f));
- let integralYtoX = (~cache, f, t: t) =>
- t |> integral(~cache) |> shapeFn(XYShape.YtoX.linear(f));
+ let integralEndY = (t: t) =>
+ t.integralSumCache |> E.O.default(t |> integral |> lastY);
+ let integralXtoY = (f, t: t) =>
+ t |> integral |> shapeFn(XYShape.XtoY.linear(f));
+ let integralYtoX = (f, t: t) =>
+ t |> integral |> shapeFn(XYShape.YtoX.linear(f));
let toContinuous = t => Some(t);
let toDiscrete = _ => None;
let normalize = (t: t): t => {
t
- |> scaleBy(~scale=1. /. integralEndY(~cache=None, t))
- |> updateKnownIntegralSum(Some(1.0));
+ |> updateIntegralCache(Some(integral(t)))
+ |> scaleBy(~scale=1. /. integralEndY(t))
+ |> updateIntegralSumCache(Some(1.0));
};
- let normalizedToContinuous = t => Some(t |> normalize);
- let normalizedToDiscrete = _ => None;
-
let mean = (t: t) => {
let indefiniteIntegralStepwise = (p, h1) => h1 *. p ** 2.0 /. 2.0;
let indefiniteIntegralLinear = (p, a, b) =>
@@ -205,23 +246,28 @@ let combineAlgebraicallyWithDiscrete =
t1: t,
t2: DistTypes.discreteShape,
) => {
- let s1 = t1 |> getShape;
- let s2 = t2.xyShape;
- let t1n = s1 |> XYShape.T.length;
- let t2n = s2 |> XYShape.T.length;
- if (t1n == 0 || t2n == 0) {
+ let t1s = t1 |> getShape;
+ let t2s = t2.xyShape; // TODO would like to use Discrete.getShape here, but current file structure doesn't allow for that
+
+ if (XYShape.T.isEmpty(t1s) || XYShape.T.isEmpty(t2s)) {
empty;
} else {
- let combinedShape =
- AlgebraicShapeCombination.combineShapesContinuousDiscrete(op, s1, s2);
+ let continuousAsLinear = switch (t1.interpolation) {
+ | `Linear => t1;
+ | `Stepwise => stepwiseToLinear(t1)
+ };
+
+ let combinedShape = AlgebraicShapeCombination.combineShapesContinuousDiscrete(op, continuousAsLinear |> getShape, t2s);
+
let combinedIntegralSum =
Common.combineIntegralSums(
(a, b) => Some(a *. b),
- t1.knownIntegralSum,
- t2.knownIntegralSum,
+ t1.integralSumCache,
+ t2.integralSumCache,
);
- // return a new Continuous distribution
- make(`Linear, combinedShape, combinedIntegralSum);
+
+ // TODO: It could make sense to automatically transform the integrals here (shift or scale)
+ make(~interpolation=t1.interpolation, ~integralSumCache=combinedIntegralSum, combinedShape)
};
};
@@ -239,10 +285,10 @@ let combineAlgebraically =
let combinedIntegralSum =
Common.combineIntegralSums(
(a, b) => Some(a *. b),
- t1.knownIntegralSum,
- t2.knownIntegralSum,
+ t1.integralSumCache,
+ t2.integralSumCache,
);
// return a new Continuous distribution
- make(`Linear, combinedShape, combinedIntegralSum);
+ make(~integralSumCache=combinedIntegralSum, combinedShape);
};
-};
\ No newline at end of file
+};
diff --git a/src/distPlus/distribution/Discrete.re b/src/distPlus/distribution/Discrete.re
index 1bea9c45..8a64a454 100644
--- a/src/distPlus/distribution/Discrete.re
+++ b/src/distPlus/distribution/Discrete.re
@@ -2,23 +2,37 @@ open Distributions;
type t = DistTypes.discreteShape;
-let make = (xyShape, knownIntegralSum): t => {xyShape, knownIntegralSum};
-let shapeMap = (fn, {xyShape, knownIntegralSum}: t): t => {
+let make = (~integralSumCache=None, ~integralCache=None, xyShape): t => {xyShape, integralSumCache, integralCache};
+let shapeMap = (fn, {xyShape, integralSumCache, integralCache}: t): t => {
xyShape: fn(xyShape),
- knownIntegralSum,
+ integralSumCache,
+ integralCache
};
let getShape = (t: t) => t.xyShape;
-let oShapeMap = (fn, {xyShape, knownIntegralSum}: t): option(t) =>
- fn(xyShape) |> E.O.fmap(make(_, knownIntegralSum));
+let oShapeMap = (fn, {xyShape, integralSumCache, integralCache}: t): option(t) =>
+ fn(xyShape) |> E.O.fmap(make(~integralSumCache, ~integralCache));
+
+let emptyIntegral: DistTypes.continuousShape = {
+ xyShape: {xs: [|neg_infinity|], ys: [|0.0|]},
+ interpolation: `Stepwise,
+ integralSumCache: Some(0.0),
+ integralCache: None,
+};
+let empty: DistTypes.discreteShape = {
+ xyShape: XYShape.T.empty,
+ integralSumCache: Some(0.0),
+ integralCache: Some(emptyIntegral),
+};
+
-let empty: t = {xyShape: XYShape.T.empty, knownIntegralSum: Some(0.0)};
let shapeFn = (fn, t: t) => t |> getShape |> fn;
let lastY = (t: t) => t |> getShape |> XYShape.T.lastY;
let combinePointwise =
(
- ~knownIntegralSumsFn,
+ ~integralSumCachesFn = (_, _) => None,
+ ~integralCachesFn: (DistTypes.continuousShape, DistTypes.continuousShape) => option(DistTypes.continuousShape) = (_, _) => None,
fn,
t1: DistTypes.discreteShape,
t2: DistTypes.discreteShape,
@@ -26,38 +40,47 @@ let combinePointwise =
: DistTypes.discreteShape => {
let combinedIntegralSum =
Common.combineIntegralSums(
- knownIntegralSumsFn,
- t1.knownIntegralSum,
- t2.knownIntegralSum,
+ integralSumCachesFn,
+ t1.integralSumCache,
+ t2.integralSumCache,
);
+ // TODO: does it ever make sense to pointwise combine the integrals here?
+ // It could be done for pointwise additions, but is that ever needed?
+
make(
+ ~integralSumCache=combinedIntegralSum,
XYShape.PointwiseCombination.combine(
- ~xsSelection=ALL_XS,
- ~xToYSelection=XYShape.XtoY.stepwiseIfAtX,
- ~fn=(a, b) => fn(E.O.default(0.0, a), E.O.default(0.0, b)), // stepwiseIfAtX returns option(float), so this fn needs to handle None
+ (+.),
+ XYShape.XtoY.discreteInterpolator,
t1.xyShape,
t2.xyShape,
),
- combinedIntegralSum,
);
};
let reduce =
- (~knownIntegralSumsFn=(_, _) => None, fn, discreteShapes)
- : DistTypes.discreteShape =>
+ (~integralSumCachesFn=(_, _) => None,
+ ~integralCachesFn=(_, _) => None,
+ fn, discreteShapes)
+ : DistTypes.discreteShape =>
discreteShapes
- |> E.A.fold_left(combinePointwise(~knownIntegralSumsFn, fn), empty);
+ |> E.A.fold_left(combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn), empty);
-let updateKnownIntegralSum = (knownIntegralSum, t: t): t => {
+let updateIntegralSumCache = (integralSumCache, t: t): t => {
...t,
- knownIntegralSum,
+ integralSumCache,
+};
+
+let updateIntegralCache = (integralCache, t: t): t => {
+ ...t,
+ integralCache,
};
/* This multiples all of the data points together and creates a new discrete distribution from the results.
Data points at the same xs get added together. It may be a good idea to downsample t1 and t2 before and/or the result after. */
let combineAlgebraically =
- (op: ExpressionTypes.algebraicOperation, t1: t, t2: t) => {
+ (op: ExpressionTypes.algebraicOperation, t1: t, t2: t): t => {
let t1s = t1 |> getShape;
let t2s = t2 |> getShape;
let t1n = t1s |> XYShape.T.length;
@@ -66,8 +89,8 @@ let combineAlgebraically =
let combinedIntegralSum =
Common.combineIntegralSums(
(s1, s2) => Some(s1 *. s2),
- t1.knownIntegralSum,
- t2.knownIntegralSum,
+ t1.integralSumCache,
+ t2.integralSumCache,
);
let fn = Operation.Algebraic.toFn(op);
@@ -87,84 +110,85 @@ let combineAlgebraically =
let combinedShape = XYShape.T.fromZippedArray(rxys);
- make(combinedShape, combinedIntegralSum);
+ make(~integralSumCache=combinedIntegralSum, combinedShape);
};
-let mapY = (~knownIntegralSumFn=previousKnownIntegralSum => None, fn, t: t) => {
- let u = E.O.bind(_, knownIntegralSumFn);
- let yMapFn = shapeMap(XYShape.T.mapY(fn));
-
- t |> yMapFn |> updateKnownIntegralSum(u(t.knownIntegralSum));
+let mapY = (~integralSumCacheFn=_ => None,
+ ~integralCacheFn=_ => None,
+ ~fn, t: t) => {
+ make(
+ ~integralSumCache=t.integralSumCache |> E.O.bind(_, integralSumCacheFn),
+ ~integralCache=t.integralCache |> E.O.bind(_, integralCacheFn),
+ t |> getShape |> XYShape.T.mapY(fn),
+ );
};
+
let scaleBy = (~scale=1.0, t: t): t => {
+ let scaledIntegralSumCache = t.integralSumCache |> E.O.fmap((*.)(scale));
+ let scaledIntegralCache = t.integralCache |> E.O.fmap(Continuous.scaleBy(~scale));
+
t
- |> mapY((r: float) => r *. scale)
- |> updateKnownIntegralSum(
- E.O.bind(t.knownIntegralSum, v => Some(scale *. v)),
- );
+ |> mapY(~fn=(r: float) => r *. scale)
+ |> updateIntegralSumCache(scaledIntegralSumCache)
+ |> updateIntegralCache(scaledIntegralCache)
};
module T =
Dist({
type t = DistTypes.discreteShape;
type integral = DistTypes.continuousShape;
- let integral = (~cache, t) =>
- if (t |> getShape |> XYShape.T.length > 0) {
- switch (cache) {
- | Some(c) => c
- | None =>
- Continuous.make(
- `Stepwise,
- XYShape.T.accumulateYs((+.), getShape(t)),
- None,
- )
- };
- } else {
- Continuous.make(
- `Stepwise,
- {xs: [|neg_infinity|], ys: [|0.0|]},
- None,
- );
+ let integral = (t) =>
+ switch (getShape(t) |> XYShape.T.isEmpty, t.integralCache) {
+ | (true, _) => emptyIntegral
+ | (false, Some(c)) => c
+ | (false, None) => {
+ let ts = getShape(t);
+ // The first xy of this integral should always be the zero, to ensure nice plotting
+ let firstX = ts |> XYShape.T.minX;
+ let prependedZeroPoint: XYShape.T.t = {xs: [|firstX -. epsilon_float|], ys: [|0.|]};
+ let integralShape =
+ ts
+ |> XYShape.T.concat(prependedZeroPoint)
+ |> XYShape.T.accumulateYs((+.));
+
+ Continuous.make(~interpolation=`Stepwise, integralShape);
+ }
};
- let integralEndY = (~cache, t: t) =>
- t.knownIntegralSum
- |> E.O.default(t |> integral(~cache) |> Continuous.lastY);
+ let integralEndY = (t: t) =>
+ t.integralSumCache
+ |> E.O.default(t |> integral |> Continuous.lastY);
let minX = shapeFn(XYShape.T.minX);
let maxX = shapeFn(XYShape.T.maxX);
let toDiscreteProbabilityMassFraction = _ => 1.0;
let mapY = mapY;
+ let updateIntegralCache = updateIntegralCache;
let toShape = (t: t): DistTypes.shape => Discrete(t);
let toContinuous = _ => None;
let toDiscrete = t => Some(t);
let normalize = (t: t): t => {
t
- |> scaleBy(~scale=1. /. integralEndY(~cache=None, t))
- |> updateKnownIntegralSum(Some(1.0));
+ |> scaleBy(~scale=1. /. integralEndY(t))
+ |> updateIntegralSumCache(Some(1.0));
};
- let normalizedToContinuous = _ => None;
- let normalizedToDiscrete = t => Some(t); // TODO: this should be normalized!
-
- let downsample = (~cache=None, i, t: t): t => {
+ let downsample = (i, t: t): t => {
// It's not clear how to downsample a set of discrete points in a meaningful way.
// The best we can do is to clip off the smallest values.
let currentLength = t |> getShape |> XYShape.T.length;
if (i < currentLength && i >= 1 && currentLength > 1) {
- let clippedShape =
- t
- |> getShape
- |> XYShape.T.zip
- |> XYShape.Zipped.sortByY
- |> Belt.Array.reverse
- |> Belt.Array.slice(_, ~offset=0, ~len=i)
- |> XYShape.Zipped.sortByX
- |> XYShape.T.fromZippedArray;
-
- make(clippedShape, None); // if someone needs the sum, they'll have to recompute it
+ t
+ |> getShape
+ |> XYShape.T.zip
+ |> XYShape.Zipped.sortByY
+ |> Belt.Array.reverse
+ |> Belt.Array.slice(_, ~offset=0, ~len=i)
+ |> XYShape.Zipped.sortByX
+ |> XYShape.T.fromZippedArray
+ |> make;
} else {
t;
};
@@ -172,17 +196,15 @@ module T =
let truncate =
(leftCutoff: option(float), rightCutoff: option(float), t: t): t => {
- let truncatedShape =
- t
- |> getShape
- |> XYShape.T.zip
- |> XYShape.Zipped.filterByX(x =>
- x >= E.O.default(neg_infinity, leftCutoff)
- || x <= E.O.default(infinity, rightCutoff)
- )
- |> XYShape.T.fromZippedArray;
-
- make(truncatedShape, None);
+ t
+ |> getShape
+ |> XYShape.T.zip
+ |> XYShape.Zipped.filterByX(x =>
+ x >= E.O.default(neg_infinity, leftCutoff)
+ && x <= E.O.default(infinity, rightCutoff)
+ )
+ |> XYShape.T.fromZippedArray
+ |> make;
};
let xToY = (f, t) =>
@@ -192,11 +214,11 @@ module T =
|> E.O.default(0.0)
|> DistTypes.MixedPoint.makeDiscrete;
- let integralXtoY = (~cache, f, t) =>
- t |> integral(~cache) |> Continuous.getShape |> XYShape.XtoY.linear(f);
+ let integralXtoY = (f, t) =>
+ t |> integral |> Continuous.getShape |> XYShape.XtoY.linear(f);
- let integralYtoX = (~cache, f, t) =>
- t |> integral(~cache) |> Continuous.getShape |> XYShape.YtoX.linear(f);
+ let integralYtoX = (f, t) =>
+ t |> integral |> Continuous.getShape |> XYShape.YtoX.linear(f);
let mean = (t: t): float => {
let s = getShape(t);
diff --git a/src/distPlus/distribution/DistPlus.re b/src/distPlus/distribution/DistPlus.re
index ff3d459b..405a0335 100644
--- a/src/distPlus/distribution/DistPlus.re
+++ b/src/distPlus/distribution/DistPlus.re
@@ -2,8 +2,7 @@ open DistTypes;
type t = DistTypes.distPlus;
-let shapeIntegral = shape =>
- Shape.T.Integral.get(~cache=None, shape);
+let shapeIntegral = shape => Shape.T.Integral.get(shape);
let make =
(
~shape,
@@ -59,7 +58,6 @@ module T =
let normalize = (t: t): t => {
let normalizedShape = t |> toShape |> Shape.T.normalize;
t |> updateShape(normalizedShape);
- // TODO: also adjust for domainIncludedProbabilityMass here.
};
let truncate = (leftCutoff, rightCutoff, t: t): t => {
@@ -71,28 +69,6 @@ module T =
t |> updateShape(truncatedShape);
};
- let normalizedToContinuous = (t: t) => {
- t
- |> toShape
- |> Shape.T.normalizedToContinuous
- |> E.O.fmap(
- Continuous.T.mapY(
- domainIncludedProbabilityMassAdjustment(t),
- ),
- );
- };
-
- let normalizedToDiscrete = (t: t) => {
- t
- |> toShape
- |> Shape.T.normalizedToDiscrete
- |> E.O.fmap(
- Discrete.T.mapY(
- domainIncludedProbabilityMassAdjustment(t),
- ),
- );
- };
-
let xToY = (f, t: t) =>
t
|> toShape
@@ -105,34 +81,36 @@ module T =
shapeFn(Shape.T.toDiscreteProbabilityMassFraction);
// This bit is kind of awkward, could probably use rethinking.
- let integral = (~cache, t: t) =>
+ let integral = (t: t) =>
updateShape(Continuous(t.integralCache), t);
- let downsample = (~cache=None, i, t): t =>
+ let updateIntegralCache = (integralCache: option(DistTypes.continuousShape), t) =>
+ update(~integralCache=E.O.default(t.integralCache, integralCache), t);
+
+ let downsample = (i, t): t =>
updateShape(t |> toShape |> Shape.T.downsample(i), t);
// todo: adjust for limit, maybe?
let mapY =
(
- ~knownIntegralSumFn=previousIntegralSum => None,
- fn,
+ ~integralSumCacheFn=previousIntegralSum => None,
+ ~integralCacheFn=previousIntegralCache => None,
+ ~fn,
{shape, _} as t: t,
)
: t =>
- Shape.T.mapY(~knownIntegralSumFn, fn, shape)
+ Shape.T.mapY(~integralSumCacheFn, ~fn, shape)
|> updateShape(_, t);
// get the total of everything
- let integralEndY = (~cache as _, t: t) => {
+ let integralEndY = (t: t) => {
Shape.T.Integral.sum(
- ~cache=Some(t.integralCache),
toShape(t),
);
};
// TODO: Fix this below, obviously. Adjust for limits
- let integralXtoY = (~cache as _, f, t: t) => {
+ let integralXtoY = (f, t: t) => {
Shape.T.Integral.xToY(
- ~cache=Some(t.integralCache),
f,
toShape(t),
)
@@ -140,8 +118,8 @@ module T =
};
// TODO: This part is broken when there is a limit, if this is supposed to be taken into account.
- let integralYtoX = (~cache as _, f, t: t) => {
- Shape.T.Integral.yToX(~cache=None, f, toShape(t));
+ let integralYtoX = (f, t: t) => {
+ Shape.T.Integral.yToX(f, toShape(t));
};
let mean = (t: t) => {
diff --git a/src/distPlus/distribution/DistPlusTime.re b/src/distPlus/distribution/DistPlusTime.re
index 064ec32c..d5b36f8e 100644
--- a/src/distPlus/distribution/DistPlusTime.re
+++ b/src/distPlus/distribution/DistPlusTime.re
@@ -23,6 +23,6 @@
include DistPlus.T.Integral;
let xToY = (f: TimeTypes.timeInVector, t: t) => {
timeInVectorToX(f, t)
- |> E.O.fmap(x => DistPlus.T.Integral.xToY(~cache=None, x, t));
+ |> E.O.fmap(x => DistPlus.T.Integral.xToY(x, t));
};
- };
\ No newline at end of file
+ };
diff --git a/src/distPlus/distribution/DistTypes.re b/src/distPlus/distribution/DistTypes.re
index 8e1b7b10..92f03ee6 100644
--- a/src/distPlus/distribution/DistTypes.re
+++ b/src/distPlus/distribution/DistTypes.re
@@ -9,25 +9,45 @@ type domain =
| RightLimited(domainLimit)
| LeftAndRightLimited(domainLimit, domainLimit);
+type distributionType = [
+ | `PDF
+ | `CDF
+];
+
type xyShape = {
xs: array(float),
ys: array(float),
};
+type interpolationStrategy = [
+ | `Stepwise
+ | `Linear
+];
+type extrapolationStrategy = [
+ | `UseZero
+ | `UseOutermostPoints
+];
+
+type interpolator = (xyShape, int, float) => float;
+
type continuousShape = {
xyShape,
- interpolation: [ | `Stepwise | `Linear],
- knownIntegralSum: option(float),
+ interpolation: interpolationStrategy,
+ integralSumCache: option(float),
+ integralCache: option(continuousShape),
};
type discreteShape = {
xyShape,
- knownIntegralSum: option(float),
+ integralSumCache: option(float),
+ integralCache: option(continuousShape),
};
type mixedShape = {
continuous: continuousShape,
discrete: discreteShape,
+ integralSumCache: option(float),
+ integralCache: option(continuousShape),
};
type shapeMonad('a, 'b, 'c) =
diff --git a/src/distPlus/distribution/Distributions.re b/src/distPlus/distribution/Distributions.re
index 00d366e0..05b57edb 100644
--- a/src/distPlus/distribution/Distributions.re
+++ b/src/distPlus/distribution/Distributions.re
@@ -4,22 +4,22 @@ module type dist = {
let minX: t => float;
let maxX: t => float;
let mapY:
- (~knownIntegralSumFn: float => option(float)=?, float => float, t) => t;
+ (~integralSumCacheFn: float => option(float)=?, ~integralCacheFn: DistTypes.continuousShape => option(DistTypes.continuousShape)=?, ~fn: float => float, t) => t;
let xToY: (float, t) => DistTypes.mixedPoint;
let toShape: t => DistTypes.shape;
let toContinuous: t => option(DistTypes.continuousShape);
let toDiscrete: t => option(DistTypes.discreteShape);
let normalize: t => t;
- let normalizedToContinuous: t => option(DistTypes.continuousShape);
- let normalizedToDiscrete: t => option(DistTypes.discreteShape);
let toDiscreteProbabilityMassFraction: t => float;
- let downsample: (~cache: option(integral)=?, int, t) => t;
+ let downsample: (int, t) => t;
let truncate: (option(float), option(float), t) => t;
- let integral: (~cache: option(integral), t) => integral;
- let integralEndY: (~cache: option(integral), t) => float;
- let integralXtoY: (~cache: option(integral), float, t) => float;
- let integralYtoX: (~cache: option(integral), float, t) => float;
+ let updateIntegralCache: (option(DistTypes.continuousShape), t) => t;
+
+ let integral: (t) => integral;
+ let integralEndY: (t) => float;
+ let integralXtoY: (float, t) => float;
+ let integralYtoX: (float, t) => float;
let mean: t => float;
let variance: t => float;
@@ -41,11 +41,11 @@ module Dist = (T: dist) => {
let toDiscrete = T.toDiscrete;
let normalize = T.normalize;
let truncate = T.truncate;
- let normalizedToContinuous = T.normalizedToContinuous;
- let normalizedToDiscrete = T.normalizedToDiscrete;
let mean = T.mean;
let variance = T.variance;
+ let updateIntegralCache = T.updateIntegralCache;
+
module Integral = {
type t = T.integral;
let get = T.integral;
@@ -59,10 +59,23 @@ module Common = {
let combineIntegralSums =
(
combineFn: (float, float) => option(float),
- t1KnownIntegralSum: option(float),
- t2KnownIntegralSum: option(float),
+ t1IntegralSumCache: option(float),
+ t2IntegralSumCache: option(float),
) => {
- switch (t1KnownIntegralSum, t2KnownIntegralSum) {
+ switch (t1IntegralSumCache, t2IntegralSumCache) {
+ | (None, _)
+ | (_, None) => None
+ | (Some(s1), Some(s2)) => combineFn(s1, s2)
+ };
+ };
+
+ let combineIntegrals =
+ (
+ combineFn: (DistTypes.continuousShape, DistTypes.continuousShape) => option(DistTypes.continuousShape),
+ t1IntegralCache: option(DistTypes.continuousShape),
+ t2IntegralCache: option(DistTypes.continuousShape),
+ ) => {
+ switch (t1IntegralCache, t2IntegralCache) {
| (None, _)
| (_, None) => None
| (Some(s1), Some(s2)) => combineFn(s1, s2)
diff --git a/src/distPlus/distribution/Mixed.re b/src/distPlus/distribution/Mixed.re
index b70da9b6..4c28d08c 100644
--- a/src/distPlus/distribution/Mixed.re
+++ b/src/distPlus/distribution/Mixed.re
@@ -1,7 +1,7 @@
open Distributions;
type t = DistTypes.mixedShape;
-let make = (~continuous, ~discrete): t => {continuous, discrete};
+let make = (~integralSumCache=None, ~integralCache=None, ~continuous, ~discrete): t => {continuous, discrete, integralSumCache, integralCache};
let totalLength = (t: t): int => {
let continuousLength =
@@ -11,29 +11,20 @@ let totalLength = (t: t): int => {
continuousLength + discreteLength;
};
-let scaleBy = (~scale=1.0, {discrete, continuous}: t): t => {
- let scaledDiscrete = Discrete.scaleBy(~scale, discrete);
- let scaledContinuous = Continuous.scaleBy(~scale, continuous);
- make(~discrete=scaledDiscrete, ~continuous=scaledContinuous);
+let scaleBy = (~scale=1.0, t: t): t => {
+ let scaledDiscrete = Discrete.scaleBy(~scale, t.discrete);
+ let scaledContinuous = Continuous.scaleBy(~scale, t.continuous);
+ let scaledIntegralCache = E.O.bind(t.integralCache, v => Some(Continuous.scaleBy(~scale, v)));
+ let scaledIntegralSumCache = E.O.bind(t.integralSumCache, s => Some(s *. scale));
+ make(~discrete=scaledDiscrete, ~continuous=scaledContinuous, ~integralSumCache=scaledIntegralSumCache, ~integralCache=scaledIntegralCache);
};
let toContinuous = ({continuous}: t) => Some(continuous);
let toDiscrete = ({discrete}: t) => Some(discrete);
-let combinePointwise = (~knownIntegralSumsFn, fn, t1: t, t2: t) => {
- let reducedDiscrete =
- [|t1, t2|]
- |> E.A.fmap(toDiscrete)
- |> E.A.O.concatSomes
- |> Discrete.reduce(~knownIntegralSumsFn, fn);
-
- let reducedContinuous =
- [|t1, t2|]
- |> E.A.fmap(toContinuous)
- |> E.A.O.concatSomes
- |> Continuous.reduce(~knownIntegralSumsFn, fn);
-
- make(~discrete=reducedDiscrete, ~continuous=reducedContinuous);
+let updateIntegralCache = (integralCache, t: t): t => {
+ ...t,
+ integralCache,
};
module T =
@@ -47,6 +38,8 @@ module T =
max(Continuous.T.maxX(continuous), Discrete.T.maxX(discrete));
let toShape = (t: t): DistTypes.shape => Mixed(t);
+ let updateIntegralCache = updateIntegralCache;
+
let toContinuous = toContinuous;
let toDiscrete = toDiscrete;
@@ -61,29 +54,35 @@ module T =
let truncatedDiscrete =
Discrete.T.truncate(leftCutoff, rightCutoff, discrete);
- make(~discrete=truncatedDiscrete, ~continuous=truncatedContinuous);
+ make(~integralSumCache=None, ~integralCache=None, ~discrete=truncatedDiscrete, ~continuous=truncatedContinuous);
};
let normalize = (t: t): t => {
+ let continuousIntegral = Continuous.T.Integral.get(t.continuous);
+ let discreteIntegral = Discrete.T.Integral.get(t.discrete);
+
+ let continuous = t.continuous |> Continuous.updateIntegralCache(Some(continuousIntegral));
+ let discrete = t.discrete |> Discrete.updateIntegralCache(Some(discreteIntegral));
+
let continuousIntegralSum =
- Continuous.T.Integral.sum(~cache=None, t.continuous);
+ Continuous.T.Integral.sum(continuous);
let discreteIntegralSum =
- Discrete.T.Integral.sum(~cache=None, t.discrete);
+ Discrete.T.Integral.sum(discrete);
let totalIntegralSum = continuousIntegralSum +. discreteIntegralSum;
let newContinuousSum = continuousIntegralSum /. totalIntegralSum;
let newDiscreteSum = discreteIntegralSum /. totalIntegralSum;
let normalizedContinuous =
- t.continuous
- |> Continuous.scaleBy(~scale=1. /. newContinuousSum)
- |> Continuous.updateKnownIntegralSum(Some(newContinuousSum));
+ continuous
+ |> Continuous.scaleBy(~scale=newContinuousSum /. continuousIntegralSum)
+ |> Continuous.updateIntegralSumCache(Some(newContinuousSum));
let normalizedDiscrete =
- t.discrete
- |> Discrete.scaleBy(~scale=1. /. newDiscreteSum)
- |> Discrete.updateKnownIntegralSum(Some(newDiscreteSum));
+ discrete
+ |> Discrete.scaleBy(~scale=newDiscreteSum /. discreteIntegralSum)
+ |> Discrete.updateIntegralSumCache(Some(newDiscreteSum));
- make(~continuous=normalizedContinuous, ~discrete=normalizedDiscrete);
+ make(~integralSumCache=Some(1.0), ~integralCache=None, ~continuous=normalizedContinuous, ~discrete=normalizedDiscrete);
};
let xToY = (x, t: t) => {
@@ -97,23 +96,22 @@ module T =
let toDiscreteProbabilityMassFraction = ({discrete, continuous}: t) => {
let discreteIntegralSum =
- Discrete.T.Integral.sum(~cache=None, discrete);
+ Discrete.T.Integral.sum(discrete);
let continuousIntegralSum =
- Continuous.T.Integral.sum(~cache=None, continuous);
+ Continuous.T.Integral.sum(continuous);
let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum;
discreteIntegralSum /. totalIntegralSum;
};
- let downsample = (~cache=None, count, {discrete, continuous}: t): t => {
+ let downsample = (count, t: t): t => {
// We will need to distribute the new xs fairly between the discrete and continuous shapes.
// The easiest way to do this is to simply go by the previous probability masses.
- // The cache really isn't helpful here, because we would need two separate caches
let discreteIntegralSum =
- Discrete.T.Integral.sum(~cache=None, discrete);
+ Discrete.T.Integral.sum(t.discrete);
let continuousIntegralSum =
- Continuous.T.Integral.sum(~cache=None, continuous);
+ Continuous.T.Integral.sum(t.continuous);
let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum;
// TODO: figure out what to do when the totalIntegralSum is zero.
@@ -123,7 +121,7 @@ module T =
int_of_float(
float_of_int(count) *. (discreteIntegralSum /. totalIntegralSum),
),
- discrete,
+ t.discrete,
);
let downsampledContinuous =
@@ -131,75 +129,71 @@ module T =
int_of_float(
float_of_int(count) *. (continuousIntegralSum /. totalIntegralSum),
),
- continuous,
+ t.continuous,
);
- {discrete: downsampledDiscrete, continuous: downsampledContinuous};
+ {...t, discrete: downsampledDiscrete, continuous: downsampledContinuous};
};
- let normalizedToContinuous = (t: t) => Some(normalize(t).continuous);
-
- let normalizedToDiscrete = ({discrete} as t: t) =>
- Some(normalize(t).discrete);
-
- let integral = (~cache, {continuous, discrete}: t) => {
- switch (cache) {
+ let integral = (t: t) => {
+ switch (t.integralCache) {
| Some(cache) => cache
| None =>
- // note: if the underlying shapes aren't normalized, then these integrals won't be either!
- let continuousIntegral =
- Continuous.T.Integral.get(~cache=None, continuous);
- let discreteIntegral = Discrete.T.Integral.get(~cache=None, discrete);
+ // note: if the underlying shapes aren't normalized, then these integrals won't be either -- but that's the way it should be.
+ let continuousIntegral = Continuous.T.Integral.get(t.continuous);
+ let discreteIntegral = Continuous.stepwiseToLinear(Discrete.T.Integral.get(t.discrete));
Continuous.make(
- `Linear,
- XYShape.PointwiseCombination.combineLinear(
- ~fn=(+.),
+ XYShape.PointwiseCombination.combine(
+ (+.),
+ XYShape.XtoY.continuousInterpolator(`Linear, `UseOutermostPoints),
Continuous.getShape(continuousIntegral),
Continuous.getShape(discreteIntegral),
),
- None,
);
};
};
- let integralEndY = (~cache, t: t) => {
- integral(~cache, t) |> Continuous.lastY;
+ let integralEndY = (t: t) => {
+ t |> integral |> Continuous.lastY;
};
- let integralXtoY = (~cache, f, t) => {
- t |> integral(~cache) |> Continuous.getShape |> XYShape.XtoY.linear(f);
+ let integralXtoY = (f, t) => {
+ t |> integral |> Continuous.getShape |> XYShape.XtoY.linear(f);
};
- let integralYtoX = (~cache, f, t) => {
- t |> integral(~cache) |> Continuous.getShape |> XYShape.YtoX.linear(f);
+ let integralYtoX = (f, t) => {
+ t |> integral |> Continuous.getShape |> XYShape.YtoX.linear(f);
};
// This pipes all ys (continuous and discrete) through fn.
- // If mapY is a linear operation, we might be able to update the knownIntegralSums as well;
+ // If mapY is a linear operation, we might be able to update the integralSumCaches as well;
// if not, they'll be set to None.
let mapY =
(
- ~knownIntegralSumFn=previousIntegralSum => None,
- fn,
- {discrete, continuous}: t,
+ ~integralSumCacheFn=previousIntegralSum => None,
+ ~integralCacheFn=previousIntegral => None,
+ ~fn,
+ t: t,
)
: t => {
- let u = E.O.bind(_, knownIntegralSumFn);
+ let yMappedDiscrete: DistTypes.discreteShape =
+ t.discrete
+ |> Discrete.T.mapY(~fn)
+ |> Discrete.updateIntegralSumCache(E.O.bind(t.discrete.integralSumCache, integralSumCacheFn))
+ |> Discrete.updateIntegralCache(E.O.bind(t.discrete.integralCache, integralCacheFn));
- let yMappedDiscrete =
- discrete
- |> Discrete.T.mapY(fn)
- |> Discrete.updateKnownIntegralSum(u(discrete.knownIntegralSum));
-
- let yMappedContinuous =
- continuous
- |> Continuous.T.mapY(fn)
- |> Continuous.updateKnownIntegralSum(u(continuous.knownIntegralSum));
+ let yMappedContinuous: DistTypes.continuousShape =
+ t.continuous
+ |> Continuous.T.mapY(~fn)
+ |> Continuous.updateIntegralSumCache(E.O.bind(t.continuous.integralSumCache, integralSumCacheFn))
+ |> Continuous.updateIntegralCache(E.O.bind(t.continuous.integralCache, integralCacheFn));
{
discrete: yMappedDiscrete,
- continuous: Continuous.T.mapY(fn, continuous),
+ continuous: yMappedContinuous,
+ integralSumCache: E.O.bind(t.integralSumCache, integralSumCacheFn),
+ integralCache: E.O.bind(t.integralCache, integralCacheFn),
};
};
@@ -208,10 +202,8 @@ module T =
let continuousMean = Continuous.T.mean(continuous);
// the combined mean is the weighted sum of the two:
- let discreteIntegralSum =
- Discrete.T.Integral.sum(~cache=None, discrete);
- let continuousIntegralSum =
- Continuous.T.Integral.sum(~cache=None, continuous);
+ let discreteIntegralSum = Discrete.T.Integral.sum(discrete);
+ let continuousIntegralSum = Continuous.T.Integral.sum(continuous);
let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum;
(
@@ -225,10 +217,8 @@ module T =
let variance = ({discrete, continuous} as t: t): float => {
// the combined mean is the weighted sum of the two:
- let discreteIntegralSum =
- Discrete.T.Integral.sum(~cache=None, discrete);
- let continuousIntegralSum =
- Continuous.T.Integral.sum(~cache=None, continuous);
+ let discreteIntegralSum = Discrete.T.Integral.sum(discrete);
+ let continuousIntegralSum = Continuous.T.Integral.sum(continuous);
let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum;
let getMeanOfSquares = ({discrete, continuous}: t) => {
@@ -279,27 +269,64 @@ let combineAlgebraically =
let ccConvResult =
Continuous.combineAlgebraically(
op,
- t1d.continuous,
- t2d.continuous,
+ t1.continuous,
+ t2.continuous,
);
let dcConvResult =
Continuous.combineAlgebraicallyWithDiscrete(
op,
- t2d.continuous,
- t1d.discrete,
+ t2.continuous,
+ t1.discrete,
);
let cdConvResult =
Continuous.combineAlgebraicallyWithDiscrete(
op,
- t1d.continuous,
- t2d.discrete,
+ t1.continuous,
+ t2.discrete,
);
let continuousConvResult =
Continuous.reduce((+.), [|ccConvResult, dcConvResult, cdConvResult|]);
// ... finally, discrete (*) discrete => discrete, obviously:
let discreteConvResult =
- Discrete.combineAlgebraically(op, t1d.discrete, t2d.discrete);
+ Discrete.combineAlgebraically(op, t1.discrete, t2.discrete);
- {discrete: discreteConvResult, continuous: continuousConvResult};
+ let combinedIntegralSum =
+ Common.combineIntegralSums(
+ (a, b) => Some(a *. b),
+ t1.integralSumCache,
+ t2.integralSumCache,
+ );
+
+ {discrete: discreteConvResult, continuous: continuousConvResult, integralSumCache: combinedIntegralSum, integralCache: None};
+};
+
+let combinePointwise = (~integralSumCachesFn = (_, _) => None, ~integralCachesFn = (_, _) => None, fn, t1: t, t2: t): t => {
+ let reducedDiscrete =
+ [|t1, t2|]
+ |> E.A.fmap(toDiscrete)
+ |> E.A.O.concatSomes
+ |> Discrete.reduce(~integralSumCachesFn, ~integralCachesFn, fn);
+
+ let reducedContinuous =
+ [|t1, t2|]
+ |> E.A.fmap(toContinuous)
+ |> E.A.O.concatSomes
+ |> Continuous.reduce(~integralSumCachesFn, ~integralCachesFn, fn);
+
+ let combinedIntegralSum =
+ Common.combineIntegralSums(
+ integralSumCachesFn,
+ t1.integralSumCache,
+ t2.integralSumCache,
+ );
+
+ let combinedIntegral =
+ Common.combineIntegrals(
+ integralCachesFn,
+ t1.integralCache,
+ t2.integralCache,
+ );
+
+ make(~integralSumCache=combinedIntegralSum, ~integralCache=combinedIntegral, ~discrete=reducedDiscrete, ~continuous=reducedContinuous);
};
diff --git a/src/distPlus/distribution/MixedShapeBuilder.re b/src/distPlus/distribution/MixedShapeBuilder.re
index 7f144ca9..97b7269f 100644
--- a/src/distPlus/distribution/MixedShapeBuilder.re
+++ b/src/distPlus/distribution/MixedShapeBuilder.re
@@ -9,8 +9,8 @@ type assumptions = {
};
let buildSimple = (~continuous: option(DistTypes.continuousShape), ~discrete: option(DistTypes.discreteShape)): option(DistTypes.shape) => {
- let continuous = continuous |> E.O.default(Continuous.make(`Linear, {xs: [||], ys: [||]}, Some(0.0)));
- let discrete = discrete |> E.O.default(Discrete.make({xs: [||], ys: [||]}, Some(0.0)));
+ let continuous = continuous |> E.O.default(Continuous.make(~integralSumCache=Some(0.0), {xs: [||], ys: [||]}));
+ let discrete = discrete |> E.O.default(Discrete.make(~integralSumCache=Some(0.0), {xs: [||], ys: [||]}));
let cLength =
continuous
|> Continuous.getShape
@@ -22,15 +22,13 @@ let buildSimple = (~continuous: option(DistTypes.continuousShape), ~discrete: op
| (0 | 1, _) => Some(Discrete(discrete))
| (_, 0) => Some(Continuous(continuous))
| (_, _) =>
- let discreteProbabilityMassFraction =
- Discrete.T.Integral.sum(~cache=None, discrete);
- let discrete = Discrete.T.normalize(discrete);
- let continuous = Continuous.T.normalize(continuous);
let mixedDist =
Mixed.make(
+ ~integralSumCache=None,
+ ~integralCache=None,
~continuous,
- ~discrete
+ ~discrete,
);
Some(Mixed(mixedDist));
};
-};
\ No newline at end of file
+};
diff --git a/src/distPlus/distribution/Shape.re b/src/distPlus/distribution/Shape.re
index 61c49947..01685c3c 100644
--- a/src/distPlus/distribution/Shape.re
+++ b/src/distPlus/distribution/Shape.re
@@ -15,42 +15,54 @@ let fmap = ((fn1, fn2, fn3), t: t): t =>
| Continuous(m) => Continuous(fn3(m))
};
+
let toMixed =
mapToAll((
m => m,
- d => Mixed.make(~discrete=d, ~continuous=Continuous.empty),
- c => Mixed.make(~discrete=Discrete.empty, ~continuous=c),
+ d => Mixed.make(~integralSumCache=d.integralSumCache, ~integralCache=d.integralCache, ~discrete=d, ~continuous=Continuous.empty),
+ c => Mixed.make(~integralSumCache=c.integralSumCache, ~integralCache=c.integralCache, ~discrete=Discrete.empty, ~continuous=c),
));
let combineAlgebraically =
(op: ExpressionTypes.algebraicOperation, t1: t, t2: t): t => {
switch (t1, t2) {
| (Continuous(m1), Continuous(m2)) =>
- DistTypes.Continuous(Continuous.combineAlgebraically(op, m1, m2))
+ Continuous.combineAlgebraically(op, m1, m2) |> Continuous.T.toShape;
+ | (Continuous(m1), Discrete(m2))
+ | (Discrete(m2), Continuous(m1)) =>
+ Continuous.combineAlgebraicallyWithDiscrete(op, m1, m2) |> Continuous.T.toShape
| (Discrete(m1), Discrete(m2)) =>
- DistTypes.Discrete(Discrete.combineAlgebraically(op, m1, m2))
+ Discrete.combineAlgebraically(op, m1, m2) |> Discrete.T.toShape
| (m1, m2) =>
- DistTypes.Mixed(
- Mixed.combineAlgebraically(op, toMixed(m1), toMixed(m2)),
+ Mixed.combineAlgebraically(
+ op,
+ toMixed(m1),
+ toMixed(m2),
)
+ |> Mixed.T.toShape
};
};
let combinePointwise =
- (~knownIntegralSumsFn=(_, _) => None, fn, t1: t, t2: t) =>
+ (~integralSumCachesFn: (float, float) => option(float) = (_, _) => None,
+ ~integralCachesFn: (DistTypes.continuousShape, DistTypes.continuousShape) => option(DistTypes.continuousShape) = (_, _) => None,
+ fn,
+ t1: t,
+ t2: t) =>
switch (t1, t2) {
| (Continuous(m1), Continuous(m2)) =>
DistTypes.Continuous(
- Continuous.combinePointwise(~knownIntegralSumsFn, fn, m1, m2),
+ Continuous.combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn, m1, m2),
)
| (Discrete(m1), Discrete(m2)) =>
DistTypes.Discrete(
- Discrete.combinePointwise(~knownIntegralSumsFn, fn, m1, m2),
+ Discrete.combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn, m1, m2),
)
| (m1, m2) =>
DistTypes.Mixed(
Mixed.combinePointwise(
- ~knownIntegralSumsFn,
+ ~integralSumCachesFn,
+ ~integralCachesFn,
fn,
toMixed(m1),
toMixed(m2),
@@ -58,15 +70,6 @@ let combinePointwise =
)
};
-// TODO: implement these functions
-let pdf = (f: float, t: t): float => {
- 0.0;
-};
-
-let inv = (f: float, t: t): float => {
- 0.0;
-};
-
module T =
Dist({
type t = DistTypes.shape;
@@ -84,7 +87,7 @@ module T =
let toContinuous = t => None;
let toDiscrete = t => None;
- let downsample = (~cache=None, i, t) =>
+ let downsample = (i, t) =>
fmap(
(
Mixed.T.downsample(i),
@@ -105,8 +108,21 @@ module T =
);
let toDiscreteProbabilityMassFraction = t => 0.0;
+
let normalize =
- fmap((Mixed.T.normalize, Discrete.T.normalize, Continuous.T.normalize));
+ fmap((
+ Mixed.T.normalize,
+ Discrete.T.normalize,
+ Continuous.T.normalize
+ ));
+
+ let updateIntegralCache = (integralCache, t: t): t =>
+ fmap((
+ Mixed.T.updateIntegralCache(integralCache),
+ Discrete.T.updateIntegralCache(integralCache),
+ Continuous.T.updateIntegralCache(integralCache),
+ ), t);
+
let toContinuous =
mapToAll((
Mixed.T.toContinuous,
@@ -127,51 +143,39 @@ module T =
Continuous.T.toDiscreteProbabilityMassFraction,
));
- let normalizedToDiscrete =
- mapToAll((
- Mixed.T.normalizedToDiscrete,
- Discrete.T.normalizedToDiscrete,
- Continuous.T.normalizedToDiscrete,
- ));
- let normalizedToContinuous =
- mapToAll((
- Mixed.T.normalizedToContinuous,
- Discrete.T.normalizedToContinuous,
- Continuous.T.normalizedToContinuous,
- ));
let minX = mapToAll((Mixed.T.minX, Discrete.T.minX, Continuous.T.minX));
- let integral = (~cache) =>
+ let integral =
mapToAll((
- Mixed.T.Integral.get(~cache=None),
- Discrete.T.Integral.get(~cache=None),
- Continuous.T.Integral.get(~cache=None),
+ Mixed.T.Integral.get,
+ Discrete.T.Integral.get,
+ Continuous.T.Integral.get,
));
- let integralEndY = (~cache) =>
+ let integralEndY =
mapToAll((
- Mixed.T.Integral.sum(~cache=None),
- Discrete.T.Integral.sum(~cache),
- Continuous.T.Integral.sum(~cache=None),
+ Mixed.T.Integral.sum,
+ Discrete.T.Integral.sum,
+ Continuous.T.Integral.sum,
));
- let integralXtoY = (~cache, f) => {
+ let integralXtoY = (f) => {
mapToAll((
- Mixed.T.Integral.xToY(~cache, f),
- Discrete.T.Integral.xToY(~cache, f),
- Continuous.T.Integral.xToY(~cache, f),
+ Mixed.T.Integral.xToY(f),
+ Discrete.T.Integral.xToY(f),
+ Continuous.T.Integral.xToY(f),
));
};
- let integralYtoX = (~cache, f) => {
+ let integralYtoX = (f) => {
mapToAll((
- Mixed.T.Integral.yToX(~cache, f),
- Discrete.T.Integral.yToX(~cache, f),
- Continuous.T.Integral.yToX(~cache, f),
+ Mixed.T.Integral.yToX(f),
+ Discrete.T.Integral.yToX(f),
+ Continuous.T.Integral.yToX(f),
));
};
let maxX = mapToAll((Mixed.T.maxX, Discrete.T.maxX, Continuous.T.maxX));
- let mapY = (~knownIntegralSumFn=previousIntegralSum => None, fn) =>
+ let mapY = (~integralSumCacheFn=previousIntegralSum => None, ~integralCacheFn=previousIntegral=>None, ~fn) =>
fmap((
- Mixed.T.mapY(~knownIntegralSumFn, fn),
- Discrete.T.mapY(~knownIntegralSumFn, fn),
- Continuous.T.mapY(~knownIntegralSumFn, fn),
+ Mixed.T.mapY(~integralSumCacheFn, ~integralCacheFn, ~fn),
+ Discrete.T.mapY(~integralSumCacheFn, ~integralCacheFn, ~fn),
+ Continuous.T.mapY(~integralSumCacheFn, ~integralCacheFn, ~fn),
));
let mean = (t: t): float =>
@@ -189,6 +193,14 @@ module T =
};
});
+let pdf = (f: float, t: t) => {
+ let mixedPoint: DistTypes.mixedPoint = T.xToY(f, t);
+ mixedPoint.continuous +. mixedPoint.discrete;
+};
+
+let inv = T.Integral.yToX;
+let cdf = T.Integral.xToY;
+
let doN = (n, fn) => {
let items = Belt.Array.make(n, 0.0);
for (x in 0 to n - 1) {
@@ -198,21 +210,24 @@ let doN = (n, fn) => {
items;
};
-let sample = (cache, t: t): float => {
+let sample = (t: t): float => {
let randomItem = Random.float(1.);
- let bar = T.Integral.yToX(~cache, randomItem, t);
+ let bar = t |> T.Integral.yToX(randomItem);
bar;
};
let sampleNRendered = (n, dist) => {
- let integralCache = T.Integral.get(~cache=None, dist);
- doN(n, () => sample(Some(integralCache), dist));
+ let integralCache = T.Integral.get(dist);
+ let distWithUpdatedIntegralCache = T.updateIntegralCache(Some(integralCache), dist);
+
+ doN(n, () => sample(distWithUpdatedIntegralCache));
};
-let operate = (distToFloatOp: ExpressionTypes.distToFloatOperation, s) =>
+let operate = (distToFloatOp: ExpressionTypes.distToFloatOperation, s): float =>
switch (distToFloatOp) {
| `Pdf(f) => pdf(f, s)
+ | `Cdf(f) => pdf(f, s)
| `Inv(f) => inv(f, s)
- | `Sample => sample(None, s)
+ | `Sample => sample(s)
| `Mean => T.mean(s)
};
diff --git a/src/distPlus/distribution/XYShape.re b/src/distPlus/distribution/XYShape.re
index 0fe1ebe0..92cffe60 100644
--- a/src/distPlus/distribution/XYShape.re
+++ b/src/distPlus/distribution/XYShape.re
@@ -19,6 +19,7 @@ module T = {
let ys = (t: t) => t.ys;
let length = (t: t) => E.A.length(t.xs);
let empty = {xs: [||], ys: [||]};
+ let isEmpty = (t: t) => length(t) == 0;
let minX = (t: t) => t |> xs |> E.A.Sorted.min |> extImp;
let maxX = (t: t) => t |> xs |> E.A.Sorted.max |> extImp;
let firstY = (t: t) => t |> ys |> E.A.first |> extImp;
@@ -32,6 +33,11 @@ module T = {
let accumulateYs = (fn, p: t) => {
fromArray((p.xs, E.A.accumulate(fn, p.ys)));
};
+ let concat = (t1: t, t2: t) => {
+ let cxs = Array.concat([t1.xs, t2.xs]);
+ let cys = Array.concat([t1.ys, t2.ys]);
+ {xs: cxs, ys: cys};
+ };
let fromZippedArray = (pairs: array((float, float))): t =>
pairs |> Belt.Array.unzip |> fromArray;
let equallyDividedXs = (t: t, newLength) => {
@@ -137,6 +143,63 @@ module XtoY = {
};
n;
};
+
+ /* Returns a between-points-interpolating function that can be used with PointwiseCombination.combine.
+ Interpolation can either be stepwise (using the value on the left) or linear. Extrapolation can be `UseZero or `UseOutermostPoints. */
+ let continuousInterpolator = (interpolation: DistTypes.interpolationStrategy, extrapolation: DistTypes.extrapolationStrategy): interpolator => {
+ switch (interpolation, extrapolation) {
+ | (`Linear, `UseZero) => (t: T.t, leftIndex: int, x: float) => {
+ if (leftIndex < 0) {
+ 0.0
+ } else if (leftIndex >= T.length(t) - 1) {
+ 0.0
+ } else {
+ let x1 = t.xs[leftIndex];
+ let x2 = t.xs[leftIndex + 1];
+ let y1 = t.ys[leftIndex];
+ let y2 = t.ys[leftIndex + 1];
+ let fraction = (x -. x1) /. (x2 -. x1);
+ y1 *. (1. -. fraction) +. y2 *. fraction;
+ };
+ }
+ | (`Linear, `UseOutermostPoints) => (t: T.t, leftIndex: int, x: float) => {
+ if (leftIndex < 0) {
+ t.ys[0];
+ } else if (leftIndex >= T.length(t) - 1) {
+ t.ys[T.length(t) - 1]
+ } else {
+ let x1 = t.xs[leftIndex];
+ let x2 = t.xs[leftIndex + 1];
+ let y1 = t.ys[leftIndex];
+ let y2 = t.ys[leftIndex + 1];
+ let fraction = (x -. x1) /. (x2 -. x1);
+ y1 *. (1. -. fraction) +. y2 *. fraction;
+ };
+ }
+ | (`Stepwise, `UseZero) => (t: T.t, leftIndex: int, x: float) => {
+ if (leftIndex < 0) {
+ 0.0
+ } else if (leftIndex >= T.length(t) - 1) {
+ 0.0
+ } else {
+ t.ys[leftIndex];
+ }
+ }
+ | (`Stepwise, `UseOutermostPoints) => (t: T.t, leftIndex: int, x: float) => {
+ if (leftIndex < 0) {
+ t.ys[0];
+ } else if (leftIndex >= T.length(t) - 1) {
+ t.ys[T.length(t) - 1]
+ } else {
+ t.ys[leftIndex];
+ }
+ }
+ }
+ };
+
+ /* Returns a between-points-interpolating function that can be used with PointwiseCombination.combine.
+ For discrete distributions, the probability density between points is zero, so we just return zero here. */
+ let discreteInterpolator: interpolator = (t: T.t, leftIndex: int, x: float) => 0.0;
};
module XsConversion = {
@@ -171,24 +234,22 @@ module Zipped = {
};
module PointwiseCombination = {
- type xsSelection =
- | ALL_XS
- | XS_EVENLY_DIVIDED(int);
- let combineLinear = [%raw {| // : (float => float => float, T.t, T.t) => T.t
+ // t1Interpolator and t2Interpolator are functions from XYShape.XtoY, e.g. linearBetweenPointsExtrapolateFlat.
+ let combine = [%raw {| // : (float => float => float, T.t, T.t, bool) => T.t
// This function combines two xyShapes by looping through both of them simultaneously.
// It always moves on to the next smallest x, whether that's in the first or second input's xs,
// and interpolates the value on the other side, thus accumulating xs and ys.
- // In this implementation (unlike in XtoY.linear above), values outside of a shape's xs are considered 0.0 (instead of firstY or lastY).
// This is written in raw JS because this can still be a bottleneck, and using refs for the i and j indices is quite painful.
- function(fn, t1, t2) {
+ function(fn, interpolator, t1, t2) {
let t1n = t1.xs.length;
let t2n = t2.xs.length;
let outX = [];
let outY = [];
let i = -1;
let j = -1;
+
while (i <= t1n - 1 && j <= t2n - 1) {
let x, ya, yb;
if (j == t2n - 1 && i < t1n - 1 ||
@@ -198,8 +259,7 @@ module PointwiseCombination = {
x = t1.xs[i];
ya = t1.ys[i];
- let bFraction = (x - (t2.xs[j] || x)) / ((t2.xs[j+1] || x) - (t2.xs[j] || 0));
- yb = (t2.ys[j] || 0) * (1-bFraction) + (t2.ys[j+1] || 0) * bFraction;
+ yb = interpolator(t2, j, x);
} else if (i == t1n - 1 && j < t2n - 1 ||
t1.xs[i+1] > t2.xs[j+1]) { // if b has to catch up to a, or if a is already done
j++;
@@ -207,8 +267,7 @@ module PointwiseCombination = {
x = t2.xs[j];
yb = t2.ys[j];
- let aFraction = (x - (t1.xs[i] || x)) / ((t1.xs[i+1] || x) - (t1.xs[i] || 0));
- ya = (t1.ys[i] || 0) * (1-aFraction) + (t1.ys[i+1] || 0) * aFraction;
+ ya = interpolator(t1, i, x);
} else if (i < t1n - 1 && j < t2n && t1.xs[i+1] === t2.xs[j+1]) { // if they happen to be equal, move both ahead
i++;
j++;
@@ -227,15 +286,16 @@ module PointwiseCombination = {
outX.push(x);
outY.push(fn(ya, yb));
}
+
return {xs: outX, ys: outY};
}
|}];
- let combine =
+ let combineEvenXs =
(
- ~xToYSelection: (float, T.t) => 'a,
- ~xsSelection=ALL_XS,
~fn,
+ ~xToYSelection,
+ sampleCount,
t1: T.t,
t2: T.t,
) => {
@@ -245,25 +305,15 @@ module PointwiseCombination = {
| (0, _) => t2
| (_, 0) => t1
| (_, _) => {
- let allXs =
- switch (xsSelection) {
- | ALL_XS => Ts.allXs([|t1, t2|])
- | XS_EVENLY_DIVIDED(sampleCount) =>
- Ts.equallyDividedXs([|t1, t2|], sampleCount)
- };
+ let allXs = Ts.equallyDividedXs([|t1, t2|], sampleCount);
- let allYs =
- allXs |> E.A.fmap(x => fn(xToYSelection(x, t1), xToYSelection(x, t2)));
+ let allYs = allXs |> E.A.fmap(x => fn(xToYSelection(x, t1), xToYSelection(x, t2)));
T.fromArrays(allXs, allYs);
}
}
};
- //let combineLinear = combine(~xToYSelection=XtoY.linear);
- let combineStepwise = combine(~xToYSelection=XtoY.stepwiseIncremental);
- let combineIfAtX = combine(~xToYSelection=XtoY.stepwiseIfAtX);
-
// TODO: I'd bet this is pretty slow. Maybe it would be faster to intersperse Xs and Ys separately.
let intersperse = (t1: T.t, t2: T.t) => {
E.A.intersperse(T.zip(t1), T.zip(t2)) |> T.fromZippedArray;
@@ -324,8 +374,31 @@ module Range = {
let derivative = mapYsBasedOnRanges(delta_y_over_delta_x);
- // TODO: It would be nicer if this the diff didn't change the first element, and also maybe if there were a more elegant way of doing this.
+ let stepwiseToLinear = ({xs, ys}: T.t): T.t => {
+ // adds points at the bottom of each step.
+ let length = E.A.length(xs);
+ let newXs: array(float) = Belt.Array.makeUninitializedUnsafe(2 * length);
+ let newYs: array(float) = Belt.Array.makeUninitializedUnsafe(2 * length);
+
+ Belt.Array.set(newXs, 0, xs[0] -. epsilon_float) |> ignore;
+ Belt.Array.set(newYs, 0, 0.) |> ignore;
+ Belt.Array.set(newXs, 1, xs[0]) |> ignore;
+ Belt.Array.set(newYs, 1, ys[0]) |> ignore;
+
+ for (i in 1 to E.A.length(xs) - 1) {
+ Belt.Array.set(newXs, i * 2, xs[i] -. epsilon_float) |> ignore;
+ Belt.Array.set(newYs, i * 2, ys[i-1]) |> ignore;
+ Belt.Array.set(newXs, i * 2 + 1, xs[i]) |> ignore;
+ Belt.Array.set(newYs, i * 2 + 1, ys[i]) |> ignore;
+ ();
+ };
+
+ {xs: newXs, ys: newYs};
+ };
+
+ // TODO: I think this isn't needed by any functions anymore.
let stepsToContinuous = t => {
+ // TODO: It would be nicer if this the diff didn't change the first element, and also maybe if there were a more elegant way of doing this.
let diff = T.xTotalRange(t) |> (r => r *. 0.00001);
let items =
switch (E.A.toRanges(Belt.Array.zip(t.xs, t.ys))) {
@@ -356,10 +429,10 @@ let pointLogScore = (prediction, answer) =>
};
let logScorePoint = (sampleCount, t1, t2) =>
- PointwiseCombination.combine(
- ~xsSelection=XS_EVENLY_DIVIDED(sampleCount),
- ~xToYSelection=XtoY.linear,
+ PointwiseCombination.combineEvenXs(
~fn=pointLogScore,
+ ~xToYSelection=XtoY.linear,
+ sampleCount,
t1,
t2,
)
diff --git a/src/distPlus/expressionTree/ExpressionTreeEvaluator.re b/src/distPlus/expressionTree/ExpressionTreeEvaluator.re
index 798ccaf6..fd6d06ee 100644
--- a/src/distPlus/expressionTree/ExpressionTreeEvaluator.re
+++ b/src/distPlus/expressionTree/ExpressionTreeEvaluator.re
@@ -20,61 +20,15 @@ module AlgebraicCombination = {
| _ => Ok(`AlgebraicCombination((operation, t1, t2)))
};
- let tryCombination = (n, algebraicOp, t1: node, t2: node) => {
- let sampleN =
- mapRenderable(Shape.sampleNRendered(n), SymbolicDist.T.sampleN(n));
- switch (sampleN(t1), sampleN(t2)) {
- | (Some(a), Some(b)) =>
- Some(
- Belt.Array.zip(a, b)
- |> E.A.fmap(((a, b)) => Operation.Algebraic.toFn(algebraicOp, a, b)),
- )
- | _ => None
- };
- };
-
- let renderIfNotRendered = (params, t) =>
- !renderable(t)
- ? switch (render(params, t)) {
- | Ok(r) => Ok(r)
- | Error(e) => Error(e)
- }
- : Ok(t);
-
- let combineAsShapes =
- (evaluationParams: evaluationParams, algebraicOp, t1: node, t2: node) => {
- let i1 = renderIfNotRendered(evaluationParams, t1);
- let i2 = renderIfNotRendered(evaluationParams, t2);
- E.R.merge(i1, i2)
- |> E.R.bind(
- _,
- ((a, b)) => {
- let samples =
- tryCombination(
- evaluationParams.samplingInputs.sampleCount,
- algebraicOp,
- a,
- b,
- );
- let shape =
- samples
- |> E.O.fmap(
- Samples.T.fromSamples(
- ~samplingInputs={
- sampleCount:
- Some(evaluationParams.samplingInputs.sampleCount),
- outputXYPoints:
- Some(evaluationParams.samplingInputs.outputXYPoints),
- kernelWidth: evaluationParams.samplingInputs.kernelWidth,
- },
- ),
- )
- |> E.O.bind(_, (r: RenderTypes.ShapeRenderer.Sampling.outputs) =>
- r.shape
- )
- |> E.O.toResult("No response");
- shape |> E.R.fmap(r => `Normalize(`RenderedDist(r)));
- },
+ let combinationByRendering =
+ (evaluationParams, algebraicOp, t1: node, t2: node)
+ : result(node, string) => {
+ E.R.merge(
+ Render.ensureIsRenderedAndGetShape(evaluationParams, t1),
+ Render.ensureIsRenderedAndGetShape(evaluationParams, t2),
+ )
+ |> E.R.fmap(((a, b)) =>
+ `RenderedDist(Shape.combineAlgebraically(algebraicOp, a, b))
);
};
@@ -92,7 +46,7 @@ module AlgebraicCombination = {
_,
fun
| `SymbolicDist(d) as t => Ok(t)
- | _ => combineAsShapes(evaluationParams, algebraicOp, t1, t2),
+ | _ => combinationByRendering(evaluationParams, algebraicOp, t1, t2),
);
};
@@ -101,16 +55,18 @@ module VerticalScaling = {
(evaluationParams: evaluationParams, scaleOp, t, scaleBy) => {
// scaleBy has to be a single float, otherwise we'll return an error.
let fn = Operation.Scale.toFn(scaleOp);
- let knownIntegralSumFn = Operation.Scale.toKnownIntegralSumFn(scaleOp);
- let renderedShape = render(evaluationParams, t);
+ let integralSumCacheFn = Operation.Scale.toIntegralSumCacheFn(scaleOp);
+ let integralCacheFn = Operation.Scale.toIntegralCacheFn(scaleOp);
+ let renderedShape = Render.render(evaluationParams, t);
switch (renderedShape, scaleBy) {
| (Ok(`RenderedDist(rs)), `SymbolicDist(`Float(sm))) =>
Ok(
`RenderedDist(
Shape.T.mapY(
- ~knownIntegralSumFn=knownIntegralSumFn(sm),
- fn(sm),
+ ~integralSumCacheFn=integralSumCacheFn(sm),
+ ~integralCacheFn=integralCacheFn(sm),
+ ~fn=fn(sm),
rs,
),
),
@@ -122,13 +78,23 @@ module VerticalScaling = {
};
module PointwiseCombination = {
- let pointwiseAdd = (evaluationParams: evaluationParams, t1, t2) => {
- switch (render(evaluationParams, t1), render(evaluationParams, t2)) {
+ let pointwiseAdd = (evaluationParams: evaluationParams, t1: t, t2: t) => {
+ switch (Render.render(evaluationParams, t1), Render.render(evaluationParams, t2)) {
| (Ok(`RenderedDist(rs1)), Ok(`RenderedDist(rs2))) =>
Ok(
`RenderedDist(
Shape.combinePointwise(
- ~knownIntegralSumsFn=(a, b) => Some(a +. b),
+ ~integralSumCachesFn=(a, b) => Some(a +. b),
+ ~integralCachesFn=
+ (a, b) =>
+ Some(
+ Continuous.combinePointwise(
+ ~distributionType=`CDF,
+ (+.),
+ a,
+ b,
+ ),
+ ),
(+.),
rs1,
rs2,
@@ -141,7 +107,7 @@ module PointwiseCombination = {
};
};
- let pointwiseMultiply = (evaluationParams: evaluationParams, t1, t2) => {
+ let pointwiseMultiply = (evaluationParams: evaluationParams, t1: t, t2: t) => {
// TODO: construct a function that we can easily sample from, to construct
// a RenderedDist. Use the xMin and xMax of the rendered shapes to tell the sampling function where to look.
Error(
@@ -150,7 +116,12 @@ module PointwiseCombination = {
};
let operationToLeaf =
- (evaluationParams: evaluationParams, pointwiseOp, t1, t2) => {
+ (
+ evaluationParams: evaluationParams,
+ pointwiseOp: pointwiseOperation,
+ t1: t,
+ t2: t,
+ ) => {
switch (pointwiseOp) {
| `Add => pointwiseAdd(evaluationParams, t1, t2)
| `Multiply => pointwiseMultiply(evaluationParams, t1, t2)
@@ -163,7 +134,7 @@ module Truncate = {
switch (leftCutoff, rightCutoff, t) {
| (None, None, t) => `Solution(t)
| (Some(lc), Some(rc), t) when lc > rc =>
- `Error("Left truncation bound must be smaller than right bound.")
+ `Error("Left truncation bound must be smaller than right truncation bound.")
| (lc, rc, `SymbolicDist(`Uniform(u))) =>
`Solution(
`SymbolicDist(`Uniform(SymbolicDist.Uniform.truncate(lc, rc, u))),
@@ -174,9 +145,9 @@ module Truncate = {
let truncateAsShape =
(evaluationParams: evaluationParams, leftCutoff, rightCutoff, t) => {
- // TODO: use named args in renderToShape; if we're lucky we can at least get the tail
+ // TODO: use named args for xMin/xMax in renderToShape; if we're lucky we can at least get the tail
// of a distribution we otherwise wouldn't get at all
- switch (render(evaluationParams, t)) {
+ switch (Render.ensureIsRendered(evaluationParams, t)) {
| Ok(`RenderedDist(rs)) =>
Ok(`RenderedDist(Shape.T.truncate(leftCutoff, rightCutoff, rs)))
| Error(e) => Error(e)
diff --git a/src/distPlus/expressionTree/ExpressionTypes.re b/src/distPlus/expressionTree/ExpressionTypes.re
index b36548e3..a76aee74 100644
--- a/src/distPlus/expressionTree/ExpressionTypes.re
+++ b/src/distPlus/expressionTree/ExpressionTypes.re
@@ -1,7 +1,13 @@
type algebraicOperation = [ | `Add | `Multiply | `Subtract | `Divide];
type pointwiseOperation = [ | `Add | `Multiply];
type scaleOperation = [ | `Multiply | `Exponentiate | `Log];
-type distToFloatOperation = [ | `Pdf(float) | `Inv(float) | `Mean | `Sample];
+type distToFloatOperation = [
+ | `Pdf(float)
+ | `Cdf(float)
+ | `Inv(float)
+ | `Mean
+ | `Sample
+];
module ExpressionTree = {
type node = [
@@ -31,26 +37,42 @@ module ExpressionTree = {
let evaluateNode = (evaluationParams: evaluationParams) =>
evaluationParams.evaluateNode(evaluationParams);
- let render = (evaluationParams: evaluationParams, r) =>
- evaluateNode(evaluationParams, `Render(r));
-
let evaluateAndRetry = (evaluationParams, fn, node) =>
node
|> evaluationParams.evaluateNode(evaluationParams)
|> E.R.bind(_, fn(evaluationParams));
- let renderable =
- fun
- | `SymbolicDist(_) => true
- | `RenderedDist(_) => true
- | _ => false;
+ module Render = {
+ type t = node;
+
+ let render = (evaluationParams: evaluationParams, r) =>
+ `Render(r) |> evaluateNode(evaluationParams);
+
+ let ensureIsRendered = (params, t) =>
+ switch (t) {
+ | `RenderedDist(_) => Ok(t)
+ | _ =>
+ switch (render(params, t)) {
+ | Ok(`RenderedDist(r)) => Ok(`RenderedDist(r))
+ | Ok(_) => Error("Did not render as requested")
+ | Error(e) => Error(e)
+ }
+ };
+
+ let ensureIsRenderedAndGetShape = (params, t) =>
+ switch (ensureIsRendered(params, t)) {
+ | Ok(`RenderedDist(r)) => Ok(r)
+ | Ok(_) => Error("Did not render as requested")
+ | Error(e) => Error(e)
+ };
+
+ let getShape = (item: node) =>
+ switch (item) {
+ | `RenderedDist(r) => Some(r)
+ | _ => None
+ };
+ };
- let mapRenderable = (renderedFn, symFn, item: node) =>
- switch (item) {
- | `SymbolicDist(s) => Some(symFn(s))
- | `RenderedDist(r) => Some(renderedFn(r))
- | _ => None
- };
};
type simplificationResult = [
diff --git a/src/distPlus/expressionTree/MathJsParser.re b/src/distPlus/expressionTree/MathJsParser.re
index bfc2e42a..c8ee68b5 100644
--- a/src/distPlus/expressionTree/MathJsParser.re
+++ b/src/distPlus/expressionTree/MathJsParser.re
@@ -237,7 +237,8 @@ module MathAdtToDistDst = {
args: array(result(ExpressionTypes.ExpressionTree.node, string)),
) => {
let toOkAlgebraic = r => Ok(`AlgebraicCombination(r));
- let toOkTrunctate = r => Ok(`Truncate(r));
+ let toOkTruncate = r => Ok(`Truncate(r));
+ let toOkFloatFromDist = r => Ok(`FloatFromDist(r))
switch (name, args) {
| ("add", [|Ok(l), Ok(r)|]) => toOkAlgebraic((`Add, l, r))
| ("add", _) => Error("Addition needs two operands")
@@ -249,11 +250,11 @@ module MathAdtToDistDst = {
| ("divide", _) => Error("Division needs two operands")
| ("pow", _) => Error("Exponentiation is not yet supported.")
| ("leftTruncate", [|Ok(d), Ok(`SymbolicDist(`Float(lc)))|]) =>
- toOkTrunctate((Some(lc), None, d))
+ toOkTruncate((Some(lc), None, d))
| ("leftTruncate", _) =>
Error("leftTruncate needs two arguments: the expression and the cutoff")
| ("rightTruncate", [|Ok(d), Ok(`SymbolicDist(`Float(rc)))|]) =>
- toOkTrunctate((None, Some(rc), d))
+ toOkTruncate((None, Some(rc), d))
| ("rightTruncate", _) =>
Error(
"rightTruncate needs two arguments: the expression and the cutoff",
@@ -266,9 +267,19 @@ module MathAdtToDistDst = {
Ok(`SymbolicDist(`Float(rc))),
|],
) =>
- toOkTrunctate((Some(lc), Some(rc), d))
+ toOkTruncate((Some(lc), Some(rc), d))
| ("truncate", _) =>
Error("truncate needs three arguments: the expression and both cutoffs")
+ | ("pdf", [|Ok(d), Ok(`SymbolicDist(`Float(v)))|]) =>
+ toOkFloatFromDist((`Pdf(v), d))
+ | ("cdf", [|Ok(d), Ok(`SymbolicDist(`Float(v)))|]) =>
+ toOkFloatFromDist((`Cdf(v), d))
+ | ("inv", [|Ok(d), Ok(`SymbolicDist(`Float(v)))|]) =>
+ toOkFloatFromDist((`Inv(v), d))
+ | ("mean", [|Ok(d)|]) =>
+ toOkFloatFromDist((`Mean, d))
+ | ("sample", [|Ok(d)|]) =>
+ toOkFloatFromDist((`Sample, d))
| _ => Error("This type not currently supported")
};
};
@@ -316,11 +327,12 @@ module MathAdtToDistDst = {
| "pow"
| "leftTruncate"
| "rightTruncate"
- | "truncate" => operationParser(name, parseArgs())
- | "mean" as n
- | "inv" as n
- | "sample" as n
- | "pdf" as n
+ | "truncate"
+ | "mean"
+ | "inv"
+ | "sample"
+ | "cdf"
+ | "pdf" => operationParser(name, parseArgs())
| n => Error(n ++ "(...) is not currently supported")
};
};
diff --git a/src/distPlus/expressionTree/Operation.re b/src/distPlus/expressionTree/Operation.re
index 33c05461..26826f3f 100644
--- a/src/distPlus/expressionTree/Operation.re
+++ b/src/distPlus/expressionTree/Operation.re
@@ -41,6 +41,7 @@ module DistToFloat = {
let format = (operation, value) =>
switch (operation) {
+ | `Cdf(f) => {j|cdf(x=$f,$value)|j}
| `Pdf(f) => {j|pdf(x=$f,$value)|j}
| `Inv(f) => {j|inv(x=$f,$value)|j}
| `Sample => "sample($value)"
@@ -63,11 +64,17 @@ module Scale = {
| `Log => {j|verticalLog($value, $scaleBy) |j}
};
- let toKnownIntegralSumFn =
+ let toIntegralSumCacheFn =
fun
| `Multiply => ((a, b) => Some(a *. b))
| `Exponentiate => ((_, _) => None)
| `Log => ((_, _) => None);
+
+ let toIntegralCacheFn =
+ fun
+ | `Multiply => ((a, b) => None) // TODO: this could probably just be multiplied out (using Continuous.scaleBy)
+ | `Exponentiate => ((_, _) => None)
+ | `Log => ((_, _) => None);
};
module T = {
diff --git a/src/distPlus/expressionTree/SamplingDistribution.re b/src/distPlus/expressionTree/SamplingDistribution.re
new file mode 100644
index 00000000..60886b40
--- /dev/null
+++ b/src/distPlus/expressionTree/SamplingDistribution.re
@@ -0,0 +1,80 @@
+open ExpressionTypes.ExpressionTree;
+
+let isSamplingDistribution: node => bool =
+ fun
+ | `SymbolicDist(_) => true
+ | `RenderedDist(_) => true
+ | _ => false;
+
+let renderIfIsNotSamplingDistribution = (params, t) =>
+ !isSamplingDistribution(t)
+ ? switch (Render.render(params, t)) {
+ | Ok(r) => Ok(r)
+ | Error(e) => Error(e)
+ }
+ : Ok(t);
+
+let map = (~renderedDistFn, ~symbolicDistFn, node: node) =>
+ node
+ |> (
+ fun
+ | `RenderedDist(r) => Some(renderedDistFn(r))
+ | `SymbolicDist(s) => Some(symbolicDistFn(s))
+ | _ => None
+ );
+
+let sampleN = n =>
+ map(
+ ~renderedDistFn=Shape.sampleNRendered(n),
+ ~symbolicDistFn=SymbolicDist.T.sampleN(n),
+ );
+
+let getCombinationSamples = (n, algebraicOp, t1: node, t2: node) => {
+ switch (sampleN(n, t1), sampleN(n, t2)) {
+ | (Some(a), Some(b)) =>
+ Some(
+ Belt.Array.zip(a, b)
+ |> E.A.fmap(((a, b)) => Operation.Algebraic.toFn(algebraicOp, a, b)),
+ )
+ | _ => None
+ };
+};
+
+let combineShapesUsingSampling =
+ (evaluationParams: evaluationParams, algebraicOp, t1: node, t2: node) => {
+ let i1 = renderIfIsNotSamplingDistribution(evaluationParams, t1);
+ let i2 = renderIfIsNotSamplingDistribution(evaluationParams, t2);
+ E.R.merge(i1, i2)
+ |> E.R.bind(
+ _,
+ ((a, b)) => {
+ let samples =
+ getCombinationSamples(
+ evaluationParams.samplingInputs.sampleCount,
+ algebraicOp,
+ a,
+ b,
+ );
+
+ // todo: This bottom part should probably be somewhere else.
+ let shape =
+ samples
+ |> E.O.fmap(
+ Samples.T.fromSamples(
+ ~samplingInputs={
+ sampleCount:
+ Some(evaluationParams.samplingInputs.sampleCount),
+ outputXYPoints:
+ Some(evaluationParams.samplingInputs.outputXYPoints),
+ kernelWidth: evaluationParams.samplingInputs.kernelWidth,
+ },
+ ),
+ )
+ |> E.O.bind(_, (r: RenderTypes.ShapeRenderer.Sampling.outputs) =>
+ r.shape
+ )
+ |> E.O.toResult("No response");
+ shape |> E.R.fmap(r => `Normalize(`RenderedDist(r)));
+ },
+ );
+};
diff --git a/src/distPlus/renderers/samplesRenderer/Samples.re b/src/distPlus/renderers/samplesRenderer/Samples.re
index 0e857c2f..ff8f6f46 100644
--- a/src/distPlus/renderers/samplesRenderer/Samples.re
+++ b/src/distPlus/renderers/samplesRenderer/Samples.re
@@ -120,7 +120,7 @@ module T = {
|> E.FloatFloatMap.fmap(r => r /. length)
|> E.FloatFloatMap.toArray
|> XYShape.T.fromZippedArray
- |> Discrete.make(_, None);
+ |> Discrete.make;
let pdf =
continuousPart |> E.A.length > 5
@@ -150,7 +150,7 @@ module T = {
~outputXYPoints=samplingInputs.outputXYPoints,
formatUnitWidth(usedUnitWidth),
)
- |> Continuous.make(`Linear, _, None)
+ |> Continuous.make
|> (r => Some((r, foo)));
}
: None;
diff --git a/src/distPlus/symbolic/SymbolicDist.re b/src/distPlus/symbolic/SymbolicDist.re
index a14609e6..efa4e49b 100644
--- a/src/distPlus/symbolic/SymbolicDist.re
+++ b/src/distPlus/symbolic/SymbolicDist.re
@@ -3,6 +3,7 @@ open SymbolicTypes;
module Exponential = {
type t = exponential;
let pdf = (x, t: t) => Jstat.exponential##pdf(x, t.rate);
+ let cdf = (x, t: t) => Jstat.exponential##cdf(x, t.rate);
let inv = (p, t: t) => Jstat.exponential##inv(p, t.rate);
let sample = (t: t) => Jstat.exponential##sample(t.rate);
let mean = (t: t) => Ok(Jstat.exponential##mean(t.rate));
@@ -12,6 +13,7 @@ module Exponential = {
module Cauchy = {
type t = cauchy;
let pdf = (x, t: t) => Jstat.cauchy##pdf(x, t.local, t.scale);
+ let cdf = (x, t: t) => Jstat.cauchy##cdf(x, t.local, t.scale);
let inv = (p, t: t) => Jstat.cauchy##inv(p, t.local, t.scale);
let sample = (t: t) => Jstat.cauchy##sample(t.local, t.scale);
let mean = (_: t) => Error("Cauchy distributions have no mean value.");
@@ -21,6 +23,7 @@ module Cauchy = {
module Triangular = {
type t = triangular;
let pdf = (x, t: t) => Jstat.triangular##pdf(x, t.low, t.high, t.medium);
+ let cdf = (x, t: t) => Jstat.triangular##cdf(x, t.low, t.high, t.medium);
let inv = (p, t: t) => Jstat.triangular##inv(p, t.low, t.high, t.medium);
let sample = (t: t) => Jstat.triangular##sample(t.low, t.high, t.medium);
let mean = (t: t) => Ok(Jstat.triangular##mean(t.low, t.high, t.medium));
@@ -30,6 +33,7 @@ module Triangular = {
module Normal = {
type t = normal;
let pdf = (x, t: t) => Jstat.normal##pdf(x, t.mean, t.stdev);
+ let cdf = (x, t: t) => Jstat.normal##cdf(x, t.mean, t.stdev);
let from90PercentCI = (low, high) => {
let mean = E.A.Floats.mean([|low, high|]);
@@ -72,6 +76,7 @@ module Normal = {
module Beta = {
type t = beta;
let pdf = (x, t: t) => Jstat.beta##pdf(x, t.alpha, t.beta);
+ let cdf = (x, t: t) => Jstat.beta##cdf(x, t.alpha, t.beta);
let inv = (p, t: t) => Jstat.beta##inv(p, t.alpha, t.beta);
let sample = (t: t) => Jstat.beta##sample(t.alpha, t.beta);
let mean = (t: t) => Ok(Jstat.beta##mean(t.alpha, t.beta));
@@ -81,6 +86,7 @@ module Beta = {
module Lognormal = {
type t = lognormal;
let pdf = (x, t: t) => Jstat.lognormal##pdf(x, t.mu, t.sigma);
+ let cdf = (x, t: t) => Jstat.lognormal##cdf(x, t.mu, t.sigma);
let inv = (p, t: t) => Jstat.lognormal##inv(p, t.mu, t.sigma);
let mean = (t: t) => Ok(Jstat.lognormal##mean(t.mu, t.sigma));
let sample = (t: t) => Jstat.lognormal##sample(t.mu, t.sigma);
@@ -126,6 +132,7 @@ module Lognormal = {
module Uniform = {
type t = uniform;
let pdf = (x, t: t) => Jstat.uniform##pdf(x, t.low, t.high);
+ let cdf = (x, t: t) => Jstat.uniform##cdf(x, t.low, t.high);
let inv = (p, t: t) => Jstat.uniform##inv(p, t.low, t.high);
let sample = (t: t) => Jstat.uniform##sample(t.low, t.high);
let mean = (t: t) => Ok(Jstat.uniform##mean(t.low, t.high));
@@ -140,6 +147,7 @@ module Uniform = {
module Float = {
type t = float;
let pdf = (x, t: t) => x == t ? 1.0 : 0.0;
+ let cdf = (x, t: t) => x >= t ? 1.0 : 0.0;
let inv = (p, t: t) => p < t ? 0.0 : 1.0;
let mean = (t: t) => Ok(t);
let sample = (t: t) => t;
@@ -162,6 +170,18 @@ module T = {
| `Float(n) => Float.pdf(x, n)
};
+ let cdf = (x, dist) =>
+ switch (dist) {
+ | `Normal(n) => Normal.cdf(x, n)
+ | `Triangular(n) => Triangular.cdf(x, n)
+ | `Exponential(n) => Exponential.cdf(x, n)
+ | `Cauchy(n) => Cauchy.cdf(x, n)
+ | `Lognormal(n) => Lognormal.cdf(x, n)
+ | `Uniform(n) => Uniform.cdf(x, n)
+ | `Beta(n) => Beta.cdf(x, n)
+ | `Float(n) => Float.cdf(x, n)
+ };
+
let inv = (x, dist) =>
switch (dist) {
| `Normal(n) => Normal.inv(x, n)
@@ -244,6 +264,7 @@ module T = {
let operate = (distToFloatOp: ExpressionTypes.distToFloatOperation, s) =>
switch (distToFloatOp) {
+ | `Cdf(f) => Ok(cdf(f, s))
| `Pdf(f) => Ok(pdf(f, s))
| `Inv(f) => Ok(inv(f, s))
| `Sample => Ok(sample(s))
@@ -295,10 +316,14 @@ module T = {
let toShape = (sampleCount, d: symbolicDist): DistTypes.shape =>
switch (d) {
| `Float(v) =>
- Discrete(Discrete.make({xs: [|v|], ys: [|1.0|]}, Some(1.0)))
+ Discrete(
+ Discrete.make(~integralSumCache=Some(1.0), {xs: [|v|], ys: [|1.0|]}),
+ )
| _ =>
let xs = interpolateXs(~xSelection=`ByWeight, d, sampleCount);
let ys = xs |> E.A.fmap(x => pdf(x, d));
- Continuous(Continuous.make(`Linear, {xs, ys}, Some(1.0)));
+ Continuous(
+ Continuous.make(~integralSumCache=Some(1.0), {xs, ys}),
+ );
};
};