.. index::

   single: moving window statistics

   single: statistics, moving window

   single: online statistics

************************

Moving Window Statistics

************************

This chapter describes routines for computing *moving window

statistics* (also called rolling statistics and running statistics),

using a window around a sample which is used to calculate various

local statistical properties of an input data stream. The window is

then slid forward by one sample to process the next data point and so on.

The functions described in this chapter are declared in the header file

:file:`gsl_movstat.h`.

Introduction

============

This chapter is concerned with calculating various statistics from

subsets of a given dataset. The main idea is to compute statistics

in the vicinity of a given data sample by defining a *window* which

includes the sample itself as well as some specified number of samples

before and after the sample in question. For a sample :math:`x_i`, we

define a window :math:`W_i^{H,J}` as

.. only:: not texinfo

   .. math:: W_i^{H,J} = \left\{ x_{i-H}, \dots, x_i, \dots, x_{i+J} \right\}

.. only:: texinfo

   ::

      W_i^{H,J} = {x_{i-H},...,x_i,...,x_{i+J}}

The parameters :math:`H` and :math:`J` are non-negative integers specifying

the number of samples to include before and after the sample :math:`x_i`.

Statistics such as the mean and standard deviation of the window :math:`W_i^{H,J}`

may be computed, and then the window is shifted forward by one sample to

focus on :math:`x_{i+1}`. The total number of samples in the window is

:math:`K = H + J + 1`. To define a symmetric window centered on :math:`x_i`,

one would set :math:`H = J = \left\lfloor K / 2 \right\rfloor`.

Handling Endpoints

==================

When processing samples near the ends of the input signal, there will not

be enough samples to fill the window :math:`W_i^{H,J}` defined above.

Therefore the user must specify how to construct the windows near the end points.

This is done by passing an input argument of type :type:`gsl_movstat_end_t`:

.. type:: gsl_movstat_end_t

   This data type specifies how to construct windows near end points and can

   be selected from the following choices:

   .. macro:: GSL_MOVSTAT_END_PADZERO

      With this option, a full window of length :math:`K` will be constructed

      by inserting zeros into the window near the signal end points. Effectively,

      the input signal is modified to

      .. only:: not texinfo

         .. math:: \tilde{x} = \{ \underbrace{0, \dots, 0}_{H \textrm{ zeros}}, x_1, x_2, \dots, x_{n-1}, x_n, \underbrace{0, \dots, 0}_{J \textrm{ zeros} } \}

      .. only:: texinfo

         ::

            x~ = {0, ..., 0, x_1, x_2, ..., x_{n-1}, x_n, 0, ..., 0}

      to ensure a well-defined window for all :math:`x_i`.

   .. macro:: GSL_MOVSTAT_END_PADVALUE

      With this option, a full window of length :math:`K` will be constructed

      by padding the window with the first and last sample in the input signal.

      Effectively, the input signal is modified to

      .. only:: not texinfo

         .. math:: \tilde{x} = \{ \underbrace{x_1, \dots, x_1}_{H}, x_1, x_2, \dots, x_{n-1}, x_n, \underbrace{x_n, \dots, x_n}_{J} \}

      .. only:: texinfo

         ::

            x~ = {x_1, ..., x_1, x_1, x_2, ..., x_{n-1}, x_n, x_n, ..., x_n}

   .. macro:: GSL_MOVSTAT_END_TRUNCATE

      With this option, no padding is performed, and the windows are simply truncated

      as the end points are approached.

.. index::

   single: moving window, allocation

Allocation for Moving Window Statistics

=======================================

.. type:: gsl_movstat_workspace

   The moving window statistical routines use a common workspace.

.. function:: gsl_movstat_workspace * gsl_movstat_alloc(const size_t K)

   This function allocates a workspace for computing symmetric, centered moving statistics with a window

   length of :math:`K` samples. In this case, :math:`H = J = \left\lfloor K/2 \right\rfloor`. The size of the workspace

   is :math:`O(7K)`.

.. function:: gsl_movstat_workspace * gsl_movstat_alloc2(const size_t H, const size_t J)

   This function allocates a workspace for computing moving statistics using a window with :math:`H`

   samples prior to the current sample, and :math:`J` samples after the current sample. The

   total window size is :math:`K = H + J + 1`. The size of the workspace is :math:`O(7K)`.

.. function:: void * gsl_movstat_free(gsl_movstat_workspace * w)

   This function frees the memory associated with :data:`w`.

.. index::

   single: moving mean

   single: rolling mean

Moving Mean

===========

