Wie kann man den Poisson-Prozess mit R abschätzen? (Oder: Wie verwende ich das NHPoisson-Paket?)

15

Ich habe eine Datenbank von Ereignissen (dh eine Variable von Daten) und zugehörigen Kovariaten.

Die Ereignisse werden durch den nicht stationären Poisson-Prozess erzeugt, wobei Parameter eine unbekannte (aber möglicherweise lineare) Funktion einiger Kovariaten sind.

Ich denke, das NHPoisson-Paket gibt es nur für diesen Zweck; aber nach 15 stunden erfolgloser recherche bin ich noch lange nicht in der lage zu wissen, wie man es benutzt.

Ich habe sogar versucht, beide Bücher zu lesen, auf die verwiesen wird: Coles, S. (2001). Eine Einführung in die statistische Modellierung von Extremwerten. Springer. Casella, G. und Berger, RL, (2002). Statistische Inferenz. Brooks / Cole.

Das einzige Beispiel in der Dokumentation von fitPP.fun scheint nicht zu meinem Setup zu passen. Ich habe keine extremen Werte! Ich habe nur nackte Ereignisse.

Kann mir bitte jemand mit einem einfachen Beispiel für die Anpassung eines Poisson-Prozesses mit dem Parameter mit einer einzelnen Kovariate und der Annahme helfen , dass ? Ich interessiere mich für die Schätzung von und . Ich biete einen zweispaltigen Datensatz mit Ereigniszeiten (zum Beispiel gemessen in Sekunden nach einer beliebigen Zeit ) und eine andere Spalte mit Werten der Kovariate ?λXλ=λ0+αXλ0αt0X

Adam Ryczkowski
quelle
Für Interessierte arbeite ich an einer Überarbeitung dieser Bibliothek, um die Benutzerfreundlichkeit zu verbessern. github.com/statwonk/NHPoisson
Statwonk

Antworten:

15

Anpassen eines stationären Poisson-Prozesses

Zuallererst ist es wichtig zu erkennen, welche Art von Eingabedaten NHPoisson benötigt.

Zuallererst benötigt NHPoisson eine Liste von Indizes von Ereignismomenten. Wenn wir das Zeitintervall und die Anzahl der Ereignisse in dem Zeitintervall aufzeichnen, müssen wir es in eine einzelne Datenspalte übersetzen und möglicherweise die Daten über das Intervall "verwischen", in dem sie aufgezeichnet wurden.

Der Einfachheit halber gehe ich davon aus, dass wir die in Sekunden gemessene Zeit verwenden und dass die "Sekunde" die natürliche Einheit von .λ

Simulieren wir Daten für einen einfachen, stationären Poisson-Prozess mit Ereignissen pro Minute:λ=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))

Das simNHP.funmacht die Simulation. Wir verwenden, um aux$posNHeine Variable mit Indizes von Momenten des simulierten Ereignisauslösens zu erhalten. Wir können sehen, dass wir ungefähr 60 * 24 = 1440 Ereignisse haben, indem wir die Länge (aux $ posNH) überprüfen.

Lassen Sie uns nun das mit :λ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

Was wir also in der Tat tun, ist, dass wir den Poisson-Prozess mit einer granularen Folge von Binomialereignissen approximieren , wobei jedes Ereignis genau eine Zeiteinheit umfasst, in Analogie zu dem Mechanismus, bei dem die Poisson-Verteilung als Grenze der Binomialverteilung im Gesetz angesehen werden kann seltener Ereignisse .

Sobald wir es verstanden haben, ist der Rest viel einfacher (zumindest für mich).

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

Anpassen eines instationären Poisson-Prozesses

NHPoisson passt für folgendes Modell:

λ=exp(PTβ)

Pλ

Nun bereiten wir einen instationären Poisson-Prozess vor.

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)

Wie zuvor aux$posNHwürden wir Indizes von Ereignissen erhalten, aber dieses Mal werden wir feststellen, dass die Intensität der Ereignisse mit der Zeit exponentiell abnimmt. Und die Rate dieser Abnahme ist ein Parameter, den wir schätzen möchten.

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

Es ist wichtig anzumerken, dass wir all.secondsals Kovariate setzen müssen, nicht lambdas. Die Exponentiation / Logaritmisierung erfolgt intern durch die fitPP.fun. Übrigens, abgesehen von den vorhergesagten Werten, erstellt die Funktion standardmäßig zwei Diagramme.

Das letzte Stück ist eine Schweizer Messerfunktion zur Modellvalidierung globalval.fun.

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

Unter anderem unterteilt die Funktion die Zeit in Intervalle, die jeweils lintlang sind, sodass grobe Diagramme erstellt werden können, die die vorhergesagte Intensität mit der beobachteten Intensität vergleichen.

Adam Ryczkowski
quelle
Tolle Erklärungen Adam, vielen Dank. Ich habe den Eindruck, wir können ein Modell nicht mit zwei Gruppen von Individuen und einer Intensität pro Gruppe kombinieren, oder?
Stéphane Laurent
@ StéphaneLaurent Sie können die Gruppenmitgliedschaft als Kovariate modellieren und - ja, Sie können Kovariaten hinzufügen. Die Intensität der Ereignisse kann für eine Gruppe unterschiedlich und für die andere unterschiedlich sein. Sowas habe ich aber noch nie gemacht.
Adam Ryczkowski
λich(t)=exp(einich+bt)bich
Adam, vielleicht war ich verwirrt. Jetzt habe ich den Eindruck, dass es kein Problem gibt. Ich komme später zurück, falls es nötig ist. Vielen Dank für Ihre Aufmerksamkeit.
Stéphane Laurent
pÖsNH,n=tichme.speinn,betein=0)ErrÖrichnfichtPP.fun(pÖsE=einuxposNH, n = time.span, beta = 0): Argument "start" fehlt, ohne Voreinstellung
vak