9ef1e74ed966fb0e1e0587de4e2dd4ba249c2963
[WebKit.git] / Websites / perf.webkit.org / public / shared / statistics.js
1 var Statistics = new (function () {
2
3     this.min = function (values) {
4         return Math.min.apply(Math, values);
5     }
6
7     this.max = function (values) {
8         return Math.max.apply(Math, values);
9     }
10
11     this.sum = function (values) {
12         return values.length ? values.reduce(function (a, b) { return a + b; }) : 0;
13     }
14
15     this.squareSum = function (values) {
16         return values.length ? values.reduce(function (sum, value) { return sum + value * value;}, 0) : 0;
17     }
18
19     // With sum and sum of squares, we can compute the sample standard deviation in O(1).
20     // See https://rniwa.com/2012-11-10/sample-standard-deviation-in-terms-of-sum-and-square-sum-of-samples/
21     this.sampleStandardDeviation = function (numberOfSamples, sum, squareSum) {
22         if (numberOfSamples < 2)
23             return 0;
24         return Math.sqrt(squareSum / (numberOfSamples - 1) - sum * sum / (numberOfSamples - 1) / numberOfSamples);
25     }
26
27     this.supportedConfidenceIntervalProbabilities = function () {
28         var supportedProbabilities = [];
29         for (var probability in tDistributionByOneSidedProbability)
30             supportedProbabilities.push(oneSidedToTwoSidedProbability(probability).toFixed(2));
31         return supportedProbabilities
32     }
33
34     // Computes the delta d s.t. (mean - d, mean + d) is the confidence interval with the specified probability in O(1).
35     this.confidenceIntervalDelta = function (probability, numberOfSamples, sum, squareSum) {
36         var oneSidedProbability = twoSidedToOneSidedProbability(probability);
37         if (!(oneSidedProbability in tDistributionByOneSidedProbability)) {
38             throw 'We only support ' + this.supportedConfidenceIntervalProbabilities().map(function (probability)
39             { return probability * 100 + '%'; } ).join(', ') + ' confidence intervals.';
40         }
41         if (numberOfSamples - 2 < 0)
42             return NaN;
43         var deltas = tDistributionByOneSidedProbability[oneSidedProbability];
44         var degreesOfFreedom = numberOfSamples - 1;
45         if (degreesOfFreedom > deltas.length)
46             throw 'We only support up to ' + deltas.length + ' degrees of freedom';
47
48         // d = c * S/sqrt(numberOfSamples) where c ~ t-distribution(degreesOfFreedom) and S is the sample standard deviation.
49         return deltas[degreesOfFreedom - 1] * this.sampleStandardDeviation(numberOfSamples, sum, squareSum) / Math.sqrt(numberOfSamples);
50     }
51
52     this.confidenceInterval = function (values, probability) {
53         var sum = this.sum(values);
54         var mean = sum / values.length;
55         var delta = this.confidenceIntervalDelta(probability || 0.95, values.length, sum, this.squareSum(values));
56         return [mean - delta, mean + delta];
57     }
58
59     // Welch's t-test (http://en.wikipedia.org/wiki/Welch%27s_t_test)
60     this.testWelchsT = function (values1, values2, probability) {
61         return this.computeWelchsT(values1, 0, values1.length, values2, 0, values2.length, probability).significantlyDifferent;
62     }
63
64     this.probabilityRangeForWelchsT = function (values1, values2) {
65         var result = this.computeWelchsT(values1, 0, values1.length, values2, 0, values2.length);
66         if (isNaN(result.t) || isNaN(result.degreesOfFreedom))
67             return {t: NaN, degreesOfFreedom:NaN, range: [null, null]};
68
69         var lowerBound = null;
70         var upperBound = null;
71         for (var probability in tDistributionByOneSidedProbability) {
72             var twoSidedProbability = oneSidedToTwoSidedProbability(probability);
73             if (result.t > tDistributionByOneSidedProbability[probability][Math.round(result.degreesOfFreedom - 1)])
74                 lowerBound = twoSidedProbability;
75             else if (lowerBound) {
76                 upperBound = twoSidedProbability;
77                 break;
78             }
79         }
80         return {t: result.t, degreesOfFreedom: result.degreesOfFreedom, range: [lowerBound, upperBound]};
81     }
82
83     this.computeWelchsT = function (values1, startIndex1, length1, values2, startIndex2, length2, probability) {
84         var stat1 = sampleMeanAndVarianceForValues(values1, startIndex1, length1);
85         var stat2 = sampleMeanAndVarianceForValues(values2, startIndex2, length2);
86         var sumOfSampleVarianceOverSampleSize = stat1.variance / stat1.size + stat2.variance / stat2.size;
87         var t = Math.abs((stat1.mean - stat2.mean) / Math.sqrt(sumOfSampleVarianceOverSampleSize));
88
89         // http://en.wikipedia.org/wiki/Welch–Satterthwaite_equation
90         var degreesOfFreedom = sumOfSampleVarianceOverSampleSize * sumOfSampleVarianceOverSampleSize
91             / (stat1.variance * stat1.variance / stat1.size / stat1.size / stat1.degreesOfFreedom
92                 + stat2.variance * stat2.variance / stat2.size / stat2.size / stat2.degreesOfFreedom);
93         var minT = tDistributionByOneSidedProbability[twoSidedToOneSidedProbability(probability || 0.8)][Math.round(degreesOfFreedom - 1)];
94         return {
95             t: t,
96             degreesOfFreedom: degreesOfFreedom,
97             significantlyDifferent: t > minT,
98         };
99     }
100
101     function sampleMeanAndVarianceForValues(values, startIndex, length) {
102         var sum = 0;
103         for (var i = 0; i < length; i++)
104             sum += values[startIndex + i];
105         var squareSum = 0;
106         for (var i = 0; i < length; i++)
107             squareSum += values[startIndex + i] * values[startIndex + i];
108         var sampleMean = sum / length;
109         // FIXME: Maybe we should be using the biased sample variance.
110         var unbiasedSampleVariance = (squareSum - sum * sum / length) / (length - 1);
111         return {
112             mean: sampleMean,
113             variance: unbiasedSampleVariance,
114             size: length,
115             degreesOfFreedom: length - 1,
116         }
117     }
118
119     this.movingAverage = function (values, backwardWindowSize, forwardWindowSize) {
120         var averages = new Array(values.length);
121         // We use naive O(n^2) algorithm for simplicy as well as to avoid accumulating round-off errors.
122         for (var i = 0; i < values.length; i++) {
123             var sum = 0;
124             var count = 0;
125             for (var j = i - backwardWindowSize; j < i + backwardWindowSize; j++) {
126                 if (j >= 0 && j < values.length) {
127                     sum += values[j];
128                     count++;
129                 }
130             }
131             averages[i] = sum / count;
132         }
133         return averages;
134     }
135
136     this.cumulativeMovingAverage = function (values) {
137         var averages = new Array(values.length);
138         var sum = 0;
139         for (var i = 0; i < values.length; i++) {
140             sum += values[i];
141             averages[i] = sum / (i + 1);
142         }
143         return averages;
144     }
145
146     this.exponentialMovingAverage = function (values, smoothingFactor) {
147         var averages = new Array(values.length);
148         var movingAverage = 0;
149         for (var i = 0; i < values.length; i++) {
150             movingAverage = smoothingFactor * values[i] + (1 - smoothingFactor) * movingAverage;
151             averages[i] = movingAverage;
152         }
153         return averages;
154     }
155
156     // The return value is the starting indices of each segment.
157     this.segmentTimeSeriesGreedyWithStudentsTTest = function (values, minLength) {
158         if (values.length < 2)
159             return [0];
160         var segments = new Array;
161         recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, 0, values.length, minLength, segments);
162         segments.push(values.length);
163         return segments;
164     }
165
166     this.debuggingSegmentation = false;
167     this.segmentTimeSeriesByMaximizingSchwarzCriterion = function (values) {
168         if (values.length < 2)
169             return [0];
170
171         // Split the time series into grids since splitIntoSegmentsUntilGoodEnough is O(n^2).
172         var gridLength = 500;
173         var totalSegmentation = [0];
174         for (var gridCount = 0; gridCount < Math.ceil(values.length / gridLength); gridCount++) {
175             var gridValues = values.slice(gridCount * gridLength, (gridCount + 1) * gridLength);
176             var segmentation = splitIntoSegmentsUntilGoodEnough(gridValues);
177
178             if (Statistics.debuggingSegmentation)
179                 console.log('grid=' + gridCount, segmentation);
180
181             for (var i = 1; i < segmentation.length - 1; i++)
182                 totalSegmentation.push(gridCount * gridLength + segmentation[i]);
183         }
184
185         if (Statistics.debuggingSegmentation)
186             console.log('Final Segmentation', totalSegmentation);
187
188         totalSegmentation.push(values.length);
189
190         return totalSegmentation;
191     }
192
193     function recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex, length, minLength, segments) {
194         var tMax = 0;
195         var argTMax = null;
196         for (var i = 1; i < length - 1; i++) {
197             var firstLength = i;
198             var secondLength = length - i;
199             if (firstLength < minLength || secondLength < minLength)
200                 continue;
201             var result = Statistics.computeWelchsT(values, startIndex, firstLength, values, startIndex + i, secondLength, 0.9);
202             if (result.significantlyDifferent && result.t > tMax) {
203                 tMax = result.t;
204                 argTMax = i;
205             }
206         }
207         if (!tMax) {
208             segments.push(startIndex);
209             return;
210         }
211         recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex, argTMax, minLength, segments);
212         recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex + argTMax, length - argTMax, minLength, segments);
213     }
214
215     var tDistributionByOneSidedProbability = {
216         0.9: [
217             3.077684, 1.885618, 1.637744, 1.533206, 1.475884, 1.439756, 1.414924, 1.396815, 1.383029, 1.372184,
218             1.363430, 1.356217, 1.350171, 1.345030, 1.340606, 1.336757, 1.333379, 1.330391, 1.327728, 1.325341,
219             1.323188, 1.321237, 1.319460, 1.317836, 1.316345, 1.314972, 1.313703, 1.312527, 1.311434, 1.310415,
220             1.309464, 1.308573, 1.307737, 1.306952, 1.306212, 1.305514, 1.304854, 1.304230, 1.303639, 1.303077,
221             1.302543, 1.302035, 1.301552, 1.301090, 1.300649, 1.300228, 1.299825, 1.299439, 1.299069, 1.298714,
222
223             1.298373, 1.298045, 1.297730, 1.297426, 1.297134, 1.296853, 1.296581, 1.296319, 1.296066, 1.295821,
224             1.295585, 1.295356, 1.295134, 1.294920, 1.294712, 1.294511, 1.294315, 1.294126, 1.293942, 1.293763,
225             1.293589, 1.293421, 1.293256, 1.293097, 1.292941, 1.292790, 1.292643, 1.292500, 1.292360, 1.292224,
226             1.292091, 1.291961, 1.291835, 1.291711, 1.291591, 1.291473, 1.291358, 1.291246, 1.291136, 1.291029,
227             1.290924, 1.290821, 1.290721, 1.290623, 1.290527, 1.290432, 1.290340, 1.290250, 1.290161, 1.290075],
228         0.95: [
229             6.313752, 2.919986, 2.353363, 2.131847, 2.015048, 1.943180, 1.894579, 1.859548, 1.833113, 1.812461,
230             1.795885, 1.782288, 1.770933, 1.761310, 1.753050, 1.745884, 1.739607, 1.734064, 1.729133, 1.724718,
231             1.720743, 1.717144, 1.713872, 1.710882, 1.708141, 1.705618, 1.703288, 1.701131, 1.699127, 1.697261,
232             1.695519, 1.693889, 1.692360, 1.690924, 1.689572, 1.688298, 1.687094, 1.685954, 1.684875, 1.683851,
233             1.682878, 1.681952, 1.681071, 1.680230, 1.679427, 1.678660, 1.677927, 1.677224, 1.676551, 1.675905,
234
235             1.675285, 1.674689, 1.674116, 1.673565, 1.673034, 1.672522, 1.672029, 1.671553, 1.671093, 1.670649,
236             1.670219, 1.669804, 1.669402, 1.669013, 1.668636, 1.668271, 1.667916, 1.667572, 1.667239, 1.666914,
237             1.666600, 1.666294, 1.665996, 1.665707, 1.665425, 1.665151, 1.664885, 1.664625, 1.664371, 1.664125,
238             1.663884, 1.663649, 1.663420, 1.663197, 1.662978, 1.662765, 1.662557, 1.662354, 1.662155, 1.661961,
239             1.661771, 1.661585, 1.661404, 1.661226, 1.661052, 1.660881, 1.660715, 1.660551, 1.660391, 1.660234],
240         0.975: [
241             12.706205, 4.302653, 3.182446, 2.776445, 2.570582, 2.446912, 2.364624, 2.306004, 2.262157, 2.228139,
242             2.200985, 2.178813, 2.160369, 2.144787, 2.131450, 2.119905, 2.109816, 2.100922, 2.093024, 2.085963,
243             2.079614, 2.073873, 2.068658, 2.063899, 2.059539, 2.055529, 2.051831, 2.048407, 2.045230, 2.042272,
244             2.039513, 2.036933, 2.034515, 2.032245, 2.030108, 2.028094, 2.026192, 2.024394, 2.022691, 2.021075,
245             2.019541, 2.018082, 2.016692, 2.015368, 2.014103, 2.012896, 2.011741, 2.010635, 2.009575, 2.008559,
246
247             2.007584, 2.006647, 2.005746, 2.004879, 2.004045, 2.003241, 2.002465, 2.001717, 2.000995, 2.000298,
248             1.999624, 1.998972, 1.998341, 1.997730, 1.997138, 1.996564, 1.996008, 1.995469, 1.994945, 1.994437,
249             1.993943, 1.993464, 1.992997, 1.992543, 1.992102, 1.991673, 1.991254, 1.990847, 1.990450, 1.990063,
250             1.989686, 1.989319, 1.988960, 1.988610, 1.988268, 1.987934, 1.987608, 1.987290, 1.986979, 1.986675,
251             1.986377, 1.986086, 1.985802, 1.985523, 1.985251, 1.984984, 1.984723, 1.984467, 1.984217, 1.983972],
252         0.99: [
253             31.820516, 6.964557, 4.540703, 3.746947, 3.364930, 3.142668, 2.997952, 2.896459, 2.821438, 2.763769,
254             2.718079, 2.680998, 2.650309, 2.624494, 2.602480, 2.583487, 2.566934, 2.552380, 2.539483, 2.527977,
255             2.517648, 2.508325, 2.499867, 2.492159, 2.485107, 2.478630, 2.472660, 2.467140, 2.462021, 2.457262,
256             2.452824, 2.448678, 2.444794, 2.441150, 2.437723, 2.434494, 2.431447, 2.428568, 2.425841, 2.423257,
257             2.420803, 2.418470, 2.416250, 2.414134, 2.412116, 2.410188, 2.408345, 2.406581, 2.404892, 2.403272,
258
259             2.401718, 2.400225, 2.398790, 2.397410, 2.396081, 2.394801, 2.393568, 2.392377, 2.391229, 2.390119,
260             2.389047, 2.388011, 2.387008, 2.386037, 2.385097, 2.384186, 2.383302, 2.382446, 2.381615, 2.380807,
261             2.380024, 2.379262, 2.378522, 2.377802, 2.377102, 2.376420, 2.375757, 2.375111, 2.374482, 2.373868,
262             2.373270, 2.372687, 2.372119, 2.371564, 2.371022, 2.370493, 2.369977, 2.369472, 2.368979, 2.368497,
263             2.368026, 2.367566, 2.367115, 2.366674, 2.366243, 2.365821, 2.365407, 2.365002, 2.364606, 2.364217]
264     };
265     function oneSidedToTwoSidedProbability(probability) { return 2 * probability - 1; }
266     function twoSidedToOneSidedProbability(probability) { return (1 - (1 - probability) / 2); }
267
268     function splitIntoSegmentsUntilGoodEnough(values) {
269         if (values.length < 2)
270             return [0, values.length];
271
272         var matrix = new SampleVarianceUpperTriangularMatrix(values);
273
274         var SchwarzCriterionBeta = Math.log1p(values.length - 1) / values.length;
275
276         var BirgeAndMassartC = 2.5; // Suggested by the authors.
277         var BirgeAndMassartPenalization = function (segmentCount) {
278             return segmentCount * (1 + BirgeAndMassartC * Math.log1p(values.length / segmentCount - 1));
279         }
280
281         var segmentation;
282         var minTotalCost = Infinity;
283         var maxK = 50;
284
285         for (var k = 1; k < maxK; k++) {
286             var start = Date.now();
287             var result = findOptimalSegmentation(values, matrix, k);
288             var cost = result.cost / values.length;
289             var penalty = SchwarzCriterionBeta * BirgeAndMassartPenalization(k);
290             if (cost + penalty < minTotalCost) {
291                 minTotalCost = cost + penalty;
292                 segmentation = result.segmentation;
293             } else
294                 maxK = Math.min(maxK, k + 3);
295             if (Statistics.debuggingSegmentation)
296                 console.log('splitIntoSegmentsUntilGoodEnough', k, Date.now() - start, cost + penalty);
297         }
298
299         return segmentation;
300     }
301
302     function findOptimalSegmentation(values, costMatrix, segmentCount) {
303         // Dynamic programming. cost[i][k] = The cost to segmenting values up to i into k segments.
304         var cost = new Array(values.length);
305         for (var i = 0; i < values.length; i++) {
306             cost[i] = new Float32Array(segmentCount + 1);
307         }
308
309         var previousNode = new Array(values.length);
310         for (var i = 0; i < values.length; i++)
311             previousNode[i] = new Array(segmentCount + 1);
312
313         cost[0] = [0]; // The cost of segmenting single value is always 0.
314         previousNode[0] = [-1];
315         for (var segmentStart = 0; segmentStart < values.length; segmentStart++) {
316             var costBySegment = cost[segmentStart];
317             for (var count = 0; count < segmentCount; count++) {
318                 if (previousNode[segmentStart][count] === undefined)
319                     continue;
320                 for (var segmentEnd = segmentStart + 1; segmentEnd < values.length; segmentEnd++) {
321                     var newCost = costBySegment[count] + costMatrix.costBetween(segmentStart, segmentEnd);
322                     if (previousNode[segmentEnd][count + 1] === undefined || newCost < cost[segmentEnd][count + 1]) {
323                         cost[segmentEnd][count + 1] = newCost;
324                         previousNode[segmentEnd][count + 1] = segmentStart;
325                     }
326                 }
327             }
328         }
329
330         if (Statistics.debuggingSegmentation) {
331             console.log('findOptimalSegmentation with k=', segmentCount);
332             for (var i = 0; i < cost.length; i++) {
333                 var t = cost[i];
334                 var s = '';
335                 for (var j = 0; j < t.length; j++) {
336                     var p = previousNode[i][j];
337                     s += '(k=' + j;
338                     if (p !== undefined)
339                         s += ' c=' + t[j] + ' p=' + p
340                     s += ')';
341                 }
342                 console.log(i, values[i], s);
343             }
344         }
345
346         var currentIndex = values.length - 1;
347         var segmentation = new Array(segmentCount);
348         segmentation[0] = values.length;
349         for (var i = 0; i < segmentCount; i++) {
350             currentIndex = previousNode[currentIndex][segmentCount - i];
351             segmentation[i + 1] = currentIndex;
352         }
353
354         return {segmentation: segmentation.reverse(), cost: cost[values.length - 1][segmentCount]};
355     }
356
357     function SampleVarianceUpperTriangularMatrix(values) {
358         // The cost of segment (i, j].
359         var costMatrix = new Array(values.length - 1);
360         for (var i = 0; i < values.length - 1; i++) {
361             var remainingValueCount = values.length - i - 1;
362             costMatrix[i] = new Float32Array(remainingValueCount);
363             var sum = values[i];
364             var squareSum = sum * sum;
365             costMatrix[i][0] = 0;
366             for (var j = i + 1; j < values.length; j++) {
367                 var currentValue = values[j];
368                 sum += currentValue;
369                 squareSum += currentValue * currentValue;
370                 var sampleSize = j - i + 1;
371                 var stdev = Statistics.sampleStandardDeviation(sampleSize, sum, squareSum);
372                 costMatrix[i][j - i - 1] = sampleSize * Math.log1p(stdev * stdev - 1);
373             }
374         }
375         this.costMatrix = costMatrix;
376     }
377
378     SampleVarianceUpperTriangularMatrix.prototype.costBetween = function(from, to) {
379         if (from >= this.costMatrix.length || from == to)
380             return 0; // The cost of the segment that starts at the last data point is 0.
381         return this.costMatrix[from][to - from - 1];
382     }
383
384     this.debuggingTestingRangeNomination = false;
385
386 })();
387
388 if (typeof module != 'undefined') {
389     for (var key in Statistics)
390         module.exports[key] = Statistics[key];
391 }