Как ограничить ширину окна взаимной корреляции в Numpy?
Я изучаю Numpy/ Scipy, исходя из опыта MATLAB. Функция xcorr в Matlab имеет необязательный аргумент "maxlag", который ограничивает диапазон задержки от –maxlag до maxlag. Это очень полезно, если вы смотрите на взаимную корреляцию между двумя очень длинными временными рядами, но интересуетесь только корреляцией в пределах определенного временного диапазона. Увеличение производительности огромно, учитывая, что кросс-корреляция невероятно дорога для вычисления.
В numpy/ scipy кажется, что есть несколько вариантов вычисления взаимной корреляции. numpy.correlate, numpy.convolve, scipy.signal.fftconvolve. Если кто-то захочет объяснить разницу между ними, я буду рад услышать, но в основном меня беспокоит то, что ни у кого из них нет функции maxlag. Это означает, что, даже если я хочу видеть только корреляции между двумя временными рядами, например, с задержкой от -100 до +100 мс, он все равно будет вычислять корреляцию для каждой задержки от -20000 до +20000 мс (что составляет длину временной ряд). Это дает 200-кратный успех! Нужно ли вручную перекодировать функцию взаимной корреляции, чтобы включить эту функцию?
6 ответов
Вот пара функций для вычисления авто- и взаимной корреляции с ограниченными лагами. Порядок умножения (и сопряжения в сложном случае) был выбран, чтобы соответствовать соответствующему поведению numpy.correlate,
import numpy as np
from numpy.lib.stride_tricks import as_strided
def _check_arg(x, xname):
x = np.asarray(x)
if x.ndim != 1:
raise ValueError('%s must be one-dimensional.' % xname)
return x
def autocorrelation(x, maxlag):
"""
Autocorrelation with a maximum number of lags.
`x` must be a one-dimensional numpy array.
This computes the same result as
numpy.correlate(x, x, mode='full')[len(x)-1:len(x)+maxlag]
The return value has length maxlag + 1.
"""
x = _check_arg(x, 'x')
p = np.pad(x.conj(), maxlag, mode='constant')
T = as_strided(p[maxlag:], shape=(maxlag+1, len(x) + maxlag),
strides=(-p.strides[0], p.strides[0]))
return T.dot(p[maxlag:].conj())
def crosscorrelation(x, y, maxlag):
"""
Cross correlation with a maximum number of lags.
`x` and `y` must be one-dimensional numpy arrays with the same length.
This computes the same result as
numpy.correlate(x, y, mode='full')[len(a)-maxlag-1:len(a)+maxlag]
The return vaue has length 2*maxlag + 1.
"""
x = _check_arg(x, 'x')
y = _check_arg(y, 'y')
py = np.pad(y.conj(), 2*maxlag, mode='constant')
T = as_strided(py[2*maxlag:], shape=(2*maxlag+1, len(y) + 2*maxlag),
strides=(-py.strides[0], py.strides[0]))
px = np.pad(x, maxlag, mode='constant')
return T.dot(px)
Например,
In [367]: x = np.array([2, 1.5, 0, 0, -1, 3, 2, -0.5])
In [368]: autocorrelation(x, 3)
Out[368]: array([ 20.5, 5. , -3.5, -1. ])
In [369]: np.correlate(x, x, mode='full')[7:11]
Out[369]: array([ 20.5, 5. , -3.5, -1. ])
In [370]: y = np.arange(8)
In [371]: crosscorrelation(x, y, 3)
Out[371]: array([ 5. , 23.5, 32. , 21. , 16. , 12.5, 9. ])
In [372]: np.correlate(x, y, mode='full')[4:11]
Out[372]: array([ 5. , 23.5, 32. , 21. , 16. , 12.5, 9. ])
(Было бы неплохо иметь такую функцию в самой NumPy.)
Пока numpy не реализует аргумент maxlag, вы можете использовать функцию ucorrelate из пикоррелатовой упаковки. ucorrelate работает на массивах NumPy и имеет maxlag ключевое слово. Он реализует корреляцию от использования цикла for и оптимизирует скорость выполнения с помощью numba.
Пример - автокорреляция с 3 временными лагами:
import numpy as np
import pycorrelate as pyc
x = np.array([2, 1.5, 0, 0, -1, 3, 2, -0.5])
c = pyc.ucorrelate(x, x, maxlag=3)
c
Результат:
Out[1]: array([20, 5, -3])
Документация Pycorrelate содержит блокнот, показывающий идеальное соответствие между pycorrelate.ucorrelate а также numpy.correlate:
matplotlib.pyplot обеспечивает синтаксис, подобный matlab, для вычисления и построения взаимной корреляции, автокорреляции и т.д.
Ты можешь использовать xcorr который позволяет определить maxlags параметр.
import matplotlib.pyplot as plt
import numpy as np
data = np.arange(0,2*np.pi,0.01)
y1 = np.sin(data)
y2 = np.cos(data)
coeff = plt.xcorr(y1,y2,maxlags=10)
print(*coeff)
[-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7
8 9 10] [ -9.81991753e-02 -8.85505028e-02 -7.88613080e-02 -6.91325329e-02
-5.93651264e-02 -4.95600447e-02 -3.97182508e-02 -2.98407146e-02
-1.99284126e-02 -9.98232812e-03 -3.45104289e-06 9.98555430e-03
1.99417667e-02 2.98641953e-02 3.97518558e-02 4.96037706e-02
5.94189688e-02 6.91964864e-02 7.89353663e-02 8.86346584e-02
9.82934198e-02] <matplotlib.collections.LineCollection object at 0x00000000074A9E80> Line2D(_line0)
Ответ @Warren Weckesser является лучшим, поскольку он использует numpy для экономии производительности (а не просто вызывает corr для каждого лага). Тем не менее, он возвращает перекрестное произведение (например, скалярный продукт между входными данными при различных задержках). Чтобы получить реальную взаимную корреляцию, я изменил его ответ с необязательнымmode аргумент, который, если установлен в 'corr', возвращает взаимную корреляцию как таковую:
def crosscorrelation(x, y, maxlag, mode='corr'):
"""
Cross correlation with a maximum number of lags.
`x` and `y` must be one-dimensional numpy arrays with the same length.
This computes the same result as
numpy.correlate(x, y, mode='full')[len(a)-maxlag-1:len(a)+maxlag]
The return vaue has length 2*maxlag + 1.
"""
py = np.pad(y.conj(), 2*maxlag, mode='constant')
T = as_strided(py[2*maxlag:], shape=(2*maxlag+1, len(y) + 2*maxlag),
strides=(-py.strides[0], py.strides[0]))
px = np.pad(x, maxlag, mode='constant')
if mode == 'dot': # get lagged dot product
return T.dot(px)
elif mode == 'corr': # gets Pearson correlation
return (T.dot(px)/px.size - (T.mean(axis=1)*px.mean())) / \
(np.std(T, axis=1) * np.std(px))
Некоторое время назад я столкнулся с той же проблемой, больше внимания уделял эффективности вычислений. Обратитесь к исходному коду функции MATLAB. xcorr.m, Я сделал простой.
import numpy as np
from scipy import signal, fftpack
import math
import time
def nextpow2(x):
if x == 0:
y = 0
else:
y = math.ceil(math.log2(x))
return y
def xcorr(x, y, maxlag):
m = max(len(x), len(y))
mx1 = min(maxlag, m - 1)
ceilLog2 = nextpow2(2 * m - 1)
m2 = 2 ** ceilLog2
X = fftpack.fft(x, m2)
Y = fftpack.fft(y, m2)
c1 = np.real(fftpack.ifft(X * np.conj(Y)))
index1 = np.arange(1, mx1+1, 1) + (m2 - mx1 -1)
index2 = np.arange(1, mx1+2, 1) - 1
c = np.hstack((c1[index1], c1[index2]))
return c
if __name__ == "__main__":
s = time.clock()
a = [1, 2, 3, 4, 5]
b = [6, 7, 8, 9, 10]
c = xcorr(a, b, 3)
e = time.clock()
print(c)
print(e-c)
В качестве примера возьмем результаты определенного прогона:
[ 29. 56. 90. 130. 110. 86. 59.]
0.0001745000000001884
сравнение с кодом MATLAB:
clear;close all;clc
tic
a = [1, 2, 3, 4, 5];
b = [6, 7, 8, 9, 10];
c = xcorr(a, b, 3)
toc
29.0000 56.0000 90.0000 130.0000 110.0000 86.0000 59.0000
时间已过 0.000279 秒。
Если бы кто-нибудь мог дать точный математический вывод по этому поводу, это было бы очень полезно.
Вот еще один ответ, полученный отсюда, кажется быстрее на полях, чем np.correlate и имеет преимущество возврата нормализованной корреляции:
def rolling_window(self, a, window):
shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
strides = a.strides + (a.strides[-1],)
return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
def xcorr(self, x,y):
N=len(x)
M=len(y)
meany=np.mean(y)
stdy=np.std(np.asarray(y))
tmp=self.rolling_window(np.asarray(x),M)
c=np.sum((y-meany)*(tmp-np.reshape(np.mean(tmp,-1),(N-M+1,1))),-1)/(M*np.std(tmp,-1)*stdy)
return c
Как я ответил здесь, /questions/24559646/ukazhite-otstavanie-v-numpycorrelate/24559656#24559656matplotlib.xcorr имеет параметр maxlags. Это на самом деле обертка numpy.correlate, поэтому нет сохранения производительности. Тем не менее он дает точно такой же результат, который дает функция взаимной корреляции Матлаба. Ниже я отредактировал код из matplotlib, чтобы он возвращал только корреляцию. Причина в том, что если мы используем matplotlib.corr как это, он также вернет сюжет. Проблема в том, что, если мы введем в качестве аргументов комплексный тип данных, мы получим предупреждение "приведение комплекса к реальному типу данных", когда matplotlib попытается нарисовать график.
<!-- language: python -->
import numpy as np
import matplotlib.pyplot as plt
def xcorr(x, y, maxlags=10):
Nx = len(x)
if Nx != len(y):
raise ValueError('x and y must be equal length')
c = np.correlate(x, y, mode=2)
if maxlags is None:
maxlags = Nx - 1
if maxlags >= Nx or maxlags < 1:
raise ValueError('maxlags must be None or strictly positive < %d' % Nx)
c = c[Nx - 1 - maxlags:Nx + maxlags]
return c
Я думаю, что нашел решение, так как столкнулся с той же проблемой:
Если у вас есть два вектора x а также y любой длины N, и требуется взаимная корреляция с окном фиксированной длины m, ты можешь сделать:
x = <some_data>
y = <some_data>
# Trim your variables
x_short = x[window:]
y_short = y[window:]
# do two xcorrelations, lagging x and y respectively
left_xcorr = np.correlate(x, y_short) #defaults to 'valid'
right_xcorr = np.correlate(x_short, y) #defaults to 'valid'
# combine the xcorrelations
# note the first value of right_xcorr is the same as the last of left_xcorr
xcorr = np.concatenate(left_xcorr, right_xcorr[1:])
Помните, что вам может потребоваться нормализовать переменные, если вы хотите ограниченную корреляцию