Initial commit

Author: rigtorp

.gitignore

@@ -0,0 +1 @@
+node_modules/

LICENSE

@@ -0,0 +1,22 @@
+
+The MIT License (MIT)
+
+Copyright (c) 2014 Erik Rigtorp <erik@rigtorp.se>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.

README.md

@@ -0,0 +1,49 @@
+statkit
+=======
+
+A statistics toolkit for javascript.
+
+Usage
+=====
+
+Install using [npm](https://npmjs.org):
+```
+npm install statkit
+```
+
+Calculate a confidence interval for a correlation using the bootstrap method:
+
+```javascript
+var sk = require("statkit");
+
+var lsat = [576, 635, 558, 578, 666, 580, 555,
+ 661, 651, 605, 653, 575, 545, 572, 594];
+var gpa = [3.39, 3.30, 2.81, 3.03, 3.44, 3.07, 3.00,
+ 3.43, 3.36, 3.13, 3.12, 2.74, 2.76, 2.88, 2.96];
+
+var corr = sk.corr(gpa, lsat);
+var ci = sk.bootci(100000, sk.corr, gpa, lsat);
+
+console.log("corr = ", corr, "ci = ", ci);
+```
+
+Functions
+=========
+
+* min(a) - Minimum
+* max(a) - Maximum
+* mean(a) - Mean
+* gmean(a) - Geometric mean
+* hmean(a) - Harmonic mean
+* var(a) - Variance
+* std(a) - Standard deviation
+* skew(a) - Skewness
+* kurt(a) - Kurtosis
+* corr(x, y) - Correlation between x and y
+* sample(a) - Sample with replacement
+* boot(nboot, bootfun, data...) - Bootstrap the bootfun statistic
+* bootci(nboot, bootfun, data...) - Calculate bootstrap confidence intervals using the normal model
+
+Credits
+=======
+(c) 2014 Erik Rigtorp <erik@rigtorp.se>. MIT License

package.json

@@ -0,0 +1,33 @@
+{
+ "name": "statkit",
+ "version": "0.1.0",
+ "description": "Statistics toolkit",
+ "main": "statkit.js",
+ "directories": {
+ "test": "test"
+ },
+ "devDependencies": {
+ "tape": "~2.13.3"
+ },
+ "scripts": {
+ "test": "tape test/*.js"
+ },
+ "repository": {
+ "type": "git",
+ "url": "git://github.com/rigtorp/statkit.git"
+ },
+ "homepage": "https://github.com/rigtorp/statkit",
+ "keywords": [
+ "statistics",
+ "resampling",
+ "bootstrap",
+ "analytics"
+ ],
+ "author": {
+ "name": "Erik Rigtorp",
+ "email": "erik@rigtorp.se",
+ "url": "http://rigtorp.se"
+ },
+ "license": "MIT",
+ "readmeFilename": "README.md"
+}

statkit.js

@@ -0,0 +1,192 @@
+"use strict"
+
+exports.min = function(a) {
+ var n = a.length;
+ var s = a[0];
+ for (var i = 1; i < n; ++i) {
+ if (a[i] < s) {
+ s = a[i];
+ }
+ }
+ return s;
+};
+
+exports.max = function(a) {
+ var n = a.length;
+ var s = a[0];
+ for (var i = 1; i < n; ++i) {
+ if (a[i] > s) {
+ s = a[i];
+ }
+ }
+ return s;
+};
+
+exports.mean = function(a) {
+ var n = a.length;
+ var s = 0.0;
+ for (var i = 0; i < n; ++i) {
+ s += a[i];
+ }
+ return s / n;
+};
+
+exports.gmean = function(a) {
+ var n = a.length;
+ var s = 0.0;
+ for (var i = 0; i < n; ++i) {
+ s += Math.log(a[i]);
+ }
+ return Math.exp(s / n);
+};
+
+exports.hmean = function(a) {
+ var n = a.length;
+ var s = 0.0;
+ for (var i = 0; i < n; ++i) {
+ s += 1.0 / a[i];
+ }
+ return n / s;
+};
+
+exports.var = function(a, m) {
+ var n = a.length;
+ var s = 0.0;
+ if (m === void 0) {
+ m = exports.mean(a);
+ }
+ for (var i = 0; i < n; ++i) {
+ var z = a[i] - m;
+ s += z * z;
+ }
+ return s / n;
+};
+
+exports.std = function(a, m) {
+ return Math.sqrt(exports.var(a, m));
+};
+
+exports.skew = function(a, m) {
+ var n = a.length;
+ var cm2 = 0.0;
+ var cm3 = 0.0;
+ if (m === void 0) {
+ m = exports.mean(a);
+ }
+ for (var i = 0; i < n; ++i) {
+ var z = a[i] - m;
+ var z2 = z * z;
+ cm2 += z2;
+ cm3 += z2 * z;
+ }
+ cm2 /= n;
+ cm3 /= n;
+ return cm3 / Math.sqrt(cm2 * cm2 * cm2);
+};
+
+exports.kurt = function(a, m) {
+ var n = a.length;
+ var cm2 = 0.0;
+ var cm4 = 0.0;
+ if (m === void 0) {
+ m = exports.mean(a);
+ }
+ for (var i = 0; i < n; ++i) {
+ var z = a[i] - m;
+ var z2 = z * z;
+ cm2 += z2;
+ cm4 += z2 * z2;
+ }
+ cm2 /= n;
+ cm4 /= n;
+ return (cm4 / (cm2 * cm2)) - 3;
+};
+
+exports.corr = function(x, y) {
+ var n = x.length;
+ var xm = exports.mean(x);
+ var ym = exports.mean(y);
+ var sxy = 0.0;
+ var sx2 = 0.0;
+ var sy2 = 0.0;
+ for (var i = 0; i < n; ++i) {
+ var xz = x[i] - xm;
+ var yz = y[i] - ym;
+ sxy += xz * yz;
+ sx2 += xz * xz;
+ sy2 += yz * yz;
+ }
+ return sxy / Math.sqrt(sx2 * sy2);
+};
+
+exports.sample = function(a) {
+ var n = a.length;
+ var s = a.slice(0);
+ for (var i = 0; i < n; ++i) {
+ s[i] = a[Math.floor(n * Math.random())];
+ }
+ return s;
+}
+
+exports.boot = function(nboot, bootfun) {
+ var data = [];
+ for (var i = 2; i < arguments.length; ++i) {
+ data[i - 2] = arguments[i];
+ }
+ var sample = [];
+ for (var i = 0; i < data.length; ++i) {
+ sample[i] = data[i].slice(0);
+ }
+ var n = data[0].length;
+ var res = Array(nboot);
+ for (var i = 0; i < nboot; ++i) {
+ for (var j = 0; j < n; ++j) {
+ var idx = Math.floor(n * Math.random());
+ for (var k = 0; k < sample.length; ++k) {
+ sample[k][j] = data[k][idx];
+ }
+ }
+ res[i] = bootfun.apply(null, sample);
+ }
+ return res;
+};
+
+exports.bootci = function(nboot, bootfun) {
+ var data = [];
+ for (var i = 2; i < arguments.length; ++i) {
+ data[i - 2] = arguments[i];
+ }
+ var v = bootfun.apply(null, data);
+ var bootstat = exports.boot.apply(null, [nboot, bootfun].concat(data));
+ var s = exports.std(bootstat);
+ return [v - 2*s, v + 2*s];
+};
+
+// http://picomath.org/javascript/erf.js.html
+exports.erf = function(x) {
+ // constants
+ var a1 = 0.254829592;
+ var a2 = -0.284496736;
+ var a3 = 1.421413741;
+ var a4 = -1.453152027;
+ var a5 = 1.061405429;
+ var p = 0.3275911;
+
+ // Save the sign of x
+ var sign = 1;
+ if (x < 0) {
+ sign = -1;
+ }
+ x = Math.abs(x);
+
+ // A&S formula 7.1.26
+ var t = 1.0/(1.0 + p*x);
+ var y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*Math.exp(-x*x);
+
+ return sign*y;
+};
+
+
+exports.normcdf = function(x) {
+ return 0.5 * (1.0 + exports.erf(x / Math.sqrt(2)));
+}

test/bootstrap.js

@@ -0,0 +1,16 @@
+var st = require("../statkit.js")
+ , tape = require("tape")
+
+tape("bootstrap", function(t) {
+
+ var lsat = [576, 635, 558, 578, 666, 580, 555, 661, 651, 605, 653, 575, 545, 572, 594];
+ var gpa = [3.39, 3.30, 2.81, 3.03, 3.44, 3.07, 3.00, 3.43, 3.36, 3.13, 3.12, 2.74, 2.76, 2.88, 2.96];
+
+ t.ok(Math.abs(st.corr(gpa, lsat) - 0.776) < 1e-3, "corr");
+
+ var ci = st.bootci(100000, st.corr, gpa, lsat);
+ t.ok(Math.abs(ci[0] - 0.509) < 1e-2, "lci");
+ t.ok(Math.abs(ci[1] - 1.043) < 1e-2, "uci");
+
+ t.end();
+})

test/test.js

@@ -0,0 +1,28 @@
+var st = require("../statkit.js")
+ , tape = require("tape")
+
+tape("stats", function(t) {
+
+ var a = [0.7004371218578347, 0.8441866428504345, 0.6765143359376254,
+ 0.7278580572480748, 0.9514579574463393, 0.012703197034767633,
+ 0.41358769878652346, 0.04881279380000003, 0.09992856132121142,
+ 0.5080663057670065];
+
+ var tol = 1e-16;
+
+ function near(a, b, msg) {
+ t.ok(Math.abs(a - b) < tol, msg)
+ }
+
+ t.equal(st.min(a), 0.01270319703476763262, "min")
+ t.equal(st.max(a), 0.95145795744633931967, "max")
+ near(st.mean(a), 0.49835526720498179998, "mean")
+ near(st.gmean(a), 0.28557127946531724039, "gmean")
+ near(st.hmean(a), 0.08325329480982600339, "hmean")
+ near(st.var(a), 0.10541207022165724472, "var")
+ near(st.std(a), 0.32467225046446029912, "std")
+ near(st.kurt(a), -1.34748749550481905501, "kurt")
+ near(st.skew(a), -0.31810224917209384277, "skew")
+
+ t.end()
+})

test/test.py

@@ -0,0 +1,27 @@
+import numpy as np
+import scipy.stats as st
+
+np.set_printoptions(precision=20)
+
+np.random.seed(666)
+a = np.random.rand(10)
+
+print a.tolist()
+
+print "min", np.array([np.min(a)])
+print "max", np.array([np.max(a)])
+print "mean", np.array([np.mean(a)])
+print "gmean", np.array([st.gmean(a)])
+print "hmean", np.array([st.hmean(a)])
+print "var", np.array([np.var(a)])
+print "std", np.array([np.std(a)])
+print "kurtosis", np.array([st.kurtosis(a)])
+print "skew", np.array([st.skew(a)])
+
+
+p = [0.0000001, 0.00001, 0.001, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65,
+ 0.75, 0.85, 0.95, 0.999, 0.99999, 0.9999999]
+
+x = st.norm.ppf(p)
+
+print x