1 (function(){science.stats = {};

     2 // Bandwidth selectors for Gaussian kernels.

     3 // Based on R's implementations in `stats.bw`.

     4 science.stats.bandwidth = {

     5 

     6   // Silverman, B. W. (1986) Density Estimation. London: Chapman and Hall.

     7   nrd0: function(x) {

     8     var hi = Math.sqrt(science.stats.variance(x));

     9     if (!(lo = Math.min(hi, science.stats.iqr(x) / 1.34)))

    10       (lo = hi) || (lo = Math.abs(x[1])) || (lo = 1);

    11     return .9 * lo * Math.pow(x.length, -.2);

    12   },

    13 

    14   // Scott, D. W. (1992) Multivariate Density Estimation: Theory, Practice, and

    15   // Visualization. Wiley.

    16   nrd: function(x) {

    17     var h = science.stats.iqr(x) / 1.34;

    18     return 1.06 * Math.min(Math.sqrt(science.stats.variance(x)), h)

    19       * Math.pow(x.length, -1/5);

    20   }

    21 };

    22 science.stats.distance = {

    23   euclidean: function(a, b) {

    24     var n = a.length,

    25         i = -1,

    26         s = 0,

    27         x;

    28     while (++i < n) {

    29       x = a[i] - b[i];

    30       s += x * x;

    31     }

    32     return Math.sqrt(s);

    33   },

    34   manhattan: function(a, b) {

    35     var n = a.length,

    36         i = -1,

    37         s = 0;

    38     while (++i < n) s += Math.abs(a[i] - b[i]);

    39     return s;

    40   },

    41   minkowski: function(p) {

    42     return function(a, b) {

    43       var n = a.length,

    44           i = -1,

    45           s = 0;

    46       while (++i < n) s += Math.pow(Math.abs(a[i] - b[i]), p);

    47       return Math.pow(s, 1 / p);

    48     };

    49   },

    50   chebyshev: function(a, b) {

    51     var n = a.length,

    52         i = -1,

    53         max = 0,

    54         x;

    55     while (++i < n) {

    56       x = Math.abs(a[i] - b[i]);

    57       if (x > max) max = x;

    58     }

    59     return max;

    60   },

    61   hamming: function(a, b) {

    62     var n = a.length,

    63         i = -1,

    64         d = 0;

    65     while (++i < n) if (a[i] !== b[i]) d++;

    66     return d;

    67   },

    68   jaccard: function(a, b) {

    69     var n = a.length,

    70         i = -1,

    71         s = 0;

    72     while (++i < n) if (a[i] === b[i]) s++;

    73     return s / n;

    74   },

    75   braycurtis: function(a, b) {

    76     var n = a.length,

    77         i = -1,

    78         s0 = 0,

    79         s1 = 0,

    80         ai,

    81         bi;

    82     while (++i < n) {

    83       ai = a[i];

    84       bi = b[i];

    85       s0 += Math.abs(ai - bi);

    86       s1 += Math.abs(ai + bi);

    87     }

    88     return s0 / s1;

    89   }

    90 };

    91 // Based on implementation in http://picomath.org/.

    92 science.stats.erf = function(x) {

    93   var a1 =  0.254829592,

    94       a2 = -0.284496736,

    95       a3 =  1.421413741,

    96       a4 = -1.453152027,

    97       a5 =  1.061405429,

    98       p  =  0.3275911;

    99 

   100   // Save the sign of x

   101   var sign = x < 0 ? -1 : 1;

   102   if (x < 0) {

   103     sign = -1;

   104     x = -x;

   105   }

   106 

   107   // A&S formula 7.1.26

   108   var t = 1 / (1 + p * x);

   109   return sign * (

   110     1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1)

   111     * t * Math.exp(-x * x));

   112 };

   113 science.stats.phi = function(x) {

   114   return .5 * (1 + science.stats.erf(x / Math.SQRT2));

   115 };

   116 // See <http://en.wikipedia.org/wiki/Kernel_(statistics)>.

   117 science.stats.kernel = {

   118   uniform: function(u) {

   119     if (u <= 1 && u >= -1) return .5;

   120     return 0;

   121   },

   122   triangular: function(u) {

   123     if (u <= 1 && u >= -1) return 1 - Math.abs(u);

   124     return 0;

   125   },

   126   epanechnikov: function(u) {

   127     if (u <= 1 && u >= -1) return .75 * (1 - u * u);

   128     return 0;

   129   },

   130   quartic: function(u) {

   131     if (u <= 1 && u >= -1) {

   132       var tmp = 1 - u * u;

   133       return (15 / 16) * tmp * tmp;

   134     }

   135     return 0;

   136   },

   137   triweight: function(u) {

   138     if (u <= 1 && u >= -1) {

   139       var tmp = 1 - u * u;

   140       return (35 / 32) * tmp * tmp * tmp;

   141     }

   142     return 0;

   143   },

   144   gaussian: function(u) {

   145     return 1 / Math.sqrt(2 * Math.PI) * Math.exp(-.5 * u * u);

   146   },

   147   cosine: function(u) {

   148     if (u <= 1 && u >= -1) return Math.PI / 4 * Math.cos(Math.PI / 2 * u);

   149     return 0;

   150   }

   151 };

   152 // http://exploringdata.net/den_trac.htm

   153 science.stats.kde = function() {

   154   var kernel = science.stats.kernel.gaussian,

   155       sample = [],

   156       bandwidth = science.stats.bandwidth.nrd;

   157 

   158   function kde(points, i) {

   159     var bw = bandwidth.call(this, sample);

   160     return points.map(function(x) {

   161       var i = -1,

   162           y = 0,

   163           n = sample.length;

   164       while (++i < n) {

   165         y += kernel((x - sample[i]) / bw);

   166       }

   167       return [x, y / bw / n];

   168     });

   169   }

   170 

   171   kde.kernel = function(x) {

   172     if (!arguments.length) return kernel;

   173     kernel = x;

   174     return kde;

   175   };

   176 

   177   kde.sample = function(x) {

   178     if (!arguments.length) return sample;

   179     sample = x;

   180     return kde;

   181   };

   182 

   183   kde.bandwidth = function(x) {

   184     if (!arguments.length) return bandwidth;

   185     bandwidth = science.functor(x);

   186     return kde;

   187   };

   188 

   189   return kde;

   190 };

   191 // Based on figue implementation by Jean-Yves Delort.

   192 // http://code.google.com/p/figue/

   193 science.stats.kmeans = function() {

   194   var distance = science.stats.distance.euclidean,

   195       maxIterations = 1000,

   196       k = 1;

   197 

   198   function kmeans(vectors) {

   199     var n = vectors.length,

   200         assignments = [],

   201         clusterSizes = [],

   202         repeat = 1,

   203         iterations = 0,

   204         centroids = science_stats_kmeansRandom(k, vectors),

   205         newCentroids,

   206         i,

   207         j,

   208         x,

   209         d,

   210         min,

   211         best;

   212 

   213     while (repeat && iterations < maxIterations) {

   214       // Assignment step.

   215       j = -1; while (++j < k) {

   216         clusterSizes[j] = 0;

   217       }

   218 

   219       i = -1; while (++i < n) {

   220         x = vectors[i];

   221         min = Infinity;

   222         j = -1; while (++j < k) {

   223           d = distance.call(this, centroids[j], x);

   224           if (d < min) {

   225             min = d;

   226             best = j;

   227           }

   228         }

   229         clusterSizes[assignments[i] = best]++;

   230       }

   231 

   232       // Update centroids step.

   233       newCentroids = [];

   234       i = -1; while (++i < n) {

   235         x = assignments[i];

   236         d = newCentroids[x];

   237         if (d == null) newCentroids[x] = vectors[i].slice();

   238         else {

   239           j = -1; while (++j < d.length) {

   240             d[j] += vectors[i][j];

   241           }

   242         }

   243       }

   244       j = -1; while (++j < k) {

   245         x = newCentroids[j];

   246         d = 1 / clusterSizes[j];

   247         i = -1; while (++i < x.length) x[i] *= d;

   248       }

   249 

   250       // Check convergence.

   251       repeat = 0;

   252       j = -1; while (++j < k) {

   253         if (!science_stats_kmeansCompare(newCentroids[j], centroids[j])) {

   254           repeat = 1;

   255           break;

   256         }

   257       }

   258       centroids = newCentroids;

   259       iterations++;

   260     }

   261     return {assignments: assignments, centroids: centroids};

   262   }

   263 

   264   kmeans.k = function(x) {

   265     if (!arguments.length) return k;

   266     k = x;

   267     return kmeans;

   268   };

   269 

   270   kmeans.distance = function(x) {

   271     if (!arguments.length) return distance;

   272     distance = x;

   273     return kmeans;

   274   };

   275 

   276   return kmeans;

   277 };

   278 

   279 function science_stats_kmeansCompare(a, b) {

   280   if (!a || !b || a.length !== b.length) return false;

   281   var n = a.length,

   282       i = -1;

   283   while (++i < n) if (a[i] !== b[i]) return false;

   284   return true;

   285 }

   286 

   287 // Returns an array of k distinct vectors randomly selected from the input

   288 // array of vectors. Returns null if k > n or if there are less than k distinct

   289 // objects in vectors.

   290 function science_stats_kmeansRandom(k, vectors) {

   291   var n = vectors.length;

   292   if (k > n) return null;

   293   

   294   var selected_vectors = [];

   295   var selected_indices = [];

   296   var tested_indices = {};

   297   var tested = 0;

   298   var selected = 0;

   299   var i,

   300       vector,

   301       select;

   302 

   303   while (selected < k) {

   304     if (tested === n) return null;

   305     

   306     var random_index = Math.floor(Math.random() * n);

   307     if (random_index in tested_indices) continue;

   308     

   309     tested_indices[random_index] = 1;

   310     tested++;

   311     vector = vectors[random_index];

   312     select = true;

   313     for (i = 0; i < selected; i++) {

   314       if (science_stats_kmeansCompare(vector, selected_vectors[i])) {

   315         select = false;

   316         break;

   317       }

   318     }

   319     if (select) {

   320       selected_vectors[selected] = vector;

   321       selected_indices[selected] = random_index;

   322       selected++;

   323     }

   324   }

   325   return selected_vectors;

   326 }

   327 science.stats.hcluster = function() {

   328   var distance = science.stats.distance.euclidean,

   329       linkage = "simple"; // simple, complete or average

   330 

   331   function hcluster(vectors) {

   332     var n = vectors.length,

   333         dMin = [],

   334         cSize = [],

   335         distMatrix = [],

   336         clusters = [],

   337         c1,

   338         c2,

   339         c1Cluster,

   340         c2Cluster,

   341         p,

   342         root,

   343         i,

   344         j;

   345 

   346     // Initialise distance matrix and vector of closest clusters.

   347     i = -1; while (++i < n) {

   348       dMin[i] = 0;

   349       distMatrix[i] = [];

   350       j = -1; while (++j < n) {

   351         distMatrix[i][j] = i === j ? Infinity : distance(vectors[i] , vectors[j]);

   352         if (distMatrix[i][dMin[i]] > distMatrix[i][j]) dMin[i] = j;

   353       }

   354     }

   355 

   356     // create leaves of the tree

   357     i = -1; while (++i < n) {

   358       clusters[i] = [];

   359       clusters[i][0] = {

   360         left: null,

   361         right: null,

   362         dist: 0,

   363         centroid: vectors[i],

   364         size: 1,

   365         depth: 0

   366       };

   367       cSize[i] = 1;

   368     }

   369 

   370     // Main loop

   371     for (p = 0; p < n-1; p++) {

   372       // find the closest pair of clusters

   373       c1 = 0;

   374       for (i = 0; i < n; i++) {

   375         if (distMatrix[i][dMin[i]] < distMatrix[c1][dMin[c1]]) c1 = i;

   376       }

   377       c2 = dMin[c1];

   378 

   379       // create node to store cluster info 

   380       c1Cluster = clusters[c1][0];

   381       c2Cluster = clusters[c2][0];

   382 

   383       newCluster = {

   384         left: c1Cluster,

   385         right: c2Cluster,

   386         dist: distMatrix[c1][c2],

   387         centroid: calculateCentroid(c1Cluster.size, c1Cluster.centroid,

   388           c2Cluster.size, c2Cluster.centroid),

   389         size: c1Cluster.size + c2Cluster.size,

   390         depth: 1 + Math.max(c1Cluster.depth, c2Cluster.depth)

   391       };

   392       clusters[c1].splice(0, 0, newCluster);

   393       cSize[c1] += cSize[c2];

   394 

   395       // overwrite row c1 with respect to the linkage type

   396       for (j = 0; j < n; j++) {

   397         switch (linkage) {

   398           case "single":

   399             if (distMatrix[c1][j] > distMatrix[c2][j])

   400               distMatrix[j][c1] = distMatrix[c1][j] = distMatrix[c2][j];

   401             break;

   402           case "complete":

   403             if (distMatrix[c1][j] < distMatrix[c2][j])

   404               distMatrix[j][c1] = distMatrix[c1][j] = distMatrix[c2][j];

   405             break;

   406           case "average":

   407             distMatrix[j][c1] = distMatrix[c1][j] = (cSize[c1] * distMatrix[c1][j] + cSize[c2] * distMatrix[c2][j]) / (cSize[c1] + cSize[j]);

   408             break;

   409         }

   410       }

   411       distMatrix[c1][c1] = Infinity;

   412 

   413       // infinity ­out old row c2 and column c2

   414       for (i = 0; i < n; i++)

   415         distMatrix[i][c2] = distMatrix[c2][i] = Infinity;

   416 

   417       // update dmin and replace ones that previous pointed to c2 to point to c1

   418       for (j = 0; j < n; j++) {

   419         if (dMin[j] == c2) dMin[j] = c1;

   420         if (distMatrix[c1][j] < distMatrix[c1][dMin[c1]]) dMin[c1] = j;

   421       }

   422 

   423       // keep track of the last added cluster

   424       root = newCluster;

   425     }

   426 

   427     return root;

   428   }

   429 

   430   hcluster.distance = function(x) {

   431     if (!arguments.length) return distance;

   432     distance = x;

   433     return hcluster;

   434   };

   435 

   436   return hcluster;

   437 };

   438 

   439 function calculateCentroid(c1Size, c1Centroid, c2Size, c2Centroid) {

   440   var newCentroid = [],

   441       newSize = c1Size + c2Size,

   442       n = c1Centroid.length,

   443       i = -1;

   444   while (++i < n) {

   445     newCentroid[i] = (c1Size * c1Centroid[i] + c2Size * c2Centroid[i]) / newSize;

   446   }

   447   return newCentroid;

   448 }

   449 science.stats.iqr = function(x) {

   450   var quartiles = science.stats.quantiles(x, [.25, .75]);

   451   return quartiles[1] - quartiles[0];

   452 };

   453 // Based on org.apache.commons.math.analysis.interpolation.LoessInterpolator

   454 // from http://commons.apache.org/math/

   455 science.stats.loess = function() {    

   456   var bandwidth = .3,

   457       robustnessIters = 2,

   458       accuracy = 1e-12;

   459 

   460   function smooth(xval, yval, weights) {

   461     var n = xval.length,

   462         i;

   463 

   464     if (n !== yval.length) throw {error: "Mismatched array lengths"};

   465     if (n == 0) throw {error: "At least one point required."};

   466 

   467     if (arguments.length < 3) {

   468       weights = [];

   469       i = -1; while (++i < n) weights[i] = 1;

   470     }

   471 

   472     science_stats_loessFiniteReal(xval);

   473     science_stats_loessFiniteReal(yval);

   474     science_stats_loessFiniteReal(weights);

   475     science_stats_loessStrictlyIncreasing(xval);

   476 

   477     if (n == 1) return [yval[0]];

   478     if (n == 2) return [yval[0], yval[1]];

   479 

   480     var bandwidthInPoints = Math.floor(bandwidth * n);

   481 

   482     if (bandwidthInPoints < 2) throw {error: "Bandwidth too small."};

   483 

   484     var res = [],

   485         residuals = [],

   486         robustnessWeights = [];

   487 

   488     // Do an initial fit and 'robustnessIters' robustness iterations.

   489     // This is equivalent to doing 'robustnessIters+1' robustness iterations

   490     // starting with all robustness weights set to 1.

   491     i = -1; while (++i < n) {

   492       res[i] = 0;

   493       residuals[i] = 0;

   494       robustnessWeights[i] = 1;

   495     }

   496 

   497     var iter = -1;

   498     while (++iter <= robustnessIters) {

   499       var bandwidthInterval = [0, bandwidthInPoints - 1];

   500       // At each x, compute a local weighted linear regression

   501       var x;

   502       i = -1; while (++i < n) {

   503         x = xval[i];

   504 

   505         // Find out the interval of source points on which

   506         // a regression is to be made.

   507         if (i > 0) {

   508           science_stats_loessUpdateBandwidthInterval(xval, weights, i, bandwidthInterval);

   509         }

   510 

   511         var ileft = bandwidthInterval[0],

   512             iright = bandwidthInterval[1];

   513 

   514         // Compute the point of the bandwidth interval that is

   515         // farthest from x

   516         var edge = (xval[i] - xval[ileft]) > (xval[iright] - xval[i]) ? ileft : iright;

   517 

   518         // Compute a least-squares linear fit weighted by

   519         // the product of robustness weights and the tricube

   520         // weight function.

   521         // See http://en.wikipedia.org/wiki/Linear_regression

   522         // (section "Univariate linear case")

   523         // and http://en.wikipedia.org/wiki/Weighted_least_squares

   524         // (section "Weighted least squares")

   525         var sumWeights = 0,

   526             sumX = 0,

   527             sumXSquared = 0,

   528             sumY = 0,

   529             sumXY = 0,

   530             denom = Math.abs(1 / (xval[edge] - x));

   531 

   532         for (var k = ileft; k <= iright; ++k) {

   533           var xk   = xval[k],

   534               yk   = yval[k],

   535               dist = k < i ? x - xk : xk - x,

   536               w    = science_stats_loessTricube(dist * denom) * robustnessWeights[k] * weights[k],

   537               xkw  = xk * w;

   538           sumWeights += w;

   539           sumX += xkw;

   540           sumXSquared += xk * xkw;

   541           sumY += yk * w;

   542           sumXY += yk * xkw;

   543         }

   544 

   545         var meanX = sumX / sumWeights,

   546             meanY = sumY / sumWeights,

   547             meanXY = sumXY / sumWeights,

   548             meanXSquared = sumXSquared / sumWeights;

   549 

   550         var beta = (Math.sqrt(Math.abs(meanXSquared - meanX * meanX)) < accuracy)

   551             ? 0 : ((meanXY - meanX * meanY) / (meanXSquared - meanX * meanX));

   552 

   553         var alpha = meanY - beta * meanX;

   554 

   555         res[i] = beta * x + alpha;

   556         residuals[i] = Math.abs(yval[i] - res[i]);

   557       }

   558 

   559       // No need to recompute the robustness weights at the last

   560       // iteration, they won't be needed anymore

   561       if (iter === robustnessIters) {

   562         break;

   563       }

   564 

   565       // Recompute the robustness weights.

   566 

   567       // Find the median residual.

   568       var sortedResiduals = residuals.slice();

   569       sortedResiduals.sort();

   570       var medianResidual = sortedResiduals[Math.floor(n / 2)];

   571 

   572       if (Math.abs(medianResidual) < accuracy)

   573         break;

   574 

   575       var arg,

   576           w;

   577       i = -1; while (++i < n) {

   578         arg = residuals[i] / (6 * medianResidual);

   579         robustnessWeights[i] = (arg >= 1) ? 0 : ((w = 1 - arg * arg) * w);

   580       }

   581     }

   582 

   583     return res;

   584   }

   585 

   586   smooth.bandwidth = function(x) {

   587     if (!arguments.length) return x;

   588     bandwidth = x;

   589     return smooth;

   590   };

   591 

   592   smooth.robustnessIterations = function(x) {

   593     if (!arguments.length) return x;

   594     robustnessIters = x;

   595     return smooth;

   596   };

   597 

   598   smooth.accuracy = function(x) {

   599     if (!arguments.length) return x;

   600     accuracy = x;

   601     return smooth;

   602   };

   603 

   604   return smooth;

   605 };

   606 

   607 function science_stats_loessFiniteReal(values) {

   608   var n = values.length,

   609       i = -1;

   610 

   611   while (++i < n) if (!isFinite(values[i])) return false;

   612 

   613   return true;

   614 }

   615 

   616 function science_stats_loessStrictlyIncreasing(xval) {

   617   var n = xval.length,

   618       i = 0;

   619 

   620   while (++i < n) if (xval[i - 1] >= xval[i]) return false;

   621 

   622   return true;

   623 }

   624 

   625 // Compute the tricube weight function.

   626 // http://en.wikipedia.org/wiki/Local_regression#Weight_function

   627 function science_stats_loessTricube(x) {

   628   return (x = 1 - x * x * x) * x * x;

   629 }

   630 

   631 // Given an index interval into xval that embraces a certain number of

   632 // points closest to xval[i-1], update the interval so that it embraces

   633 // the same number of points closest to xval[i], ignoring zero weights.

   634 function science_stats_loessUpdateBandwidthInterval(

   635   xval, weights, i, bandwidthInterval) {

   636 

   637   var left = bandwidthInterval[0],

   638       right = bandwidthInterval[1];

   639 

   640   // The right edge should be adjusted if the next point to the right

   641   // is closer to xval[i] than the leftmost point of the current interval

   642   var nextRight = science_stats_loessNextNonzero(weights, right);

   643   if ((nextRight < xval.length) && (xval[nextRight] - xval[i]) < (xval[i] - xval[left])) {

   644     var nextLeft = science_stats_loessNextNonzero(weights, left);

   645     bandwidthInterval[0] = nextLeft;

   646     bandwidthInterval[1] = nextRight;

   647   }

   648 }

   649 

   650 function science_stats_loessNextNonzero(weights, i) {

   651   var j = i + 1;

   652   while (j < weights.length && weights[j] === 0) j++;

   653   return j;

   654 }

   655 // Welford's algorithm.

   656 science.stats.mean = function(x) {

   657   var n = x.length;

   658   if (n === 0) return NaN;

   659   var m = 0,

   660       i = -1;

   661   while (++i < n) m += (x[i] - m) / (i + 1);

   662   return m;

   663 };

   664 science.stats.median = function(x) {

   665   return science.stats.quantiles(x, [.5])[0];

   666 };

   667 science.stats.mode = function(x) {

   668   x = x.slice().sort(science.ascending);

   669   var mode,

   670       n = x.length,

   671       i = -1,

   672       l = i,

   673       last = null,

   674       max = 0,

   675       tmp,

   676       v;

   677   while (++i < n) {

   678     if ((v = x[i]) !== last) {

   679       if ((tmp = i - l) > max) {

   680         max = tmp;

   681         mode = last;

   682       }

   683       last = v;

   684       l = i;

   685     }

   686   }

   687   return mode;

   688 };

   689 // Uses R's quantile algorithm type=7.

   690 science.stats.quantiles = function(d, quantiles) {

   691   d = d.slice().sort(science.ascending);

   692   var n_1 = d.length - 1;

   693   return quantiles.map(function(q) {

   694     if (q === 0) return d[0];

   695     else if (q === 1) return d[n_1];

   696 

   697     var index = 1 + q * n_1,

   698         lo = Math.floor(index),

   699         h = index - lo,

   700         a = d[lo - 1];

   701 

   702     return h === 0 ? a : a + h * (d[lo] - a);

   703   });

   704 };

   705 // Unbiased estimate of a sample's variance.

   706 // Also known as the sample variance, where the denominator is n - 1.

   707 science.stats.variance = function(x) {

   708   var n = x.length;

   709   if (n < 1) return NaN;

   710   if (n === 1) return 0;

   711   var mean = science.stats.mean(x),

   712       i = -1,

   713       s = 0;

   714   while (++i < n) {

   715     var v = x[i] - mean;

   716     s += v * v;

   717   }

   718   return s / (n - 1);

   719 };

   720 })()