Pitanje:
Razdoblje otkrivanja generičke vremenske serije
gianluca
2010-08-04 05:32:13 UTC
view on stackexchange narkive permalink

Ovaj je post nastavak drugog posta koji se odnosi na generičku metodu za otkrivanje odstupanja u vremenskim serijama. U osnovi, u ovom trenutku zanima me robustan način otkrivanja periodičnosti / sezonalnosti generičku vremensku seriju na koju utječe puno buke. S gledišta programera, želio bih jednostavno sučelje kao što je:

unsigned int Discover_period (vector<double> v);

Gdje je v niz koji sadrži uzorke, a povratna vrijednost je period signala. Glavna poanta je da, opet, ne mogu pretpostaviti vezano uz analizirao sam signal.Već sam probao pristup zasnovan na autokorelaciji signala (otkrivanje vrhova korelograma), ali nije robustan kao što bih želio.

Jeste li probali xts :: periodicity?
Sedam odgovori:
#1
+53
Rob Hyndman
2010-08-04 10:41:03 UTC
view on stackexchange narkive permalink

Ako stvarno nemate pojma koja je periodičnost, vjerojatno je najbolji pristup pronaći frekvenciju koja odgovara maksimumu spektralne gustoće. Međutim, trend će utjecati na spektar na niskim frekvencijama, pa prvo morate ukloniti seriju. Sljedeća funkcija R trebala bi obaviti posao za većinu serija. Daleko je od savršenog, ali testirao sam ga na nekoliko desetaka primjera i čini mi se da djeluje u redu. Vratit će 1 za podatke koji nemaju jaku periodičnost, a u suprotnom trajanje razdoblja.

Ažuriranje: Verzija 2 funkcije. Ovo je mnogo brže i čini se da je robusnije.

  find.freq <- funkcija (x) {n <- length (x) spec <- spec.ar (c (x), plot = FALSE) if (max (spec $ spec) >10) # Proizvoljni prag odabran metodom pokušaja i pogrešaka. {period <- krug (1 / spec $ freq [koji.max (spec $ spec)]]) if (period == Inf) # Pronađi sljedeći lokalni maksimum {j <- koji (razlika (spec $ spec) >0) if ( dužina (j) >0) {nextmax <- j [1] + which.max (spec. $ spec [j [1]: 500]) razdoblje <- okruglo (1 / spec $ freq [nextmax])} ostalo razdoblje <- 1}} else period <- 1 povratak (period)}  
Hvala vam. Ponovno ću pokušati s ovim pristupom što je prije moguće i ovdje ću napisati konačne rezultate.
Vaša je ideja prilično dobra, ali u mom slučaju ne uspijeva otkriti periodičnost doista jednostavne (i ne tako bučne) vremenske serije poput http://dl.dropbox.com/u/540394/chart.png. Uz moj "empirijski" pristup (zasnovan na autokorelaciji), jednostavni algoritam koji sam napisao vraća točno razdoblje od 1008 (imajući uzorak svakih 10 minuta, to znači 1008/24/6 = 7, dakle tjednu periodičnost). Moji glavni problemi su: 1) Presporo se konvergira (zahtijeva puno povijesnih podataka) i potreban mi je reaktivni mrežni pristup; 2) Neučinkovit je s gledišta korištenja memorije; 3) Nije robustan u svi;
Hvala vam. Nažalost, ovo još uvijek ne funkcionira kako bih očekivao. Za istu vremensku seriju prethodnog komentara vraća se 166, što je samo djelomično u pravu (s mog stajališta je očiglednije tjedno razdoblje zanimljivije). A pomoću vrlo bučne vremenske serije, poput ove http://dl.dropbox.com/u/540394/chart2.png (analiza prozora TCP prijemnika), funkcija vraća 10, dok bih očekivao 1 (mogu ' ne vidim nikakvu očitu periodičnost). BTW Znam da će biti jako teško pronaći ono što tražim, jer imam posla s previše različitim signalima.
166 nije loša procjena od 168. Ako znate da se podaci promatraju svaki sat tjedno, zašto onda uopće procjenjivati ​​učestalost?
Budući da moram analizirati puno vremenskih serija (pretpostavimo 100 mrežnih mjernih podataka), a samo neki od njih imaju tjednu periodičnost. U svakom slučaju, pretpostavljam da ću u svojoj implementaciji koristiti algoritam sličan vašoj funkciji i ručno ću razlikovati tjednu periodičnost. Doista hvala na podršci, zaista cijenim (i nastavite s dobrim radom s bibliotekom predviđanja :-))
Testirao sam ovu funkciju na jednostavnom primjeru: x = c (58.89446, 37.31097, 53.99865, 26.13904, 34.74298) i y = ts (rep_len (x, 15 * dužina (x)). Koristeći gornju definiciju, očekivao sam da ću dobiti15 kao find.freq (y) (ili nešto blisko), ali dobivam 3. Što mi ovdje nedostaje?
Zašto ovo ne uključiti u neki paket?Puno je zadataka kad je periodičnost nepoznata ..
Poboljšana verzija nalazi se u paketu predviđanja pod nazivom `findfrequency`
#2
+10
Rich
2010-08-10 23:41:11 UTC
view on stackexchange narkive permalink

Ako očekujete da će postupak biti stacionaran - periodičnost / sezonalnost neće se mijenjati s vremenom - tada bi nešto poput Chi-kvadrat periodograma (vidi npr. Sokolove i Bushell, 1978) mogao biti dobar izbor. Obično se koristi u analizi cirkadijskih podataka koji mogu sadržavati izuzetno velike količine buke, ali očekuje se da će imati vrlo stabilnu periodičnost.

Ovaj pristup ne pretpostavlja oblik valnog oblika (osim onog dosljedan je iz ciklusa u ciklus), ali zahtijeva da bilo koji šum bude konstantne srednje vrijednosti i da nije u korelaciji sa signalom.

  chisq.pd <- funkcija (x, min. razdoblje, maks. razdoblje, alfa) {N <- duljine (x) varijance = NULLperiods = seq (min.period, max.period) list row = NULLfor (lc u razdobljima) {ncol = lc nrow = floor (N / ncol) row row = c ( rowlist, nrow) x.trunc = x [1: (ncol * nrow)] x.reshape = t (array (x.trunc, c (ncol, nrow))) variance = c (variance, var (colMeans (x. preoblikovati)))} Qp = (popis redova * razdoblja * varijance) / var (x) df = razdoblja - 1pvali = 1-pchisq (Qp, df) pass.periods = razdoblja [pvals<alpha] pass.pvals = pvali [pvals<alpha] # return (cbind (pass.periods, pass.pvals)) return (cbind (razdoblja [pvali == min (pvali)], pvali [pvali == min (pvali)]))} x = cos ((2 * pi / 37) * (1: 1000)) + rnorm (1000) chisq.pd (x, 2, 72, .05)  

Posljednja dva retka samo su primjer koji pokazuje da može identificirati razdoblje čiste trigonometrijske funkcije, čak i s puno aditivnog šuma.

Kao što je napisano, posljednji argument ( alpha ) u pozivu je suvišan, funkcija jednostavno vraća 'najbolje' razdoblje koje može pronaći; raskomentirajte prvu izjavu return i komentirajte drugu kako bi se ona vratila na popis svih razdoblja značajnih na razini alpha .

Ova funkcija ne vrši bilo kakvu provjeru zdravstvenog stanja kako bi bila sigurna da ste unijeli razdoblja koja se mogu identificirati, niti (može li) raditi s razlomljenim razdobljima, niti je ugrađena bilo kakva kontrola višestruke usporedbe odlučili ste pogledati više razdoblja. Ali osim toga, trebao bi biti razumno robustan.

Izgleda zanimljivo, ali ne razumijem rezultate, ne govori mi odakle započinje razdoblje i većina vrijednosti 1.
#3
+4
Wesley Burr
2010-08-06 07:48:10 UTC
view on stackexchange narkive permalink

Možda biste željeli jasnije definirati što želite (sebi, ako ne ovdje). Ako je ono što tražite statistički najznačajnije stacionarno razdoblje sadržano u vašim bučnim podacima, u osnovi postoje dva puta:

1) izračunajte robusnu procjenu autokorelacije i uzmite maksimalni koeficijent
2) izračunajte robusnu procjenu spektralne gustoće snage i uzmite maksimum spektra

