Add more functions

Author: rigtorp

README.md

@@ -30,19 +30,23 @@ 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
+* `min(a)` - [Minimum](http://en.wikipedia.org/wiki/Minimum)
+* `max(a)` - [Maximum](http://en.wikipedia.org/wiki/Maximum)
+* `mean(a)` - [Mean](http://en.wikipedia.org/wiki/Mean)
+* `gmean(a)` - [Geometric mean](http://en.wikipedia.org/wiki/Mean)
+* `hmean(a)` - [Harmonic mean](http://en.wikipedia.org/wiki/Mean)
+* `var(a)` - [Variance](http://en.wikipedia.org/wiki/Variance)
+* `std(a)` - [Standard deviation](http://en.wikipedia.org/wiki/Standard_deviation)
+* `skew(a)` - [Skewness](http://en.wikipedia.org/wiki/Skewness)
+* `kurt(a)` - [Kurtosis](http://en.wikipedia.org/wiki/Kurtosis)
+* `corr(x, y)` - [Correlation](http://en.wikipedia.org/wiki/Correlation) between x and y
+* `entropy(p)` - [Entropy](http://en.wikipedia.org/wiki/Entropy_(information_theory))
+* `kldiv(p, q)` - [Kullback–Leibler divergence](http://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence)
+* `shuffle(a)` - [Shuffle](http://en.wikipedia.org/wiki/Random_permutation) using the [Fisher–Yates shuffle](http://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle)
+* `sample(a)` - [Sample](http://en.wikipedia.org/wiki/Sampling_(statistics)) with replacement
+* `boot(nboot, bootfun, data...)` - [Bootstrap](http://en.wikipedia.org/wiki/Bootstrapping_(statistics)) the bootfun statistic
+* `bootci(nboot, bootfun, data...)` - Calculate bootstrap confidence intervals using the normal model
+* `randn()` - Draw random sample from the [standard normal distribution](http://en.wikipedia.org/wiki/Normal_distribution) using the [Marsaglia polar method](http://en.wikipedia.org/wiki/Marsaglia_polar_method)
 
 Credits
 =======

examples/bootstrap.js

@@ -0,0 +1,11 @@
+var sk = require("../statkit.js");
+
+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);

statkit.js

@@ -119,6 +119,38 @@ exports.corr = function(x, y) {
  return sxy / Math.sqrt(sx2 * sy2);
 };
 
+exports.entropy = function(p) {
+ var n = p.length;
+ var e = 0.0;
+ for (var i = 0; i < n; ++i) {
+ e -= p[i] * Math.log(p[i]);
+ }
+ return e;
+}
+
+exports.kldiv = function(p, q) {
+ var n = p.length;
+ var e = 0.0;
+ var ce = 0.0;
+ for (var i = 0; i < n; ++i) {
+ e -= p[i] * Math.log(p[i]);
+ ce -= p[i] * Math.log(q[i]);
+ }
+ return ce - e;
+}
+
+// http://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle
+exports.shuffle = function(a) {
+ var n = a.length;
+ for (var i = n - 1; i > 0; i--) {
+ var j = Math.random() * i | 0; // 0 ≤ j < i
+ var t = array[j];
+ array[j] = array[i];
+ array[i] = t;
+ }
+ return a;
+}
+
 exports.sample = function(a) {
  var n = a.length;
  var s = a.slice(0);
@@ -162,6 +194,19 @@ exports.bootci = function(nboot, bootfun) {
  return [v - 2*s, v + 2*s];
 };
 
+// http://en.wikipedia.org/wiki/Marsaglia_polar_method
+// TODO: implement http://en.wikipedia.org/wiki/Ziggurat_algorithm
+exports.randn = function() {
+ var u, v, s;
+ do {
+ u = Math.rand() * 2 - 1;
+ v = Math.rand() * 2 - 1;
+ s = u * u + v * v;
+ } while (s >= 1 || s == 0);
+ s = Math.sqrt(-2 * Math.log(s) / s);
+ return u * s;
+};
+
 // http://picomath.org/javascript/erf.js.html
 exports.erf = function(x) {
  // constants

test/test.py

@@ -18,6 +18,9 @@
 print "kurtosis", np.array([st.kurtosis(a)])
 print "skew", np.array([st.skew(a)])
 
+p = np.ones(5) / 5;
+print "entropy", np.array([st.entropy(p)])
+
 
 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]