From 7ab81249bf95d350916f1b5a26135ef83da05003 Mon Sep 17 00:00:00 2001
From: Ozzie Gooen <ozzieagooen@gmail.com>
Date: Fri, 3 Apr 2020 20:41:37 +0100
Subject: [PATCH] Kde with auto kernel width

---
 __tests__/Bandwidth__Test.re |  4 ++--
 package.json                 |  1 -
 src/Samples.re               |  3 +--
 src/symbolic/Jstat.re        |  6 ++++--
 src/utility/Bandwidth.re     | 30 ++++++++++++++++++++++++++++++
 src/utility/Science.re       | 28 ----------------------------
 6 files changed, 37 insertions(+), 35 deletions(-)
 create mode 100644 src/utility/Bandwidth.re
 delete mode 100644 src/utility/Science.re

diff --git a/__tests__/Bandwidth__Test.re b/__tests__/Bandwidth__Test.re
index 013402fa..54cd111d 100644
--- a/__tests__/Bandwidth__Test.re
+++ b/__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);
   });
 });
\ No newline at end of file
diff --git a/package.json b/package.json
index d483b362..dda61fae 100644
--- a/package.json
+++ b/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": {
diff --git a/src/Samples.re b/src/Samples.re
index d7f3b5d2..360f8279 100644
--- a/src/Samples.re
+++ b/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);
   };
 
diff --git a/src/symbolic/Jstat.re b/src/symbolic/Jstat.re
index 4c361bd2..5f1c6c51 100644
--- a/src/symbolic/Jstat.re
+++ b/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";
\ No newline at end of file
+external quantiles: (array(float), array(float)) => array(float) = "quantiles";
+[@bs.module "jstat"]
+external percentile: (array(float), float, bool) => float = "percentile";
\ No newline at end of file
diff --git a/src/utility/Bandwidth.re b/src/utility/Bandwidth.re
new file mode 100644
index 00000000..d2315639
--- /dev/null
+++ b/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);
+};
\ No newline at end of file
diff --git a/src/utility/Science.re b/src/utility/Science.re
deleted file mode 100644
index b755ad5b..00000000
--- a/src/utility/Science.re
+++ /dev/null
@@ -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);
-};
\ No newline at end of file