Testowanie znaczenia pików w gęstości widmowej

Czasami używamy wykresu gęstości widmowej do analizy okresowości w szeregach czasowych. Zwykle analizujemy fabułę poprzez kontrolę wzrokową, a następnie próbujemy wyciągnąć wnioski na temat okresowości. Ale czy statystycy opracowali jakiś test, aby sprawdzić, czy jakiekolwiek skoki na wykresie różnią się statystycznie od białego szumu? Czy R-eksperci opracowali jakiś pakiet do analizy gęstości widmowej i do przeprowadzenia tego rodzaju testu? Świetnie, jeśli ktoś może pomóc.

Pozdrawiam
P.

r time-series hypothesis-testing

— Pantera
źródło

Naciskany przez @Wesley, usunąłem moje szybkie przemyślenia na temat funkcji autokorelacji i periodogramu (być może rzeczywiście jest on guru analizy domen częstotliwości, ale osobiście nie sądzę, że Bartlett, pracując z autokorelacją w dziedzinie czasu), ale nadal myślę, że mój druga sugestia bootspecdensmoże być pomocna.

— Dmitrij Celov

Opieram swoje przypuszczenie na reakcji ludzi na „czym jest autokorelacja?” w wyglądzie literatury, w którym prawie wszystkie przypadki, w których stosowana jest autokorelacja, są standardową, obliczoną w dziedzinie czasu autokorelacją Barletta. I niestety jest źle! :) Doceniam sugestię bootspecdensDmitrija; czekam na sprawdzenie.

— Wesley Burr

Należy pamiętać, że szacowanie widm mocy za pomocą periodogramu nie jest zalecane, a faktycznie jest złą praktyką od ~ 1896 r. Jest to niespójny estymator dla mniej niż milionów próbek danych (a nawet wtedy ...) i ogólnie stronniczy. Dokładnie to samo dotyczy stosowania standardowych oszacowań autokorelacji (tj. Bartletta), ponieważ są to pary transformacji Fouriera. Pod warunkiem, że korzystasz ze spójnego narzędzia do szacowania, dostępne są pewne opcje.

Najlepszym z nich jest oszacowanie widma mocy w wielu oknach (lub stożkach). W takim przypadku, korzystając ze współczynników każdego okna z częstotliwością będącą przedmiotem zainteresowania, można obliczyć statystykę harmonicznych F na podstawie hipotezy zerowej białego szumu. Jest to doskonałe narzędzie do wykrywania elementów linii w hałasie i jest wysoce zalecane. Jest to domyślny wybór w środowisku przetwarzającym sygnały do wykrywania okresowości hałasu przy założeniu stacjonarności.

Możesz uzyskać dostęp zarówno do metody multitaperwielopunktowej oceny widma, jak i powiązanego testu F za pośrednictwem pakietu w R (dostępnego przez CRAN). Dokumentacja dostarczona z pakietem powinna wystarczyć, aby zacząć; test F jest prostą opcją w wywołaniu funkcji spec.mtm.

Pierwotne odniesienie, które definiuje obie te techniki i podaje dla nich algorytmy, to: Spectrum Estimation and Harmonic Analysis , DJ Thomson, Proceedings of the IEEE, vol. 70, str. 1055–1096, 1982.

Oto przykład użycia dołączonego zestawu danych z multitaperpakietem.

require(multitaper);
data(willamette);
resSpec <- spec.mtm(willamette, k=10, nw=5.0, nFFT = "default",
                    centreWithSlepians = TRUE, Ftest = TRUE,
                    jackknife = FALSE, maxAdaptiveIterations = 100,
                    plot = TRUE, na.action = na.fail)

Parametry, o których powinieneś wiedzieć, to k i nw : jest to liczba okien (ustawiona na 10 powyżej) i iloczyn przepustowości czasowej (5.0 powyżej). Możesz łatwo pozostawić te wartości quasi-domyślne dla większości aplikacji. Komenda centreWithSlepians usuwa rzetelne oszacowanie średniej szeregów czasowych za pomocą rzutowania na okna Slepian - jest to również zalecane, ponieważ pozostawienie wartości średniej powoduje wytworzenie dużej mocy przy niskich częstotliwościach.

Poleciłbym również wykreślić wyjście widma z „spec.mtm” w skali logarytmicznej, ponieważ znacznie to oczyszcza. Jeśli potrzebujesz więcej informacji, po prostu opublikuj i chętnie je przekażę.

— Wesley Burr
źródło

Burrowi, Silvie i Celovowi - wielkie dzięki za interesujące odpowiedzi i sugestie. Nie mogę się doczekać, aby przetestować te estymatory. Z pozdrowieniami

— Pantera

(+1) tej nocy dokładnie przemyślałem twoje sugestie i zdecydowałem, że domena czasu rzeczywiście jest ostatnią rzeczą (ze względu na skrócenie opóźnienia i słabe właściwości w małych próbkach), aby spróbować wyszukać zachowania rowerowe. Osobiście niepokoją mnie założenia do statystyki F i właściwości małego rozmiaru próby sugerowanego schematu. No i prawdopodobnie dobrze jest zacząć osobne pytanie dotyczące optymalnego wyboru okna, ponieważ rzeczywiście jest ich wiele.

— Dmitrij Celov,

Rzeczywiście istnieje wiele możliwości wyboru okien, chociaż dwie najczęstsze to Sekwencje Kuliste Dyskretne Prolate (lub Slepian ) i sinusoidy zwężają się. Jeśli szukasz maksymalnego stężenia energii w lokalnej szerokości pasma, okazało się, że Slepianie są optymalni i faktycznie są wynikiem integralnej postaci równania gęstości widmowej (szczegółowe informacje znajdują się w artykule, o którym wspominałem). Jeśli chodzi o statystyki F, z pewnością istnieją pewne problemy ze stopniami swobody, ale ogólnie działają one całkiem dobrze, z dostępnym ~ 2k-2 dof.

— Wesley Burr

Wygładzony periodogram wykorzystuje również stożek, pozwala na FFT, książka Davida Stoffera uczy również, jak obliczać przedziały ufności. multitaperWydaje się, że w tym pakiecie zastosowano bardziej zaawansowane techniki zmniejszania i obliczania przedziału ufności. Ale myślę, że pomysł był taki sam, zdaniem Davida Stoffera. To jedyna rzecz, jaką mogłem pomyśleć o tym nauczaniu waniliowego perydogoramu, nadal ma sens.

— stucash

ok, więc jesteś jednym z autorów tego pakietu i użyłeś bardzo silnych słów przeciwko periodogramowi. Mam nadzieję, że pewnego dnia możesz wrócić z większą ilością dowodów. Powszechne zalety i wady Periodogramu są dobrze znane, podobnie jak jego wybuchowa wariancja, dlatego nie jest to dobry spójny estymator widma, ale wygładzony periodogram nie jest wcale taki zły, nie tak zły, jak tutaj powiedziałeś.

— stucash

Próbowaliśmy próbę rozwiązania tego problemu przez transformaty falkowej testu widmowej oparte niedawno w tym artykule . Zasadniczo należy wziąć pod uwagę rozkład rzędnych periodogramu, podobnie jak w artykule Fishera, wspomnianym we wcześniejszych odpowiedziach. Kolejny artykuł Koena jest taki . Niedawno opublikowaliśmy pakiet R. hwwntest .

— Delyan Savchev
źródło

Savchev, dziękuję bardzo za komentarz i referencje. Z niecierpliwością czekam na testowanie twojego pakietu R.

— Pantera,

$f(\omega_k)$

Więcej informacji na temat testu można znaleźć w MB Priestley, Spectral Analysis and Time Series , Academic Press, Londyn, 1981, strona 406.

W R pakiet GeneCycle zawiera funkcję fisher.g.test():

library(GeneCycle)
?fisher.g.test

Mam nadzieję że to pomoże.

— Washington S. Silva
źródło

to świetnie, ale test g pakietu opiera się na własnej funkcji periodogramu, która ma bardzo ograniczone opcje obliczania widm mocy ...

— stucash