pandas cum扩展之cumvar：增量方差算法

pandas cum扩展之cumvar：增量方差算法

累计方差计算问题

在处理时间序列问题时，通常会遇到累计计算问题，即：对于每个时间$t=1,2,......$，计算从开始到时间$t$的样本$x_0,x_1,...,x_t$的统计量值。pandas中有现成的cumsum、cumprod、cummax和cummin分别用于计算累计求和、累计乘积、累计最大值和最小值。不过没有cumvar函数用来计算累计方差。当然可以通过简单循环来计算每一个时间$t$的累计方差，不过在数据量大的情况下，这样效率太低了。

增量方差算法

有变量$X$的一组长度为$M$的历史样本（记为$H$）：
$$h_1,h_2,h_3,...,h_M$$
其均值为：
$$\overline{H}=\frac{1}{M}\sum _{i=1}^M{h}_i$$
方差为：
$${\delta }_H^2(d)=\frac{1}{M-d}\sum _{i=1}^M(h_i-\overline{H}{)}^2$$
其中$d$为自由度，一般计算样本方差时$d=1$（无偏估计），计算总体方差时$d=0$（有偏估计）。

现在有一组长度为$N$的增量样本（记为$A$）：
$$a_1,a_2,a_3,...,a_N$$
其均值和方差分别为：
$$\begin{gathered} \overline{A}=\frac{1}{N}\sum _{j=1}^N{a}_j \\ {\delta }_A^2(d)=\frac{1}{N-d}\sum _{j=1}^N(a_j-\overline{A}{)}^2 \end{gathered}$$
我们要通过上面给出的两组样本分别的均值和方差公式来计算两组样本合并在一起之后的全样本：
$$h_1,h_2,h_3,...,h_M,a_1,a_2,a_3,...,a_N$$
的均值和方差。

全样本均值为：
$$\begin{array}{rl}\overline{X}& =\frac{1}{M+N}\left[\sum _{i=1}^M{h}_i+\sum _{j=1}^N{a}_j\right]\\ & \\ & =\frac{M\overline{H}+N\overline{A}}{M+N}\end{array}$$
全样本方差为：
$$\begin{array}{rl}{\delta }_X^2(d)& =\frac{1}{M+N-d}\left[\sum _{i=1}^M(h_i-\overline{X}{)}^2+\sum _{j=1}^N(a_j-\overline{X}{)}^2\right]\\ & \\ & =\frac{1}{M+N-d}\left[\sum _{i=1}^M{\left((h_i-\overline{H})-(\overline{X}-\overline{H})\right)}^2+\sum _{j=1}^N{\left((a_j-\overline{A})-(\overline{X}-\overline{A})\right)}^2\right]\\ & \\ & =\frac{1}{M+N-d}[\sum _{i=1}^M\left((h_i-\overline{H}{)}^2-2(h_i-\overline{H})(\overline{X}-\overline{H})+(\overline{X}-\overline{H}{)}^2\right)\\ & \\ & +\sum _{j=1}^N\left((a_j-\overline{A}{)}^2-2(a_j-\overline{A})(\overline{X}-\overline{A})+(\overline{X}-\overline{A}{)}^2\right)]\\ & \\ & =\frac{1}{M+N-d}\left[(M-d){\delta }_H^2(d)+M(\overline{X}-\overline{H}{)}^2+(N-d){\delta }_A^2(d)+N(\overline{X}-\overline{A}{)}^2\right]\end{array}$$

Python实现

通过Python实现增量方差算法。

以10万个样本作为测试，迭代算法用时13.93秒，增量算法用时3.42秒。

以100万个样本作为测试，迭代算法用时81.68秒，增量算法用时6.81秒。

迭代算法用时呈指数递增，增量算法使用是线性递增的，效率提升十分明显。

# -*- coding: utf-8 -*-

import time
import numpy as np


def cumvar_iter(series, ddof=1):
    '''累计方差计算——迭代'''
    
    cumvar = np.nan * np.zeros(len(series),)
    for k in range(len(series)):
        cumvar[k] = np.var(series[:k+1], ddof=ddof)
        
    return cumvar


def cumvar_delta(series, ddof=1):
    '''累计方差计算——增量算法'''
    
    def delta_var(n0, mean0, var0, n1, mean1, var1, ddof=1):
        '''
        增量方差算法
        '''
        n = n0+n1
        if n0 <= ddof or n1 <= ddof or n <= ddof: # 样本量必须大于自由度
            return np.nan
        fm = n - ddof
        mean = (n0 * mean0 + n1 * mean1) / n
        fz1 = (n0-ddof) * var0 + n0 * (mean - mean0) ** 2
        fz2 = (n1-ddof) * var1 + n1 * (mean - mean1) ** 2
        var = (fz1 + fz2) / fm
        return var
    
    # 累计均值
    cummean = np.cumsum(series) / np.arange(1, len(series)+1)
    
    # 累计方差
    cumvar = np.nan * np.ones(len(series),)
    if ddof == 0:
        cumvar[0] = 0
    else:
        for k in range(ddof, ddof+ddof+1):
            cumvar[k] = np.var(series[:k+1], ddof=ddof)
    for k in range(ddof+ddof+1, len(series)):
        var0, mean0, n0 = cumvar[k-ddof-1], cummean[k-ddof-1], k-ddof
        var1 = np.var(series[k-ddof:k+1], ddof=ddof)
        mean1 = np.mean(series[k-ddof:k+1])
        # 增量方差
        cumvar[k] = delta_var(n0, mean0, var0, ddof+1, mean1, var1,
                              ddof=ddof)
    
    return cumvar


if __name__ == '__main__':
    
    # 生成一个长度为十万的测试序列
    series = np.random.randint(10, 1000, (100000,))
    
    start_time = time.time()
    cumvar1 = cumvar_iter(series, ddof=1)
    print(f'迭代算法用时: {round(time.time()-start_time, 6)}s.')
    
    start_time = time.time()
    cumvar2 = cumvar_delta(series, ddof=1)
    print(f'增量算法用时: {round(time.time()-start_time, 6)}s.')
    
	# 结果：
    # 迭代算法用时: 13.930764s.
	# 增量算法用时: 3.415858s.

参考：

​ 算法之美 之 小小方差增量算法带来的大大收益

欢迎关注公众号：一本正经d胡说