# How can I estimate the sliding window standard deviation of a stream?

I am processing a stream of database records. At current levels, about 250 million records are added per week, but this will increase. I wish to compute the 90-day sliding window standard deviation of values taken from these records. I may later decide that I need multiple windows of different lengths if I need to track seasonality, so 365-day windows are also a possibility.

The traditional approach is to store one full windows' worth of data. As data drops off the end, its contribution is subtracted. As new values are added, their contribution is added. I am acquainted with Welford's method as described in Knuth (Art of Computer Programming, Vol 2, page 232, 3rd edition) which guards against the catastrophic loss of precision found in the obvious algorithm.

However, storing the entire history is expensive. I wish to use an approximate algorithm that uses less memory/table space to maintain a sketch of the data from which an estimate may be derived.

I have already found suitable sketches for Max and Min by adapting the PartitionGreedy algorithm from the paper "Competitive analysis of aggregate max in windowed streaming" by Luca Becchetti and Elias Koutsoupias. I am using exponential smoothing with a suitable alpha value keyed to my window length to estimate the Mean and a similar technique for the Sum. Nevertheless, the standard deviation eludes me.

I have studied Count-Min Sketch and its generalization to windowed problems called ECM Sketch ("Sketch-based Querying of Distributed Sliding-Window Data Streams" by Papapetrou, Garofalakis and Deligiannakis), but they seem like overkill and are very complicated for this bear of little brain.

I also tried the range heuristics of SD = (Max-Min) / 4 or SD = (Max-Min) / 6 but the accuracy was horrible. Anything that is reasonably accurate (±5% would be great, ±10% would be tolerable) would be appreciated.

UPDATE:

I found an utterly simple, reasonably accurate, low memory footprint algorithm for estimating quantiles which I can use in place of standard deviation. I was astonished when I read the paper.

For a layman's overview blog post which links to a cool interactive simulator, see http://research.neustar.biz/2013/09/16/sketch-of-the-day-frugal-streaming/

For the paper itself, see "Frugal Streaming for Estimating Quantiles: One (or two) memory suffices" by Qiang Ma, S. Muthukrishnan and Mark Sandler at this URL: http://arxiv.org/pdf/1407.1121v1.pdf

• You shouldn't need the entire history to update the variance; just a history of the width of your window (which can be "wrapped around"; you just need to keep track of the current ends of the window) Commented Aug 19, 2015 at 0:43
• The solution I am seeking should require less data even than the length of the window. I found such an algorithm and updated my question. Commented Aug 19, 2015 at 12:57

You might be able to adapt a technique dating from the dark ages when people calculated standard deviations with hand-operated calculators, so they kept running tallies of both the sums of the observations and of the sums of the squares of the observations. Quoting from the Wikipedia page on mean squared error

the "corrected sample variance" [is]:'

$$S^2_{n-1} = \frac{1}{n-1}\sum_{i=1}^n\left(X_i-\overline{X}\,\right)^2 =\frac{1}{n-1}\left(\sum_{i=1}^n X_i^2-n\overline{X}^2\right)$$

With this formula you get the sample variance directly from the running sum and the running sum of squares.

So if you have an appropriate way to keep a smoothed sum, and it's also appropriate for a smoothed sum of squares, then your problem might be solved.

• This will suffer from exactly the sort of catastrophic cancellation mentioned in the question. Commented Aug 19, 2015 at 0:39

I implemented the frugal streaming algorithm and made a small enhancement that improved the convergence and reduced the error. I am well satisfied with the results for quantiles: less than 5% error, 95% of the time. Using this to compute the first and third quartiles and other code that estimates the max and min, I then estimated the standard deviation using formula 13 in "Estimating the sample mean and standard deviation from the sample size, median, range and/or interquartile range" by Wan, Wang and Tong in BMC Medical Research Methodology (2014).

