acf

mdtools.statistics.acf(x, axis=None, dt=1, dtau=1, tau_max=None, center=True, unbiased=False)[source]

Calculate the autocorrelation function of an array.

Parameters:
  • x (array_like) – The input array.

  • axis (int or None, optional) – The axis along which to calculate the ACF. By default, the flattened input array is used.

  • dt (int or None, optional) – Distance between restarting times (actually restarting indices) \(t_0\). If None is given, dt is set to the current lag time \(\tau\) which results in independent sampling windows.

  • dtau (int, optional) – Distance between evaluated lag times (actually lag indices) \(\tau\).

  • tau_max (int or None, optional) – Maximum lag time \(\tau\) upon which to compute the ACF. Because the first lag time is 0, the last evaluated lag time is actually tau_max - 1. By default, tau_max is set to x.size is axis is None and to x.shape[axis] otherwise. tau_max must lie within the interval [0, len(x)].

  • center (bool, optional) – If True, center the input array around its mean, i.e. subtract the sample mean \(\bar{X}\) when calculating the variance and covariance.

  • unbiased (bool, optional) – If True, the covariance \(\text{Cov}[X_{t_0}, X_{t_0 + \tau}]\) is normed by the actual number of sample points \(t_0\) (which depends on the lag time \(\tau\)). If False, the covariance is for all lag times normed by the number of sampling points at \(t_0 = 0\) (which is equal to the length of the input array x). Note that setting unbiased to True might result in values of the ACF that lie outside the interval [-1, 1], especially for lag times larger than len(x) // 2 if center is True.

Returns:

ac (numpy.ndarray) – Autocorrelation function of x. If axis is None, the output array has shape (n,) where n is int(np.ceil(tau_max / dtau)). If a specific axis is given, ac has the same shape as x except along the given axis where the length is given by n.

Raises:

ValueError : – If the first element of ac is not unity or one or more elements of ac fall outside the interval [-1, 1].

See also

mdtools.statistics.acf_np()

Different implementation of the ACF using numpy.correlate()

mdtools.statistics.acf_se()

Calculate the standard errors of an autocorrelation function

mdtools.statistics.acf_confint()

Calculate the confidence intervals of an autocorrelation function

Notes

This function computes the autocorrelation function (ACF) of the input array x assuming that x describes a weakly stationary stochastic process. That means, the mean does not vary with time, the variance is finite for all times and the autocovariance does not depend on the starting time \(t_0\) but only on the lag time \(\tau\). In this case the ACF is given by:[1],[2]

\[C(\tau) = \frac{\text{Cov}[X_{t_0}, X_{t_0 + \tau}]}{\text{Var}[X]} = \frac{\langle (X_{t_0} - \mu) \cdot (X_{t_0 + \tau} - \mu) \rangle}{\langle (X - \mu)^2 \rangle}\]

The ACF strictly fulfills \(C(\tau) \in [-1, 1] \text{ } \forall \text{ } \tau \in \mathbb{R}\). Numerically, the ACF is calculated using the following formula:[3],[4]

\[C_\tau = \frac{\frac{1}{T^*} \sum_{t_0=0}^{T-\tau} (X_{t_0} - \bar{X}) (X_{t_0 + \tau} - \bar{X})}{\frac{1}{T} \sum_{t_0=0}^T (X_{t_0} - \bar{X})^2} = \frac{\frac{1}{T^*} \sum_{t_0=0}^{T-\tau} (X_{t_0} - \bar{X}) (X_{t_0 + \tau} - \bar{X})}{C_0}\]

The increment of the sums in the formula of \(C_\tau\) is given by dt. \(\bar{X} = \frac{1}{T} \sum_{t=1}^T X_t\) is the sample mean. For the sample mean, the increment of the sum is always 1, independent of the choice of dt. Note that the sample mean is only subtracted when center is True. T is always the length of the input array x, independent of the choice of tau_max. \(T^*\) is either \(T - \tau\) if unbiased is True or simply \(T\) if unbiased is False.

For the sake of completeness, the ACF of a non-stationary stochastic process is given by:[1],[2]

\[C(t_1, t_2) = \frac{\text{Cov}[X_{t_1}, X_{t_2}]}{\sqrt{\text{Var}[X_{t_1}]} \cdot \sqrt{\text{Var}[X_{t_2}]}} = \frac{\langle (X_{t_1} - \mu_{t_1}) \cdot (X_{t_2} - \mu_{t_2}) \rangle}{\sqrt{\langle (X_{t_1} - \mu_{t_1})^2 \rangle} \cdot \sqrt{\langle (X_{t_2} - \mu_{t_2})^2 \rangle}}\]

References

Examples

>>> a = np.array([-2,  2, -2,  2, -2])
>>> mdt.stats.acf(a)
array([ 1.        , -0.8       ,  0.56666667, -0.4       ,  0.13333333])
>>> mdt.stats.acf(a, dt=None)
array([ 1.        , -0.8       ,  0.26666667, -0.2       ,  0.13333333])
>>> mdt.stats.acf(a, unbiased=True)
array([ 1.        , -1.        ,  0.94444444, -1.        ,  0.66666667])
>>> mdt.stats.acf(a, unbiased=True, dt=None)
array([ 1.        , -1.        ,  0.66666667, -1.        ,  0.66666667])
>>> mdt.stats.acf(a, center=False)
array([ 1. , -0.8,  0.6, -0.4,  0.2])
>>> mdt.stats.acf(a, center=False, dtau=2)
array([1. , 0.6, 0.2])
>>> mdt.stats.acf(a, center=False, tau_max=len(a)//2)
array([ 1. , -0.8])
>>> mdt.stats.acf(a, center=False, dt=None)
array([ 1. , -0.8,  0.4, -0.2,  0.2])
>>> mdt.stats.acf(a, center=False, unbiased=True)
array([ 1., -1.,  1., -1.,  1.])
>>> mdt.stats.acf(a, center=False, unbiased=True, dt=None)
array([ 1., -1.,  1., -1.,  1.])
>>> a = np.array([-2, -2, -2,  2,  2])
>>> mdt.stats.acf(a)
array([ 1.        ,  0.36666667, -0.26666667, -0.4       , -0.2       ])
>>> mdt.stats.acf(a, dt=None)
array([ 1.        ,  0.36666667, -0.06666667, -0.2       , -0.2       ])
>>> mdt.stats.acf(a, unbiased=True)
array([ 1.        ,  0.45833333, -0.44444444, -1.        , -1.        ])
>>> mdt.stats.acf(a, unbiased=True, dt=None)
array([ 1.        ,  0.45833333, -0.16666667, -1.        , -1.        ])
>>> mdt.stats.acf(a, center=False)
array([ 1. ,  0.4, -0.2, -0.4, -0.2])
>>> mdt.stats.acf(a, center=False, dt=None)
array([ 1. ,  0.4,  0. , -0.2, -0.2])
>>> mdt.stats.acf(a, center=False, unbiased=True)
array([ 1.        ,  0.5       , -0.33333333, -1.        , -1.        ])
>>> mdt.stats.acf(a, center=False, unbiased=True, dt=None)
array([ 1. ,  0.5,  0. , -1. , -1. ])
>>> a = np.array([[-2,  2, -2,  2, -2],
...               [-2,  2, -2,  2, -2],
...               [-2, -2, -2,  2,  2],
...               [-2, -2, -2,  2,  2]])
>>> mdt.stats.acf(a, center=False)
array([ 1.  , -0.15,  0.  , -0.25, -0.2 ,  0.55,  0.  ,  0.15, -0.1 ,
       -0.25,  0.1 ,  0.05,  0.1 ,  0.05, -0.2 ,  0.05,  0.  ,  0.05,
        0.  , -0.05])
>>> mdt.stats.acf(a, center=False, axis=0)
array([[ 1.  ,  1.  ,  1.  ,  1.  ,  1.  ],
       [ 0.75,  0.25,  0.75,  0.75,  0.25],
       [ 0.5 , -0.5 ,  0.5 ,  0.5 , -0.5 ],
       [ 0.25, -0.25,  0.25,  0.25, -0.25]])
>>> mdt.stats.acf(a, center=False, axis=0, unbiased=True)
array([[ 1.        ,  1.        ,  1.        ,  1.        ,  1.        ],
       [ 1.        ,  0.33333333,  1.        ,  1.        ,  0.33333333],
       [ 1.        , -1.        ,  1.        ,  1.        , -1.        ],
       [ 1.        , -1.        ,  1.        ,  1.        , -1.        ]])
>>> mdt.stats.acf(a, center=False, axis=1)
array([[ 1. , -0.8,  0.6, -0.4,  0.2],
       [ 1. , -0.8,  0.6, -0.4,  0.2],
       [ 1. ,  0.4, -0.2, -0.4, -0.2],
       [ 1. ,  0.4, -0.2, -0.4, -0.2]])
>>> mdt.stats.acf(a, center=False, axis=1, unbiased=True)
array([[ 1.        , -1.        ,  1.        , -1.        ,  1.        ],
       [ 1.        , -1.        ,  1.        , -1.        ,  1.        ],
       [ 1.        ,  0.5       , -0.33333333, -1.        , -1.        ],
       [ 1.        ,  0.5       , -0.33333333, -1.        , -1.        ]])
>>> a = np.array([[[-2,  2, -2,  2, -2],
...                [-2,  2, -2,  2, -2]],
...
...               [[-2, -2, -2,  2,  2],
...                [-2, -2, -2,  2,  2]]])
>>> mdt.stats.acf(a, center=False)
array([ 1.  , -0.15,  0.  , -0.25, -0.2 ,  0.55,  0.  ,  0.15, -0.1 ,
       -0.25,  0.1 ,  0.05,  0.1 ,  0.05, -0.2 ,  0.05,  0.  ,  0.05,
        0.  , -0.05])
>>> mdt.stats.acf(a, center=False, axis=0)
array([[[ 1. ,  1. ,  1. ,  1. ,  1. ],
        [ 1. ,  1. ,  1. ,  1. ,  1. ]],

       [[ 0.5, -0.5,  0.5,  0.5, -0.5],
        [ 0.5, -0.5,  0.5,  0.5, -0.5]]])
>>> mdt.stats.acf(a, center=False, axis=0, unbiased=True)
array([[[ 1.,  1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.,  1.]],

       [[ 1., -1.,  1.,  1., -1.],
        [ 1., -1.,  1.,  1., -1.]]])
>>> mdt.stats.acf(a, center=False, axis=1)
array([[[1. , 1. , 1. , 1. , 1. ],
        [0.5, 0.5, 0.5, 0.5, 0.5]],

       [[1. , 1. , 1. , 1. , 1. ],
        [0.5, 0.5, 0.5, 0.5, 0.5]]])
>>> mdt.stats.acf(a, center=False, axis=1, unbiased=True)
array([[[1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1.]],

       [[1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1.]]])
>>> mdt.stats.acf(a, center=False, axis=2)
array([[[ 1. , -0.8,  0.6, -0.4,  0.2],
        [ 1. , -0.8,  0.6, -0.4,  0.2]],

       [[ 1. ,  0.4, -0.2, -0.4, -0.2],
        [ 1. ,  0.4, -0.2, -0.4, -0.2]]])
>>> mdt.stats.acf(a, center=False, axis=2, unbiased=True)
array([[[ 1.        , -1.        ,  1.        , -1.        ,
          1.        ],
        [ 1.        , -1.        ,  1.        , -1.        ,
          1.        ]],

       [[ 1.        ,  0.5       , -0.33333333, -1.        ,
         -1.        ],
        [ 1.        ,  0.5       , -0.33333333, -1.        ,
         -1.        ]]])