The moving window mean calculates the mean of the values of each window :math:`W_i^{H,J}`.

.. only:: not texinfo

   .. math:: \hat{\mu}_i = \frac{1}{\left| W_i^{H,J} \right|} \sum_{x_m \in W_i^{H,J}} x_m

.. only:: texinfo

   ::

      \hat{\mu}_i = 1/| W_i^{H,J} | \sum_{x_m \in W_i^{H,J}} x_m

Here, :math:`\left| W_i^{H,J} \right|` represents the number of elements in the window

:math:`W_i^{H,J}`. This will normally be :math:`K`, unless the :macro:`GSL_MOVSTAT_END_TRUNCATE`

option is selected, in which case it could be less than :math:`K` near the signal end points.

.. function:: int gsl_movstat_mean(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function computes the moving window mean of the input vector :data:`x`, storing

   the output in :data:`y`. The parameter :data:`endtype` specifies how windows near

   the ends of the input should be handled. It is allowed to have :data:`x` = :data:`y`

   for an in-place moving mean.

.. index::

   single: moving variance

   single: moving standard deviation

   single: rolling variance

   single: rolling standard deviation

Moving Variance and Standard Deviation

======================================

The moving window variance calculates the *sample variance* of the values of each window :math:`W_i^{H,J}`,

defined by

.. only:: not texinfo

   .. math:: \hat{\sigma}_i^2 = \frac{1}{\left( \left| W_i^{H,J} \right| - 1 \right)} \sum_{x_m \in W_i^{H,J}} \left( x_m - \hat{\mu}_i \right)^2

.. only:: texinfo

   ::

      \hat{\sigma}_i^2 = 1/(|W_i^{H,J}| - 1) \sum_{x_m \in W_i^{H,J}} ( x_m - \hat{\mu}_i )^2

where :math:`\hat{\mu}_i` is the mean of :math:`W_i^{H,J}` defined above. The standard deviation :math:`\hat{\sigma}_i`

is the square root of the variance.

.. function:: int gsl_movstat_variance(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function computes the moving window variance of the input vector :data:`x`, storing

   the output in :data:`y`. The parameter :data:`endtype` specifies how windows near

   the ends of the input should be handled. It is allowed to have :data:`x` = :data:`y`

   for an in-place moving variance.

.. function:: int gsl_movstat_sd(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function computes the moving window standard deviation of the input vector :data:`x`, storing

   the output in :data:`y`. The parameter :data:`endtype` specifies how windows near

   the ends of the input should be handled. It is allowed to have :data:`x` = :data:`y`

   for an in-place moving standard deviation.

.. index::

   single: moving minimum

   single: moving maximum

   single: rolling minimum

   single: rolling maximum

Moving Minimum and Maximum

==========================

The moving minimum/maximum calculates the minimum and maximum values of

each window :math:`W_i^{H,J}`.

.. only:: not texinfo

    .. math::

       y_i^{min} &= \min \left( W_i^{H,J} \right) \\

       y_i^{max} &= \max \left( W_i^{H,J} \right)

.. only:: texinfo

   ::

      y_i^{min} = \min W_i^{H,J}

      y_i^{max} = \max W_i^{H,J}

.. function:: int gsl_movstat_min(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function computes the moving minimum of the input vector :data:`x`, storing

   the result in :data:`y`. The parameter :data:`endtype` specifies how windows near

   the ends of the input should be handled.

   It is allowed to have :data:`x` = :data:`y` for an in-place moving minimum.

.. function:: int gsl_movstat_max(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function computes the moving maximum of the input vector :data:`x`, storing

   the result in :data:`y`. The parameter :data:`endtype` specifies how windows near

   the ends of the input should be handled.

   It is allowed to have :data:`x` = :data:`y` for an in-place moving maximum.

.. function:: int gsl_movstat_minmax(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y_min, gsl_vector * y_max, gsl_movstat_workspace * w)

   This function computes the moving minimum and maximum of the input vector :data:`x`, storing

   the window minimums in :data:`y_min` and the window maximums in :data:`y_max`.

   The parameter :data:`endtype` specifies how windows near the ends of the input should be handled.

.. index::

   single: moving sum

   single: rolling sum

Moving Sum

==========

The moving window sum calculates the sum of the values of each window :math:`W_i^{H,J}`.

.. math:: y_i = \sum_{x_m \in W_i^{H,J}} x_m

.. function:: int gsl_movstat_sum(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function computes the moving window sum of the input vector :data:`x`, storing

   the output in :data:`y`. The parameter :data:`endtype` specifies how windows near

   the ends of the input should be handled. It is allowed to have :data:`x` = :data:`y`

   for an in-place moving sum.

.. index::

   single: moving median

   single: rolling median

Moving Median

=============

The moving median calculates the median of the window :math:`W_i^{H,J}` for

each sample :math:`x_i`:

.. only:: not texinfo

    .. math:: y_i = \textrm{median} \left( W_i^{H,J} \right)

.. only:: texinfo

   ::

      y_i = median(W_i^{H,J})

.. function:: int gsl_movstat_median(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function computes the moving median of the input vector :data:`x`, storing

   the output in :data:`y`. The parameter :data:`endtype` specifies how windows near

   the ends of the input should be handled. It is allowed for

   :data:`x` = :data:`y` for an in-place moving window median.

Robust Scale Estimation

=======================

A common problem in statistics is to quantify the dispersion (also known as the variability, scatter, and spread) of

a set of data. Often this is done by calculating the variance or standard deviation. However these statistics

are strongly influenced by outliers, and can often provide erroneous results when even a small number of outliers

are present.

Several useful statistics have emerged to provide robust estimates of scale which are not as susceptible to data outliers.

A few of these statistical scale estimators are described below.

.. index::

   single: moving median absolute deviation

   single: rolling median absolute deviation

Moving MAD

----------

The median absolute deviation (MAD) for the window :math:`W_i^{H,J}` is defined

to be the median of the absolute deviations from the window's median:

.. only:: not texinfo

    .. math:: MAD_i = 1.4826 \times \textrm{median} \left( \left| W_i^{H,J} - \textrm{median} \left( W_i^{H,J} \right) \right| \right)

.. only:: texinfo

   ::

      MAD_i = 1.4826 * median[ |W_i^{H,J} - median(W_i^{H,J})| ]

The factor of :math:`1.4826` makes the MAD an unbiased estimator of the standard deviation

for Gaussian data. The MAD has an efficiency of 37%.  See :ref:`here <sec_mad-statistic>` for more information.

.. function:: int gsl_movstat_mad0(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xmedian, gsl_vector * xmad, gsl_movstat_workspace * w)

              int gsl_movstat_mad(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xmedian, gsl_vector * xmad, gsl_movstat_workspace * w)

   These functions compute the moving MAD of the input vector :data:`x` and store the result

   in :data:`xmad`. The medians of each window :math:`W_i^{H,J}` are stored in :data:`xmedian`

   on output. The inputs :data:`x`, :data:`xmedian`, and :data:`xmad` must all be the same length.

   The parameter :data:`endtype` specifies how windows near the ends of the input should be handled.

   The function :code:`mad0` does not include the scale factor of :math:`1.4826`, while the

   function :code:`mad` does include this factor.

.. index::

   single: moving quantile range

   single: rolling quantile range

Moving QQR

----------

The q-quantile range (QQR) is the difference between the :math:`(1-q)` and :math:`q` quantiles

of a set of data,

.. math:: QQR = Q_{1-q} - Q_q

The case :math:`q = 0.25` corresponds to the well-known *interquartile range (IQR)*, which

is the difference between the 75th and 25th percentiles of a set of data. The QQR is

a *trimmed estimator*, the main idea being to discard the largest and smallest values in

a data window and compute a scale estimate from the remaining middle values. In the case

of the IQR, the largest and smallest 25% of the data are discarded and the scale is

estimated from the remaining (middle) 50%.

.. function:: int gsl_movstat_qqr(const gsl_movstat_end_t endtype, const gsl_vector * x, const double q, gsl_vector * xqqr, gsl_movstat_workspace * w)

   This function computes the moving QQR of the input vector :data:`x` and stores the q-quantile ranges

   of each window :math:`W_i^{H,J}` in :data:`xqqr`. The quantile parameter :data:`q` must be between

   :math:`0` and :math:`0.5`. The input :math:`q = 0.25` corresponds to the IQR.

   The inputs :data:`x` and :data:`xqqr` must be the same length.

   The parameter :data:`endtype` specifies how windows near the ends of the input should be handled.

Moving :math:`S_n`

------------------

The :math:`S_n` statistic proposed by Croux and Rousseeuw is based on pairwise differences between

all samples in the window. It has an efficiency of 58%, significantly higher than the MAD.

See :ref:`here <sec_Sn-statistic>` for more information.

.. function:: int gsl_movstat_Sn(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xscale, gsl_movstat_workspace * w)

   This function computes the moving :math:`S_n` of the input vector :data:`x` and stores the output

   in :data:`xscale`. The inputs :data:`x` and :data:`xscale` must be the same length.

   The parameter :data:`endtype` specifies how windows near the ends of the input should be handled.

   It is allowed for :data:`x` = :data:`xscale` for an in-place moving window :math:`S_n`.

Moving :math:`Q_n`

------------------

The :math:`Q_n` statistic proposed by Croux and Rousseeuw is loosely based on the Hodges-Lehmann location

estimator. It has a relatively high efficiency of 82%. See :ref:`here <sec_Qn-statistic>` for more information.

.. function:: int gsl_movstat_Qn(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xscale, gsl_movstat_workspace * w)

   This function computes the moving :math:`Q_n` of the input vector :data:`x` and stores the output

   in :data:`xscale`. The inputs :data:`x` and :data:`xscale` must be the same length.

   The parameter :data:`endtype` specifies how windows near the ends of the input should be handled.

   It is allowed for :data:`x` = :data:`xscale` for an in-place moving window :math:`Q_n`.

.. index::

   single: moving window accumulators

   single: rolling window accumulators

User-defined Moving Statistics

==============================

GSL offers an interface for users to define their own moving window statistics

functions, without needing to implement the edge-handling and accumulator

machinery. This can be done by explicitly constructing the windows

:math:`W_i^{H,J}` for a given input signal (:func:`gsl_movstat_fill`), or by calculating a user-defined

function for each window automatically. In order to apply a user-defined

function to each window, users must define a variable of type

:type:`gsl_movstat_function` to pass into :func:`gsl_movstat_apply`.

This structure is defined as follows.

.. type:: gsl_movstat_function

   Structure specifying user-defined moving window statistical function::

     typedef struct

     {

       double (* function) (const size_t n, double x[], void * params);

       void * params;

     } gsl_movstat_function;

   This structure contains a pointer to the user-defined function as well

   as possible parameters to pass to the function.

   .. member:: double (* function) (const size_t n, double x[], void * params)

      This function returns the user-defined statistic of the array :data:`x`

      of length :data:`n`. User-specified parameters are passed in via :data:`params`.

      It is allowed to modify the array :data:`x`.

   .. member:: void * params

      User-specified parameters to be passed into the function.

.. function:: int gsl_movstat_apply(const gsl_movstat_end_t endtype, const gsl_movstat_function * F, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)

   This function applies the user-defined moving window statistic specified in :data:`F`

   to the input vector :data:`x`, storing the output in :data:`y`.

   The parameter :data:`endtype` specifies how windows near the ends of the input should be handled.

   It is allowed for :data:`x` = :data:`y` for an in-place moving window calculation.

.. function:: size_t gsl_movstat_fill(const gsl_movstat_end_t endtype, const gsl_vector * x, const size_t idx, const size_t H, const size_t J, double * window)

   This function explicitly constructs the sliding window for the input vector :data:`x` which

   is centered on the sample :data:`idx`. On output, the array :data:`window` will contain

   :math:`W_{idx}^{H,J}`. The number of samples to the left and right

   of the sample :data:`idx` are specified by :data:`H` and :data:`J` respectively.

   The parameter :data:`endtype` specifies how windows near the ends of the input should be handled.

   The function returns the size of the window.

Accumulators

============

Many of the algorithms of this chapter are based on an accumulator design, which

process the input vector one sample at a time, updating calculations of the

desired statistic for the current window. Each accumulator is stored in the

following structure:

.. type:: gsl_movstat_accum

   Structure specifying accumulator for moving window statistics::

     typedef struct

     {

       size_t (* size) (const size_t n);

       int (* init) (const size_t n, void * vstate);

       int (* insert) (const double x, void * vstate);

       int (* delete) (void * vstate);

       int (* get) (void * params, double * result, const void * vstate);

     } gsl_movstat_accum;

   The structure contains function pointers responsible for performing

   different tasks for the accumulator.

   .. member:: size_t (* size) (const size_t n)

        This function returns the size of the workspace (in bytes) needed by the accumulator

        for a moving window of length :data:`n`.

   .. member:: int (* init) (const size_t n, void * vstate)

        This function initializes the workspace :data:`vstate` for a moving window of length :data:`n`.

   .. member:: int (* insert) (const double x, void * vstate)

        This function inserts a single sample :data:`x` into the accumulator, updating internal

        calculations of the desired statistic. If the accumulator is full (i.e. :math:`n` samples

        have already been inserted), then the oldest sample is deleted from the accumulator.

   .. member:: int (* delete) (void * vstate)

        This function deletes the oldest sample from the accumulator, updating internal

        calculations of the desired statistic.

   .. member:: int (* get) (void * params, double * result, const void * vstate)

        This function stores the desired statistic for the current window in

        :data:`result`. The input :data:`params` specifies optional parameters

        for calculating the statistic.

The following accumulators of type :type:`gsl_movstat_accum` are defined by GSL to perform moving window statistics

calculations.

.. var:: gsl_movstat_accum_min

         gsl_movstat_accum_max

         gsl_movstat_accum_minmax

   These accumulators calculate moving window minimum/maximums efficiently, using

   the algorithm of D. Lemire.

.. var:: gsl_movstat_accum_mean

         gsl_movstat_accum_sd

         gsl_movstat_accum_variance

   These accumulators calculate the moving window mean, standard deviation, and variance,

   using the algorithm of B. P. Welford.

.. var:: gsl_movstat_median

   This accumulator calculates the moving window median using the min/max heap algorithm

   of Härdle and Steiger.

.. var:: gsl_movstat_accum_Sn

         gsl_movstat_accum_Qn

   These accumulators calculate the moving window :math:`S_n` and :math:`Q_n` statistics

   developed by Croux and Rousseeuw.

.. var:: gsl_movstat_accum_sum

   This accumulator calculates the moving window sum.

.. var:: gsl_movstat_accum_qqr

   This accumulator calculates the moving window q-quantile range.

Examples

========

Example 1

---------

The following example program computes the moving mean, minimum and maximum of a noisy

sinusoid signal of length :math:`N = 500` with a symmetric moving window of size :math:`K = 11`.

.. _fig_movstat1:

.. figure:: /images/movstat1.png

   :scale: 60%

   Original signal time series (gray) with moving mean (green), moving minimum (blue),

   and moving maximum (orange).

The program is given below.

.. include:: examples/movstat1.c

   :code:

Example 2: Robust Scale

-----------------------

The following example program analyzes a time series of length :math:`N = 1000` composed

of Gaussian random variates with zero mean whose standard deviation changes in a piecewise constant fashion

as shown in the table below.

============ ==============

Sample Range :math:`\sigma`

============ ==============

1-200        1.0

201-450      5.0

451-600      1.0

601-850      3.0

851-1000     5.0

============ ==============

Additionally, about 1% of the samples are perturbed to represent outliers by adding

:math:`\pm 15` to the random Gaussian variate.

The program calculates the moving statistics MAD, IQR, :math:`S_n`, :math:`Q_n`, and

the standard deviation using a symmetric moving window of length :math:`K = 41`. The results are shown in

:numref:`fig_movstat2`.

.. _fig_movstat2:

.. figure:: /images/movstat2.png

   :scale: 60%

   Top: time series of piecewise constant variance. Bottom: scale estimates using a moving

   window; the true sigma value is in light blue, MAD in green, IQR in red, :math:`S_n` in yellow, and

   :math:`Q_n` in dark blue. The moving standard deviation is shown in gray.

The robust statistics follow the true standard deviation piecewise changes well, without being

influenced by the outliers. The moving standard deviation (gray curve) is heavily influenced by

the presence of the outliers. The program is given below.

.. include:: examples/movstat2.c

   :code:

Example 3: User-defined Moving Window

-------------------------------------

This example program illustrates how a user can define their own moving window function to apply

to an input vector. It constructs a random noisy time series of length :math:`N = 1000` with

some outliers added. Then it applies a moving window trimmed mean to the time series with

trim parameter :math:`\alpha = 0.1`. The length of the moving window is :math:`K = 11`, so

the smallest and largest sample of each window is discarded prior to computing the mean.

The results are shown in :numref:`fig_movstat3`.

.. _fig_movstat3:

.. figure:: /images/movstat3.png

   :scale: 60%

   Noisy time series data (black) with moving window trimmed mean (red)

The program is given below.

.. include:: examples/movstat3.c

   :code:

References and Further Reading

==============================

The following publications are relevant to the algorithms described

in this chapter,

* W.Hardle and W. Steiger, *Optimal Median Smoothing*, Appl. Statist., 44 (2), 1995.

* D. Lemire, *Streaming Maximum-Minimum Filter Using No More than Three Comparisons per Element*,

  Nordic Journal of Computing, 13 (4), 2006 (https://arxiv.org/abs/cs/0610046).

* B. P. Welford, *Note on a method for calculating corrected sums of squares and products*,

  Technometrics, 4 (3), 1962.