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