Problem br. 2 je taj što ćete za bilo koju bučnu vremensku seriju dobiti veliku količinu snage na niskim frekvencijama, što čini teško je razlikovati. Postoje neke tehnike za rješavanje ovog problema (tj. Prebijeliti, a zatim procijeniti PSD), ali ako je istinito razdoblje iz vaših podataka dovoljno dugo, automatsko otkrivanje bit će nedovoljno.

Vaša je najbolja oklada vjerojatno za provedbu robusne rutine autokorelacije kakvu možete pronaći u poglavlju 8.6, 8.7 u Robust Statistics - Theory and Methods autora Maronna, Martin i Yohai. Pretraživanje Googlea za "robusnim durbin-levinsonom" također će donijeti neke rezultate.

Ako tražite samo jednostavan odgovor, nisam siguran da postoji. Otkrivanje razdoblja u vremenskim serijama može biti komplicirano, a traženje automatizirane rutine koja može izvoditi magiju može biti previše.

Hvala vam na dragocjenim informacijama, sigurno ću pogledati tu knjigu.
#4
+4
babelproofreader
2010-08-10 22:29:28 UTC
view on stackexchange narkive permalink

Mogli biste upotrijebiti Hilbertovu transformaciju iz DSP teorije za mjerenje trenutne učestalosti vaših podataka. Web stranica http://ta-lib.org/ ima otvoreni kod za mjerenje dominantnog razdoblja ciklusa financijskih podataka; relevantna funkcija naziva se HT_DCPERIOD; možda ćete moći koristiti ovo ili prilagoditi kod svojim potrebama.

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

