带有时间序列的python递归向量化 [英] python recursive vectorization with timeseries
问题描述
我有一个时间序列，需要对其进行递归处理才能获得时间序列结果(res).这是我的示例代码:
I have a Timeseries (s) which need to be processed recursively to get a timeseries result (res). Here is my sample code:
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是常数.它可以正常工作，但是我想使用向量化，并且尝试了以下方法:
where c1,c2,c3 are constants. It works properly but I'd like to use vectorization and I tried with:
res[2:]=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]
但我收到"ValueError:操作数无法与形状(1016)(1018)一起广播"
如果我尝试
but I get "ValueError: operands could not be broadcast together with shapes (1016) (1018) "
if I try with
res=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]
没有给出任何错误，但是我没有得到正确的结果，因为在进行计算之前必须先初始化res [0]和res [1].
有没有一种通过矢量化处理它的方法?
任何帮助将不胜感激，谢谢！
doesn't give any error, but I don't get a correct result, because res[0] and res[1] have to be initialized before the calculation will take place.
Is there a way to process it with vectorization?
Any help will be appreciated, thanks!
推荐答案
此表达式
res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]
表示res是带有输入s的线性滤波器(或ARMA进程)的输出.几个库具有用于对此进行计算的功能.这是使用scipy函数的方法 scipy.signal.lfilter .
says that res is the output of a linear filter (or ARMA process) with input s. Several libraries have functions for computing this. Here's how you can use the scipy function scipy.signal.lfilter.
通过检查递归关系，我们得到了滤波器传递函数的分子(b)和分母(a)的系数:
From inspection of the recurrence relation, we get the coefficients of the numerator (b) and denominator (a) of the filter's transfer function:
b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])
我们还需要一个适当的初始条件，以便lfilter处理res[:2] == [0, k].为此，我们使用 scipy.signal.lfiltic :
We'll also need an appropriate initial condition for lfilter to handle res[:2] == [0, k]. For this, we use scipy.signal.lfiltic:
zi = lfiltic(b, a, [k, 0], x=s[1::-1])
在最简单的情况下，将这样调用lfilter:
In the simplest case, one would call lfilter like this:
y = lfilter(b, a, s)
在初始条件为zi的情况下，我们使用:
With an initial condition zi, we use:
y, zo = lfilter(b, a, s, zi=zi)
但是，要完全匹配问题中提供的计算，我们需要输出y以[0, k]开头.因此，我们将分配一个数组y，用[0, k]初始化前两个元素，然后将lfilter的输出分配给y[2:]:
However, to exactly match the calculation provided in the question, we need the output y to start with [0, k]. So we'll allocate an array y, initialize the first two elements with [0, k], and assign the output of lfilter to y[2:]:
y = np.empty_like(s)
y[:2] = [0, k]
y[2:], zo = lfilter(b, a, s[2:], zi=zi)
这是带有原始循环和lfilter的完整脚本:
Here's a complete script with the original loop and with 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参数，因此设置初始条件需要循环.以下脚本显示了一个示例.
lfilter accepts an axis argument, so you can filter an array of signals with a single call. lfiltic does not have an axis argument, so setting up the initial conditions requires a loop. The following script shows an example.
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递归向量化的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
