Как я могу использовать numpy.correlate для автокорреляции?

Мне нужно сделать автокорреляцию набора чисел, что, как я понимаю, просто корреляция набора с самим собой.

Я пробовал использовать функцию корреляции numpy, но я не верю результату, так как он почти всегда дает вектор, где первое число не самым большим, как должно быть.

Итак, это действительно два вопроса:

1. Что именно numpy.correlate делает?
2. Как я могу использовать его (или что-то еще) для автокорреляции?

Чтобы ответить на ваш первый вопрос, numpy.correlate(a, v, mode)выполняет свертку aс обратным vи выдает результаты, обрезанные указанным режимом. Определение свертки , С (т) = Е _{-∞ <я <∞ я} v _{т + я} где -∞ <т <∞, позволяет сделать результаты от -∞ до ∞, но вы , очевидно , не может хранить бесконечно долго массив. Таким образом, он должен быть обрезан, и именно здесь вступает в силу режим. Есть 3 разных режима: полный, одинаковый и действительный:

• «полный» режим возвращается результаты для каждого , tгде оба aи vимеют некоторое перекрытие.
• «такой же» режим возвращает результат той же длины, что и самый короткий вектор ( aили v).
• «действительный» возвращает режим результаты только когда aи vполностью перекрывают друг друга. В документации для numpy.convolveболее подробно описаны режимы.

Что касается вашего второго вопроса, я думаю, что numpy.correlate он дает вам автокорреляцию, он просто дает вам немного больше. Автокорреляция используется, чтобы определить, насколько сигнал или функция похож на сам себя в определенную разницу во времени. При разнице во времени, равной 0, автокорреляция должна быть максимальной, потому что сигнал идентичен самому себе, поэтому вы ожидали, что первый элемент в массиве результатов автокорреляции будет максимальным. Однако корреляция не начинается при разнице во времени, равной 0. Она начинается при отрицательной разнице во времени, закрывается до 0, а затем становится положительной. То есть вы ожидали:

автокорреляция (a) = ∑ _{-∞ <i <∞} a _i v _{t + i,} где 0 <= t <∞

Но у вас было:

автокорреляция (a) = ∑ _{-∞ <i <∞} a _i v _{t + i,} где -∞ <t <∞

Вам нужно взять последнюю половину результата корреляции, и это должна быть автокорреляция, которую вы ищете. Простая функция python для этого:

def autocorr(x):
    result = numpy.correlate(x, x, mode='full')
    return result[result.size/2:]

Вам, конечно, потребуется проверка ошибок, чтобы убедиться, что xэто действительно одномерный массив. Кроме того, это объяснение, вероятно, не самое математически строгое. Я бросал бесконечности, потому что определение свертки использует их, но это не обязательно применимо к автокорреляции. Таким образом, теоретическая часть этого объяснения может быть немного шаткой, но, надеюсь, практические результаты будут полезными. Эти страницы, посвященные автокорреляции, очень полезны и могут дать вам гораздо лучшую теоретическую базу, если вы не против пробираться через нотации и сложные концепции.

Автокорреляция бывает двух версий: статистической и свертки. Они оба делают то же самое, за исключением небольшой детали: статистическая версия нормализована для интервала [-1,1]. Вот пример того, как вы делаете статистический:

def acf(x, length=20):
    return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:])[0,1]  \
        for i in range(1, length)])

Используйте numpy.corrcoefфункцию вместо, numpy.correlateчтобы вычислить статистическую корреляцию для запаздывания t:

def autocorr(x, t=1):
    return numpy.corrcoef(numpy.array([x[:-t], x[t:]]))

Думаю, есть 2 вещи, которые добавляют путаницу в эту тему:

1. определение статистической и обработки сигналов: как указывали другие, в статистике мы нормализуем автокорреляцию до [-1,1].
2. частичное или неполное среднее значение / дисперсия: когда временной ряд сдвигается с задержкой> 0, их размер перекрытия всегда будет <исходной длины. Используем ли мы среднее и стандартное значение исходного (неполного) или всегда вычисляем новое среднее значение и стандартное значение, используя постоянно меняющееся перекрытие (частичное), имеет значение. (Вероятно, есть формальный термин для этого, но сейчас я буду использовать «частичный»).

Я создал 5 функций, которые вычисляют автокорреляцию 1d-массива с частичным отличием от неполного. Некоторые используют формулы из статистики, некоторые используют коррелят в смысле обработки сигналов, что также можно сделать с помощью БПФ. Но все результаты являются автокорреляциями в определении статистики , поэтому они иллюстрируют, как они связаны друг с другом. Код ниже:

import numpy
import matplotlib.pyplot as plt

