← Back to branch summary

~smagoun/whoopsie/whoopsie-lp1017637

« back to all changes in this revision

Viewing changes to backend/stats/static/js/d3/lib/science/science.stats.js

  • Committer: Evan Dandrea
  • Date: 2012-05-09 05:53:45 UTC
  • Revision ID: evan.dandrea@canonical.com-20120509055345-z2j41tmcbf4as5uf
The backend now lives in lp:daisy and the website (errors.ubuntu.com) now lives in lp:errors.

Show diffs side-by-side

added added

removed removed

Lines of Context:
backend/stats/static/js/d3/lib/science/science.stats.js
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
 
})()
 
 
b'\\ No newline at end of file'