Como estimar o processo de Poisson usando R? (Ou: como usar o pacote NHPoisson?)

15

Eu tenho um banco de dados de eventos (ou seja, uma variável de datas) e covariáveis ​​associadas.

Os eventos são gerados pelo processo Poisson não estacionário, com o parâmetro sendo uma função desconhecida (mas possivelmente linear) de algumas covariáveis.

Eu acho que o pacote NHPoisson existe apenas para esse fim; mas, após 15 horas de pesquisa malsucedida, ainda não estou nem perto de saber como usá-la.

Heck, eu até tentei ler os dois livros referenciados: Coles, S. (2001). Uma introdução à modelagem estatística de valores extremos. Springer. Casella, G. e Berger, RL, (2002). Inferência estatística. Brooks / Cole.

O único exemplo na documentação do fitPP.fun não parece se encaixar na minha configuração; Eu não tenho valores extremos! Eu só tenho eventos simples.

Alguém pode me ajudar com um exemplo simples de ajuste de um processo de Poisson com o parâmetro com covariável único X e suposição de que λ = λ 0 + α X ? Estou interessado na estimativa de λ 0 e α . Eu forneço um conjunto de dados de duas colunas com tempos de eventos (digamos, medidos em segundos após algum tempo arbitrário t 0 ) e outra coluna com valores da covariável X ?λXλ=λ0+αXλ0αt0X

Adam Ryczkowski
fonte
Para os interessados, estou trabalhando na reescrita desta biblioteca para melhorar a usabilidade. github.com/statwonk/NHPoisson
Statwonk

Respostas:

15

Ajustando um processo estacionário de Poisson

Antes de tudo, é importante perceber que tipo de dados de entrada o NHPoisson precisa.

Em primeiro lugar, o NHPoisson precisa de uma lista de índices de momentos do evento. Se registrarmos o intervalo de tempo e o número de eventos no intervalo de tempo, então devemos traduzi-lo em uma única coluna de datas, possivelmente "manchando" as datas no intervalo em que estão gravadas.

Pela simplicidade, assumirei que usamos o tempo medido em segundos e que o "segundo" é a unidade natural de .λ

Vamos simular dados para um processo Poisson simples e estacionário, que possui eventos por minuto:λ=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))

O simNHP.funfaz a simulação. Usamos get aux$posNH, uma variável com índices de momentos de disparo simulado de eventos. Podemos ver que temos aproximadamente 60 * 24 = 1440 eventos, verificando `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

Então, o que fazemos de fato é que aproximamos o processo de Poisson com uma sequência granular de eventos binomiais, cada evento abrange exatamente uma unidade de tempo, em analogia ao mecanismo pelo qual a distribuição de Poisson pode ser vista como um limite da distribuição binomial na lei de eventos raros .

Uma vez que entendemos, o resto é muito mais simples (pelo menos para mim).

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

Ajustando um processo Poisson não estacionário

NHPoisson se encaixa no seguinte modelo:

λ=exp(PTβ)

Pλ

Agora vamos preparar o processo Poisson não estacionário.

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)

Assim como antes, aux$posNHnos daria índices de eventos, mas desta vez perceberemos que a intensidade dos eventos diminui exponencialmente com o tempo. E a taxa dessa diminuição é um 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)

É importante notar que precisamos colocar all.secondscomo covariável, não lambdas. A exponenciação / logaritmização é feita internamente pelo fitPP.fun. Aliás, além dos valores previstos, a função cria dois gráficos por padrão.

A última peça é uma função de faca suíça para validação de modelo globalval.fun,.

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

Entre outras coisas, a função divide o tempo em intervalos, cada uma das lintamostras, para que seja possível criar gráficos brutos que comparem a intensidade prevista com a intensidade observada.

Adam Ryczkowski
fonte
Ótimas explicações Adam, muito obrigado. Tenho a impressão de que não podemos encaixar um modelo com dois grupos de indivíduos e uma intensidade por grupo, estou certo?
Stéphane Laurent
@ StéphaneLaurent Você pode modelar a associação ao grupo como covariável e - sim, você pode adicionar covariáveis. Pode haver intensidade de eventos diferentes para um grupo e diferente para o outro. Eu nunca fiz algo assim, no entanto.
Adam Ryczkowski 24/03
Obrigado Adam. Eu vejo como colocar um "intercepto" por grupo, por exemplo, uma intensidade que tem formaλEu(t)=exp(umaEu+bt), mas não vejo como colocar uma inclinação bEupor grupo.
Stéphane Laurent
Adam, talvez eu estivesse confuso. Agora estou com a impressão de que não há problema. Volto mais tarde em caso de necessidade. Muito obrigado pela atenção.
Stéphane Laurent
resposta muito boa! btw,> out <-fitPP.fun (posE = auxposNH,n=tEume.spuman,betuma=0 0)ErrorEunfEutPP.fvocên(posE=umavocêxposNH, n = time.span, beta = 0): o argumento "start" está ausente, sem padrão
vak 14/02/15