Рекурсивная векторизация python с временными рядами

Question

У меня есть Timeseries (ы), которые нужно обрабатывать рекурсивно, чтобы получить результат таймсерия (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] 

не дает какой-либо ошибки, но я не получаю правильный результат, так как разреш [0] и разрешением [1] должны быть инициализирована до расчета будет иметь место. Есть ли способ обработать его векторизацией?
Любая помощь будет оценена, спасибо!

Answer 1

Это выражение

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

говорит, что res является выходом линейного фильтра (или АРМА процесса) с входным сигналом 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() 

Участок: