python递归矢量化与时间序列

Question

我有一个时间序列需要递归处理以获得时间序列结果(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].有没有办法用矢量化处理它？
任何帮助将不胜感激,谢谢!

Answer 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.

通过检查递归关系,我们得到滤波器传递函数的numerator(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()

情节: