Merge pull request #440 from quantified-uncertainty/fromsamples

`fromSamples`

packages/squiggle-lang/__tests__/TS/SampleSet_test.ts

@@ -1,5 +1,5 @@
 import { Distribution } from "../../src/js/index";
-import { expectErrorToBeBounded, failDefault } from "./TestHelpers";
+import { expectErrorToBeBounded, failDefault, testRun } from "./TestHelpers";
 import * as fc from "fast-check";
 
 // Beware: float64Array makes it appear in an infinite loop.
@@ -212,3 +212,18 @@ describe("mean is mean", () => {
     );
   });
 });
+
+describe("fromSamples function", () => {
+  test.skip("gives a mean near the mean of the input", () => {
+    fc.assert(
+      fc.property(arrayGen(), (xs_) => {
+        let xs = Array.from(xs_);
+        let xsString = xs.toString();
+        let squiggleString = `x = fromSamples([${xsString}]); mean(x)`;
+        let squiggleResult = testRun(squiggleString);
+        let mean = xs.reduce((a, b) => a + b, 0.0) / xs.length;
+        expect(squiggleResult.value).toBeCloseTo(mean, 4);
+      })
+    );
+  });
+});

packages/squiggle-lang/src/rescript/Distributions/DistributionOperation/DistributionOperation.res

@@ -189,6 +189,11 @@ let rec run = (~env, functionCallInfo: functionCallInfo): outputType => {
     ->GenericDist.mixture(~scaleMultiplyFn=scaleMultiply, ~pointwiseAddFn=pointwiseAdd)
     ->E.R2.fmap(r => Dist(r))
     ->OutputLocal.fromResult
+  | FromSamples(xs) => xs
+    ->SampleSetDist.make
+    ->E.R2.errMap(x => DistributionTypes.SampleSetError(x))
+    ->E.R2.fmap(x => x->DistributionTypes.SampleSet->Dist)
+    ->OutputLocal.fromResult
   }
 }
 
