我有一个时间序列 (s),需要对其进行递归处理以获得时间序列结果 (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] 必须在计算发生之前进行初始化。
有没有办法用向量化处理?
任何帮助将不胜感激,谢谢!
最佳答案
这个表达式
res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]
表示 res
是具有输入 s
的线性滤波器(或 ARMA 过程)的输出。有几个库有计算这个的函数。以下是如何使用 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()
剧情:
关于带有时间序列的python递归矢量化,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21336794/