Jak oszacować proces Poissona za pomocą R? (Lub: jak korzystać z pakietu NHPoisson?)

15

Mam bazę danych zdarzeń (tj. Zmienną dat) i powiązanych zmiennych towarzyszących.

Zdarzenia są generowane przez niestacjonarny proces Poissona z parametrem będącym nieznaną (ale być może liniową) funkcją niektórych zmiennych towarzyszących.

Myślę, że pakiet NHPoisson istnieje właśnie w tym celu; ale po 15 godzinach nieudanych badań wciąż nie jestem w pobliżu, aby dowiedzieć się, jak z niego korzystać.

Heck, nawet próbowałem czytać obie przywoływane książki: Coles, S. (2001). Wprowadzenie do modelowania statystycznego wartości ekstremalnych. Skoczek. Casella, G. i Berger, RL, (2002). Wnioskowanie statystyczne. Brooks / Cole.

Jeden przykład w dokumentacji fitPP.fun wydaje się nie pasować do mojej konfiguracji; Nie mam ekstremalnych wartości! Mam tylko nagie wydarzenia.

Czy ktoś może mi pomóc z prostym przykładem dopasowania procesu Poissona za pomocą parametru do pojedynczej zmiennej objaśniającej X i założenia, że λ = λ 0 + α X ? Interesuje mnie oszacowanie λ 0 i . Podaję dwukolumnowy zestaw danych z czasami zdarzeń (powiedzmy, mierzonymi w sekundach po pewnym dowolnym czasie ) i kolejną kolumnę z wartościami współzmiennej ?λXλ=λ0+αXλ0t 0 Xαt0X

Adam Ryczkowski
źródło
Dla zainteresowanych pracuję nad przepisem tej biblioteki, aby poprawić użyteczność. github.com/statwonk/NHPoisson
Statwonk 13.04.16

Odpowiedzi:

15

Dopasowanie stacjonarnego procesu Poissona

Przede wszystkim należy zdać sobie sprawę, jakiego rodzaju danych wejściowych potrzebuje NHPoisson.

Przede wszystkim NHPoisson potrzebuje listy wskaźników momentów zdarzeń. Jeśli rejestrujemy przedział czasowy i liczbę zdarzeń w danym przedziale czasowym, wówczas musimy go przetłumaczyć na pojedynczą kolumnę dat, prawdopodobnie „rozmazując” daty w przedziale, w którym są rejestrowane.

λ

λ=1

lambda=1/60 #1 event per minute
time.span=60*60*24 #24 hours, with time granularity one second

aux<-simNHP.fun(rep(lambda,time.span))

To simNHP.funsprawia, że ​​symulacja. Używamy get aux$posNH, zmiennej ze wskaźnikami momentów symulowanego odpalenia zdarzenia. Widzimy, że z grubsza mamy 60 * 24 = 1440 zdarzeń, sprawdzając `length (aux $ posNH).

λfitPP.fun

out<-fitPP.fun(posE=aux$posNH,n=time.span,start=list(b0=0)) # b0=0 is our guess at initial value for optimization, which is internally made with `nlminb` function

λ>0fitPP

Więc co robimy w rzeczywistości jest to, że możemy w przybliżeniu proces Poissona z granulowanym sekwencji dwumiennych zdarzeń, każdy przęsła zdarzeń dokładnie jedna jednostka czasu, w sposób analogiczny do mechanizmu, w którym rozkład Poissona może być postrzegane jako ograniczenie rozkładu dwumianowego w prawo rzadkich wydarzeń .

Kiedy to zrozumiemy, reszta jest znacznie prostsza (przynajmniej dla mnie).

λbetaexp(coef(out)[1])NHPoissonλλ

Dopasowanie niestacjonarnego procesu Poissona

NHPoisson pasuje do następującego modelu:

λ=exp(P.T.β)

P.λ

Teraz przygotujmy niestacjonarny proces Poissona.

time.span=60*60*24 #24 hours, with time granularity one second
all.seconds<-seq(1,time.span,length.out=time.span)
lambdas=0.05*exp(-0.0001*all.seconds) #we can't model a linear regression with NHPoisson. It must have the form with exp.
aux<-simNHP.fun(lambdas)

Tak jak poprzednio, aux$posNHdałby nam wskaźniki zdarzeń, ale tym razem zauważymy, że intensywność wydarzeń maleje wykładniczo z czasem. Szybkość tego zmniejszania się jest parametrem, który chcemy oszacować.

out<-fitPP.fun(tind=TRUE,covariates=cbind(all.seconds),
        posE=aux$posNH,
        start=list(b0=0,b1=0),modSim=TRUE)

Ważne jest, aby pamiętać, że musimy all.secondspodawać jako zmienną towarzyszącą, a nie lambdas. Potęgowanie / logarytmizacja odbywa się wewnętrznie przez fitPP.fun. BTW, oprócz przewidywanych wartości, funkcja domyślnie tworzy dwa wykresy.

Ostatni kawałek jest funkcją szwajcarski nóż do walidacji modelu globalval.fun.

aux<-globalval.fun(obFPP=out,lint=2000,
        covariates=cbind(all.seconds),typeI='Disjoint',
        typeRes='Raw',typeResLV='Raw',resqqplot=FALSE)

Między innymi funkcja dzieli czas na przedziały, każda lintpróbka jest długa, więc możliwe jest tworzenie surowych wykresów, które porównują przewidywane natężenie z obserwowanym natężeniem.

Adam Ryczkowski
źródło
Świetne wyjaśnienia Adam, dziękuję bardzo. Mam wrażenie, że nie możemy dopasować modelu z dwiema grupami osób i jedną intensywnością na grupę, prawda?
Stéphane Laurent
@ StéphaneLaurent Możesz modelować członkostwo w grupie jako zmienną towarzyszącą i - tak, możesz dodać zmienne towarzyszące. Intensywność wydarzeń może być różna dla jednej grupy i inna dla drugiej. Jednak nigdy czegoś takiego nie zrobiłem.
Adam Ryczkowski
λja(t)=exp(zaja+bt)bja
Adam, może byłem zdezorientowany. Teraz mam wrażenie, że nie ma problemu. Wrócę później w razie potrzeby. Dziękuję bardzo za uwagę.
Stéphane Laurent
posN.H.,n=tjammi.spzan,bmitza=0)mirrorjanfajatP.P..faun(posmi=zauxposNH, n = time.span, beta = 0): brakuje argumentu „start”, bez domyślnego ustawienia
vak