Add unit tests for config.json and statistics.js
[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 + forwardWindowSize; 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 = values[0];
149         averages[0] = movingAverage;
150         for (var i = 1; i < values.length; i++) {
151             movingAverage = smoothingFactor * values[i] + (1 - smoothingFactor) * movingAverage;
152             averages[i] = movingAverage;
153         }
154         return averages;
155     }
156
157     // The return value is the starting indices of each segment.
158     this.segmentTimeSeriesGreedyWithStudentsTTest = function (values, minLength) {
159         if (values.length < 2)
160             return [0];
161         var segments = new Array;
162         recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, 0, values.length, minLength, segments);
163         segments.push(values.length);
164         return segments;
165     }
166
167     this.debuggingSegmentation = false;
168     this.segmentTimeSeriesByMaximizingSchwarzCriterion = function (values) {
169         if (values.length < 2)
170             return [0];
171
172         // Split the time series into grids since splitIntoSegmentsUntilGoodEnough is O(n^2).
173         var gridLength = 500;
174         var totalSegmentation = [0];
175         for (var gridCount = 0; gridCount < Math.ceil(values.length / gridLength); gridCount++) {
176             var gridValues = values.slice(gridCount * gridLength, (gridCount + 1) * gridLength);
177             var segmentation = splitIntoSegmentsUntilGoodEnough(gridValues);
178
179             if (Statistics.debuggingSegmentation)
180                 console.log('grid=' + gridCount, segmentation);
181
182             for (var i = 1; i < segmentation.length - 1; i++)
183                 totalSegmentation.push(gridCount * gridLength + segmentation[i]);
184         }
185
186         if (Statistics.debuggingSegmentation)
187             console.log('Final Segmentation', totalSegmentation);
188
189         totalSegmentation.push(values.length);
190
191         return totalSegmentation;
192     }
193
194     function recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex, length, minLength, segments) {
195         var tMax = 0;
196         var argTMax = null;
197         for (var i = 1; i < length - 1; i++) {
198             var firstLength = i;
199             var secondLength = length - i;
200             if (firstLength < minLength || secondLength < minLength)
201                 continue;
202             var result = Statistics.computeWelchsT(values, startIndex, firstLength, values, startIndex + i, secondLength, 0.9);
203             if (result.significantlyDifferent && result.t > tMax) {
204                 tMax = result.t;
205                 argTMax = i;
206             }
207         }
208         if (!tMax) {
209             segments.push(startIndex);
210             return;
211         }
212         recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex, argTMax, minLength, segments);
213         recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex + argTMax, length - argTMax, minLength, segments);
214     }
215
216     var tDistributionByOneSidedProbability = {
217         0.9: [
218             3.077684, 1.885618, 1.637744, 1.533206, 1.475884, 1.439756, 1.414924, 1.396815, 1.383029, 1.372184,
219             1.363430, 1.356217, 1.350171, 1.345030, 1.340606, 1.336757, 1.333379, 1.330391, 1.327728, 1.325341,
220             1.323188, 1.321237, 1.319460, 1.317836, 1.316345, 1.314972, 1.313703, 1.312527, 1.311434, 1.310415,
221             1.309464, 1.308573, 1.307737, 1.306952, 1.306212, 1.305514, 1.304854, 1.304230, 1.303639, 1.303077,
222             1.302543, 1.302035, 1.301552, 1.301090, 1.300649, 1.300228, 1.299825, 1.299439, 1.299069, 1.298714,
223
224             1.298373, 1.298045, 1.297730, 1.297426, 1.297134, 1.296853, 1.296581, 1.296319, 1.296066, 1.295821,
225             1.295585, 1.295356, 1.295134, 1.294920, 1.294712, 1.294511, 1.294315, 1.294126, 1.293942, 1.293763,
226             1.293589, 1.293421, 1.293256, 1.293097, 1.292941, 1.292790, 1.292643, 1.292500, 1.292360, 1.292224,
227             1.292091, 1.291961, 1.291835, 1.291711, 1.291591, 1.291473, 1.291358, 1.291246, 1.291136, 1.291029,
228             1.290924, 1.290821, 1.290721, 1.290623, 1.290527, 1.290432, 1.290340, 1.290250, 1.290161, 1.290075],
229         0.95: [
230             6.313752, 2.919986, 2.353363, 2.131847, 2.015048, 1.943180, 1.894579, 1.859548, 1.833113, 1.812461,
231             1.795885, 1.782288, 1.770933, 1.761310, 1.753050, 1.745884, 1.739607, 1.734064, 1.729133, 1.724718,
232             1.720743, 1.717144, 1.713872, 1.710882, 1.708141, 1.705618, 1.703288, 1.701131, 1.699127, 1.697261,
233             1.695519, 1.693889, 1.692360, 1.690924, 1.689572, 1.688298, 1.687094, 1.685954, 1.684875, 1.683851,
234             1.682878, 1.681952, 1.681071, 1.680230, 1.679427, 1.678660, 1.677927, 1.677224, 1.676551, 1.675905,
235
236             1.675285, 1.674689, 1.674116, 1.673565, 1.673034, 1.672522, 1.672029, 1.671553, 1.671093, 1.670649,
237             1.670219, 1.669804, 1.669402, 1.669013, 1.668636, 1.668271, 1.667916, 1.667572, 1.667239, 1.666914,
238             1.666600, 1.666294, 1.665996, 1.665707, 1.665425, 1.665151, 1.664885, 1.664625, 1.664371, 1.664125,
239             1.663884, 1.663649, 1.663420, 1.663197, 1.662978, 1.662765, 1.662557, 1.662354, 1.662155, 1.661961,
240             1.661771, 1.661585, 1.661404, 1.661226, 1.661052, 1.660881, 1.660715, 1.660551, 1.660391, 1.660234],
241         0.975: [
242             12.706205, 4.302653, 3.182446, 2.776445, 2.570582, 2.446912, 2.364624, 2.306004, 2.262157, 2.228139,
243             2.200985, 2.178813, 2.160369, 2.144787, 2.131450, 2.119905, 2.109816, 2.100922, 2.093024, 2.085963,
244             2.079614, 2.073873, 2.068658, 2.063899, 2.059539, 2.055529, 2.051831, 2.048407, 2.045230, 2.042272,
245             2.039513, 2.036933, 2.034515, 2.032245, 2.030108, 2.028094, 2.026192, 2.024394, 2.022691, 2.021075,
246             2.019541, 2.018082, 2.016692, 2.015368, 2.014103, 2.012896, 2.011741, 2.010635, 2.009575, 2.008559,
247
248             2.007584, 2.006647, 2.005746, 2.004879, 2.004045, 2.003241, 2.002465, 2.001717, 2.000995, 2.000298,
249             1.999624, 1.998972, 1.998341, 1.997730, 1.997138, 1.996564, 1.996008, 1.995469, 1.994945, 1.994437,
250             1.993943, 1.993464, 1.992997, 1.992543, 1.992102, 1.991673, 1.991254, 1.990847, 1.990450, 1.990063,
251             1.989686, 1.989319, 1.988960, 1.988610, 1.988268, 1.987934, 1.987608, 1.987290, 1.986979, 1.986675,
252             1.986377, 1.986086, 1.985802, 1.985523, 1.985251, 1.984984, 1.984723, 1.984467, 1.984217, 1.983972],
253         0.99: [
254             31.820516, 6.964557, 4.540703, 3.746947, 3.364930, 3.142668, 2.997952, 2.896459, 2.821438, 2.763769,
255             2.718079, 2.680998, 2.650309, 2.624494, 2.602480, 2.583487, 2.566934, 2.552380, 2.539483, 2.527977,
256             2.517648, 2.508325, 2.499867, 2.492159, 2.485107, 2.478630, 2.472660, 2.467140, 2.462021, 2.457262,
257             2.452824, 2.448678, 2.444794, 2.441150, 2.437723, 2.434494, 2.431447, 2.428568, 2.425841, 2.423257,
258             2.420803, 2.418470, 2.416250, 2.414134, 2.412116, 2.410188, 2.408345, 2.406581, 2.404892, 2.403272,
259
260             2.401718, 2.400225, 2.398790, 2.397410, 2.396081, 2.394801, 2.393568, 2.392377, 2.391229, 2.390119,
261             2.389047, 2.388011, 2.387008, 2.386037, 2.385097, 2.384186, 2.383302, 2.382446, 2.381615, 2.380807,
262             2.380024, 2.379262, 2.378522, 2.377802, 2.377102, 2.376420, 2.375757, 2.375111, 2.374482, 2.373868,
263             2.373270, 2.372687, 2.372119, 2.371564, 2.371022, 2.370493, 2.369977, 2.369472, 2.368979, 2.368497,
264             2.368026, 2.367566, 2.367115, 2.366674, 2.366243, 2.365821, 2.365407, 2.365002, 2.364606, 2.364217]
265     };
266     function oneSidedToTwoSidedProbability(probability) { return 2 * probability - 1; }
267     function twoSidedToOneSidedProbability(probability) { return (1 - (1 - probability) / 2); }
268
269     function splitIntoSegmentsUntilGoodEnough(values) {
270         if (values.length < 2)
271             return [0, values.length];
272
273         var matrix = new SampleVarianceUpperTriangularMatrix(values);
274
275         var SchwarzCriterionBeta = Math.log1p(values.length - 1) / values.length;
276
277         var BirgeAndMassartC = 2.5; // Suggested by the authors.
278         var BirgeAndMassartPenalization = function (segmentCount) {
279             return segmentCount * (1 + BirgeAndMassartC * Math.log1p(values.length / segmentCount - 1));
280         }
281
282         var segmentation;
283         var minTotalCost = Infinity;
284         var maxK = Math.min(50, values.length);
285
286         for (var k = 1; k < maxK; k++) {
287             var start = Date.now();
288             var result = findOptimalSegmentation(values, matrix, k);
289             var cost = result.cost / values.length;
290             var penalty = SchwarzCriterionBeta * BirgeAndMassartPenalization(k);
291             if (cost + penalty < minTotalCost) {
292                 minTotalCost = cost + penalty;
293                 segmentation = result.segmentation;
294             } else
295                 maxK = Math.min(maxK, k + 3);
296             if (Statistics.debuggingSegmentation)
297                 console.log('splitIntoSegmentsUntilGoodEnough', k, Date.now() - start, cost + penalty);
298         }
299
300         return segmentation;
301     }
302
303     function findOptimalSegmentation(values, costMatrix, segmentCount) {
304         // Dynamic programming. cost[i][k] = The cost to segmenting values up to i into k segments.
305         var cost = new Array(values.length);
306         for (var segmentEnd = 0; segmentEnd < values.length; segmentEnd++)
307             cost[segmentEnd] = new Float32Array(segmentCount + 1);
308
309         // previousNode[i][k] = The start of the last segment in an optimal segmentation that ends at i with k segments.
310         var previousNode = new Array(values.length);
311         for (var i = 0; i < values.length; i++)
312             previousNode[i] = new Array(segmentCount + 1);
313
314         cost[0] = [0]; // The cost of segmenting single value is always 0.
315         previousNode[0] = [-1];
316         for (var segmentStart = 0; segmentStart < values.length; segmentStart++) {
317             var costOfOptimalSegmentationThatEndAtCurrentStart = cost[segmentStart];
318             for (var k = 0; k < segmentCount; k++) {
319                 var noSegmentationOfLenghtKEndsAtCurrentStart = previousNode[segmentStart][k] === undefined;
320                 if (noSegmentationOfLenghtKEndsAtCurrentStart)
321                     continue;
322                 for (var segmentEnd = segmentStart + 1; segmentEnd < values.length; segmentEnd++) {
323                     var costOfOptimalSegmentationOfLengthK = costOfOptimalSegmentationThatEndAtCurrentStart[k];
324                     var costOfCurrentSegment = costMatrix.costBetween(segmentStart, segmentEnd);
325                     var totalCost = costOfOptimalSegmentationOfLengthK + costOfCurrentSegment;
326                     if (previousNode[segmentEnd][k + 1] === undefined || totalCost < cost[segmentEnd][k + 1]) {
327                         cost[segmentEnd][k + 1] = totalCost;
328                         previousNode[segmentEnd][k + 1] = segmentStart;
329                     }
330                 }
331             }
332         }
333
334         if (Statistics.debuggingSegmentation) {
335             console.log('findOptimalSegmentation with', segmentCount, 'segments');
336             for (var end = 0; end < values.length; end++) {
337                 for (var k = 0; k <= segmentCount; k++) {
338                     var start = previousNode[end][k];
339                     if (start === undefined)
340                         continue;
341                     console.log(`C(segment=[${start}, ${end + 1}], segmentCount=${k})=${cost[end][k]}`);
342                 }
343             }
344         }
345
346         var segmentEnd = values.length - 1;
347         var segmentation = new Array(segmentCount + 1);
348         segmentation[segmentCount] = values.length;
349         for (var k = segmentCount; k > 0; k--) {
350             segmentEnd = previousNode[segmentEnd][k];
351             segmentation[k - 1] = segmentEnd;
352         }
353         var costOfOptimalSegmentation = cost[values.length - 1][segmentCount];
354
355         if (Statistics.debuggingSegmentation)
356             console.log('Optimal segmentation:', segmentation, 'with cost =', costOfOptimalSegmentation);
357
358         return {segmentation: segmentation, cost: costOfOptimalSegmentation};
359     }
360
361     function SampleVarianceUpperTriangularMatrix(values) {
362         // The cost of segment (i, j].
363         var costMatrix = new Array(values.length - 1);
364         for (var i = 0; i < values.length - 1; i++) {
365             var remainingValueCount = values.length - i - 1;
366             costMatrix[i] = new Float32Array(remainingValueCount);
367             var sum = values[i];
368             var squareSum = sum * sum;
369             costMatrix[i][0] = 0;
370             for (var j = i + 1; j < values.length; j++) {
371                 var currentValue = values[j];
372                 sum += currentValue;
373                 squareSum += currentValue * currentValue;
374                 var sampleSize = j - i + 1;
375                 var stdev = Statistics.sampleStandardDeviation(sampleSize, sum, squareSum);
376                 costMatrix[i][j - i - 1] = sampleSize * Math.log1p(stdev * stdev - 1);
377             }
378         }
379         this.costMatrix = costMatrix;
380     }
381
382     SampleVarianceUpperTriangularMatrix.prototype.costBetween = function(from, to) {
383         if (from >= this.costMatrix.length || from == to)
384             return 0; // The cost of the segment that starts at the last data point is 0.
385         return this.costMatrix[from][to - from - 1];
386     }
387
388     this.debuggingTestingRangeNomination = false;
389
390 })();
391
392 if (typeof module != 'undefined')
393     module.exports = Statistics;