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

Ответ 1

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

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

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

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

Но у вас было:

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

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

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

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

Ответ 2

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

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

Ответ 3

Автокорреляция поставляется в двух версиях: статистическая и свертка. Они оба делают то же самое, за исключением небольшой детализации: статистическая версия нормирована на интервал [-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)])

Ответ 4

Как только я столкнулся с той же проблемой, я хотел бы поделиться с вами несколькими строками кода. На самом деле в настоящее время имеется несколько довольно похожих сообщений об автокорреляции в 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()

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

Ответ 5

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

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

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

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

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

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

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

Код для этого:

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)

Ответ 6

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

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

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

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 (красная линия) отличается, потому что он не заполнял временные ряды нулями перед выполнением FFT, поэтому он циклический FFT. Я не могу подробно объяснить, почему, то, что я узнал из других.

Ответ 7

Я использую 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)

Ответ 8

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

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

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

Ответ 9

Я вычислительный биолог, и когда мне нужно было вычислить авто /np.correlate корреляции между парами временных рядов случайных процессов, я понял, что 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

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

Ответ 10

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

import numpy as np

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

Ответ 11

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

Сложность времени 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:]

Ответ 12

Построить статистическую автокорреляцию с данными 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)

Ответ 13

Альтернатива 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= для вашего конкретного приложения. В нижней части страницы есть ссылка на статистику этой функции.