Vraag:
Periodedetectie van een generieke tijdreeks
gianluca
2010-08-04 05:32:13 UTC
view on stackexchange narkive permalink

Dit bericht is de voortzetting van een ander bericht met betrekking tot een generieke methode voor het opsporen van uitbijter in tijdreeksen. Op dit punt ben ik in wezen geïnteresseerd in een robuuste manier om de periodiciteit / seizoensgebondenheid van een generieke tijdreeks die wordt beïnvloed door veel ruis.Vanuit het oogpunt van een ontwikkelaar zou ik graag een eenvoudige interface willen, zoals:

unsigned int discovery_period (vector<double> v);

Waar v de array is die de samples bevat, en de geretourneerde waarde de periode van het signaal. Het belangrijkste punt is dat ik, nogmaals, geen aanname kan doen met betrekking tot de geanalyseerd signaal. ik heb al een benadering geprobeerd op basis van de signaalautocorrelatie (het detecteren van de pieken van een correlogram), maar het is niet robuust zoals ik zou willen.

Heb je xts :: periodicity geprobeerd?
Zeven antwoorden:
#1
+53
Rob Hyndman
2010-08-04 10:41:03 UTC
view on stackexchange narkive permalink

Als je echt geen idee hebt wat de periodiciteit is, is de beste benadering waarschijnlijk om de frequentie te vinden die overeenkomt met het maximum van de spectrale dichtheid. Het spectrum bij lage frequenties wordt echter beïnvloed door de trend, dus u moet eerst de reeks afbreken. De volgende R-functie zou het werk moeten doen voor de meeste series. Het is verre van perfect, maar ik heb het op een paar dozijn voorbeelden getest en het lijkt goed te werken. Het retourneert 1 voor gegevens die geen sterke periodiciteit hebben, en anders de lengte van de periode.

