Numpy.cealrate를 사용하여 자기 상관을 수행하려면 어떻게해야합니까?

https://stackoverflow.com/questions/643699

22-07-2019
|

문제

나는 숫자 세트의 자동 상관을 수행해야합니다. 숫자는 세트 자체와의 상관 관계 일뿐입니다.

Numpy의 상관 관계 기능을 사용하여 시도했지만 결과는 거의 항상 벡터를 제공하기 때문에 결과를 믿지 않습니다. ~ 아니다 가장 큰 것입니다.

따라서이 질문은 실제로 두 가지 질문입니다.

정확히 무엇입니까 numpy.correlate 행위?
자동 상관을 수행하기 위해 어떻게 (또는 다른)를 사용할 수 있습니까?

해결책

첫 번째 질문에 답하기 위해 numpy.correlate(a, v, mode) 컨볼 루션을 수행하고 있습니다 a 반대와 함께 v 지정된 모드로 결과를 제공합니다. 그만큼 컨볼 루션의 정의, c (t) = ∑ _{-∞ <i <∞} ㅏ_나V_t+i 여기서 -∞ <t <∞는 -∞에서 ∞에서 ∞에서 결과를 허용하지만, 무한히 긴 배열을 저장할 수는 없습니다. 따라서 클리핑을해야하며 모드가 들어오는 곳입니다. 전체, 동일 및 유효한 3 가지 모드가 있습니다.

"전체"모드는 모든 결과를 반환합니다 t 둘 다 a 그리고 v 겹침이 있습니다.
"동일한"모드는 가장 짧은 벡터와 같은 길이의 결과를 반환합니다 (a 또는 v).
"유효한"모드는 다음과 같은 결과를 반환합니다 a 그리고 v 서로 완전히 겹칩니다. 그만큼 선적 서류 비치 ~을 위한 numpy.convolve 모드에 대한 자세한 내용을 제공합니다.

두 번째 질문은 생각합니다 numpy.correlate ~이다 당신에게 자기 상관을 주면, 그것은 당신에게도 조금 더 줄 것입니다. 자기 상관은 특정 시차에서 신호 또는 기능이 얼마나 유사한지를 찾는 데 사용됩니다. 0의 시간 차이에서, 신호가 자체와 동일하기 때문에 자동 상관이 가장 높아야하므로 자기 상관 결과 배열의 첫 번째 요소가 가장 큰 것으로 예상했다. 그러나 상관 관계는 0의 시차에서 시작되지 않습니다. 음의 시차에서 시작하여 0으로 닫힌 다음 양수가됩니다. 즉, 당신은 다음을 기대하고있었습니다.

자기 상관 (a) = ∑ _{-∞ <i <∞} ㅏ_나V_t+i 여기서 0 <= t <∞

하지만 당신이 얻은 것은 다음과 같습니다.

자기 상관 (a) = ∑ _{-∞ <i <∞} ㅏ_나V_t+i 여기서 -∞ <t <∞

당신이해야 할 일은 상관 관계 결과의 마지막 절반을 취하는 것입니다. 그리고 그것은 당신이 찾고있는 자기 상관이어야합니다. 그렇게하는 간단한 파이썬 기능은 다음과 같습니다.

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

물론 당신은 x 실제로 1D 배열입니다. 또한이 설명은 아마도 수학적으로 가장 엄격하지 않을 것입니다. 컨볼 루션의 정의가 그것들을 사용하기 때문에 무한대를 던지고 있지만, 반드시 자기 상관에 적용되는 것은 아닙니다. 따라서이 설명의 이론적 부분은 약간 원래가 될 수 있지만 실제 결과가 도움이되기를 바랍니다. 이것들 페이지 자기 상관은 매우 도움이되며 표기법과 무거운 개념을 넘어서는 것을 신경 쓰지 않으면 훨씬 더 나은 이론적 배경을 줄 수 있습니다.

다른 팁

자동 상관은 통계 및 컨볼 루션이라는 두 가지 버전으로 제공됩니다. 약간의 세부 사항을 제외하고는 모두 동일하게 수행합니다. 통계 버전은 간격으로 정규화됩니다 [-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) 이것은 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]로 정상화합니다.
부분 대 비 당면 평균/분산 : Timeseries가 LAG> 0으로 이동하면 오버랩 크기는 항상 <원래 길이입니다. 우리는 원본의 평균과 STD를 사용하거나 항상 변화하는 오버랩 (부분적)을 사용하여 새로운 평균 및 성병을 계산하는 것이 차이를 만듭니다. (아마도 이것에 대한 공식적인 용어가 있지만 지금은 "부분적"을 사용할 것입니다).

부분 대 비평가 구분과 함께 1D 배열의 자동 상관을 계산하는 5 가지 기능을 만들었습니다. 일부는 통계에서 공식을 사용하며 일부 사용은 신호 처리 의미에서 상관 관계가 있으며 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 줄이 모두 보라색으로 겹치기 때문에 5 줄이 모두 보이지 않습니다. 오버랩은 모두 비평가 자동 상관입니다. 신호 처리 방법의 계산이 있기 때문입니다 (np.correlate, fft) 각 오버랩마다 다른 평균/std를 계산하지 마십시오.

또한 fft, no padding, non-partial (빨간 선) 결과는 FFT를 수행하기 전에 타임시를 0으로 채우지 않았기 때문에 원형 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*로그 (n)

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

다음은 정규화되고 편견이없는 버전입니다. n*로그 (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]

A. Levy가 제공 한 방법은 작품이지만 PC에서 테스트했는데 시간 복잡성이있는 것 같습니다. n*n

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

나는 OP의 질문에 대한 진정한 답변이 Numpy.ceallate 문서에서 발췌 한 내용에 간결하게 포함된다고 생각합니다.

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

이는 '모드'정의가 없을 때 Numpy.cealrate 함수가 두 개의 입력 인수에 대해 동일한 벡터 (예 : 자기 상관을 수행 할 때)가 주어지면 스칼라를 반환한다는 것을 의미합니다.

팬더가없는 간단한 솔루션 :

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)

나는 계산 생물 학자이며, 시계열의 확률 론적 프로세스 커플 사이의 자동/교차 상관을 계산해야 할 때 나는 그것을 깨달았습니다. 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

이 자동/교차 상관 의이 인스턴스를 다루는 이전 답변 중 어느 것도 나에게 보이지 않는 것 같습니다.

Numpy.cealrate에 대한 대안이 있습니다 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에서 중지하는 것이지만 이는 다음과 같이 조정할 수 있습니다. nlag= 특정 응용 프로그램의 옵션. 페이지 하단에 인용이 있습니다. 함수의 통계.

라이센스 : CC-BY-SA ~와 함께 속성

제휴하지 않습니다 StackOverflow