Add time series segmentation algorithms as moving averages
authorrniwa@webkit.org <rniwa@webkit.org@268f45cc-cd09-0410-ab3c-d52691b4dbfc>
Fri, 3 Apr 2015 20:46:11 +0000 (20:46 +0000)
committerrniwa@webkit.org <rniwa@webkit.org@268f45cc-cd09-0410-ab3c-d52691b4dbfc>
Fri, 3 Apr 2015 20:46:11 +0000 (20:46 +0000)
https://bugs.webkit.org/show_bug.cgi?id=143362

Reviewed by Chris Dumez.

This patch implements two preliminary time series segmentation algorithms as moving averages.

Recursive t-test: Compute Welch's t-statistic at each point in a given segment of the time series.
If Welch's t-test implicates a statistically significance difference, then split the segment into two
sub segments with the maximum t-statistic (i.e. the point at which if split would yield the highest
probability that two segments do not share the same "underlying" mean in classical / frequentist sense).
We repeat this process recursively. See [1] for the evaluation of this particular algorithm.

Schwarz criterion: Use Schwarz or Bayesian information criterion to heuristically find the optimal
segmentation. Intuitively, the problem of finding the best segmentation comes down to minimizing the
residual sum of squares in each segment as in linear regressions. That is, for a given segment with
values y_1 through y_n with mean y_avg, we want to minimize the sum of (y_i - y_avg)^2 over i = 1
through i = n. However, we also don't want to split every data point into a separate segment so we need
to account the "cost" of introducing new segments. We use a cost function that's loosely based on two
models discussed in [2] for simplicity. We will tune this cost function further in the future.

The problem of finding the best segmentation then reduces to a search problem. Unfortunately, our problem
space is exponential with respect to the size of the time series since we could split at each data point.
We workaround this problem by first splitting the time series into a manageable smaller grids, and only
considering segmentation of a fixed size (i.e. the number of segments is constant). Since time series
tend to contain a lot more data points than segments, this strategy finds the optimal solution without
exploring much of the problem space.

Finding the optimal segmentation of a fixed size is, itself, another search problem that is equivalent to
finding the shortest path of a fixed length in DAG. Here, we use dynamic programming with a matrix of size
n by n where n is the length of the time series (grid). Each entry in this matrix at (i, k) stores
the minimum cost of segmenting data points 1 through i using k segments. We start our search at i = 1.
Clearly C(1, 0) = 0 (note the actual code uses 0-based index). In i-th iteration, we compute the cost
S(i, j) of each segment starting at i and ending at another point j after i and update C(j, k + 1) by
min( C(j, k + 1), C(i, k) + S(i, j) ) for all values of j above i.

[1] Kensuke Fukuda, H. Eugene Stanley, and Luis A. Nunes Amaral, "Heuristic segmentation of
a nonstationary time series", Physical Review E 69, 021108 (2004)

