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

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 ?$\lambda$$X$$\lambda = \lambda_0 + \alpha \cdot X$$\lambda_0$t 0$\alpha$$t_0$$X$

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.

$\lambda$

$\lambda=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).

$\lambda$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

$\lambda>0$fitPP

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).

$\lambda$betaexp(coef(out)[1])NHPoisson$\lambda$$\lambda$

Dopasowanie niestacjonarnego procesu Poissona

NHPoisson pasuje do następującego modelu:

$\lambda = \exp(\vec{P}^T \cdot \vec{\beta} )$

$\vec{P}$$\lambda$

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.