Here is the C# code I wrote for quantile computations. I use a 3rd party library called FastRandom. You can replace FastRandom with the C# Random class and keep the same method calls without any problems.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Algorithms
{
/// <summary>
/// Maintain a running estimate of a quantile over a stream with very small memory requirements
/// using the algorithm frugal_2u found in:
/// http://arxiv.org/pdf/1407.1121v1.pdf
/// "Frugal Streaming for Estimating Quantiles: One (or two) memory suffices" by Ma, Muthukrishnan and Sandler (2014).
///
/// One can, for instance, track the median value of a stream of data, or the 68th percentile, or the third decile.
/// This estimate follows recent values of data; it is not an estimate over all time.
/// Thus if the quantile you are measuring changes, this will adapt and track the new value.
///
/// Caveat: The published algorithm uses integers. While this implementation uses doubles, the quantile values cannot
/// be resolved any finer than one, the minimum step size. To resolve to finer values would require small
/// changes to this algorithm and much testing to decide how to balance convergence speed with accuracy.
///
/// Usage:
///
///   // Let's track the median, which has quantile = 0.5.
///   var seed = 100; // Educated guess for the median.
///   var estimator = new FrugalQuantile(seed, 0.5, FrugalQuantile.LinearStepAdjuster);
///   IEnumerable data = ... your data ...;
///   foreach (var item in data) {
///       // Do something with estimate...
///   }
///
/// Author: Paul A. Chernoch
/// </summary>
public class FrugalQuantile
{
#region Standard functions you can use for StepAdjuster.

/// <summary>
/// Best step adjuster found so far because it converges fast without overshooting.
/// Every time the step grows by an amount that increases by one:
///    1, 2, 4, 7, 11, 16, 22, 29...
/// </summary>
public static Func<double, double> LinearStepAdjuster = oldStep => oldStep + 1;

/// <summary>
/// Step adjuster used in the published paper, which is good, but not as good as LinearStepAdjuster.
/// Every time the step increases by one:
///    1, 2, 3, 4, 5, 6...
/// </summary>
public static Func<double, double> ConstantStepAdjuster = oldStep => 1;

#endregion

#region Input parameters

/// <summary>
/// Quantile whose estimate will be maintained.
/// If 0.5, the median will be estimated.
/// If 0.75, the third quartile will be estimated.
/// Id 0.2, the second decile will be estimated.
/// etc...
/// </summary>
public double Quantile { get; set; }

/// <summary>
/// Function to dynamically adjust the step size based on the previous step size.
///
/// NOTE: Best function found so far:
///    StepAdjuster = step => step + 1;
/// </summary>
public Func<double, double> StepAdjuster { get; set; }

#endregion

#region Output parameters

/// <summary>
/// The running estimate of the value found at the given quantile.
///
/// This is the value returned by the most recent call to Add.
/// </summary>
public double Estimate { get; set; }

#endregion

#region Internal state

/// <summary>
/// Amount to add to or subtract from the current estimate, depending on whether our estimate is too low or too high.
///
/// As the algorithm proceeds, this is adjusted up and down to improve convergence.
/// </summary>
private double Step { get; set; }

/// <summary>
/// Tracks whether the previous adjustment was to increase the Estimate or decrease it.
///
/// If +1, the Estimate increased.
/// If -1, the Estimate decreased.
/// This should always have the value +1 or -1.
/// </summary>
private SByte Sign { get; set; }

/// <summary>
/// Random number generator.
///
/// Note: One could refactor to use the C# Random class instead. I prefer FastRandom.
/// </summary>
private FastRandom Rand { get; set; }

#endregion

#region Constructors

/// <summary>
/// Create a FrugalQuantile to track a running estimate of a quantile value.
/// </summary>
/// <param name="seed">Initial estimate for the quantile.
/// A good initial estimate permits more rapid convergence.</param>
/// <param name="quantile">Quantile to estimate, in the exclusive range [0,1].
/// The default is 0.5, the median.
/// </param>
/// <param name="stepAdjuster">Function that can update the step size to improve the rate of convergence.
/// Its parameter is the previous step size.
/// The default lambda for this parameter is good, but there are better functions, like this one:
///     stepAdjuster = step => step + 1
/// Researching the function best for your data is recommended.
/// </param>
public FrugalQuantile(double seed, double quantile = 0.5, Func<double,double> stepAdjuster = null)
{
if (quantile <= 0 || quantile >= 1)
throw new ArgumentOutOfRangeException("quantile", "Must be between zero and one, exclusive.");
Quantile = quantile;
Estimate = seed;
Step = 1;
Sign = 1;
// Default lambda for StepAdjuster shown below always return a step change of 1.
// This default is per the published algorithm but testing shows a different function works much better:
Rand = new FastRandom();
}

#endregion

/// <summary>
/// Update the quantile Estimate to reflect the latest value arriving from the stream and return that estimate.
/// </summary>
/// <param name="item">Data Item arriving from the stream.
/// Note: This algorithm was designed for use on non-negative integers. Its accuracy or suitability
/// for negative values is not guaranteed.
/// </param>
/// <returns>The new Estimate.</returns>
{
// This is implemented to resemble as close as possible the pseudo code for function frugal_2u
// http://research.neustar.biz/2013/09/16/sketch-of-the-day-frugal-streaming/
var m = Estimate;
var q = Quantile;
var random = Rand.NextDouble();
if (item > m && random > 1 - q) {
// Increment the step size if and only if the estimate keeps moving in
// the same direction. Step size is incremented by the result of applying
// the specified step function to the previous step size.
Step += (Sign > 0 ? 1 : -1) * f(Step);
// Increment the estimate by step size if step is positive. Otherwise,
// increment the step size by one.
m += Step > 0 ? Step : 1;
// Mark that the estimate increased this step
Sign = 1;
// If the estimate overshot the item in the stream, pull the estimate back
// and re-adjust the step size.
if (m > item) {
Step += (item - m);
m = item;
}
}
else if (item < m && random > q) {
// If the item is less than the stream, follow all of the same steps as
// above, with signs reversed.
Step += (Sign < 0 ? 1 : -1) * f(Step);
m -= Step > 0 ? Step : 1;
Sign = -1;
if (m < item) {
Step += (m - item);
m = item;
}
}
// Damp down the step size to avoid oscillation.
if ((m - item) * Sign < 0 && Step > 1)
Step = 1;

Estimate = m;
return Estimate;
}

}
}