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

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

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

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

  1. Что именно делает numpy.correlate?
  2. Как я могу использовать его (или что-то еще) для автокорреляции?
81 голос | спросил Ben 13 MarpmFri, 13 Mar 2009 20:07:53 +03002009-03-13T20:07:53+03:0008 2009, 20:07:53

12 ответов


0

Чтобы ответить на ваш первый вопрос, numpy.correlate(a, v, mode) выполняет свертку a с обратным v и выдачей результатов, обрезанных в указанном режиме. определение свертки , C (t) = ∑ -∞ <я <∞ a i v t + i , где -∞ <т <∞, позволяет получать результаты от -∞ до ∞, но вы, очевидно, не можете хранить бесконечно длинный массив. Таким образом, он должен быть обрезан, и именно здесь режим входит. Есть 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 , где -∞ <т <∞

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

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

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

ответил A. Levy 24 MaramTue, 24 Mar 2009 09:09:02 +03002009-03-24T09:09:02+03:0009 2009, 09:09:02
0

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

def autocorr(x, t=1):
    return numpy.corrcoef(numpy.array([x[0:len(x)-t], x[t:len(x)]]))
ответил Ramon J Romero y Vigil 18 FebruaryEurope/MoscowbTue, 18 Feb 2014 23:10:37 +0400000000pmTue, 18 Feb 2014 23:10:37 +040014 2014, 23:10:37
0

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

def acf(x, length=20):
    return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:]) \
        for i in range(1, length)])
ответил jonathf 2 32011vEurope/Moscow11bEurope/MoscowWed, 02 Nov 2011 17:27:37 +0400 2011, 17:27:37
0

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

ответил maschu 13 J0000006Europe/Moscow 2013, 18:52:10
0

Ваш вопрос 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)
ответил Ruggero 20 +03002016-10-20T15:46:02+03:00312016bEurope/MoscowThu, 20 Oct 2016 15:46:02 +0300 2016, 15:46:02
0

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

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

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

ответил Jason 4 J000000Wednesday18 2018, 10:32:40
0

Я использую 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)
ответил litepresence 22 TueEurope/Moscow2015-12-22T08:06:48+03:00Europe/Moscow12bEurope/MoscowTue, 22 Dec 2015 08:06:48 +0300 2015, 08:06:48
0

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

import numpy as np

def auto_corrcoef(x):
   return np.corrcoef(x[1:-1], x[2:])[0,1]
ответил dignitas 23 J000000Monday18 2018, 16:22:38
0

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

ответил dbanas 2 FebruaryEurope/MoscowbMon, 02 Feb 2015 00:17:42 +0300000000amMon, 02 Feb 2015 00:17:42 +030015 2015, 00:17:42
0

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

Действительно, в np.correlate, по-видимому, отсутствует усреднение по всем возможным парам временных точек на расстоянии \ tau.

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

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

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

ответил Orso 12 PMpThu, 12 Apr 2018 18:25:30 +030025Thursday 2018, 18:25:30
0

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

Сложность времени 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:]
ответил wwwjjj 17 thEurope/Moscowp30Europe/Moscow09bEurope/MoscowMon, 17 Sep 2018 09:20:55 +0300 2018, 09:20:55
0

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

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)
ответил Antonio Catalano 8 +03002018-10-08T15:42:32+03:00312018bEurope/MoscowMon, 08 Oct 2018 15:42:32 +0300 2018, 15:42:32

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132