def autocorr1(x,lags):
    '''numpy.corrcoef, partial'''

    corr=[1. if l==0 else numpy.corrcoef(x[l:],x[:-l])[0][1] for l in lags]
    return numpy.array(corr)

def autocorr2(x,lags):
    '''manualy compute, non partial'''

    mean=numpy.mean(x)
    var=numpy.var(x)
    xp=x-mean
    corr=[1. if l==0 else numpy.sum(xp[l:]*xp[:-l])/len(x)/var for l in lags]

    return numpy.array(corr)

def autocorr3(x,lags):
    '''fft, pad 0s, non partial'''

    n=len(x)
    # pad 0s to 2n-1
    ext_size=2*n-1
    # nearest power of 2
    fsize=2**numpy.ceil(numpy.log2(ext_size)).astype('int')

    xp=x-numpy.mean(x)
    var=numpy.var(x)

    # do fft and ifft
    cf=numpy.fft.fft(xp,fsize)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real
    corr=corr/var/n

    return corr[:len(lags)]

def autocorr4(x,lags):
    '''fft, don't pad 0s, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean

    cf=numpy.fft.fft(xp)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real/var/len(x)

    return corr[:len(lags)]

def autocorr5(x,lags):
    '''numpy.correlate, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean
    corr=numpy.correlate(xp,xp,'full')[len(x)-1:]/var/len(x)

    return corr[:len(lags)]


if __name__=='__main__':

    y=[28,28,26,19,16,24,26,24,24,29,29,27,31,26,38,23,13,14,28,19,19,\
            17,22,2,4,5,7,8,14,14,23]
    y=numpy.array(y).astype('float')

    lags=range(15)
    fig,ax=plt.subplots()

    for funcii, labelii in zip([autocorr1, autocorr2, autocorr3, autocorr4,
        autocorr5], ['np.corrcoef, partial', 'manual, non-partial',
            'fft, pad 0s, non-partial', 'fft, no padding, non-partial',
            'np.correlate, non-partial']):

        cii=funcii(y,lags)
        print(labelii)
        print(cii)
        ax.plot(lags,cii,label=labelii)

    ax.set_xlabel('lag')
    ax.set_ylabel('correlation coefficient')
    ax.legend()
    plt.show()

Вот результат:

Мы не видим все 5 линий, потому что 3 из них перекрываются (фиолетовые). Перекрытия - это все неполные автокорреляции. Это связано с тем, что вычисления на основе методов обработки сигналов ( np.correlate, FFT) не вычисляют другое среднее / стандартное значение для каждого перекрытия.

Также обратите внимание, что результат fft, no padding, non-partial(красная линия) отличается, потому что он не заполняет временные ряды нулями перед выполнением БПФ, поэтому это круговое БПФ. Я не могу подробно объяснить почему, это то, что я узнал из других источников.

Поскольку я только что столкнулся с той же проблемой, я хотел бы поделиться с вами несколькими строками кода. Фактически, к настоящему времени существует несколько довольно похожих сообщений об автокорреляции в stackoverflow. Если вы определяете автокорреляцию как a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2)[это определение, данное в функции IDL a_correlate, и оно согласуется с тем, что я вижу в ответе 2 на вопрос № 12269834 ], то, похоже, следующее дает правильные результаты:

import numpy as np
import matplotlib.pyplot as plt

# generate some data
x = np.arange(0.,6.12,0.01)
y = np.sin(x)
# y = np.random.uniform(size=300)
yunbiased = y-np.mean(y)
ynorm = np.sum(yunbiased**2)
acor = np.correlate(yunbiased, yunbiased, "same")/ynorm
# use only second half
acor = acor[len(acor)/2:]

plt.plot(acor)
plt.show()

Как видите, я тестировал это с кривой греха и равномерным случайным распределением, и оба результата выглядят так, как я ожидал их. Обратите внимание, что я использовал mode="same"вместо того, mode="full"как другие.

Ваш вопрос 1 уже подробно обсуждался в нескольких отличных ответах здесь.

Я решил поделиться с вами несколькими строками кода, которые позволят вам вычислить автокорреляцию сигнала, основываясь только на математических свойствах автокорреляции. То есть автокорреляцию можно вычислить следующим образом:

1. вычесть среднее значение из сигнала и получить несмещенный сигнал

2. вычислить преобразование Фурье несмещенного сигнала

3. вычислить спектральную плотность мощности сигнала, взяв квадратную норму каждого значения преобразования Фурье несмещенного сигнала

4. вычислить обратное преобразование Фурье спектральной плотности мощности

5. нормализовать обратное преобразование Фурье спектральной плотности мощности на сумму квадратов несмещенного сигнала и взять только половину результирующего вектора

