Kde with auto kernel width

__tests__/Bandwidth__Test.re

@@ -4,10 +4,10 @@ open Expect;
 describe("Bandwidth", () => {
   test("nrd0()", () => {
     let data = [|1., 4., 3., 2.|];
-    expect(Science.nrd0(data)) |> toEqual(0.7635139420854616);
+    expect(Bandwidth.nrd0(data)) |> toEqual(0.7625801874014622);
   });
   test("nrd()", () => {
     let data = [|1., 4., 3., 2.|];
-    expect(Science.nrd(data)) |> toEqual(0.899249754011766);
+    expect(Bandwidth.nrd(data)) |> toEqual(0.8981499984950554);
   });
 });

package.json

@@ -60,7 +60,6 @@
     "react-use": "^13.27.0",
     "reason-react": ">=0.7.0",
     "reschema": "1.3.0",
-    "science": "^1.9.3",
     "tailwindcss": "1.2.0"
   },
   "alias": {

src/Samples.re

@@ -89,11 +89,10 @@ module T = {
   };
 
   let kde = (~samples, ~outputXYPoints) => {
-    let width = Science.nrd0(samples);
+    let width = Bandwidth.nrd0(samples);
     let xyPointRange = E.A.Sorted.range(samples) |> E.O.default(0.0);
     let xyPointWidth = xyPointRange /. float_of_int(outputXYPoints);
     let kernelWidth = int_of_float(Jstat.max([|(width /. xyPointWidth), 1.0 |]));
-    Js.log4(samples, width, xyPointWidth, kernelWidth);
     KDE.normalSampling(samples, outputXYPoints, kernelWidth);
   };
 

src/symbolic/Jstat.re

@@ -97,6 +97,8 @@ type binomial = {
 [@bs.module "jstat"] external deviation: array(float) => float = "deviation";
 [@bs.module "jstat"] external stdev: array(float) => float = "stdev";
 [@bs.module "jstat"]
-external quantiles: (array(float), array(float)) => float = "quantiles";
+external quartiles: (array(float)) => array(float) = "quartiles";
 [@bs.module "jstat"]
-external percentile: (array(float), array(float)) => float = "percentile";
+external quantiles: (array(float), array(float)) => array(float) = "quantiles";
+[@bs.module "jstat"]
+external percentile: (array(float), float, bool) => float = "percentile";

src/utility/Bandwidth.re

@@ -0,0 +1,30 @@
+//The math here was taken from https://github.com/jasondavies/science.js/blob/master/src/stats/bandwidth.js
+
+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);
+};
+
+// 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 e = Js_math.abs_float(x[1]);
+  let lo' =
+    switch (lo, hi, e) {
+    | (lo, _, _) when !Js.Float.isNaN(lo) => lo
+    | (_, hi, _) when !Js.Float.isNaN(hi) => hi
+    | (_, _, e) when !Js.Float.isNaN(e) => e
+    | _ => 1.0
+    };
+  0.9 *. lo' *. Js.Math.pow_float(~base=len(x), ~exp=-0.2);
+};
+
+// Scott, D. W. (1992) Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley.
+let nrd = x => {
+  let h = iqr(x) /. 1.34;
+  1.06
+  *. Js.Math.min_float(Js.Math.sqrt(Jstat.variance(x)), h)
+  *. Js.Math.pow_float(~base=len(x), ~exp=(-1.0) /. 5.0);
+};

src/utility/Science.re

@@ -1,28 +0,0 @@
-[@bs.val] [@bs.module "science"] [@bs.scope "stats"]
-external variance: array(float) => float = "variance";
-
-[@bs.val] [@bs.module "science"] [@bs.scope "stats"]
-external iqr: array(float) => float = "iqr";
-
-let len = x => E.A.length(x) |> float_of_int;
-
-let nrd0 = x => {
-  let hi = Js_math.sqrt(variance(x));
-  let lo = Js_math.minMany_float([|hi, iqr(x) /. 1.34|]);
-  let e = Js_math.abs_float(x[1]);
-  let lo' =
-    switch (lo, hi, e) {
-    | (lo, hi, e) when !Js.Float.isNaN(lo) => lo
-    | (lo, hi, e) when !Js.Float.isNaN(hi) => hi
-    | (lo, hi, e) when !Js.Float.isNaN(e) => e
-    | _ => 1.0
-    };
-  0.9 *. lo' *. Js_math.pow_float(len(x), -0.2);
-};
-
-let nrd = x => {
-  let h = iqr(x) /. 1.34;
-  1.06
-  *. Js_math.min_float(Js_math.sqrt(variance(x)), h)
-  *. Js_math.pow_float(len(x), (-1.0) /. 5.0);
-};