¿Cómo estimar el proceso de Poisson utilizando R? (O: ¿cómo utilizar el paquete NHPoisson?)

Question

Tengo una base de datos de eventos (es decir, una variable de fechas) y covariables asociadas.

Los sucesos son generados por el proceso de Poisson no estacionario con un parámetro que es una función desconocida (pero posiblemente lineal) de algunas covariables.

Creo que el paquete NHPoisson existe justo para este propósito; pero después de 15 horas de investigación infructuosa todavía estoy lejos de saber cómo usarlo.

Diablos, incluso intenté leer los dos libros referenciados: Coles, S. (2001). An introduction to statistical modelling of extreme values. Springer. Casella, G. y Berger, R.L., (2002). Statistical inference. Brooks/Cole.

El único ejemplo que aparece en la documentación de fitPP.fun no parece ajustarse a mi configuración; ¡no tengo valores extremos! Sólo tengo eventos desnudos.

¿Puede alguien ayudarme con un ejemplo sencillo de ajuste de un proceso de Poisson con parámetro $\lambda$ con una única covariable $X$ y suposición, que $\lambda = \lambda_0 + \alpha \cdot X$ ? Estoy interesado en la estimación de $\lambda_0$ y $\alpha$ . Proporciono un conjunto de datos de dos columnas con tiempos de eventos (digamos, medidos en segundos después de algún tiempo arbitrario $t_0$ ) y otra columna con los valores de la covariable $X$ ?

Answer 1

Ajuste de un proceso de Poisson estacionario

En primer lugar, es importante saber qué tipo de datos de entrada necesita NHPoisson.

Ante todo, NHPoisson necesita una lista de índices de los momentos del acontecimiento. Si registramos el intervalo de tiempo y el número de eventos en el intervalo de tiempo, entonces debemos traducirlo en una sola columna de fechas, posiblemente "emborronando" las fechas sobre el intervalo en el que se registran.

Para simplificar, supondré que utilizamos el tiempo medido en segundos, y que el "segundo" es la unidad natural de $\lambda$ .

Vamos a simular los datos de un proceso de Poisson simple y estacionario, que tiene $\lambda=1$ eventos por minuto:

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

En simNHP.fun realiza la simulación. Utilizamos para obtener aux$posNH una variable con índices de momentos de disparo de eventos simulados. Podemos ver que tenemos aproximadamente 60*24 = 1440 eventos, comprobando `length(aux$posNH).

Ahora hagamos ingeniería inversa del $\lambda$ con 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

Como la función sólo toma índices de eventos, necesita también una medida de cuántos índices posibles hay. Y esta es una parte muy confusa porque en el verdadero proceso de Poisson es posible tener número infinito de acontecimientos posibles (si sólo $\lambda>0$ ). Pero desde la perspectiva de la fitPP tenemos que elegir alguna unidad de tiempo lo suficientemente pequeña. La elegimos tan pequeña que podamos suponer como máximo un suceso por unidad de tiempo.

Así que lo que hacemos de hecho es que aproximado el proceso de Poisson con secuencia granular de eventos binomiales, cada evento abarca exactamente una unidad de tiempo, en analogía con el mecanismo en el que la distribución de Poisson puede verse como un límite de la distribución binomial en la ley de las rarezas .

Una vez que lo entendemos, el resto es mucho más sencillo (al menos para mí).

Para obtener la aproximación de nuestro $\lambda$ necesito tomar el exponente del parámetro ajustado beta , exp(coef(out)[1]) . Y aquí viene otro dato importante, que falta en el manual : con NHPoisson nos ajustamos a la logaritmo de $\lambda$ no a $\lambda$ sí mismo.

Ajuste de un proceso de Poisson no estacionario

NHPoisson se ajusta al siguiente modelo:

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

es decir, se ajusta a la combinación lineal de los parámetros $\vec{P}$ (denominadas covariables) a la logaritmo de la $\lambda$ .

Preparemos ahora un proceso de Poisson no estacionario.

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)

Igual que antes, aux$posNH nos daría índices de sucesos, pero esta vez notaremos, que la intensidad de los sucesos disminuye exponencialmente con el tiempo. Y la tasa de esta disminución es un parámetro que queremos estimar.

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

Es importante señalar que tenemos que poner all.seconds como covariable, no lambdas . La exponenciación/logaritmización se realiza internamente por el programa fitPP.fun . BTW, aparte de los valores predichos, la función hace dos gráficos por defecto.

La última pieza es una función swiss-knife para la validación del modelo, globalval.fun .

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

Entre otras cosas, la función divide el tiempo en intervalos, cada uno de los cuales lint de muestras, por lo que es posible crear gráficos rudimentarios que comparen la intensidad prevista con la observada.