[2] Marc Lavielle, Gilles Teyssi`ere, "Detection of Multiple Change–Points in Multivariate Time Series"
Lithuanian Mathematical Journal, vol 46, 2006

* public/v2/index.html: Show the optional description for the chosen moving average strategy.
* public/v2/js/statistics.js:
(Statistics.testWelchsT):
(Statistics.computeWelchsT): Extracted from testWelchsT. Generalized to take the offset and the length
of each value array between which Welch's t-statistic is computed. This generalization helps the
Schwarz criterion segmentation algorithm avoid splitting values array O(n^2) times.
(.sampleMeanAndVarianceForValues): Ditto for the generalization.
(.recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent): Added. Implements recursive t-test.
(.splitIntoSegmentsUntilGoodEnough): Added. Implements Schwarz criterion.
(.findOptimalSegmentation): Added. Implements the algorithm to find the optimal segmentation of a fixed
segment count.
(.SampleVarianceUpperTriangularMatrix): Added. Stores S(i, j) used by findOptimalSegmentation.
(.SampleVarianceUpperTriangularMatrix.prototype.costBetween): Added.

git-svn-id: https://svn.webkit.org/repository/webkit/trunk@182330 268f45cc-cd09-0410-ab3c-d52691b4dbfc

Websites/perf.webkit.org/ChangeLog
Websites/perf.webkit.org/public/v2/index.html
Websites/perf.webkit.org/public/v2/js/statistics.js

index bb80e12..608d015 100644 (file)
@@ -1,5 +1,63 @@
 2015-04-03  Ryosuke Niwa  <rniwa@webkit.org>
 
+        Add time series segmentation algorithms as moving averages
+        https://bugs.webkit.org/show_bug.cgi?id=143362
+
+        Reviewed by Chris Dumez.
+
+        This patch implements two preliminary time series segmentation algorithms as moving averages.
+
+        Recursive t-test: Compute Welch's t-statistic at each point in a given segment of the time series.
+        If Welch's t-test implicates a statistically significance difference, then split the segment into two
+        sub segments with the maximum t-statistic (i.e. the point at which if split would yield the highest
+        probability that two segments do not share the same "underlying" mean in classical / frequentist sense).
+        We repeat this process recursively. See [1] for the evaluation of this particular algorithm.
+
+        Schwarz criterion: Use Schwarz or Bayesian information criterion to heuristically find the optimal
+        segmentation. Intuitively, the problem of finding the best segmentation comes down to minimizing the
+        residual sum of squares in each segment as in linear regressions. That is, for a given segment with
+        values y_1 through y_n with mean y_avg, we want to minimize the sum of (y_i - y_avg)^2 over i = 1
+        through i = n. However, we also don't want to split every data point into a separate segment so we need
+        to account the "cost" of introducing new segments. We use a cost function that's loosely based on two
+        models discussed in [2] for simplicity. We will tune this cost function further in the future.
+
+        The problem of finding the best segmentation then reduces to a search problem. Unfortunately, our problem
+        space is exponential with respect to the size of the time series since we could split at each data point.
+        We workaround this problem by first splitting the time series into a manageable smaller grids, and only
+        considering segmentation of a fixed size (i.e. the number of segments is constant). Since time series
+        tend to contain a lot more data points than segments, this strategy finds the optimal solution without
+        exploring much of the problem space.
+
+        Finding the optimal segmentation of a fixed size is, itself, another search problem that is equivalent to
+        finding the shortest path of a fixed length in DAG. Here, we use dynamic programming with a matrix of size
+        n by n where n is the length of the time series (grid). Each entry in this matrix at (i, k) stores
+        the minimum cost of segmenting data points 1 through i using k segments. We start our search at i = 1.
+        Clearly C(1, 0) = 0 (note the actual code uses 0-based index). In i-th iteration, we compute the cost
+        S(i, j) of each segment starting at i and ending at another point j after i and update C(j, k + 1) by
+        min( C(j, k + 1), C(i, k) + S(i, j) ) for all values of j above i.
+
+        [1] Kensuke Fukuda, H. Eugene Stanley, and Luis A. Nunes Amaral, "Heuristic segmentation of
+        a nonstationary time series", Physical Review E 69, 021108 (2004)
+
+        [2] Marc Lavielle, Gilles Teyssi`ere, "Detection of Multiple Change–Points in Multivariate Time Series"
+        Lithuanian Mathematical Journal, vol 46, 2006
+
+        * public/v2/index.html: Show the optional description for the chosen moving average strategy.
+        * public/v2/js/statistics.js:
+        (Statistics.testWelchsT):
+        (Statistics.computeWelchsT): Extracted from testWelchsT. Generalized to take the offset and the length
+        of each value array between which Welch's t-statistic is computed. This generalization helps the
+        Schwarz criterion segmentation algorithm avoid splitting values array O(n^2) times.
+        (.sampleMeanAndVarianceForValues): Ditto for the generalization.
+        (.recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent): Added. Implements recursive t-test.
+        (.splitIntoSegmentsUntilGoodEnough): Added. Implements Schwarz criterion.
+        (.findOptimalSegmentation): Added. Implements the algorithm to find the optimal segmentation of a fixed
+        segment count.
+        (.SampleVarianceUpperTriangularMatrix): Added. Stores S(i, j) used by findOptimalSegmentation.
+        (.SampleVarianceUpperTriangularMatrix.prototype.costBetween): Added.
+
+2015-04-03  Ryosuke Niwa  <rniwa@webkit.org>
+
         REGRESSION: Perf dashboard sometimes fails to update zooming level
         https://bugs.webkit.org/show_bug.cgi?id=143359
 
