Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Update science.js to 1.7.0: fix science.stats.kde.

  • Loading branch information...
commit 3c0cdb2285e8ce82df2f0bc6955800d423ac4bbb 1 parent f90bc4f
@jasondavies jasondavies authored
View
222 lib/science/science.js
@@ -1,5 +1,225 @@
-(function(){science = {version: "0.0.1"}; // semver
+(function(){science = {version: "1.7.0"}; // semver
+science.ascending = function(a, b) {
+ return a - b;
+};
+// Euler's constant.
+science.EULER = .5772156649015329;
+// Compute exp(x) - 1 accurately for small x.
+science.expm1 = function(x) {
+ return (x < 1e-5 && x > -1e-5) ? x + .5 * x * x : Math.exp(x) - 1;
+};
science.functor = function(v) {
return typeof v === "function" ? v : function() { return v; };
};
+// Based on:
+// http://www.johndcook.com/blog/2010/06/02/whats-so-hard-about-finding-a-hypotenuse/
+science.hypot = function(x, y) {
+ x = Math.abs(x);
+ y = Math.abs(y);
+ var max,
+ min;
+ if (x > y) { max = x; min = y; }
+ else { max = y; min = x; }
+ var r = min / max;
+ return max * Math.sqrt(1 + r * r);
+};
+science.quadratic = function() {
+ var complex = false;
+
+ function quadratic(a, b, c) {
+ var d = b * b - 4 * a * c;
+ if (d > 0) {
+ d = Math.sqrt(d) / (2 * a);
+ return complex
+ ? [{r: -b - d, i: 0}, {r: -b + d, i: 0}]
+ : [-b - d, -b + d];
+ } else if (d === 0) {
+ d = -b / (2 * a);
+ return complex ? [{r: d, i: 0}] : [d];
+ } else {
+ if (complex) {
+ d = Math.sqrt(-d) / (2 * a);
+ return [
+ {r: -b, i: -d},
+ {r: -b, i: d}
+ ];
+ }
+ return [];
+ }
+ }
+
+ quadratic.complex = function(x) {
+ if (!arguments.length) return complex;
+ complex = x;
+ return quadratic;
+ };
+
+ return quadratic;
+};
+// Constructs a multi-dimensional array filled with zeroes.
+science.zeroes = function(n) {
+ var i = -1,
+ a = [];
+ if (arguments.length === 1)
+ while (++i < n)
+ a[i] = 0;
+ else
+ while (++i < n)
+ a[i] = science.zeroes.apply(
+ this, Array.prototype.slice.call(arguments, 1));
+ return a;
+};
+science.vector = {};
+science.vector.cross = function(a, b) {
+ // TODO how to handle non-3D vectors?
+ // TODO handle 7D vectors?
+ return [
+ a[1] * b[2] - a[2] * b[1],
+ a[2] * b[0] - a[0] * b[2],
+ a[0] * b[1] - a[1] * b[0]
+ ];
+};
+science.vector.dot = function(a, b) {
+ var s = 0,
+ i = -1,
+ n = Math.min(a.length, b.length);
+ while (++i < n) s += a[i] * b[i];
+ return s;
+};
+science.vector.length = function(p) {
+ return Math.sqrt(science.vector.dot(p, p));
+};
+science.vector.normalize = function(p) {
+ var length = science.vector.length(p);
+ return p.map(function(d) { return d / length; });
+};
+// 4x4 matrix determinant.
+science.vector.determinant = function(matrix) {
+ var m = matrix[0].concat(matrix[1]).concat(matrix[2]).concat(matrix[3]);
+ return (
+ m[12] * m[9] * m[6] * m[3] - m[8] * m[13] * m[6] * m[3] -
+ m[12] * m[5] * m[10] * m[3] + m[4] * m[13] * m[10] * m[3] +
+ m[8] * m[5] * m[14] * m[3] - m[4] * m[9] * m[14] * m[3] -
+ m[12] * m[9] * m[2] * m[7] + m[8] * m[13] * m[2] * m[7] +
+ m[12] * m[1] * m[10] * m[7] - m[0] * m[13] * m[10] * m[7] -
+ m[8] * m[1] * m[14] * m[7] + m[0] * m[9] * m[14] * m[7] +
+ m[12] * m[5] * m[2] * m[11] - m[4] * m[13] * m[2] * m[11] -
+ m[12] * m[1] * m[6] * m[11] + m[0] * m[13] * m[6] * m[11] +
+ m[4] * m[1] * m[14] * m[11] - m[0] * m[5] * m[14] * m[11] -
+ m[8] * m[5] * m[2] * m[15] + m[4] * m[9] * m[2] * m[15] +
+ m[8] * m[1] * m[6] * m[15] - m[0] * m[9] * m[6] * m[15] -
+ m[4] * m[1] * m[10] * m[15] + m[0] * m[5] * m[10] * m[15]);
+};
+// Performs in-place Gauss-Jordan elimination.
+//
+// Based on Jarno Elonen's Python version (public domain):
+// http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html
+science.vector.gaussjordan = function(m, eps) {
+ if (!eps) eps = 1e-10;
+
+ var h = m.length,
+ w = m[0].length,
+ y = -1,
+ y2,
+ x;
+
+ while (++y < h) {
+ var maxrow = y;
+
+ // Find max pivot.
+ y2 = y; while (++y2 < h) {
+ if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y]))
+ maxrow = y2;
+ }
+
+ // Swap.
+ var tmp = m[y];
+ m[y] = m[maxrow];
+ m[maxrow] = tmp;
+
+ // Singular?
+ if (Math.abs(m[y][y]) <= eps) return false;
+
+ // Eliminate column y.
+ y2 = y; while (++y2 < h) {
+ var c = m[y2][y] / m[y][y];
+ x = y - 1; while (++x < w) {
+ m[y2][x] -= m[y][x] * c;
+ }
+ }
+ }
+
+ // Backsubstitute.
+ y = h; while (--y >= 0) {
+ var c = m[y][y];
+ y2 = -1; while (++y2 < y) {
+ x = w; while (--x >= y) {
+ m[y2][x] -= m[y][x] * m[y2][y] / c;
+ }
+ }
+ m[y][y] /= c;
+ // Normalize row y.
+ x = h - 1; while (++x < w) {
+ m[y][x] /= c;
+ }
+ }
+ return true;
+};
+// Find matrix inverse using Gauss-Jordan.
+science.vector.inverse = function(m) {
+ var n = m.length
+ i = -1;
+
+ // Check if the matrix is square.
+ if (n !== m[0].length) return;
+
+ // Augment with identity matrix I to get AI.
+ m = m.map(function(row, i) {
+ var identity = new Array(n),
+ j = -1;
+ while (++j < n) identity[j] = i === j ? 1 : 0;
+ return row.concat(identity);
+ });
+
+ // Compute IA^-1.
+ science.vector.gaussjordan(m);
+
+ // Remove identity matrix I to get A^-1.
+ while (++i < n) {
+ m[i] = m[i].slice(n);
+ }
+
+ return m;
+};
+science.vector.multiply = function(a, b) {
+ var m = a.length,
+ n = b[0].length,
+ p = b.length,
+ i = -1,
+ j,
+ k;
+ if (p !== a[0].length) throw {"error": "columns(a) != rows(b); " + a[0].length + " != " + p};
+ var ab = new Array(m);
+ while (++i < m) {
+ ab[i] = new Array(n);
+ j = -1; while(++j < n) {
+ var s = 0;
+ k = -1; while (++k < p) s += a[i][k] * b[k][j];
+ ab[i][j] = s;
+ }
+ }
+ return ab;
+};
+science.vector.transpose = function(a) {
+ var m = a.length,
+ n = a[0].length,
+ i = -1,
+ j,
+ b = new Array(n);
+ while (++i < n) {
+ b[i] = new Array(m);
+ j = -1; while (++j < m) b[i][j] = a[j][i];
+ }
+ return b;
+};
})()
View
27 lib/science/science.lin.js
@@ -0,0 +1,27 @@
+(function(){science.lin = {};
+/**
+ * Solves tridiagonal systems of linear equations.
+ *
+ * Source: http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
+ *
+ * @param {number[]} a
+ * @param {number[]} b
+ * @param {number[]} c
+ * @param {number[]} d
+ * @param {number[]} x
+ * @param {number} n
+ */
+science.lin.tridag = function(a, b, c, d, x, n) {
+ var i,
+ m;
+ for (i = 1; i < n; i++) {
+ m = a[i] / b[i - 1];
+ b[i] -= m * c[i - 1];
+ d[i] -= m * d[i - 1];
+ }
+ x[n - 1] = d[n - 1] / b[n - 1];
+ for (i = n - 2; i >= 0; i--) {
+ x[i] = (d[i] - c[i] * x[i + 1]) / b[i];
+ }
+};
+})()
View
1  lib/science/science.lin.min.js
@@ -0,0 +1 @@
+(function(){science.lin={},science.lin.tridag=function(a,b,c,d,e,f){var g,h;for(g=1;g<f;g++)h=a[g]/b[g-1],b[g]-=h*c[g-1],d[g]-=h*d[g-1];e[f-1]=d[f-1]/b[f-1];for(g=f-2;g>=0;g--)e[g]=(d[g]-c[g]*e[g+1])/b[g]}})()
View
2  lib/science/science.min.js
@@ -1 +1 @@
-(function(){science={version:"0.0.1"},science.functor=function(a){return typeof a=="function"?a:function(){return a}}})()
+(function(){science={version:"1.7.0"},science.ascending=function(a,b){return a-b},science.EULER=.5772156649015329,science.expm1=function(a){return a<1e-5&&a>-0.00001?a+.5*a*a:Math.exp(a)-1},science.functor=function(a){return typeof a=="function"?a:function(){return a}},science.hypot=function(a,b){a=Math.abs(a),b=Math.abs(b);var c,d;a>b?(c=a,d=b):(c=b,d=a);var e=d/c;return c*Math.sqrt(1+e*e)},science.quadratic=function(){function b(b,c,d){var e=c*c-4*b*d;if(e>0){e=Math.sqrt(e)/(2*b);return a?[{r:-c-e,i:0},{r:-c+e,i:0}]:[-c-e,-c+e]}if(e===0){e=-c/(2*b);return a?[{r:e,i:0}]:[e]}if(a){e=Math.sqrt(-e)/(2*b);return[{r:-c,i:-e},{r:-c,i:e}]}return[]}var a=!1;b.complex=function(c){if(!arguments.length)return a;a=c;return b};return b},science.zeroes=function(a){var b=-1,c=[];if(arguments.length===1)while(++b<a)c[b]=0;else while(++b<a)c[b]=science.zeroes.apply(this,Array.prototype.slice.call(arguments,1));return c},science.vector={},science.vector.cross=function(a,b){return[a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]},science.vector.dot=function(a,b){var c=0,d=-1,e=Math.min(a.length,b.length);while(++d<e)c+=a[d]*b[d];return c},science.vector.length=function(a){return Math.sqrt(science.vector.dot(a,a))},science.vector.normalize=function(a){var b=science.vector.length(a);return a.map(function(a){return a/b})},science.vector.determinant=function(a){var b=a[0].concat(a[1]).concat(a[2]).concat(a[3]);return b[12]*b[9]*b[6]*b[3]-b[8]*b[13]*b[6]*b[3]-b[12]*b[5]*b[10]*b[3]+b[4]*b[13]*b[10]*b[3]+b[8]*b[5]*b[14]*b[3]-b[4]*b[9]*b[14]*b[3]-b[12]*b[9]*b[2]*b[7]+b[8]*b[13]*b[2]*b[7]+b[12]*b[1]*b[10]*b[7]-b[0]*b[13]*b[10]*b[7]-b[8]*b[1]*b[14]*b[7]+b[0]*b[9]*b[14]*b[7]+b[12]*b[5]*b[2]*b[11]-b[4]*b[13]*b[2]*b[11]-b[12]*b[1]*b[6]*b[11]+b[0]*b[13]*b[6]*b[11]+b[4]*b[1]*b[14]*b[11]-b[0]*b[5]*b[14]*b[11]-b[8]*b[5]*b[2]*b[15]+b[4]*b[9]*b[2]*b[15]+b[8]*b[1]*b[6]*b[15]-b[0]*b[9]*b[6]*b[15]-b[4]*b[1]*b[10]*b[15]+b[0]*b[5]*b[10]*b[15]},science.vector.gaussjordan=function(a,b){b||(b=1e-10);var c=a.length,d=a[0].length,e=-1,f,g;while(++e<c){var h=e;f=e;while(++f<c)Math.abs(a[f][e])>Math.abs(a[h][e])&&(h=f);var i=a[e];a[e]=a[h],a[h]=i;if(Math.abs(a[e][e])<=b)return!1;f=e;while(++f<c){var j=a[f][e]/a[e][e];g=e-1;while(++g<d)a[f][g]-=a[e][g]*j}}e=c;while(--e>=0){var j=a[e][e];f=-1;while(++f<e){g=d;while(--g>=e)a[f][g]-=a[e][g]*a[f][e]/j}a[e][e]/=j,g=c-1;while(++g<d)a[e][g]/=j}return!0},science.vector.inverse=function(a){var b=a.length;i=-1;if(b===a[0].length){a=a.map(function(a,c){var d=Array(b),e=-1;while(++e<b)d[e]=c===e?1:0;return a.concat(d)}),science.vector.gaussjordan(a);while(++i<b)a[i]=a[i].slice(b);return a}},science.vector.multiply=function(a,b){var c=a.length,d=b[0].length,e=b.length,f=-1,g,h;if(e!==a[0].length)throw{error:"columns(a) != rows(b); "+a[0].length+" != "+e};var i=Array(c);while(++f<c){i[f]=Array(d),g=-1;while(++g<d){var j=0;h=-1;while(++h<e)j+=a[f][h]*b[h][g];i[f][g]=j}}return i},science.vector.transpose=function(a){var b=a.length,c=a[0].length,d=-1,e,f=Array(c);while(++d<c){f[d]=Array(b),e=-1;while(++e<b)f[d][e]=a[e][d]}return f}})()
View
554 lib/science/science.stats.js
@@ -19,6 +19,100 @@ science.stats.bandwidth = {
* Math.pow(x.length, -1/5);
}
};
+science.stats.distance = {
+ euclidean: function(a, b) {
+ var n = a.length,
+ i = -1,
+ s = 0,
+ x;
+ while (++i < n) {
+ x = a[i] - b[i];
+ s += x * x;
+ }
+ return Math.sqrt(s);
+ },
+ manhattan: function(a, b) {
+ var n = a.length,
+ i = -1,
+ s = 0;
+ while (++i < n) s += Math.abs(a[i] - b[i]);
+ return s;
+ },
+ minkowski: function(p) {
+ return function(a, b) {
+ var n = a.length,
+ i = -1,
+ s = 0;
+ while (++i < n) s += Math.pow(Math.abs(a[i] - b[i]), p);
+ return Math.pow(s, 1 / p);
+ };
+ },
+ chebyshev: function(a, b) {
+ var n = a.length,
+ i = -1,
+ max = 0,
+ x;
+ while (++i < n) {
+ x = Math.abs(a[i] - b[i]);
+ if (x > max) max = x;
+ }
+ return max;
+ },
+ hamming: function(a, b) {
+ var n = a.length,
+ i = -1,
+ d = 0;
+ while (++i < n) if (a[i] !== b[i]) d++;
+ return d;
+ },
+ jaccard: function(a, b) {
+ var n = a.length,
+ i = -1,
+ s = 0;
+ while (++i < n) if (a[i] === b[i]) s++;
+ return s / n;
+ },
+ braycurtis: function(a, b) {
+ var n = a.length,
+ i = -1,
+ s0 = 0,
+ s1 = 0,
+ ai,
+ bi;
+ while (++i < n) {
+ ai = a[i];
+ bi = b[i];
+ s0 += Math.abs(ai - bi);
+ s1 += Math.abs(ai + bi);
+ }
+ return s0 / s1;
+ }
+};
+// Based on implementation in http://picomath.org/.
+science.stats.erf = function(x) {
+ var a1 = 0.254829592,
+ a2 = -0.284496736,
+ a3 = 1.421413741,
+ a4 = -1.453152027,
+ a5 = 1.061405429,
+ p = 0.3275911;
+
+ // Save the sign of x
+ var sign = x < 0 ? -1 : 1;
+ if (x < 0) {
+ sign = -1;
+ x = -x;
+ }
+
+ // A&S formula 7.1.26
+ var t = 1 / (1 + p * x);
+ return sign * (
+ 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1)
+ * t * Math.exp(-x * x));
+};
+science.stats.phi = function(x) {
+ return .5 * (1 + science.stats.erf(x / Math.SQRT2));
+};
// See <http://en.wikipedia.org/wiki/Kernel_(statistics)>.
science.stats.kernel = {
uniform: function(u) {
@@ -94,10 +188,470 @@ science.stats.kde = function() {
return kde;
};
+// Based on figue implementation by Jean-Yves Delort.
+// http://code.google.com/p/figue/
+science.stats.kmeans = function() {
+ var distance = science.stats.distance.euclidean,
+ maxIterations = 1000,
+ k = 1;
+
+ function kmeans(vectors) {
+ var n = vectors.length,
+ assignments = [],
+ clusterSizes = [],
+ repeat = 1,
+ iterations = 0,
+ centroids = science_stats_kmeansRandom(k, vectors),
+ newCentroids,
+ i,
+ j,
+ x,
+ d,
+ min,
+ best;
+
+ while (repeat && iterations < maxIterations) {
+ // Assignment step.
+ j = -1; while (++j < k) {
+ clusterSizes[j] = 0;
+ }
+
+ i = -1; while (++i < n) {
+ x = vectors[i];
+ min = Infinity;
+ j = -1; while (++j < k) {
+ d = distance.call(this, centroids[j], x);
+ if (d < min) {
+ min = d;
+ best = j;
+ }
+ }
+ clusterSizes[assignments[i] = best]++;
+ }
+
+ // Update centroids step.
+ newCentroids = [];
+ i = -1; while (++i < n) {
+ x = assignments[i];
+ d = newCentroids[x];
+ if (d == null) newCentroids[x] = vectors[i].slice();
+ else {
+ j = -1; while (++j < d.length) {
+ d[j] += vectors[i][j];
+ }
+ }
+ }
+ j = -1; while (++j < k) {
+ x = newCentroids[j];
+ d = 1 / clusterSizes[j];
+ i = -1; while (++i < x.length) x[i] *= d;
+ }
+
+ // Check convergence.
+ repeat = 0;
+ j = -1; while (++j < k) {
+ if (!science_stats_kmeansCompare(newCentroids[j], centroids[j])) {
+ repeat = 1;
+ break;
+ }
+ }
+ centroids = newCentroids;
+ iterations++;
+ }
+ return {assignments: assignments, centroids: centroids};
+ }
+
+ kmeans.k = function(x) {
+ if (!arguments.length) return k;
+ k = x;
+ return kmeans;
+ };
+
+ kmeans.distance = function(x) {
+ if (!arguments.length) return distance;
+ distance = x;
+ return kmeans;
+ };
+
+ return kmeans;
+};
+
+function science_stats_kmeansCompare(a, b) {
+ if (!a || !b || a.length !== b.length) return false;
+ var n = a.length,
+ i = -1;
+ while (++i < n) if (a[i] !== b[i]) return false;
+ return true;
+}
+
+// Returns an array of k distinct vectors randomly selected from the input
+// array of vectors. Returns null if k > n or if there are less than k distinct
+// objects in vectors.
+function science_stats_kmeansRandom(k, vectors) {
+ var n = vectors.length;
+ if (k > n) return null;
+
+ var selected_vectors = [];
+ var selected_indices = [];
+ var tested_indices = {};
+ var tested = 0;
+ var selected = 0;
+ var i,
+ vector,
+ select;
+
+ while (selected < k) {
+ if (tested === n) return null;
+
+ var random_index = Math.floor(Math.random() * n);
+ if (random_index in tested_indices) continue;
+
+ tested_indices[random_index] = 1;
+ tested++;
+ vector = vectors[random_index];
+ select = true;
+ for (i = 0; i < selected; i++) {
+ if (science_stats_kmeansCompare(vector, selected_vectors[i])) {
+ select = false;
+ break;
+ }
+ }
+ if (select) {
+ selected_vectors[selected] = vector;
+ selected_indices[selected] = random_index;
+ selected++;
+ }
+ }
+ return selected_vectors;
+}
+science.stats.hcluster = function() {
+ var distance = science.stats.distance.euclidean,
+ linkage = "simple"; // simple, complete or average
+
+ function hcluster(vectors) {
+ var n = vectors.length,
+ dMin = [],
+ cSize = [],
+ distMatrix = [],
+ clusters = [],
+ c1,
+ c2,
+ c1Cluster,
+ c2Cluster,
+ p,
+ root,
+ i,
+ j;
+
+ // Initialise distance matrix and vector of closest clusters.
+ i = -1; while (++i < n) {
+ dMin[i] = 0;
+ distMatrix[i] = [];
+ j = -1; while (++j < n) {
+ distMatrix[i][j] = i === j ? Infinity : distance(vectors[i] , vectors[j]);
+ if (distMatrix[i][dMin[i]] > distMatrix[i][j]) dMin[i] = j;
+ }
+ }
+
+ // create leaves of the tree
+ i = -1; while (++i < n) {
+ clusters[i] = [];
+ clusters[i][0] = {
+ left: null,
+ right: null,
+ dist: 0,
+ centroid: vectors[i],
+ size: 1,
+ depth: 0
+ };
+ cSize[i] = 1;
+ }
+
+ // Main loop
+ for (p = 0; p < n-1; p++) {
+ // find the closest pair of clusters
+ c1 = 0;
+ for (i = 0; i < n; i++) {
+ if (distMatrix[i][dMin[i]] < distMatrix[c1][dMin[c1]]) c1 = i;
+ }
+ c2 = dMin[c1];
+
+ // create node to store cluster info
+ c1Cluster = clusters[c1][0];
+ c2Cluster = clusters[c2][0];
+
+ newCluster = {
+ left: c1Cluster,
+ right: c2Cluster,
+ dist: distMatrix[c1][c2],
+ centroid: calculateCentroid(c1Cluster.size, c1Cluster.centroid,
+ c2Cluster.size, c2Cluster.centroid),
+ size: c1Cluster.size + c2Cluster.size,
+ depth: 1 + Math.max(c1Cluster.depth, c2Cluster.depth)
+ };
+ clusters[c1].splice(0, 0, newCluster);
+ cSize[c1] += cSize[c2];
+
+ // overwrite row c1 with respect to the linkage type
+ for (j = 0; j < n; j++) {
+ switch (linkage) {
+ case "single":
+ if (distMatrix[c1][j] > distMatrix[c2][j])
+ distMatrix[j][c1] = distMatrix[c1][j] = distMatrix[c2][j];
+ break;
+ case "complete":
+ if (distMatrix[c1][j] < distMatrix[c2][j])
+ distMatrix[j][c1] = distMatrix[c1][j] = distMatrix[c2][j];
+ break;
+ case "average":
+ distMatrix[j][c1] = distMatrix[c1][j] = (cSize[c1] * distMatrix[c1][j] + cSize[c2] * distMatrix[c2][j]) / (cSize[c1] + cSize[j]);
+ break;
+ }
+ }
+ distMatrix[c1][c1] = Infinity;
+
+ // infinity ­out old row c2 and column c2
+ for (i = 0; i < n; i++)
+ distMatrix[i][c2] = distMatrix[c2][i] = Infinity;
+
+ // update dmin and replace ones that previous pointed to c2 to point to c1
+ for (j = 0; j < n; j++) {
+ if (dMin[j] == c2) dMin[j] = c1;
+ if (distMatrix[c1][j] < distMatrix[c1][dMin[c1]]) dMin[c1] = j;
+ }
+
+ // keep track of the last added cluster
+ root = newCluster;
+ }
+
+ return root;
+ }
+
+ hcluster.distance = function(x) {
+ if (!arguments.length) return distance;
+ distance = x;
+ return hcluster;
+ };
+
+ return hcluster;
+};
+
+function calculateCentroid(c1Size, c1Centroid, c2Size, c2Centroid) {
+ var newCentroid = [],
+ newSize = c1Size + c2Size,
+ n = c1Centroid.length,
+ i = -1;
+ while (++i < n) {
+ newCentroid[i] = (c1Size * c1Centroid[i] + c2Size * c2Centroid[i]) / newSize;
+ }
+ return newCentroid;
+}
science.stats.iqr = function(x) {
var quartiles = science.stats.quantiles(x, [.25, .75]);
return quartiles[1] - quartiles[0];
};
+// Based on org.apache.commons.math.analysis.interpolation.LoessInterpolator
+// from http://commons.apache.org/math/
+science.stats.loess = function() {
+ var bandwidth = .3,
+ robustnessIters = 2,
+ accuracy = 1e-12;
+
+ function smooth(xval, yval, weights) {
+ var n = xval.length,
+ i;
+
+ if (n !== yval.length) throw {error: "Mismatched array lengths"};
+ if (n == 0) throw {error: "At least one point required."};
+
+ if (arguments.length < 3) {
+ weights = [];
+ i = -1; while (++i < n) weights[i] = 1;
+ }
+
+ science_stats_loessFiniteReal(xval);
+ science_stats_loessFiniteReal(yval);
+ science_stats_loessFiniteReal(weights);
+ science_stats_loessStrictlyIncreasing(xval);
+
+ if (n == 1) return [yval[0]];
+ if (n == 2) return [yval[0], yval[1]];
+
+ var bandwidthInPoints = Math.floor(bandwidth * n);
+
+ if (bandwidthInPoints < 2) throw {error: "Bandwidth too small."};
+
+ var res = [],
+ residuals = [],
+ robustnessWeights = [];
+
+ // Do an initial fit and 'robustnessIters' robustness iterations.
+ // This is equivalent to doing 'robustnessIters+1' robustness iterations
+ // starting with all robustness weights set to 1.
+ i = -1; while (++i < n) {
+ res[i] = 0;
+ residuals[i] = 0;
+ robustnessWeights[i] = 1;
+ }
+
+ var iter = -1;
+ while (++iter <= robustnessIters) {
+ var bandwidthInterval = [0, bandwidthInPoints - 1];
+ // At each x, compute a local weighted linear regression
+ var x;
+ i = -1; while (++i < n) {
+ x = xval[i];
+
+ // Find out the interval of source points on which
+ // a regression is to be made.
+ if (i > 0) {
+ science_stats_loessUpdateBandwidthInterval(xval, weights, i, bandwidthInterval);
+ }
+
+ var ileft = bandwidthInterval[0],
+ iright = bandwidthInterval[1];
+
+ // Compute the point of the bandwidth interval that is
+ // farthest from x
+ var edge = (xval[i] - xval[ileft]) > (xval[iright] - xval[i]) ? ileft : iright;
+
+ // Compute a least-squares linear fit weighted by
+ // the product of robustness weights and the tricube
+ // weight function.
+ // See http://en.wikipedia.org/wiki/Linear_regression
+ // (section "Univariate linear case")
+ // and http://en.wikipedia.org/wiki/Weighted_least_squares
+ // (section "Weighted least squares")
+ var sumWeights = 0,
+ sumX = 0,
+ sumXSquared = 0,
+ sumY = 0,
+ sumXY = 0,
+ denom = Math.abs(1 / (xval[edge] - x));
+
+ for (var k = ileft; k <= iright; ++k) {
+ var xk = xval[k],
+ yk = yval[k],
+ dist = k < i ? x - xk : xk - x,
+ w = science_stats_loessTricube(dist * denom) * robustnessWeights[k] * weights[k],
+ xkw = xk * w;
+ sumWeights += w;
+ sumX += xkw;
+ sumXSquared += xk * xkw;
+ sumY += yk * w;
+ sumXY += yk * xkw;
+ }
+
+ var meanX = sumX / sumWeights,
+ meanY = sumY / sumWeights,
+ meanXY = sumXY / sumWeights,
+ meanXSquared = sumXSquared / sumWeights;
+
+ var beta = (Math.sqrt(Math.abs(meanXSquared - meanX * meanX)) < accuracy)
+ ? 0 : ((meanXY - meanX * meanY) / (meanXSquared - meanX * meanX));
+
+ var alpha = meanY - beta * meanX;
+
+ res[i] = beta * x + alpha;
+ residuals[i] = Math.abs(yval[i] - res[i]);
+ }
+
+ // No need to recompute the robustness weights at the last
+ // iteration, they won't be needed anymore
+ if (iter === robustnessIters) {
+ break;
+ }
+
+ // Recompute the robustness weights.
+
+ // Find the median residual.
+ var sortedResiduals = residuals.slice();
+ sortedResiduals.sort();
+ var medianResidual = sortedResiduals[Math.floor(n / 2)];
+
+ if (Math.abs(medianResidual) < accuracy)
+ break;
+
+ var arg,
+ w;
+ i = -1; while (++i < n) {
+ arg = residuals[i] / (6 * medianResidual);
+ robustnessWeights[i] = (arg >= 1) ? 0 : ((w = 1 - arg * arg) * w);
+ }
+ }
+
+ return res;
+ }
+
+ smooth.bandwidth = function(x) {
+ if (!arguments.length) return x;
+ bandwidth = x;
+ return smooth;
+ };
+
+ smooth.robustnessIterations = function(x) {
+ if (!arguments.length) return x;
+ robustnessIters = x;
+ return smooth;
+ };
+
+ smooth.accuracy = function(x) {
+ if (!arguments.length) return x;
+ accuracy = x;
+ return smooth;
+ };
+
+ return smooth;
+};
+
+function science_stats_loessFiniteReal(values) {
+ var n = values.length,
+ i = -1;
+
+ while (++i < n) if (!isFinite(values[i])) return false;
+
+ return true;
+}
+
+function science_stats_loessStrictlyIncreasing(xval) {
+ var n = xval.length,
+ i = 0;
+
+ while (++i < n) if (xval[i - 1] >= xval[i]) return false;
+
+ return true;
+}
+
+// Compute the tricube weight function.
+// http://en.wikipedia.org/wiki/Local_regression#Weight_function
+function science_stats_loessTricube(x) {
+ return (x = 1 - x * x * x) * x * x;
+}
+
+// Given an index interval into xval that embraces a certain number of
+// points closest to xval[i-1], update the interval so that it embraces
+// the same number of points closest to xval[i], ignoring zero weights.
+function science_stats_loessUpdateBandwidthInterval(
+ xval, weights, i, bandwidthInterval) {
+
+ var left = bandwidthInterval[0],
+ right = bandwidthInterval[1];
+
+ // The right edge should be adjusted if the next point to the right
+ // is closer to xval[i] than the leftmost point of the current interval
+ var nextRight = science_stats_loessNextNonzero(weights, right);
+ if ((nextRight < xval.length) && (xval[nextRight] - xval[i]) < (xval[i] - xval[left])) {
+ var nextLeft = science_stats_loessNextNonzero(weights, left);
+ bandwidthInterval[0] = nextLeft;
+ bandwidthInterval[1] = nextRight;
+ }
+}
+
+function science_stats_loessNextNonzero(weights, i) {
+ var j = i + 1;
+ while (j < weights.length && weights[j] === 0) j++;
+ return j;
+}
// Welford's algorithm.
science.stats.mean = function(x) {
var n = x.length;
View
2  lib/science/science.stats.min.js
@@ -1 +1 @@
-(function(){science.stats={},science.stats.bandwidth={nrd0:function(a){var b=Math.sqrt(science.stats.variance(a));(lo=Math.min(b,science.stats.iqr(a)/1.34))||(lo=b)||(lo=Math.abs(a[1]))||(lo=1);return.9*lo*Math.pow(a.length,-0.2)},nrd:function(a){var b=science.stats.iqr(a)/1.34;return 1.06*Math.min(Math.sqrt(science.stats.variance(a)),b)*Math.pow(a.length,-0.2)}},science.stats.kernel={uniform:function(a){if(a<=1&&a>=-1)return.5;return 0},triangular:function(a){if(a<=1&&a>=-1)return 1-Math.abs(a);return 0},epanechnikov:function(a){if(a<=1&&a>=-1)return.75*(1-a*a);return 0},quartic:function(a){if(a<=1&&a>=-1){var b=1-a*a;return.9375*b*b}return 0},triweight:function(a){if(a<=1&&a>=-1){var b=1-a*a;return 35/32*b*b*b}return 0},gaussian:function(a){return 1/Math.sqrt(2*Math.PI)*Math.exp(-0.5*a*a)},cosine:function(a){if(a<=1&&a>=-1)return Math.PI/4*Math.cos(Math.PI/2*a);return 0}},science.stats.kde=function(){function d(d,e){var f=c.call(this,b);return d.map(function(c){var d=-1,e=0,g=b.length;while(++d<g)e+=a((c-b[d])/f);return[c,e/f/g]})}var a=science.stats.kernel.gaussian,b=[],c=science.stats.bandwidth.nrd;d.kernel=function(b){if(!arguments.length)return a;a=b;return d},d.sample=function(a){if(!arguments.length)return b;b=a;return d},d.bandwidth=function(a){if(!arguments.length)return c;c=science.functor(a);return d};return d},science.stats.iqr=function(a){var b=science.stats.quantiles(a,[.25,.75]);return b[1]-b[0]},science.stats.mean=function(a){var b=a.length;if(b===0)return NaN;var c=0,d=-1;while(++d<b)c+=(a[d]-c)/(d+1);return c},science.stats.median=function(a){return science.stats.quantiles(a,[.5])[0]},science.stats.mode=function(a){a=a.slice().sort(science.ascending);var b,c=a.length,d=-1,e=d,f=null,g=0,h,i;while(++d<c)(i=a[d])!==f&&((h=d-e)>g&&(g=h,b=f),f=i,e=d);return b},science.stats.quantiles=function(a,b){a=a.slice().sort(science.ascending);var c=a.length-1;return b.map(function(b){if(b===0)return a[0];if(b===1)return a[c];var e=1+b*c,f=Math.floor(e),g=e-f,h=a[f-1];return g===0?h:h+g*(a[f]-h)})},science.stats.variance=function(a){var b=a.length;if(b<1)return NaN;if(b===1)return 0;var c=science.stats.mean(a),d=-1,e=0;while(++d<b){var f=a[d]-c;e+=f*f}return e/(b-1)}})()
+(function(){function h(a,b){var c=b+1;while(c<a.length&&a[c]===0)c++;return c}function g(a,b,c,d){var e=d[0],f=d[1],g=h(b,f);if(g<a.length&&a[g]-a[c]<a[c]-a[e]){var i=h(b,e);d[0]=i,d[1]=g}}function f(a){return(a=1-a*a*a)*a*a}function e(a){var b=a.length,c=0;while(++c<b)if(a[c-1]>=a[c])return!1;return!0}function d(a){var b=a.length,c=-1;while(++c<b)if(!isFinite(a[c]))return!1;return!0}function c(a,b,c,d){var e=[],f=a+c,g=b.length,h=-1;while(++h<g)e[h]=(a*b[h]+c*d[h])/f;return e}function b(b,c){var d=c.length;if(b>d)return null;var e=[],f=[],g={},h=0,i=0,j,k,l;while(i<b){if(h===d)return null;var m=Math.floor(Math.random()*d);if(m in g)continue;g[m]=1,h++,k=c[m],l=!0;for(j=0;j<i;j++)if(a(k,e[j])){l=!1;break}l&&(e[i]=k,f[i]=m,i++)}return e}function a(a,b){if(!a||!b||a.length!==b.length)return!1;var c=a.length,d=-1;while(++d<c)if(a[d]!==b[d])return!1;return!0}science.stats={},science.stats.bandwidth={nrd0:function(a){var b=Math.sqrt(science.stats.variance(a));(lo=Math.min(b,science.stats.iqr(a)/1.34))||(lo=b)||(lo=Math.abs(a[1]))||(lo=1);return.9*lo*Math.pow(a.length,-0.2)},nrd:function(a){var b=science.stats.iqr(a)/1.34;return 1.06*Math.min(Math.sqrt(science.stats.variance(a)),b)*Math.pow(a.length,-0.2)}},science.stats.distance={euclidean:function(a,b){var c=a.length,d=-1,e=0,f;while(++d<c)f=a[d]-b[d],e+=f*f;return Math.sqrt(e)},manhattan:function(a,b){var c=a.length,d=-1,e=0;while(++d<c)e+=Math.abs(a[d]-b[d]);return e},minkowski:function(a){return function(b,c){var d=b.length,e=-1,f=0;while(++e<d)f+=Math.pow(Math.abs(b[e]-c[e]),a);return Math.pow(f,1/a)}},chebyshev:function(a,b){var c=a.length,d=-1,e=0,f;while(++d<c)f=Math.abs(a[d]-b[d]),f>e&&(e=f);return e},hamming:function(a,b){var c=a.length,d=-1,e=0;while(++d<c)a[d]!==b[d]&&e++;return e},jaccard:function(a,b){var c=a.length,d=-1,e=0;while(++d<c)a[d]===b[d]&&e++;return e/c},braycurtis:function(a,b){var c=a.length,d=-1,e=0,f=0,g,h;while(++d<c)g=a[d],h=b[d],e+=Math.abs(g-h),f+=Math.abs(g+h);return e/f}},science.stats.erf=function(a){var b=.254829592,c=-0.284496736,d=1.421413741,e=-1.453152027,f=1.061405429,g=.3275911,h=a<0?-1:1;a<0&&(h=-1,a=-a);var i=1/(1+g*a);return h*(1-((((f*i+e)*i+d)*i+c)*i+b)*i*Math.exp(-a*a))},science.stats.phi=function(a){return.5*(1+science.stats.erf(a/Math.SQRT2))},science.stats.kernel={uniform:function(a){return a<=1&&a>=-1?.5:0},triangular:function(a){return a<=1&&a>=-1?1-Math.abs(a):0},epanechnikov:function(a){return a<=1&&a>=-1?.75*(1-a*a):0},quartic:function(a){if(a<=1&&a>=-1){var b=1-a*a;return.9375*b*b}return 0},triweight:function(a){if(a<=1&&a>=-1){var b=1-a*a;return 35/32*b*b*b}return 0},gaussian:function(a){return 1/Math.sqrt(2*Math.PI)*Math.exp(-0.5*a*a)},cosine:function(a){return a<=1&&a>=-1?Math.PI/4*Math.cos(Math.PI/2*a):0}},science.stats.kde=function(){function d(d,e){var f=c.call(this,b);return d.map(function(c){var d=-1,e=0,g=b.length;while(++d<g)e+=a((c-b[d])/f);return[c,e/f/g]})}var a=science.stats.kernel.gaussian,b=[],c=science.stats.bandwidth.nrd;d.kernel=function(b){if(!arguments.length)return a;a=b;return d},d.sample=function(a){if(!arguments.length)return b;b=a;return d},d.bandwidth=function(a){if(!arguments.length)return c;c=science.functor(a);return d};return d},science.stats.kmeans=function(){function f(f){var g=f.length,h=[],i=[],j=1,l=0,m=b(e,f),n,o,p,q,r,s,t;while(j&&l<d){p=-1;while(++p<e)i[p]=0;o=-1;while(++o<g){q=f[o],s=Infinity,p=-1;while(++p<e)r=c.call(this,m[p],q),r<s&&(s=r,t=p);i[h[o]=t]++}n=[],o=-1;while(++o<g){q=h[o],r=n[q];if(r==null)n[q]=f[o].slice();else{p=-1;while(++p<r.length)r[p]+=f[o][p]}}p=-1;while(++p<e){q=n[p],r=1/i[p],o=-1;while(++o<q.length)q[o]*=r}j=0,p=-1;while(++p<e)if(!a(n[p],m[p])){j=1;break}m=n,l++}return{assignments:h,centroids:m}}var c=science.stats.distance.euclidean,d=1e3,e=1;f.k=function(a){if(!arguments.length)return e;e=a;return f},f.distance=function(a){if(!arguments.length)return c;c=a;return f};return f},science.stats.hcluster=function(){function d(d){var e=d.length,f=[],g=[],h=[],i=[],j,k,l,m,n,o,p,q;p=-1;while(++p<e){f[p]=0,h[p]=[],q=-1;while(++q<e)h[p][q]=p===q?Infinity:a(d[p],d[q]),h[p][f[p]]>h[p][q]&&(f[p]=q)}p=-1;while(++p<e)i[p]=[],i[p][0]={left:null,right:null,dist:0,centroid:d[p],size:1,depth:0},g[p]=1;for(n=0;n<e-1;n++){j=0;for(p=0;p<e;p++)h[p][f[p]]<h[j][f[j]]&&(j=p);k=f[j],l=i[j][0],m=i[k][0],newCluster={left:l,right:m,dist:h[j][k],centroid:c(l.size,l.centroid,m.size,m.centroid),size:l.size+m.size,depth:1+Math.max(l.depth,m.depth)},i[j].splice(0,0,newCluster),g[j]+=g[k];for(q=0;q<e;q++)switch(b){case"single":h[j][q]>h[k][q]&&(h[q][j]=h[j][q]=h[k][q]);break;case"complete":h[j][q]<h[k][q]&&(h[q][j]=h[j][q]=h[k][q]);break;case"average":h[q][j]=h[j][q]=(g[j]*h[j][q]+g[k]*h[k][q])/(g[j]+g[q])}h[j][j]=Infinity;for(p=0;p<e;p++)h[p][k]=h[k][p]=Infinity;for(q=0;q<e;q++)f[q]==k&&(f[q]=j),h[j][q]<h[j][f[j]]&&(f[j]=q);o=newCluster}return o}var a=science.stats.distance.euclidean,b="simple";d.distance=function(b){if(!arguments.length)return a;a=b;return d};return d},science.stats.iqr=function(a){var b=science.stats.quantiles(a,[.25,.75]);return b[1]-b[0]},science.stats.loess=function(){function h(h,i,j){var k=h.length,l;if(k!==i.length)throw{error:"Mismatched array lengths"};if(k==0)throw{error:"At least one point required."};if(arguments.length<3){j=[],l=-1;while(++l<k)j[l]=1}d(h),d(i),d(j),e(h);if(k==1)return[i[0]];if(k==2)return[i[0],i[1]];var m=Math.floor(a*k);if(m<2)throw{error:"Bandwidth too small."};var n=[],o=[],p=[];l=-1;while(++l<k)n[l]=0,o[l]=0,p[l]=1;var q=-1;while(++q<=b){var r=[0,m-1],s;l=-1;while(++l<k){s=h[l],l>0&&g(h,j,l,r);var t=r[0],u=r[1],v=h[l]-h[t]>h[u]-h[l]?t:u,w=0,x=0,y=0,z=0,A=0,B=Math.abs(1/(h[v]-s));for(var C=t;C<=u;++C){var D=h[C],E=i[C],F=C<l?s-D:D-s,G=f(F*B)*p[C]*j[C],H=D*G;w+=G,x+=H,y+=D*H,z+=E*G,A+=E*H}var I=x/w,J=z/w,K=A/w,L=y/w,M=Math.sqrt(Math.abs(L-I*I))<c?0:(K-I*J)/(L-I*I),N=J-M*I;n[l]=M*s+N,o[l]=Math.abs(i[l]-n[l])}if(q===b)break;var O=o.slice();O.sort();var P=O[Math.floor(k/2)];if(Math.abs(P)<c)break;var Q,G;l=-1;while(++l<k)Q=o[l]/(6*P),p[l]=Q>=1?0:(G=1-Q*Q)*G}return n}var a=.3,b=2,c=1e-12;h.bandwidth=function(b){if(!arguments.length)return b;a=b;return h},h.robustnessIterations=function(a){if(!arguments.length)return a;b=a;return h},h.accuracy=function(a){if(!arguments.length)return a;c=a;return h};return h},science.stats.mean=function(a){var b=a.length;if(b===0)return NaN;var c=0,d=-1;while(++d<b)c+=(a[d]-c)/(d+1);return c},science.stats.median=function(a){return science.stats.quantiles(a,[.5])[0]},science.stats.mode=function(a){a=a.slice().sort(science.ascending);var b,c=a.length,d=-1,e=d,f=null,g=0,h,i;while(++d<c)(i=a[d])!==f&&((h=d-e)>g&&(g=h,b=f),f=i,e=d);return b},science.stats.quantiles=function(a,b){a=a.slice().sort(science.ascending);var c=a.length-1;return b.map(function(b){if(b===0)return a[0];if(b===1)return a[c];var d=1+b*c,e=Math.floor(d),f=d-e,g=a[e-1];return f===0?g:g+f*(a[e]-g)})},science.stats.variance=function(a){var b=a.length;if(b<1)return NaN;if(b===1)return 0;var c=science.stats.mean(a),d=-1,e=0;while(++d<b){var f=a[d]-c;e+=f*f}return e/(b-1)}})()
Please sign in to comment.
Something went wrong with that request. Please try again.