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

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

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

  • Преобразование Pandas SparseDataframe в Scipy sparse csc_matrix
  • Анализ основных компонентов в Python
  • Создает ли anaconda отдельную переменную PYTHONPATH для каждой новой среды?
  • Поиск корреляционной матрицы
  • Монтаж и построение логнормального
  • pdist для тензора аана
  • Как выбрать элементы по строке из массива NumPy?
  • Эффективный расчет расстояния между N точками и ссылкой в ​​numpy / scipy
  • 4 Solutions collect form web for “Как ограничить ширину окна корреляции в Numpy?”

    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) 

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

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

    Если у вас есть два вектора 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:]) 

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

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

    Легенда не появляется в участке сложенной площади Матплотлиба

    Пропустить первую строку (поле) в цикле с использованием файла CSV?

    Matplotlib imshow, динамическая пересборка на основе масштабирования

    Python: отображение Dict of Dicts с использованием дерева пользовательского интерфейса для ключей и любого другого виджета для значений

    Можно ли выполнять побитовые операции над строкой в ​​Python?

    Какова цель унарного оператора + (pos) в Python?

    Связывание GSL (или другой библиотеки) статически в общей библиотеке

    Numpy где () на двумерной матрице

    Метакласса Python и базовый класс объекта

    Проверьте, существует ли отношение OneToOne в Django

    Заменить не-ASCII-символы одним пространством

    Как установить pip в новую установку python

    Django- Как отображать сообщения, отправленные между пользователями

    получение элементов списка из request.POST в django / python

    Преобразование строки в список бит и наоборот

    Python - лучший язык программирования в мире.