如何使用 NumPy / SciPy 计算移动/运行/滚动任意函数（例如峰度和偏度）

Question

我正在处理时间序列数据。要从数据中获取特征，我必须计算移动均值、中值、众数、斜率、峰度、偏度等。我熟悉 scipy.stat，它提供了一种简单的方法来计算这些量以进行直接计算。但是对于移动/跑步部分，我已经探索了整个互联网，一无所获。

使用numpy 可以很容易地计算出移动均值、中位数和众数。不幸的是，没有用于计算峰度和偏度的内置函数。 如果有人可以提供帮助，如何用 scipy 计算移动峰度和偏度？非常感谢

Answer 1

Pandas 提供了一个DataFrame.rolling() 方法，可以结合其Rolling.apply() 方法（即df.rolling().apply()）将任意函数应用于指定的滚动窗口。

---

如果您正在寻找基于 NumPy 的解决方案，您可以使用FlyingCircus Numeric（免责声明：我是它的主要作者）。

在那里，您可以找到以下内容：

1. flyingcircus_numeric.running_apply(): 可以对一维数组应用任何函数并支持权重，但速度很慢；
2. flyingcircus_numeric.moving_apply()：可以将任何支持axis: int 参数的函数应用于一维数组并支持权重，而且速度很快（但很耗内存）；
3. flyingcircus_numeric.rolling_apply_nd()：可以将支持 axis: int|Sequence[int] 参数的任何函数应用于任何 ND 数组，并且速度快（并且内存效率高），但不支持权重。

根据您的要求，我建议使用rolling_apply_nd()，例如：

import numpy as np
import scipy as sp
import flyingcircus_numeric as fcn

import scipy.stats


NUM = 30
arr = np.arange(NUM)

window = 4
new_arr = fcn.rolling_apply_nd(arr, window, func=sp.stats.kurtosis)
print(new_arr)
# [-1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36
#  -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36 -1.36
#  -1.36 -1.36 -1.36]

当然，请随意查看source code，它是开源的 (GPL)。

---

编辑

只是为了感受一下我们所说的那种速度，以下是 FlyingCircus 中实施的解决方案的基准：

一般方法flyingcircus_numeric.running_apply() 比flyingcircus_numeric.rolling_apply_nd() 或flyingcircus_numeric.moving_apply() 慢几个数量级，第一个大约是。比第二个快一个数量级。 这显示了一般性或支持加权的速度价格。

以上图表是使用here的脚本和以下代码获得的：

import scipy as sp
import flyingcircus_numeric as fcn 
import scipy.stats


WINDOW = 4
FUNC = sp.stats.kurtosis


def my_rolling_apply_nd(arr, window=WINDOW, func=FUNC):
    return fcn.rolling_apply_nd(arr, window, func=FUNC)


def my_moving_apply(arr, window=WINDOW, func=FUNC):
    return fcn.moving_apply(arr, window, func)


def my_running_apply(arr, window=WINDOW, func=FUNC):
    return fcn.running_apply(arr, window, func)


def equal_output(a, b):
    return np.all(np.isclose(a, b))


input_sizes = (5, 10, 50, 100, 500, 1000, 5000, 10000, 50000, 100000)
funcs = my_rolling_apply_nd, my_moving_apply, my_running_apply

runtimes, input_sizes, labels, results = benchmark(
    funcs, gen_input=np.random.random, equal_output=equal_output,
    input_sizes=input_sizes)

plot_benchmarks(runtimes, input_sizes, labels, units='s')
plot_benchmarks(runtimes, input_sizes, labels, units='ms', zoom_fastest=8)

---

（已编辑以反映FlyingCircus 的一些重构）

Answer 2

玩了之后，我想出了一个纯粹基于numpy 和scipy 的解决方案。当然是使用scipy.statskurtosis和skew。

import numpy as np
from scipy.stats import kurtosis, skew

# Window size
N = 4

# Some random data
m = np.array([2, 3, 10, 11, 0, 4, 8, 2, 5, 9])

# Running Kurtosis
def runningKurt(x, N):
    # Initilize placeholder array
    y = np.zeros((len(x) - (N - 1),))
    for i in range(len(x) - (N - 1)):

         y[i] = kurtosis(x[i:(i + N)])

    return y

# Running Kurtosis

def runningSkew(x, N):
    # Initilize placeholder array
    y = np.zeros((len(x) - (N - 1),))
    for i in range(len(x) - (N - 1)):

         y[i] = skew(x[i:(i + N)])

    return y

kurt = runningKurt(m, N)
print("kurtosis : ", kurt)
# kurtosis :  [-1.93940828 -1.77879935 -1.61464214 -1.40236694 -1.15428571 -1.07626667 -1.42666667]


skw = runningSkew(m, N)
print("skew : ", skw)
# skew :  [ 0.         -0.1354179  -0.26356495 -0.13814702  0.43465076  0.32331615 -0.36514837]