Мне нужно сделать автокорреляцию набора чисел, который, как я понимаю, является просто корреляцией набора с самим собой.
Я попробовал это, используя функцию корреляции numpy, но я не верю результату, так как он почти всегда дает вектор, где первое число не самое большое, как и должно быть.
Итак, этот вопрос на самом деле два вопроса:
numpy.correlate?
Чтобы ответить на ваш первый вопрос, 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-м массивом. Кроме того, это объяснение, вероятно, не является наиболее математически строгим. Я обманывал бесконечности, потому что определение свертки использует их, но это не обязательно относится к автокорреляции. Таким образом, теоретическая часть этого объяснения может быть слегка отвратительной, но, надеюсь, практические результаты полезны. Эти страницы по автокорреляции очень полезны и могут дать вам гораздо лучше теоретического фона, если вы не против пробираться через нотацию и тяжелые понятия.
Используйте функцию numpy.corrcoef вместо numpy.correlate для расчета статистической корреляции для задержки t:
def autocorr(x, t=1):
return numpy.corrcoef(numpy.array([x[:-t], x[t:]]))
Автокорреляция поставляется в двух версиях: статистическая и свертка. Они оба делают то же самое, за исключением небольшой детализации: статистическая версия нормирована на интервал [-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)])
Как только я столкнулся с той же проблемой, я хотел бы поделиться с вами несколькими строками кода. На самом деле в настоящее время имеется несколько довольно похожих сообщений об автокорреляции в 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", как это делали другие.
Ваш вопрос 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)
Я думаю, что есть две вещи, которые добавляют путаницу к этой теме:
Я создал 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. Я не могу подробно объяснить, почему, то, что я узнал из других.
Я использую 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)
Я думаю, что реальный ответ на вопрос 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 вернет скаляр, если ему задан один и тот же вектор для двух входных аргументов (т.е. - когда используется для выполнения автокорреляции).
Я вычислительный биолог, и когда мне нужно было вычислить авто /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
Мне кажется, что ни один из предыдущих ответов не охватывает этот случай авто/взаимной корреляции: надеюсь, этот ответ может быть полезен для кого-то, работающего над случайными процессами, такими как я.
Простое решение без панд:
import numpy as np
def auto_corrcoef(x):
return np.corrcoef(x[1:-1], x[2:])[0,1]
Использование преобразования Фурье и теорема о свертке
Сложность времени 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:]
Построить статистическую автокорреляцию с данными 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)
Альтернатива 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= для вашего конкретного приложения. В нижней части страницы есть ссылка на статистику этой функции.