Drugačiji pristup mogao bi biti dekompozicija empirijskog načina.Paket R naziva se EMD koji je razvio izumitelj metode:

  require (EMD) ndata <- 3000 tt2 <- seq (0, 9, length =ndata) xt2 <- sin (pi * tt2) + sin (2 * pi * tt2) + sin (6 * pi * tt2) + 0,5 * tt2 pokušaj <- emd (xt2, tt2, border = "val") ### Ucrtavanje MMF-ovog nominalnog (mfrow = c (pokušaj $ nimf + 1, 1), mar = c (2,1,2,1)) rangeimf <- raspona (try $ imf) for (i in 1: try $ nimf) {plot (tt2, pokušajte $ imf [, i], type = "l", xlab = "", ylab = "", ylim = rangeimf, main = paste (i, "-ti MMF", sep = ""));abline (h = 0)} plot (tt2, pokušajte $ ostatak, xlab = "", ylab = "", main = "ostatak", type = "l", osi = FALSE);box ()  

Metoda je s dobrim razlogom označena kao "Empirijska" i postoji rizik da se funkcije unutarnjeg načina rada (pojedinačne komponente aditiva) pomiješaju.S druge strane, metoda je vrlo intuitivna i može biti korisna za brzi vizualni pregled cikličnosti.

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

U odnosu na gornji post Roba Hyndmana https://stats.stackexchange.com/a/1214/70282

Funkcija find.freq sjajno djeluje. Na dnevnom skupu podataka koji koristim ispravno je utvrdio da je učestalost 7.

Kad sam ga isprobao samo u tjednima, spomenuo je da je učestalost 23, što je izuzetno blizu 21.42857 = 29,6 * 5/7 što je prosječni broj radnih dana u mjesecu. (Ili obratno, 23 * 7/5 je 32.)

Osvrćući se na svoje dnevne podatke, eksperimentirao sam s naslućivanjem prvog razdoblja, usrednjavajući po tome, a zatim pronalazeći sljedeće razdoblje itd. Vidi dolje:

 find.freq.all = function (x) {f = find.freq (x); frekvencije = c (f); while (f> 1) {start = 1; # također pokušajte start = f; x = razdoblje.primjena (x, seq (početak, dužina (x), f), srednja vrijednost); f = pronađi.freq (x); frekvencije = c (freqs, f); } if (length (freqs) == 1) {return (freqs); } for (i u 2: dužina (freqs)) {freqs [i] = freqs [i] * freqs [i-1]; } freqs [1: (length (freqs) -1)];} find.freq.all (dailyts) # using daily data 

Gore navedeno daje (7,28) ili (7,35) ovisno na ako seq započinje s 1 ili f. (Vidi komentar gore.)

Što bi značilo da sezonska razdoblja za msts (...) trebaju biti (7,28) ili (7,35).

Logika čini se osjetljivim na početne uvjete s obzirom na osjetljivost parametara algoritma. Srednja vrijednost 28 i 35 je 31,5 što je približno prosječnoj duljini mjeseca.

Pretpostavljam da sam ponovno izmislio kotač, kako se zove ovaj algoritam? Postoji li negdje bolja implementacija u R-u?

Kasnije sam pokrenuo gornji kôd pokušavajući sve pokrete od 1 do 7, a za drugi sam dobio 35,35,28,28,28,28,28 razdoblje. Prosjek iznosi 30, što je prosječan broj dana u mjesecu. Zanimljivo ...

Imate li kakvih misli ili komentara?

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

Također se može koristiti Ljung-Boxov test kako bi se utvrdilo koja sezonska razlika doseže najbolju stacionarnost. Radio sam na drugoj temi i ovo sam zapravo koristio u iste svrhe. Isprobajte različita razdoblja, poput 3 do 24 za mjesečne podatke. I testirajte svakog od njih na Ljung-Boxu i pohranite rezultate Chi-Square-a. I odaberite razdoblje s najmanjom vrijednošću hi-kvadrat.

Evo jednostavnog koda za to.

  minval0 <- 5000 # dodijelite veliki broj kako biste bili sigurni da su Chi vrijednosti manje minindex0 <- 0periyot <- 0for (i za 3:24) {# pronađi optimalno razdoblje po Qtestovima u odnosu na izvorne podatke d0D1 <- diff (a, lag = i) #store results Qtest_d0D1 [[i]] <- Box.test (d0D1, lag = 20, type = "Ljung-Box") #store statistika Chi-Square sira0 [i] <- Qtest_d0D1 [[i]] [1]} # skretanje popisa na podatkovni okvir, a zatim matrixdatam0 <- data.frame (matrica (unlist (sira0), nrow = length (Qtest_d0D1) -2, byrow = T)) datamtrx0 <- as.matrix (datam0 []) # get min value's indexminindex0 <- koji (datamtrx0 == min (datamtrx0), arr. ind = F) periyot <- minindex0 + 2  


Ova pitanja su automatski prevedena s engleskog jezika.Izvorni sadržaj dostupan je na stackexchange-u, što zahvaljujemo na cc by-sa 2.0 licenci pod kojom se distribuira.
Loading...