@@ -229,6 +234,7 @@ module Constructors = {
   let isNormalized = (~env, dist) => C.isNormalized(dist)->run(~env)->toBoolR
   let toPointSet = (~env, dist) => C.toPointSet(dist)->run(~env)->toDistR
   let toSampleSet = (~env, dist, n) => C.toSampleSet(dist, n)->run(~env)->toDistR
+  let fromSamples = (~env, xs) => C.fromSamples(xs)->run(~env)->toDistR
   let truncate = (~env, dist, leftCutoff, rightCutoff) =>
     C.truncate(dist, leftCutoff, rightCutoff)->run(~env)->toDistR
   let inspect = (~env, dist) => C.inspect(dist)->run(~env)->toDistR

packages/squiggle-lang/src/rescript/Distributions/DistributionOperation/DistributionOperation.resi

@@ -61,6 +61,8 @@ module Constructors: {
   @genType
   let toSampleSet: (~env: env, genericDist, int) => result<genericDist, error>
   @genType
+  let fromSamples: (~env: env, SampleSetDist.t) => result<genericDist, error>
+  @genType
   let truncate: (~env: env, genericDist, option<float>, option<float>) => result<genericDist, error>
   @genType
   let inspect: (~env: env, genericDist) => result<genericDist, error>

packages/squiggle-lang/src/rescript/Distributions/DistributionTypes.res

@@ -11,7 +11,7 @@ type error =
   | NotYetImplemented
   | Unreachable
   | DistributionVerticalShiftIsInvalid
-  | TooFewSamples
+  | SampleSetError(SampleSetDist.sampleSetError)
   | ArgumentError(string)
   | OperationError(Operation.Error.t)
   | PointSetConversionError(SampleSetDist.pointsetConversionError)
@@ -35,7 +35,8 @@ module Error = {
     | DistributionVerticalShiftIsInvalid => "Distribution Vertical Shift is Invalid"
     | ArgumentError(s) => `Argument Error ${s}`
     | LogarithmOfDistributionError(s) => `Logarithm of input error: ${s}`
-    | TooFewSamples => "Too Few Samples"
+    | SampleSetError(TooFewSamples) => "Too Few Samples"
+    | SampleSetError(NonNumericInput(err)) => `Found a non-number in input: ${err}`
     | OperationError(err) => Operation.Error.toString(err)
     | PointSetConversionError(err) => SampleSetDist.pointsetConversionErrorToString(err)
     | SparklineError(err) => PointSetTypes.sparklineErrorToString(err)
@@ -47,10 +48,7 @@ module Error = {
   let resultStringToResultError: result<'a, string> => result<'a, error> = n =>
     n->E.R2.errMap(r => r->fromString)
 
-  let sampleErrorToDistErr = (err: SampleSetDist.sampleSetError): error =>
-    switch err {
-    | TooFewSamples => TooFewSamples
-    }
+  let sampleErrorToDistErr = (err: SampleSetDist.sampleSetError): error => SampleSetError(err)
 }
 
 @genType
@@ -99,6 +97,7 @@ module DistributionOperation = {
   type genericFunctionCallInfo =
     | FromDist(fromDist, genericDist)
     | FromFloat(fromDist, float)
+    | FromSamples(array<float>)
     | Mixture(array<(genericDist, float)>)
 
   let distCallToString = (distFunction: fromDist): string =>
@@ -124,6 +123,7 @@ module DistributionOperation = {
     switch d {
     | FromDist(f, _) | FromFloat(f, _) => distCallToString(f)
     | Mixture(_) => `mixture`
+    | FromSamples(_) => `fromSamples`
     }
 }
 module Constructors = {
@@ -140,6 +140,7 @@ module Constructors = {
     let isNormalized = (dist): t => FromDist(ToBool(IsNormalized), dist)
     let toPointSet = (dist): t => FromDist(ToDist(ToPointSet), dist)
     let toSampleSet = (dist, r): t => FromDist(ToDist(ToSampleSet(r)), dist)
+    let fromSamples = (xs): t => FromSamples(xs)
     let truncate = (dist, left, right): t => FromDist(ToDist(Truncate(left, right)), dist)
     let inspect = (dist): t => FromDist(ToDist(Inspect), dist)
     let toString = (dist): t => FromDist(ToString(ToString), dist)

packages/squiggle-lang/src/rescript/Distributions/SampleSetDist/SampleSetDist.res

@@ -1,11 +1,12 @@
 @genType
 module Error = {
   @genType
-  type sampleSetError = TooFewSamples
+  type sampleSetError = TooFewSamples | NonNumericInput(string)
 
   let sampleSetErrorToString = (err: sampleSetError): string =>
     switch err {
     | TooFewSamples => "Too few samples when constructing sample set"
+    | NonNumericInput(err) => `Found a non-number in input: ${err}`
     }
 
   @genType

packages/squiggle-lang/src/rescript/Distributions/SampleSetDist/SampleSetDist_Bandwidth.res

@@ -1,27 +1,30 @@
 //The math here was taken from https://github.com/jasondavies/science.js/blob/master/src/stats/SampleSetDist_Bandwidth.js
-
+let {iqr_percentile, nrd0_lo_denominator, one, nrd0_coef, nrd_coef, nrd_fractionalPower} = module(
+  MagicNumbers.SampleSetBandwidth
+)
 let len = x => E.A.length(x) |> float_of_int
 
-let iqr = x => Jstat.percentile(x, 0.75, true) -. Jstat.percentile(x, 0.25, true)
+let iqr = x =>
+  Jstat.percentile(x, iqr_percentile, true) -. Jstat.percentile(x, 1.0 -. iqr_percentile, true)
 
 // Silverman, B. W. (1986) Density Estimation. London: Chapman and Hall.
 let nrd0 = x => {
   let hi = Js_math.sqrt(Jstat.variance(x))
-  let lo = Js_math.minMany_float([hi, iqr(x) /. 1.34])
+  let lo = Js_math.minMany_float([hi, iqr(x) /. nrd0_lo_denominator])
   let e = Js_math.abs_float(x[1])
   let lo' = switch (lo, hi, e) {
   | (lo, _, _) if !Js.Float.isNaN(lo) => lo
   | (_, hi, _) if !Js.Float.isNaN(hi) => hi
   | (_, _, e) if !Js.Float.isNaN(e) => e
-  | _ => 1.0
+  | _ => one
   }
-  0.9 *. lo' *. Js.Math.pow_float(~base=len(x), ~exp=-0.2)
+  nrd0_coef *. lo' *. Js.Math.pow_float(~base=len(x), ~exp=nrd_fractionalPower)
 }
 
 // Scott, D. W. (1992) Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley.
 let nrd = x => {
-  let h = iqr(x) /. 1.34
-  1.06 *.
+  let h = iqr(x) /. nrd0_lo_denominator
+  nrd_coef *.
   Js.Math.min_float(Js.Math.sqrt(Jstat.variance(x)), h) *.
-  Js.Math.pow_float(~base=len(x), ~exp=-1.0 /. 5.0)
+  Js.Math.pow_float(~base=len(x), ~exp=nrd_fractionalPower)
 }

packages/squiggle-lang/src/rescript/MagicNumbers.res

@@ -35,3 +35,16 @@ module ToPointSet = {
  */
   let minDiscreteToKeep = samples => max(20, E.A.length(samples) / 50)
 }
+
+module SampleSetBandwidth = {
+  // Silverman, B. W. (1986) Density Estimation. London: Chapman and Hall.
+  // Scott, D. W. (1992) Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley.
+  let iqr_percentile = 0.75
+  let iqr_percentile_complement = 1.0 -. iqr_percentile
+  let nrd0_lo_denominator = 1.34
+  let one = 1.0
+  let nrd0_coef = 0.9
+
+  let nrd_coef = 1.06
+  let nrd_fractionalPower = -0.2
+}

packages/squiggle-lang/src/rescript/ReducerInterface/ReducerInterface_GenericDistribution.res

@@ -218,6 +218,14 @@ let dispatchToGenericOutput = (call: ExpressionValue.functionCall): option<
     Helpers.toDistFn(ToSampleSet(Belt.Int.fromFloat(float)), dist)
   | ("toSampleSet", [EvDistribution(dist)]) =>
     Helpers.toDistFn(ToSampleSet(MagicNumbers.Environment.defaultSampleCount), dist)
+  | ("fromSamples", [EvArray(inputArray)]) => {
+      let _wrapInputErrors = x => SampleSetDist.NonNumericInput(x)
+      let parsedArray = Helpers.parseNumberArray(inputArray)->E.R2.errMap(_wrapInputErrors)
+      switch parsedArray {
+      | Ok(array) => runGenericOperation(FromSamples(array))
+      | Error(e) => GenDistError(SampleSetError(e))
+      }->Some
+    }
   | ("inspect", [EvDistribution(dist)]) => Helpers.toDistFn(Inspect, dist)
   | ("truncateLeft", [EvDistribution(dist), EvNumber(float)]) =>
     Helpers.toDistFn(Truncate(Some(float), None), dist)

packages/squiggle-lang/src/rescript/Utility/E.res

@@ -289,6 +289,13 @@ module R = {
     | Ok(r) => r->Ok
     | Error(x) => x->f->Error
     }
+
+  //I'm not sure what to call this.
+  let unify = (a: result<'a, 'b>, c: 'b => 'a): 'a =>
+    switch a {
+    | Ok(x) => x
+    | Error(x) => c(x)
+    }
 }
 
 module R2 = {
@@ -307,6 +314,8 @@ module R2 = {
     | Ok(x) => x->Ok
     | Error(x) => x->f->Error
     }
+
+  let toExn = (a, b) => R.toExn(b, a)
 }
 
 let safe_fn_of_string = (fn, s: string): option<'a> =>

packages/website/docs/Features/Functions.mdx

@@ -98,6 +98,19 @@ bound `a`, mode `b` and upper bound `c`.
 
 Squiggle, when the context is right, automatically casts a float to a constant distribution.
 
+## `fromSamples`
+
+The last distribution constructor takes an array of samples and constructs a sample set distribution.
+
+<SquiggleEditor initialSquiggleString="fromSamples([1,2,3,4,6,5,5,5])" />
+
+#### Validity
+
+For `fromSamples(xs)`,
+
+- `xs.length > 5`
+- Strictly every element of `xs` must be a number.
+
 ## Operating on distributions
 
 Here are the ways we combine distributions.
@@ -315,6 +328,16 @@ Or `PointSet` format
 
 <SquiggleEditor initialSquiggleString="toPointSet(normal(5, 10))" />
 
+### `toSampleSet` has two signatures
+
+Above, we saw the unary `toSampleSet`, which uses an internal hardcoded number of samples. If you'd like to provide the number of samples, it has a binary signature as well (floored)
+
+<SquiggleEditor initialSquiggleString="toSampleSet(0.1 to 1, 100.1)" />
+
+#### Validity
+
+- Second argument to `toSampleSet` must be a number.
+
 ## Normalization
 
 Some distribution operations (like horizontal shift) return an unnormalized distriibution.
@@ -333,18 +356,6 @@ We provide a predicate `isNormalized`, for when we have simple control flow
 
 - Input to `isNormalized` must be a dist
 
-## Convert any distribution to a sample set distribution
-
-`toSampleSet` has two signatures
-
-It is unary when you use an internal hardcoded number of samples
-
-<SquiggleEditor initialSquiggleString="toSampleSet(0.1 to 1)" />
-
-And binary when you provide a number of samples (floored)
-
-<SquiggleEditor initialSquiggleString="toSampleSet(0.1 to 1, 100)" />
-
 ## `inspect`
 
 You may like to debug by right clicking your browser and using the _inspect_ functionality on the webpage, and viewing the _console_ tab. Then, wrap your squiggle output with `inspect` to log an internal representation.