Update: versie 2 van de functie. Dit is veel sneller en lijkt robuuster te zijn.

  find.freq <- functie (x) {n <- lengte (x) specificatie <- spec.ar (c (x), plot = FALSE) if (max (spec $ spec) >10) # Willekeurige drempel gekozen door vallen en opstaan. {period <- round (1 / spec $ freq [which.max (spec $ spec)]) if (period == Inf) # Vind het volgende lokale maximum {j <- which (diff (spec $ spec) >0) if ( lengte (j) >0) {nextmax <- j [1] + which.max (spec $ spec [j [1]: 500]) periode <- round (1 / spec $ freq [nextmax])} anders periode <- 1}} else period <- 1 return (period)}  
Dank je. Nogmaals, ik zal deze aanpak zo snel mogelijk proberen en zal hier de definitieve resultaten schrijven.
Uw idee is redelijk goed, maar in mijn geval slaagt het er niet in de periodiciteit te detecteren van een heel eenvoudige (en niet zo luidruchtige) tijdserie zoals http://dl.dropbox.com/u/540394/chart.png. Met mijn "empirische" benadering (gebaseerd op de autocorrelatie), geeft het eenvoudige algoritme dat ik schreef een exacte periode van 1008 terug (met elke 10 minuten een steekproef, dit betekent 1008/24/6 = 7, dus een wekelijkse periodiciteit). Mijn belangrijkste problemen zijn: 1) het is te traag om samen te komen (het vereist veel historische gegevens) en ik heb een reactieve, online benadering nodig; 2) het is inefficiënt als de hel vanuit het oogpunt van geheugengebruik; 3) het is niet robuust alle;
Dank je. Helaas werkt dit nog steeds niet zoals ik zou verwachten. Voor dezelfde tijdreeks van de vorige opmerking wordt 166 geretourneerd, wat slechts gedeeltelijk juist is (naar mijn mening is de duidelijke wekelijkse periode interessanter). En met behulp van een zeer luidruchtige tijdreeks, zoals deze http://dl.dropbox.com/u/540394/chart2.png (een vensteranalyse van een TCP-ontvanger), retourneert de functie 10, terwijl ik 1 zou verwachten (ik kan ' zie geen duidelijke periodiciteit). Trouwens, ik weet dat het erg moeilijk zal zijn om te vinden wat ik zoek, aangezien ik met te verschillende signalen te maken heb.
166 is geen slechte schatting van 168. Als u weet dat de gegevens elk uur met een wekelijks patroon worden waargenomen, waarom zou u dan de frequentie überhaupt schatten?
Omdat ik veel tijdreeksen moet analyseren (stel dat er 100 netwerkstatistieken zijn), en slechts enkele daarvan hebben een wekelijkse periodiciteit. In elk geval denk ik dat ik in mijn implementatie een algoritme zal gebruiken dat lijkt op jouw functie, en ik zal handmatig de wekelijkse periodiciteit onderscheiden. Echt bedankt voor je steun, ik waardeer het (en ga zo door met het goede werk met de prognosebibliotheek :-))
Ik heb deze functie getest met een eenvoudig voorbeeld: x = c (58.89446, 37.31097, 53.99865, 26.13904, 34.74298) en y = ts (rep_len (x, 15 * lengte (x)). Met behulp van de bovenstaande definitie verwachtte ik15 als find.freq (y) (of iets dichtbij), maar ik krijg er 3. Wat mis ik hier?
Waarom zou je dit niet in een pakket opnemen?Er zijn veel taken waarvan de periodiciteit onbekend is.
Een verbeterde versie zit in het voorspellingspakket als `vindfrequentie`
#2
+10
Rich
2010-08-10 23:41:11 UTC
view on stackexchange narkive permalink

Als u verwacht dat het proces stationair verloopt - de periodiciteit / seizoensinvloeden zullen in de loop van de tijd niet veranderen - dan kan zoiets als een Chi-kwadraat-periodogram (zie bijv.Sokolove en Bushell, 1978) een goede keuze zijn. Het wordt vaak gebruikt bij de analyse van circadiane gegevens die extreem grote hoeveelheden ruis kunnen bevatten, maar naar verwachting zeer stabiele periodiciteiten hebben.

Deze benadering gaat niet uit van de vorm van de golfvorm (behalve dat het is consistent van cyclus tot cyclus), maar vereist wel dat elke ruis constant gemiddeld is en niet gecorreleerd is met het signaal.

  chisq.pd <- functie (x, min. periode, max. periode, alfa) {N <- lengte (x) varianties = NULLperiods = seq (min. periode, max. periode) rowlist = NULLfor (lc in periodes) {ncol = lc nrow = floor (N / ncol) rowlist = c ( rowlist, nrow) x.trunc = x [1: (ncol * nrow)] x.reshape = t (matrix (x.trunc, c (ncol, nrow))) varianties = c (varianties, var (colMeans (x. reshape)))} Qp = (rowlist * periodes * varianties) / var (x) df = periodes - 1pvals = 1-pchisq (Qp, df) pass.periods = periodes [pvals<alpha] pass.pvals = pvals [pvals<alpha] # return (cbind (pass.periods, pass.pvals)) return (cbind (periodes [pvals == min (pvals)], pvals [pvals == min (pvals)]))} x = cos ((2 * pi / 37) * (1: 1000)) + rnorm (1000) chisq.pd (x, 2, 72, .05)  

De laatste twee regels zijn slechts een voorbeeld en laten zien dat het de periode van een pure trigonometrische functie kan identificeren, zelfs met veel additieve ruis.

Zoals geschreven, is het laatste argument ( alpha ) in de aanroep overbodig, de functie geeft simpelweg de 'beste' periode terug die ze kan vinden; verwijder commentaar op de eerste return -instructie en maak commentaar op de tweede om een ​​lijst te laten zien van alle periodes die significant zijn op het niveau alpha .

Deze functie doet geen enkele vorm van gezondheidscontrole om er zeker van te zijn dat je identificeerbare periodes hebt ingevoerd, noch werkt het (kan het) werken met fractionele periodes, noch is er enige vorm van meervoudige vergelijkingscontrole ingebouwd als je besluit om naar meerdere periodes te kijken. Maar verder moet het redelijk robuust zijn.

Ziet er interessant uit, maar ik begrijp de output niet, het vertelt me niet waar de periode begint, en de meeste p-waarden van 1.
#3
+4
Wesley Burr
2010-08-06 07:48:10 UTC
view on stackexchange narkive permalink

Misschien wilt u duidelijker definiëren wat u wilt (voor uzelf, zo niet hier). Als u op zoek bent naar de meest statistisch significante stationaire periode in uw gegevens met veel ruis, zijn er in wezen twee routes die u kunt nemen:

1) bereken een robuuste autocorrelatieschatting en neem de maximale coëfficiënt
2) bereken een robuuste schatting van de spectrale vermogensdichtheid en neem het maximale uit het spectrum