Код для этого следующий:

def autocorrelation (x) :
    """
    Compute the autocorrelation of the signal, based on the properties of the
    power spectral density of the signal.
    """
    xp = x-np.mean(x)
    f = np.fft.fft(xp)
    p = np.array([np.real(v)**2+np.imag(v)**2 for v in f])
    pi = np.fft.ifft(p)
    return np.real(pi)[:x.size/2]/np.sum(xp**2)

Я компьютерный биолог, и когда мне пришлось вычислить авто / кросс-корреляции между парами временных рядов стохастических процессов, я понял, что np.correlateне выполняет нужную мне работу.

В самом деле, что , кажется, отсутствует np.correlateявляется усреднение по всем возможным парам точек времени на расстоянии т.

Вот как я определил функцию, выполняющую то, что мне нужно:

def autocross(x, y):
    c = np.correlate(x, y, "same")
    v = [c[i]/( len(x)-abs( i - (len(x)/2)  ) ) for i in range(len(c))]
    return v

Мне кажется, что ни один из предыдущих ответов не касается этого случая авто / взаимной корреляции: надеюсь, этот ответ может быть полезен кому-то, кто работает со стохастическими процессами, такими как я.

Я использую talib.CORREL для такой автокорреляции, я подозреваю, что вы можете сделать то же самое с другими пакетами:

def autocorrelate(x, period):

    # x is a deep indicator array 
    # period of sample and slices of comparison

    # oldest data (period of input array) may be nan; remove it
    x = x[-np.count_nonzero(~np.isnan(x)):]
    # subtract mean to normalize indicator
    x -= np.mean(x)
    # isolate the recent sample to be autocorrelated
    sample = x[-period:]
    # create slices of indicator data
    correls = []
    for n in range((len(x)-1), period, -1):
        alpha = period + n
        slices = (x[-alpha:])[:period]
        # compare each slice to the recent sample
        correls.append(ta.CORREL(slices, sample, period)[-1])
    # fill in zeros for sample overlap period of recent correlations    
    for n in range(period,0,-1):
        correls.append(0)
    # oldest data (autocorrelation period) will be nan; remove it
    correls = np.array(correls[-np.count_nonzero(~np.isnan(correls)):])      

    return correls

# CORRELATION OF BEST FIT
# the highest value correlation    
max_value = np.max(correls)
# index of the best correlation
max_index = np.argmax(correls)

Используя преобразование Фурье и теорему о свертке

Сложность времени составляет N * log (N)

def autocorr1(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    return r2[:len(x)//2]

Вот нормализованная и беспристрастная версия, это тоже N * log (N)

def autocorr2(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    c=(r2/x.shape-np.mean(x)**2)/np.std(x)**2
    return c[:len(x)//2]

Метод, предложенный А. Леви, работает, но я тестировал его на своем ПК, его временная сложность составляет N * N

def autocorr(x):
    result = numpy.correlate(x, x, mode='full')
    return result[result.size/2:]

Альтернатива numpy.correlate доступна в statsmodels.tsa.stattools.acf () . Это дает непрерывно убывающую функцию автокорреляции, подобную описанной OP. Реализовать это довольно просто:

from statsmodels.tsa import stattools
# x = 1-D array
# Yield normalized autocorrelation function of number lags
autocorr = stattools.acf( x )

# Get autocorrelation coefficient at lag = 1
autocorr_coeff = autocorr[1]

Поведение по умолчанию - остановка на 40 nlags, но это можно настроить с помощью nlag=параметра для вашего конкретного приложения. Внизу страницы есть ссылка на статистику функции .

Я думаю, что настоящий ответ на вопрос OP кратко содержится в этом отрывке из документации Numpy.correlate:

mode : {'valid', 'same', 'full'}, optional
    Refer to the `convolve` docstring.  Note that the default
    is `valid`, unlike `convolve`, which uses `full`.

Это означает, что при использовании без определения «режима» функция Numpy.correlate вернет скаляр, если ему задан один и тот же вектор для двух входных аргументов (то есть - когда используется для выполнения автокорреляции).

Простое решение без панд:

import numpy as np

def auto_corrcoef(x):
   return np.corrcoef(x[1:-1], x[2:])[0,1]

Постройте статистическую автокорреляцию с учетом серии результатов pandas datatime:

import matplotlib.pyplot as plt

def plot_autocorr(returns, lags):
    autocorrelation = []
    for lag in range(lags+1):
        corr_lag = returns.corr(returns.shift(-lag)) 
        autocorrelation.append(corr_lag)
    plt.plot(range(lags+1), autocorrelation, '--o')
    plt.xticks(range(lags+1))
    return np.array(autocorrelation)