带有时间序列的python递归向量化

Posted 2023-03-12

【中文标题】带有时间序列的python递归向量化

【英文标题】：python recursive vectorization with timeseries

【发布时间】：2014-02-15 16:34:17

【问题描述】：

我有一个时间序列，需要递归处理以获得时间序列结果 (res)。这是我的示例代码：

res=s.copy()*0  
res[1]=k # k is a constant  
for i in range(2,len(s)):  
    res[i]=c1*(s[i]+s[i-1])/2 +c2*res[i-1]+c3*res[i-2]

其中 c1,c2,c3 是常量。它工作正常，但我想使用矢量化，我尝试过：

res[2:]=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]  

但我得到“ValueError：操作数无法与形状一起广播 (1016) (1018)” 如果我尝试使用

res=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]  

没有给出任何错误，但我没有得到正确的结果，因为 res[0] 和 res[1] 必须在计算发生之前初始化。 有没有办法用矢量化处理它？ 任何帮助将不胜感激，谢谢！

【问题讨论】：

用你用于常量的值更新它。如果你让它可以复制粘贴运行，你可能会得到更多帮助。另外 - 你知道 numpy 吗？

这是因为您想加快速度吗？ （如果你稍微重写一下，我猜记忆化可以加快这样的事情？）

对于我喜欢使用 numexpr (code.google.com/p/numexpr) 的此类问题，通常可以毫不费力地显着提高速度。

第一行最好写成res = np.zeros_like(s)

【参考方案1】：

这个表达式

    res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]

说res 是线性滤波器（或ARMA 进程）的输出，输入为s。有几个库有计算这个的函数。以下是如何使用 scipy 函数scipy.signal.lfilter。

通过检查递推关系，我们得到滤波器传递函数的分子（b）和分母（a）的系数：

b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

我们还需要lfilter 的适当初始条件来处理res[:2] == [0, k]。为此，我们使用scipy.signal.lfiltic：

zi = lfiltic(b, a, [k, 0], x=s[1::-1])

在最简单的情况下，可以像这样调用lfilter：

y = lfilter(b, a, s)

使用初始条件zi，我们使用：

y, zo = lfilter(b, a, s, zi=zi)

但是，为了与问题中提供的计算完全匹配，我们需要输出 y 以 [0, k] 开头。所以我们将分配一个数组y，用[0, k]初始化前两个元素，并将lfilter的输出分配给y[2:]：

y = np.empty_like(s)
y[:2] = [0, k]
y[2:], zo = lfilter(b, a, s[2:], zi=zi)

这是带有原始循环和lfilter 的完整脚本：

import numpy as np
from scipy.signal import lfilter, lfiltic


c1 = 0.125
c2 = 0.5
c3 = 0.25

np.random.seed(123)
s = np.random.rand(8)
k = 3.0

# Original version (edited lightly)

res = np.zeros_like(s)
res[1] = k  # k is a constant  
for i in range(2, len(s)):  
    res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]


# Using scipy.signal.lfilter

# Coefficients of the filter's transfer function.
b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

# Create the initial condition of the filter such that
#     y[:2] == [0, k]
zi = lfiltic(b, a, [k, 0], x=s[1::-1])

y = np.empty_like(s)
y[:2] = [0, k]
y[2:], zo = lfilter(b, a, s[2:], zi=zi)

np.set_printoptions(precision=5)
print "res:", res
print "y:  ", y

输出是：

res: [ 0.       3.       1.53206  1.56467  1.24477  1.08496  0.94142  0.84605]
y:   [ 0.       3.       1.53206  1.56467  1.24477  1.08496  0.94142  0.84605]

---

lfilter 接受 axis 参数，因此您可以通过一次调用过滤信号数组。 lfiltic 没有 axis 参数，因此设置初始条件需要一个循环。以下脚本显示了一个示例。

import numpy as np
from scipy.signal import lfilter, lfiltic
import matplotlib.pyplot as plt


# Parameters
c1 = 0.2
c2 = 1.1
c3 = -0.5
k = 1

# Create an array of signals for the demonstration.
np.random.seed(123)
nsamples = 50
nsignals = 4
s = np.random.randn(nsamples, nsignals)

# Coefficients of the filter's transfer function.
b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

# Create the initial condition of the filter for each signal
# such that
#     y[:2] == [0, k]
# We need a loop here, because lfiltic is not vectorized.
zi = np.empty((2, nsignals))
for i in range(nsignals):
    zi[:, i] = lfiltic(b, a, [k, 0], x=s[1::-1, i])

# Create the filtered signals.
y = np.empty_like(s)
y[:2, :] = np.array([0, k]).reshape(-1, 1)
y[2:, :], zo = lfilter(b, a, s[2:], zi=zi, axis=0)

# Plot the filtered signals.
plt.plot(y, linewidth=2, alpha=0.6)
ptp = y.ptp()
plt.ylim(y.min() - 0.05*ptp, y.max() + 0.05*ptp)
plt.grid(True)
plt.show()

剧情：

【讨论】：

非常感谢，您的解决方案比我的更高效、更快捷！有没有办法一次将它应用于信号数据帧（每个数据帧列一个信号），而无需迭代数据帧列并每次都调用 lfilter？

我更新了我的答案以展示如何处理一组信号。