Het probleem met # 2 is dat je voor elke lawaaierige tijdreeks een grote hoeveelheid stroom krijgt in lage frequenties, waardoor het is moeilijk te onderscheiden. Er zijn enkele technieken om dit probleem op te lossen (bijv. Voor wit maken, dan de PSD schatten), maar als de werkelijke periode van uw gegevens lang genoeg is, zal automatische detectie twijfelachtig zijn.

Waarschijnlijk is uw beste gok om een ​​robuuste autocorrelatieroutine te implementeren, zoals te vinden in hoofdstuk 8.6, 8.7 in Robuuste Statistieken - Theorie en Methoden door Maronna, Martin en Yohai. Zoeken op Google naar "robuuste durbin-levinson" zal ook enkele resultaten opleveren.

Als u alleen op zoek bent naar een eenvoudig antwoord, weet ik niet zeker of er een bestaat. Periodedetectie in tijdreeksen kan ingewikkeld zijn, en het kan te veel zijn om een ​​geautomatiseerde routine te vragen die magie kan uitvoeren.

Bedankt voor je kostbare informatie, ik zal zeker naar dat boek kijken.
#4
+4
babelproofreader
2010-08-10 22:29:28 UTC
view on stackexchange narkive permalink

Je zou de Hilbert Transformation from DSP-theorie kunnen gebruiken om de momentane frequentie van je gegevens te meten. De site http://ta-lib.org/ heeft open source-code voor het meten van de dominante cyclusperiode van financiële gegevens; de relevante functie heet HT_DCPERIOD; u kunt dit wellicht gebruiken of de code aanpassen aan uw doeleinden.

#5
+3
Fabrizio Maccallini
2016-12-29 19:15:47 UTC
view on stackexchange narkive permalink

Een andere benadering zou Empirical Mode Decomposition kunnen zijn.Het R-pakket heet EMD en is ontwikkeld door de uitvinder van de methode:

  vereisen (EMD) ndata <- 3000 tt2 <- seq (0, 9, length =ndata) xt2 <- sin (pi * tt2) + sin (2 * pi * tt2) + sin (6 * pi * tt2) + 0,5 * tt2 probeer <- emd (xt2, tt2, boundary = "wave") ### Het uitzetten van de par van het IMF (mfrow = c (probeer $ nimf + 1, 1), mar = c (2,1,2,1)) rangeimf <- bereik (probeer $ imf) voor (i in 1: probeer $ nimf) {plot (tt2, probeer $ imf [, i], type = "l", xlab = "", ylab = "", ylim = rangeimf, main = paste (i, "-th IMF", sep = ""));abline (h = 0)} plot (tt2, probeer $ residu, xlab = "", ylab = "", main = "residu", type = "l", axes = FALSE);box ()  

De methode werd om een goede reden 'Empirisch' genoemd en het risico bestaat dat de intrinsieke modusfuncties (de individuele additieve componenten) door elkaar worden gehaald.Aan de andere kant is de methode erg intuïtief en kan ze nuttig zijn voor een snelle visuele inspectie van cycliciteit.

#6
  0
Chris
2015-05-02 20:24:14 UTC
view on stackexchange narkive permalink

