Comment estimer le processus de Poisson en utilisant R? (Ou: comment utiliser le package NHPoisson?)

15

J'ai une base de données d'événements (c'est-à-dire une variable de dates) et des covariables associées.

Les événements sont générés par le processus de Poisson non stationnaire, le paramètre étant une fonction inconnue (mais peut-être linéaire) de certaines covariables.

Je pense que le package NHPoisson existe uniquement à cette fin; mais après 15 heures de recherches infructueuses, je suis encore loin de savoir comment l'utiliser.

Heck, j'ai même essayé de lire les deux ouvrages référencés: Coles, S. (2001). Une introduction à la modélisation statistique des valeurs extrêmes. Springer. Casella, G. et Berger, RL, (2002). Inférence statistique. Brooks / Cole.

Le seul exemple dans la documentation de fitPP.fun ne semble pas correspondre à ma configuration; Je n'ai pas de valeurs extrêmes! J'ai juste des événements nus.

Quelqu'un peut-il m'aider avec un exemple simple d'ajustement d'un processus de Poisson avec le paramètre avec une seule covariable , et en supposant que ? Je suis intéressé par l'estimation de et . Je fournis un ensemble de données à deux colonnes avec des temps d'événements (disons, mesurés en secondes après un certain temps arbitraire ) et une autre colonne avec des valeurs de la covariable ?λXλ=λ0+αXλ0αt0X

Adam Ryczkowski
la source
Pour ceux qui sont intéressés, je travaille sur une réécriture de cette bibliothèque pour améliorer la convivialité. github.com/statwonk/NHPoisson
Statwonk

Réponses:

15

Montage d'un processus de Poisson stationnaire

Tout d'abord, il est important de savoir quel type de données d'entrée NHPoisson a besoin.

Avant tout, NHPoisson a besoin d'une liste d' indices de moments d'événements. Si nous enregistrons l'intervalle de temps et le nombre d'événements dans l'intervalle de temps, alors je dois le traduire en une seule colonne de dates, éventuellement "étalant" les dates sur l'intervalle sur lequel elles sont enregistrées.

Pour la simplicité, je suppose que nous utilisons du temps mesuré en secondes et que la "seconde" est l'unité naturelle de .λ

Simulons les données d'un processus de Poisson simple et stationnaire, qui a événements par 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))

Le simNHP.funfait la simulation. Nous utilisons pour obtenir aux$posNH, une variable avec des indices de moments de déclenchement d'événements simulés. Nous pouvons voir que nous avons environ 60 * 24 = 1440 événements, en vérifiant `length (aux $ posNH).

Maintenant, rétroconcevons le avec :λ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

Étant donné que la fonction ne prend que les indices d'événements, elle a également besoin d'une mesure du nombre d'indices possibles. Et c'est une partie très déroutante , car dans le vrai processus de Poisson, il est possible d'avoir un nombre infini d'événements possibles (si seulement ). Mais du point de vue de la, nous devons choisir une unité de temps suffisamment petite. Nous le choisissons si petit que nous pouvons supposer au maximum un événement par unité de temps.λ>0fitPP

Donc, ce que nous faisons en fait, c'est que nous approchons le processus de Poisson avec une séquence granulaire d'événements binomiaux, chaque événement s'étend exactement sur une unité de temps, par analogie avec le mécanisme dans lequel la distribution de Poisson peut être considérée comme une limite de la distribution binomiale dans la loi d'événements rares .

Une fois que nous l'avons compris, le reste est beaucoup plus simple (du moins pour moi).

Pour obtenir l'approximation de notre me dois prendre exposant du paramètre ajusté , . Et voici une autre information importante, qui manque dans le manuel : avec nous nous adaptons au logarithme de , pas à lui-même.λbetaexp(coef(out)[1])NHPoissonλλ

Montage d'un processus de Poisson non stationnaire

NHPoisson s'adapte au modèle suivant:

λ=exp(PTβ)

c'est-à-dire qu'il adapte la combinaison linéaire des paramètres (appelés covariables) au logarithme du .Pλ

Préparons maintenant le processus de Poisson non stationnaire.

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)

Tout comme auparavant, aux$posNHnous donnerait des indices d'événements, mais cette fois nous remarquerons que l'intensité des événements diminue exponentiellement avec le temps. Et le taux de cette diminution est un paramètre que nous voulons estimer.

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

Il est important de noter que nous devons mettre en all.secondstant que covariable, non lambdas. L'exponentiation / logaritmisation se fait en interne par le fitPP.fun. BTW, en dehors des valeurs prédites, la fonction crée deux graphiques par défaut.

La dernière pièce est une fonction couteau suisse pour la validation du modèle, globalval.fun.

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

Entre autres choses, la fonction divise le temps en intervalles, chacun étant lintlong, il est donc possible de créer des graphiques bruts qui comparent l'intensité prédite avec l'intensité observée.

Adam Ryczkowski
la source
Excellentes explications Adam, merci beaucoup. J'ai l'impression que nous ne pouvons pas adapter un modèle avec deux groupes d'individus et une intensité par groupe, ai-je raison?
Stéphane Laurent
@ StéphaneLaurent Vous pouvez modéliser l'appartenance à un groupe comme covariable, et - oui, vous pouvez ajouter des covariables. L'intensité des événements peut être différente pour un groupe et différente pour l'autre. Mais je n'ai jamais rien fait de tel.
Adam Ryczkowski
Merci Adam. Je vois comment mettre une "interception" par groupe, par exemple une intensité ayant la forme , mais je ne vois pas comment mettre une pente par groupe. λje(t)=exp(uneje+bt)bje
Stéphane Laurent
Adam, j'étais peut-être confus. Maintenant, j'ai l'impression qu'il n'y a pas de problème. Je reviendrai plus tard en cas de besoin. Merci beaucoup pour votre attention.
Stéphane Laurent
très bonne réponse! btw,> out <-fitPP.fun (posE = aux posNH, n = time.span, beta = 0): argument "start "est manquant, sans défautposNH,n=tjeme.spunen,betune=0)ErrorjenFjetPP.Fun(posE=uneuX
vak