index 4ecf3bb..c7cd030 100755 (executable)
                     optionValuePath='content'
                     optionLabelPath='content.label'
                     selection=chosenMovingAverageStrategy}}</label>
+                    {{#if chosenMovingAverageStrategy.description}}
+                        <p class="description">{{chosenMovingAverageStrategy.description}}</p>
+                    {{/if}}
                 {{#each chosenMovingAverageStrategy.parameterList}}
                     <label>{{label}}: {{input type="number" value=value min=min max=max step=step}}</label>
                 {{/each}}
index da30e1d..a752240 100755 (executable)
@@ -58,32 +58,64 @@ var Statistics = new (function () {
 
     // Welch's t-test (http://en.wikipedia.org/wiki/Welch%27s_t_test)
     this.testWelchsT = function (values1, values2, probability) {
-        var stat1 = sampleMeanAndVarianceForValues(values1);
-        var stat2 = sampleMeanAndVarianceForValues(values2);
+        return this.computeWelchsT(values1, 0, values1.length, values2, 0, values2.length, probability).significantlyDifferent;
+    }
+
+    this.computeWelchsT = function (values1, startIndex1, length1, values2, startIndex2, length2, probability) {
+        var stat1 = sampleMeanAndVarianceForValues(values1, startIndex1, length1);
+        var stat2 = sampleMeanAndVarianceForValues(values2, startIndex2, length2);
         var sumOfSampleVarianceOverSampleSize = stat1.variance / stat1.size + stat2.variance / stat2.size;
-        var t = (stat1.mean - stat2.mean) / Math.sqrt(sumOfSampleVarianceOverSampleSize);
+        var t = Math.abs((stat1.mean - stat2.mean) / Math.sqrt(sumOfSampleVarianceOverSampleSize));
 
         // http://en.wikipedia.org/wiki/Welch–Satterthwaite_equation
         var degreesOfFreedom = sumOfSampleVarianceOverSampleSize * sumOfSampleVarianceOverSampleSize
             / (stat1.variance * stat1.variance / stat1.size / stat1.size / stat1.degreesOfFreedom
                 + stat2.variance * stat2.variance / stat2.size / stat2.size / stat2.degreesOfFreedom);
-
-        // They're different beyond the confidence interval of the specified probability.
-        return Math.abs(t) > tDistributionQuantiles[probability || 0.9][Math.round(degreesOfFreedom - 1)];
+        return {
+            t: t,
+            degreesOfFreedom: degreesOfFreedom,
+            significantlyDifferent: t > tDistributionQuantiles[probability || 0.9][Math.round(degreesOfFreedom - 1)],
+        };
     }
 
-    function sampleMeanAndVarianceForValues(values) {
-        var sum = Statistics.sum(values);
-        var squareSum = Statistics.squareSum(values);
-        var sampleMean = sum / values.length;
+    function sampleMeanAndVarianceForValues(values, startIndex, length) {
+        var sum = 0;
+        for (var i = 0; i < length; i++)
+            sum += values[startIndex + i];
+        var squareSum = 0;
+        for (var i = 0; i < length; i++)
+            squareSum += values[startIndex + i] * values[startIndex + i];
+        var sampleMean = sum / length;
         // FIXME: Maybe we should be using the biased sample variance.
-        var unbiasedSampleVariance = (squareSum - sum * sum / values.length) / (values.length - 1);
+        var unbiasedSampleVariance = (squareSum - sum * sum / length) / (length - 1);
         return {
             mean: sampleMean,
             variance: unbiasedSampleVariance,
-            size: values.length,
-            degreesOfFreedom: values.length - 1,
+            size: length,
+            degreesOfFreedom: length - 1,
+        }
+    }
+
+    function recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex, length, minLength, segments) {
+        var tMax = 0;
+        var argTMax = null;
+        for (var i = 1; i < length - 1; i++) {
+            var firstLength = i;
+            var secondLength = length - i;
+            if (firstLength < minLength || secondLength < minLength)
+                continue;
+            var result = Statistics.computeWelchsT(values, startIndex, firstLength, values, startIndex + i, secondLength, 0.9);
+            if (result.significantlyDifferent && result.t > tMax) {
+                tMax = result.t;
+                argTMax = i;
+            }
+        }
+        if (!tMax) {
+            segments.push(values.slice(startIndex, startIndex + length));
+            return;
         }
+        recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex, argTMax, minLength, segments);
+        recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, startIndex + argTMax, length - argTMax, minLength, segments);
     }
 
     // One-sided t-distribution.
@@ -193,8 +225,190 @@ var Statistics = new (function () {
                 return averages;
             }
         },
+        {
+            id: 4,
+            label: 'Segmentation: Recursive t-test',
+            description: "Recursively split values into two segments if Welch's t-test detects a statistically significant difference.",
+            parameterList: [{label: "Minimum segment length", value: 20, min: 5}],
+            execute: function (minLength, values) {
+                if (values.length < 2)
+                    return null;
+
+                var averages = new Array(values.length);
+                var segments = new Array;
+                recursivelySplitIntoTwoSegmentsAtMaxTIfSignificantlyDifferent(values, 0, values.length, minLength, segments);
+                var averageIndex = 0;
+                for (var j = 0; j < segments.length; j++) {
+                    var values = segments[j];
+                    var mean = Statistics.sum(values) / values.length;
+                    for (var i = 0; i < values.length; i++)
+                        averages[averageIndex++] = mean;
+                }
+
+                return averages;
+            }
+        },
+        {
+            id: 5,
+            label: 'Segmentation: Schwarz criterion',
+            description: 'Adaptive algorithm that maximizes the Schwarz criterion (BIC).',
+            // Based on Detection of Multiple Change–Points in Multivariate Time Series by Marc Lavielle (July 2006).
+            execute: function (values) {
+                if (values.length < 2)
+                    return null;
+
+                var averages = new Array(values.length);
+                var averageIndex = 0;
+
+                // Split the time series into grids since splitIntoSegmentsUntilGoodEnough is O(n^2).
+                var gridLength = 500;
+                var totalSegmentation = [0];
+                for (var gridCount = 0; gridCount < Math.ceil(values.length / gridLength); gridCount++) {
+                    var gridValues = values.slice(gridCount * gridLength, (gridCount + 1) * gridLength);
+                    var segmentation = splitIntoSegmentsUntilGoodEnough(gridValues);
+
+                    if (Statistics.debuggingSegmentation)
+                        console.log('grid=' + gridCount, segmentation);
+
+                    for (var i = 1; i < segmentation.length - 1; i++)
+                        totalSegmentation.push(gridCount * gridLength + segmentation[i]);
+                }
+
+                if (Statistics.debuggingSegmentation)
+                    console.log('Final Segmentation', totalSegmentation);
+
+                totalSegmentation.push(values.length);
+
+                for (var i = 1; i < totalSegmentation.length; i++) {
+                    var segment = values.slice(totalSegmentation[i - 1], totalSegmentation[i]);
+                    var average = Statistics.sum(segment) / segment.length;
+                    for (var j = 0; j < segment.length; j++)
+                        averages[averageIndex++] = average;
+                }
+
+                return averages;
+            }
+        },
     ];
 
+    this.debuggingSegmentation = false;
+
+    function splitIntoSegmentsUntilGoodEnough(values) {
+        if (values.length < 2)
+            return [0, values.length];
+
+        var matrix = new SampleVarianceUpperTriangularMatrix(values);
+
+        var SchwarzCriterionBeta = Math.log1p(values.length - 1) / values.length;
+
+        var BirgeAndMassartC = 2.5; // Suggested by the authors.
+        var BirgeAndMassartPenalization = function (segmentCount) {
+            return segmentCount * (1 + BirgeAndMassartC * Math.log1p(values.length / segmentCount - 1));
+        }
+
+        var segmentation;
+        var minTotalCost = Infinity;
+        var maxK = 50;
+
+        for (var k = 1; k < maxK; k++) {
+            var start = Date.now();
+            var result = findOptimalSegmentation(values, matrix, k);
+            var cost = result.cost / values.length;
+            var penalty = SchwarzCriterionBeta * BirgeAndMassartPenalization(k);
+            if (cost + penalty < minTotalCost) {
+                minTotalCost = cost + penalty;
+                segmentation = result.segmentation;
+            } else
+                maxK = Math.min(maxK, k + 3);
+            if (Statistics.debuggingSegmentation)
+                console.log('splitIntoSegmentsUntilGoodEnough', k, Date.now() - start, cost + penalty);
+        }
+
+        return segmentation;
+    }
+
+    function findOptimalSegmentation(values, costMatrix, segmentCount) {
+        // Dynamic programming. cost[i][k] = The cost to segmenting values up to i into k segments.
+        var cost = new Array(values.length);
+        for (var i = 0; i < values.length; i++) {
+            cost[i] = new Float32Array(segmentCount + 1);
+        }
+
+        var previousNode = new Array(values.length);
+        for (var i = 0; i < values.length; i++)
+            previousNode[i] = new Array(segmentCount + 1);
+
+        cost[0] = [0]; // The cost of segmenting single value is always 0.
+        previousNode[0] = [-1];
+        for (var segmentStart = 0; segmentStart < values.length; segmentStart++) {
+            var costBySegment = cost[segmentStart];
+            for (var count = 0; count < segmentCount; count++) {
+                if (previousNode[segmentStart][count] === undefined)
+                    continue;
+                for (var segmentEnd = segmentStart + 1; segmentEnd < values.length; segmentEnd++) {
+                    var newCost = costBySegment[count] + costMatrix.costBetween(segmentStart, segmentEnd);
+                    if (previousNode[segmentEnd][count + 1] === undefined || newCost < cost[segmentEnd][count + 1]) {
+                        cost[segmentEnd][count + 1] = newCost;
+                        previousNode[segmentEnd][count + 1] = segmentStart;
+                    }
+                }
+            }
+        }
+
+        if (Statistics.debuggingSegmentation) {
+            console.log('findOptimalSegmentation with k=', segmentCount);
+            for (var i = 0; i < cost.length; i++) {
+                var t = cost[i];
+                var s = '';
+                for (var j = 0; j < t.length; j++) {
+                    var p = previousNode[i][j];
+                    s += '(k=' + j;
+                    if (p !== undefined)
+                        s += ' c=' + t[j] + ' p=' + p
+                    s += ')';
+                }
+                console.log(i, values[i], s);
+            }
+        }
+
+        var currentIndex = values.length - 1;
+        var segmentation = new Array(segmentCount);
+        segmentation[0] = values.length;
+        for (var i = 0; i < segmentCount; i++) {
+            currentIndex = previousNode[currentIndex][segmentCount - i];
+            segmentation[i + 1] = currentIndex;
+        }
+
+        return {segmentation: segmentation.reverse(), cost: cost[values.length - 1][segmentCount]};
+    }
+
+    function SampleVarianceUpperTriangularMatrix(values) {
+        // The cost of segment (i, j].
+        var costMatrix = new Array(values.length - 1);
+        for (var i = 0; i < values.length - 1; i++) {
+            var remainingValueCount = values.length - i - 1;
+            costMatrix[i] = new Float32Array(remainingValueCount);
+            var sum = values[i];
+            var squareSum = sum * sum;
+            costMatrix[i][0] = 0;
+            for (var j = i + 1; j < values.length; j++) {
+                var currentValue = values[j];
+                sum += currentValue;
+                squareSum += currentValue * currentValue;
+                var sampleSize = j - i + 1;
+                var stdev = Statistics.sampleStandardDeviation(sampleSize, sum, squareSum);
+                costMatrix[i][j - i - 1] = sampleSize * Math.log1p(stdev * stdev - 1);
+            }
+        }
+        this.costMatrix = costMatrix;
+    }
+
+    SampleVarianceUpperTriangularMatrix.prototype.costBetween = function(from, to) {
+        if (from >= this.costMatrix.length || from == to)
+            return 0; // The cost of the segment that starts at the last data point is 0.
+        return this.costMatrix[from][to - from - 1];
+    }
+
     this.EnvelopingStrategies = [
         {
             id: 100,