Met betrekking tot het bericht van Rob Hyndman hierboven https://stats.stackexchange.com/a/1214/70282

De functie find.freq werkt uitstekend. Op de dagelijkse dataset die ik gebruik, is de frequentie correct berekend op 7.

Toen ik het alleen op de weekdagen probeerde, vermeldde het dat de frequentie 23 is, wat opmerkelijk dicht bij 21 ligt.42857 = 29,6 * 5/7 wat het gemiddelde aantal werkdagen in een maand is. (Of omgekeerd 23 * 7/5 is 32.)

Terugkijkend naar mijn dagelijkse gegevens, experimenteerde ik met het voorgevoel van het nemen van de eerste periode, daarmee een gemiddelde nemen en dan de volgende periode zoeken, enz. Zie hieronder:

 find.freq.all = functie (x) {f = find.freq (x); freqs = c (f); while (f> 1) {start = 1; # probeer ook start = f; x = period.apply (x, seq (start, lengte (x), f), gemiddelde); f = find.freq (x); freqs = c (freqs, f); } if (length (freqs) == 1) {return (freqs); } voor (i in 2: lengte (freqs)) {freqs [i] = freqs [i] * freqs [i-1]; } freqs [1: (length (freqs) -1)];} find.freq.all (dailyts) # using daily data 

Het bovenstaande geeft (7,28) of (7,35) afhankelijk aan als de seq begint met 1 of f. (Zie opmerking hierboven.)

Wat zou impliceren dat de seizoensperioden voor msts (...) (7,28) of (7,35) zouden moeten zijn.

De logica lijkt gevoelig voor beginvoorwaarden gezien de gevoeligheid van de algoritme parameters. Het gemiddelde van 28 en 35 is 31,5, wat dicht bij de gemiddelde lengte van een maand ligt.

Ik vermoed dat ik het wiel opnieuw heb uitgevonden, wat is de naam van dit algoritme? Is er ergens een betere implementatie in R?

Later heb ik de bovenstaande code uitgevoerd om alle starts van 1 tot en met 7 te proberen en ik kreeg 35,35,28,28,28,28,28 voor de tweede periode. Het gemiddelde komt uit op 30, het gemiddelde aantal dagen in een maand. Interessant ...

Eventuele gedachten of opmerkingen?

#7
  0
ali
2016-09-27 17:10:59 UTC
view on stackexchange narkive permalink

Men kan ook de Ljung-Box-test gebruiken om erachter te komen welk seizoensverschil de beste stationariteit bereikt. Ik werkte aan een ander onderwerp en ik gebruikte dit eigenlijk voor dezelfde doeleinden. Probeer verschillende periodes, zoals 3 tot 24 voor maandelijkse gegevens. En test ze allemaal met Ljung-Box en bewaar de Chi-Square-resultaten. En kies de periode met de laagste chikwadraatwaarde.

Hier is een eenvoudige code om dat te doen.

  minval0 <- 5000 # wijs een groot getal toe om er zeker van te zijn dat de Chi-waarden klein zijnerminindex0 <- 0periyot <- 0for (i in 3:24) {#vind optimale periode door Qtests over originele gegevens d0D1 <- diff (a, lag = i) #opslagresultaten Qtest_d0D1 [[i]] <- Box.test (d0D1, lag = 20, type = "Ljung-Box") #store Chi-Square statistieken sira0 [i] <- Qtest_d0D1 [[i]] [1]} # draai lijst om naar een dataframe, dan matrixdatam0 <- data.frame (matrix (unlist (sira0), nrow = lengte (Qtest_d0D1) -2, byrow = T)) datamtrx0 <- as.matrix (datam0 []) # haal min waarde's indexminindex0 <- welke (datamtrx0 == min (datamtrx0), arr. ind = F) periyot <- minindex0 + 2  


Deze Q&A is automatisch vertaald vanuit de Engelse taal.De originele inhoud is beschikbaar op stackexchange, waarvoor we bedanken voor de cc by-sa 2.0-licentie waaronder het wordt gedistribueerd.
Loading...