Как ограничить ширину окна корреляции в Numpy?

Question

Я изучаю numpy/scipy, исходя из фона MATLAB. xcorr function in Matlab имеет необязательный аргумент «maxlag», который ограничивает диапазон от -maxlag до maxlag. Это очень полезно, если вы смотрите на кросс-корреляцию между двумя очень длинными временными рядами, но интересуетесь только корреляцией в течение определенного временного интервала. Производительность огромна, учитывая, что кросс-корреляция невероятно дорога для вычисления.

В numpy/scipy кажется, что существует несколько вариантов вычисления взаимной корреляции. numpy.correlate, numpy.convolve, scipy.signal.fftconvolve. Если кто-то хочет объяснить разницу между ними, я буду рад услышать, но в основном меня беспокоит то, что ни у кого из них нет функции maxlag. Это означает, что даже если я только хочу видеть корреляции между двумя временными рядами с задержками между -100 и +100 мс, например, он все равно вычислит корреляцию для каждого лага между -20000 и +20000 мс (это длина временные ряды). Это дает 200-кратное повышение производительности! Нужно ли вручную перекодировать функцию взаимной корреляции, чтобы включить эту функцию?

Answer 1

matplotlib.pyplot обеспечивает MATLAB подобный синтаксис для computating и черчения перекрестной корреляции, автокорреляции и т.д.

Вы можете использовать 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) 

Answer 2

Я думаю, я нашел решение, как я столкнулся с той же проблемой:

Если у вас есть два вектора x и y любой длины N, и хотите кросс-корреляции с окном фиксированный Len 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:]) 

Помните, что вы, возможно, потребуется normalise переменные, если вы хотите, ограниченная корреляция

Answer 3

Вот еще один ответ, найденный от here, кажется, быстрее, на полях, чем 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   

Answer 4

Вот несколько функций для вычисления авто- и кросс-корреляции с ограниченной запаздывает. Порядок умножения (и сопряжения в сложном случае) был выбран в соответствии с соответствующим поведением 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.)

Answer 5

Пока NumPy не реализует maxlag аргумент, вы можете использовать функцию ucorrelate из pycorrelate package. 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 содержит a notebook показывает идеальный матч между pycorrelate.ucorrelate и numpy.correlate:

Answer 6

как я здесь ответил, https://stackoverflow.com/a/47897581/5122657 matplotlib.xcorr имеет параметр maxlags. На самом деле это обертка numpy.correlate, поэтому нет экономии производительности. Тем не менее он дает точно такой же результат, который дает кросс-корреляционная функция Матлаба. Ниже я редактировал код из matplotlib, чтобы он возвращал только корреляцию. Причина в том, что, если мы используем matplotlib.corr, как он есть, он также вернет участок. Проблема в том, что если мы поместим в нее сложные типы данных, мы получим предупреждение «casting complex to real datatype», когда 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