Ik moet autocorrelatie van een reeks getallen doen, wat, zoals ik het begrijp, gewoon de correlatie is van de reeks met zichzelf.
Ik heb het geprobeerd met de correlaatfunctie van numpy, maar ik geloof het resultaat niet, omdat het bijna altijd een vector geeft waarbij het eerste getal niethet grootste is, zoals het zou moeten zijn .
Dus deze vraag bestaat eigenlijk uit twee vragen:
- Wat is
numpy.correlateaan het doen? - Hoe kan ik het (of iets anders) gebruiken om autocorrelatie te doen?
Antwoord 1, autoriteit 100%
Om je eerste vraag te beantwoorden: numpy.correlate(a, v, mode)voert de convolutie uit van amet het omgekeerde van ven het geven van de resultaten geknipt door de gespecificeerde modus. De definitie van convolutie, C(t)=? -? < ik < ? aivt+iwaar -? < t < ?, zorgt voor resultaten van -? naar ?, maar je kunt natuurlijk geen oneindig lange array opslaan. Het moet dus worden geknipt, en dat is waar de modus binnenkomt. Er zijn 3 verschillende modi: vol, hetzelfde, & geldig:
- “volledige” modus geeft resultaten voor elke
twaarbij zowelaalsvenige overlap hebben. - “dezelfde” modus geeft een resultaat met dezelfde lengte als de kortste vector (
aofv). - ‘valid’-modus geeft alleen resultaten als
aenvelkaar volledig overlappen. De documentatievoornumpy.convolvegeeft meer details over de modi.
Voor je tweede vraag, ik denk dat numpy.correlateje de autocorrelatie geeft, het geeft je ook iets meer. De autocorrelatie wordt gebruikt om te bepalen hoe vergelijkbaar een signaal of functie is met zichzelf op een bepaald tijdsverschil. Bij een tijdsverschil van 0 zou de autocorrelatie het hoogst moeten zijn omdat het signaal identiek is aan zichzelf, dus je verwachtte dat het eerste element in de autocorrelatieresultaatarray het grootst zou zijn. De correlatie begint echter niet bij een tijdsverschil van 0. Het begint bij een negatief tijdsverschil, sluit bij 0 en wordt dan positief. Dat wil zeggen, u verwachtte:
autocorrelatie(a) = ? -? < ik < ? aivt+iwaarbij 0 <= t < ?
Maar wat je kreeg was:
autocorrelatie(a) = ? -? < ik < ? aivt+iwaar -? < t < ?
Wat u moet doen, is de laatste helft van uw correlatieresultaat nemen, en dat zou de autocorrelatie moeten zijn waarnaar u op zoek bent. Een eenvoudige python-functie om dat te doen zou zijn:
def autocorr(x):
result = numpy.correlate(x, x, mode='full')
return result[result.size/2:]
Je zult natuurlijk een foutcontrole nodig hebben om er zeker van te zijn dat xecht een 1-d array is. Ook is deze verklaring waarschijnlijk niet de meest wiskundig rigoureuze. Ik heb oneindigheden rondgegooid omdat de definitie van convolutie ze gebruikt, maar dat geldt niet noodzakelijk voor autocorrelatie. Dus het theoretische gedeelte van deze uitleg is misschien wat wankel, maar hopelijk zijn de praktische resultaten nuttig. Dezepagina’sover autocorrelatie zijn erg nuttig en kunnen je een veel betere theoretische achtergrond geven als je het niet erg vindt om door de notatie en zware concepten te bladeren.
Antwoord 2, autoriteit 25%
Autocorrelatie is er in twee versies: statistisch en convolutie. Ze doen allebei hetzelfde, behalve een klein detail: de statistische versie is genormaliseerd op het interval [-1,1]. Hier is een voorbeeld van hoe je de statistische doet:
def acf(x, length=20):
return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:])[0,1] \
for i in range(1, length)])
Antwoord 3, autoriteit 19%
Gebruik de numpy.corrcoef-functie in plaats van numpy.correlateom de statistische correlatie voor een vertraging van t te berekenen:
def autocorr(x, t=1):
return numpy.corrcoef(numpy.array([x[:-t], x[t:]]))
Antwoord 4, autoriteit 17%
Ik denk dat er twee dingen zijn die dit onderwerp verwarrend maken:
- statistisch vs. signaalverwerkingsdefinitie: zoals anderen hebben opgemerkt, normaliseren we in statistieken autocorrelatie in [-1,1].
- gedeeltelijk vs. niet-partieel gemiddelde/variantie: wanneer de tijdreeks verschuift met een vertraging>0, zal hun overlapgrootte altijd < originele lengte. Gebruiken we het gemiddelde en standaard van het origineel (niet-gedeeltelijk), of berekenen we altijd een nieuw gemiddelde en std met behulp van de steeds veranderende overlap (gedeeltelijk) maakt een verschil. (Hier is waarschijnlijk een formele term voor, maar ik ga nu “gedeeltelijk” gebruiken).
Ik heb 5 functies gemaakt die de automatische correlatie van een 1d-array berekenen, met gedeeltelijke vs. onpartijdige onderscheidingen. Sommige gebruiken formules uit statistieken, andere gebruiken correleren in de zin van signaalverwerking, wat ook via FFT kan. Maar alle resultaten zijn autocorrelaties in de statistieken-definitie, dus ze illustreren hoe ze aan elkaar zijn gekoppeld. Code hieronder:
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()
Hier is het outputcijfer:
We zien niet alle 5 regels omdat 3 ervan elkaar overlappen (bij de paarse). De overlappingen zijn allemaal niet-gedeeltelijke autocorrelaties. Dit komt omdat berekeningen van de signaalverwerkingsmethoden (np.correlate, FFT) geen verschillende gemiddelde/std berekenen voor elke overlap.
Houd er ook rekening mee dat het resultaat fft, no padding, non-partial(rode lijn) anders is, omdat het de tijdreeks niet met nullen opvulde voordat FFT werd uitgevoerd, dus het is circulaire FFT. Ik kan niet in detail uitleggen waarom, dat heb ik elders geleerd.
Antwoord 5, autoriteit 10%
Omdat ik net hetzelfde probleem tegenkwam, wil ik graag een paar regels code met je delen. In feite zijn er inmiddels verschillende redelijk vergelijkbare berichten over autocorrelatie in stackoverflow. Als je de autocorrelatie definieert als a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2)[dit is de definitie gegeven in IDL’s a_correlate functie en het komt overeen met wat ik zie in antwoord 2 van vraag #12269834], dan lijkt het volgende de juiste resultaten te geven:
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()
Zoals je ziet heb ik dit getest met een sin-curve en een uniforme willekeurige verdeling, en beide resultaten zien eruit zoals ik ze zou verwachten. Merk op dat ik mode="same"gebruikte in plaats van mode="full"zoals de anderen deden.
Antwoord 6, autoriteit 9%
Uw vraag 1 is hier al uitgebreid besproken in verschillende uitstekende antwoorden.
Ik dacht om een paar regels code met je te delen waarmee je de autocorrelatie van een signaal kunt berekenen op basis van alleen de wiskundige eigenschappen van de autocorrelatie. Dat wil zeggen, de autocorrelatie kan op de volgende manier worden berekend:
-
trek het gemiddelde van het signaal af en verkrijg een onbevooroordeeld signaal
-
bereken de Fourier-transformatie van het zuivere signaal
-
bereken de spectrale vermogensdichtheid van het signaal door de vierkante norm te nemen van elke waarde van de Fourier-transformatie van het zuivere signaal
-
bereken de inverse Fourier-transformatie van de spectrale vermogensdichtheid
-
normaliseer de inverse Fourier-transformatie van de spectrale vermogensdichtheid met de som van de kwadraten van het zuivere signaal, en neem slechts de helft van de resulterende vector
De code om dit te doen is de volgende:
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)
Antwoord 7, autoriteit 6%
Een alternatief voor numpy.correlate is beschikbaar in statsmodels.tsa .stattools.acf(). Dit levert een continu afnemende autocorrelatiefunctie op zoals beschreven door OP. Het implementeren ervan is vrij eenvoudig:
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]
Het standaardgedrag is om te stoppen bij 40 nlags, maar dit kan worden aangepast met de optie nlag=voor uw specifieke toepassing. Er is een citaat onderaan de pagina voor de statistieken achter de functie.
Antwoord 8, autoriteit 2%
Ik ben een computerbioloog, en toen ik de auto-/kruiscorrelaties tussen paren van tijdreeksen van stochastische processen moest berekenen, realiseerde ik me dat np.correlateniet het werk deed dat ik nodig had .
Inderdaad, wat lijkt te ontbreken in np.correlateis het gemiddelde over alle mogelijke paren tijdpuntenop afstand ?? .
Hier is hoe ik een functie definieerde die deed wat ik nodig had:
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
Het lijkt me dat geen van de eerdere antwoorden betrekking heeft op dit geval van auto/cross-correlatie: ik hoop dat dit antwoord nuttig kan zijn voor iemand die aan stochastische processen werkt zoals ik.
Antwoord 9
Ik gebruik talib.CORREL voor autocorrelatie zoals deze, ik vermoed dat je hetzelfde zou kunnen doen met andere pakketten:
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)
Antwoord 10
Fouriertransformatie en de convolutiestelling gebruiken
De tijdscomplexiteit is N*log(N)
def autocorr1(x):
r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
return r2[:len(x)//2]
Hier is een genormaliseerde en onbevooroordeelde versie, het is ook 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]
De methode van A. Levy werkt, maar ik heb het op mijn pc getest, de complexiteit van de tijd lijkt N*N
te zijn
def autocorr(x):
result = numpy.correlate(x, x, mode='full')
return result[result.size/2:]
Antwoord 11
Ik denk dat het echte antwoord op de vraag van de OP beknopt is opgenomen in dit fragment uit de Numpy.correlate-documentatie:
mode : {'valid', 'same', 'full'}, optional
Refer to the `convolve` docstring. Note that the default
is `valid`, unlike `convolve`, which uses `full`.
Dit houdt in dat, wanneer gebruikt zonder ‘mode’-definitie, de Numpy.correlate-functie een scalair zal retourneren, wanneer dezelfde vector wordt gegeven voor zijn twee invoerargumenten (d.w.z. – wanneer gebruikt om autocorrelatie uit te voeren).
Antwoord 12
Een eenvoudige oplossing zonder panda’s:
import numpy as np
def auto_corrcoef(x):
return np.corrcoef(x[1:-1], x[2:])[0,1]
Antwoord 13
Plot de statistische autocorrelatie gegeven een panda’s datatime Reeks van rendementen:
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)