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

https://stackoverflow.com/questions/643699

22-07-2019
|

Вопрос

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

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

Итак, этот вопрос на самом деле состоит из двух вопросов:

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

Решение

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

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

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

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

Но то, что ты получил, было:

автокорреляция (a) = ∑ _{-∞ < i < ∞} a_iv_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:]]))

Поскольку я столкнулся с той же проблемой, я хотел бы поделиться с вами несколькими строками кода. На самом деле на данный момент есть несколько довольно похожих сообщений об автокорреляции в stackoverflow. Если вы определяете автокорреляцию как a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2) [это определение дано в функции a_correlate IDL и оно согласуется с тем, что я вижу в ответе 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", как и другие.

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

<Ол>

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

частичное против Неполное среднее значение / дисперсия: когда временной ряд сдвигается с задержкой > 0, их размер перекрытия всегда будет < оригинальная длина Используем ли мы среднее и стандартное отклонение оригинала (не частичное), или всегда вычисляем новое среднее значение, и стандартное отклонение, используя постоянно меняющееся перекрытие (частичное), имеет значение. (Вероятно, есть формальный термин для этого, но я собираюсь использовать & Quot; частичный & Quot; пока).

Я создал 5 функций, которые вычисляют автокорреляцию массива 1d с частичным v.s. не частичные различия. Некоторые используют формулу из статистики, некоторые используют коррелят в смысле обработки сигнала, что также может быть сделано через 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. Я не могу подробно объяснить, почему, это то, что я узнал из других источников.

Ваш вопрос 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)

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

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

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

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

import numpy as np

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

Постройте статистическую автокорреляцию с учетом панды дата-время. Серия возвратов:

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)

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

Действительно, в <=>, похоже, отсутствует усреднение по всем возможным парам временных точек на расстоянии & # 120591;.

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

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

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

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

Лицензировано под: CC-BY-SA с атрибуция

Не связан с StackOverflow