Como saber se um dado segue uma distribuição de Poisson em R?

Sou estudante de graduação e tenho um projeto para a minha aula de probabilidade. Basicamente, tenho um conjunto de dados sobre os furacões que afetaram meu país por uma série de anos.

No meu Livro de Probabilidades, (Probabilidade e Estatística com R), há um exemplo (não completo) de como verificar se os dados seguem uma distribuição de Poisson, eles começam a tentar provar que esses 3 critérios são seguidos: (Do meu livro, página 120 (critério) página 122-123 exemplo)

1- O número de resultados em intervalos sem sobreposição é independente. Em outras palavras, o número de resultados no intervalo de tempo (0, t] é independente do número de resultados no intervalo de tempo (t, t + h], h> 0

2- A probabilidade de dois ou mais resultados em um intervalo suficientemente curto é praticamente zero. Em outras palavras, desde que h seja suficientemente pequeno, a probabilidade de obter dois ou mais resultados no intervalo (t, t + h] é insignificante em comparação com a probabilidade de obter um ou zero resultados no mesmo intervalo de tempo.

3- A probabilidade de exatamente um resultado em um intervalo suficientemente curto ou em uma região pequena é proporcional à duração do intervalo ou região. Em outras palavras, a probabilidade de um resultado em um intervalo de comprimento h é lambda * h.

Mas o critério 3 é deixado "como um exercício".

A- Alguém pode me dizer se existe uma maneira mais "fácil" de ver se meu conjunto de dados segue uma distribuição de Poisson?

B- Alguém pode me explicar os critérios 1 e 3 com algum tipo de exemplo (se for com R, fantástico)?

Obrigado!

Nota: desculpe pela postagem longa. Além disso, eu tenho que converter os dados para que eu tenha uma tabela como:

  number of hurricanes       | 0 | 1 | 2  etc.
  -----------------------------------------
  total years that have      |   |   |
  that number of hurricanes  |   |   |

Há um número infinito de maneiras para uma distribuição ser ligeiramente diferente de uma distribuição de Poisson; você não pode identificar que um conjunto de dados é extraído de uma distribuição Poisson. O que você pode fazer é procurar inconsistência com o que você deve ver com um Poisson, mas uma falta de inconsistência óbvia não o torna Poisson.

No entanto, o que você está falando lá, verificando esses três critérios, não é verificar se os dados provêm de uma distribuição de Poisson por meios estatísticos (por exemplo, olhando dados), mas avaliando se o processo em que os dados são gerados satisfaz a condições de um processo de Poisson; se todas as condições mantiveram ou quase mantiveram (e isso é uma consideração do processo de geração de dados), você pode ter algo de ou muito perto de um processo de Poisson, que, por sua vez, seria uma maneira de obter dados extraídos de algo próximo a um Distribuição de veneno.

Mas as condições não se sustentam de várias maneiras ... e o mais longe de ser verdade é o número 3. Não há nenhuma razão específica nessa base para afirmar um processo de Poisson, embora as violações não sejam tão ruins que os dados resultantes estejam longe. de Poisson.

Então, voltamos aos argumentos estatísticos resultantes do exame dos próprios dados. Como os dados mostram que a distribuição era Poisson, e não algo parecido?

Como mencionado no início, o que você pode fazer é verificar se os dados não são obviamente inconsistentes com a distribuição subjacente sendo Poisson, mas isso não indica que eles foram retirados de um Poisson (você já pode ter certeza de que eles são não).

Você pode fazer essa verificação através de testes de qualidade de ajuste.

O qui-quadrado mencionado é um deles, mas eu não recomendaria o teste do qui-quadrado para essa situação **; tem baixo poder contra desvios interessantes. Se seu objetivo é ter um bom poder, você não conseguirá assim (se não se importa com o poder, por que você testaria?). Seu principal valor é a simplicidade e possui valor pedagógico; Fora isso, não é competitivo como um teste de qualidade.

** Adicionado em edição posterior: agora que está claro que isso é lição de casa, as chances de você fazer um teste qui-quadrado para verificar os dados não são inconsistentes com um Poisson aumentam bastante. Veja meu exemplo de teste de ajuste do qui-quadrado, feito abaixo do primeiro gráfico de Poissonness

As pessoas costumam fazer esses testes pelo motivo errado (por exemplo, porque querem dizer 'portanto, não há problema em fazer outra coisa estatística com os dados que pressupõem que os dados são Poisson'). A verdadeira questão é "quão errado isso poderia acontecer?" ... e a qualidade dos testes de ajuste não ajuda muito nessa questão. Muitas vezes, a resposta para essa pergunta é, na melhor das hipóteses, independente do tamanho da amostra (/ quase independente) - e, em alguns casos, uma com consequências que tendem a desaparecer com o tamanho da amostra ... enquanto um teste de qualidade do ajuste é inútil com o tamanho da amostra. amostras pequenas (onde o risco de violações de suposições costuma ser maior).

Se você deve testar uma distribuição Poisson, existem algumas alternativas razoáveis. Um seria fazer algo semelhante a um teste de Anderson-Darling, com base na estatística AD, mas usando uma distribuição simulada sob o nulo (para explicar os problemas gêmeos de uma distribuição discreta e que você deve estimar parâmetros).

Uma alternativa mais simples pode ser um Teste Suave para adequação - estes são uma coleção de testes projetados para distribuições individuais, modelando os dados usando uma família de polinômios ortogonais em relação à função de probabilidade no nulo. Alternativas de ordem baixa (ou seja, interessantes) são testadas testando se os coeficientes dos polinômios acima da base são diferentes de zero e, geralmente, eles podem lidar com a estimativa de parâmetros, omitindo os termos de ordem mais baixa do teste. Existe um teste para o Poisson. Eu posso desenterrar uma referência, se você precisar.

$n(1-r^2)$$\log(x_k)+\log(k!)$$k$

Aqui está um exemplo desse cálculo (e gráfico), feito em R:

y=rpois(100,5)
n=length(y)
(x=table(y))
y
 0  1  2  3  4  5  6  7  8  9 10 
 1  2  7 15 19 25 14  7  5  1  4 

k=as.numeric(names(x))
plot(k,log(x)+lfactorial(k))

Aqui está a estatística que sugeri que poderia ser usada para um teste de ajuste de um Poisson:

n*(1-cor(k,log(x)+lfactorial(k))^2)
[1] 1.0599

Obviamente, para calcular o valor-p, você também precisará simular a distribuição da estatística de teste sob o valor nulo (e eu não discuti como se pode lidar com contagens zero dentro do intervalo de valores). Isso deve produzir um teste razoavelmente poderoso. Existem vários outros testes alternativos.

Aqui está um exemplo de como fazer um gráfico de Poissonness em uma amostra de tamanho 50 de uma distribuição geométrica (p = 0,3):

Como você vê, ele exibe uma clara 'torção', indicando não linearidade

As referências para o gráfico de Poissonness seriam:

David C. Hoaglin (1980),
"A Poissonness Plot",
The American Statistician
Vol. 34, No. 3 (ago.), Pp. 146-149

e

Hoaglin, D. J. e Tukey (1985),
"9. Verificar a Forma de Discrete Distributions",
Explorando tabelas de dados, Trends e formas ,
(Hoaglin, Mosteller & Tukey Eds)
John Wiley & Sons

A segunda referência contém um ajuste no gráfico para pequenas contagens; você provavelmente gostaria de incorporá-lo (mas não tenho a referência em mãos).

Exemplo de realização de um teste de ajuste de qualidade qui-quadrado:

Além de executar a qualidade do ajuste do qui-quadrado, da maneira que normalmente se espera que seja feito em muitas classes (embora não da maneira que eu faria):

1: começando com seus dados (que considerarei os dados gerados aleatoriamente em 'y' acima, gere a tabela de contagens:

(x=table(y))
y
 0  1  2  3  4  5  6  7  8  9 10 
 1  2  7 15 19 25 14  7  5  1  4 

2: calcule o valor esperado em cada célula, assumindo um Poisson ajustado por ML:

 (expec=dpois(0:10,lambda=mean(y))*length(y))
 [1]  0.7907054  3.8270142  9.2613743 14.9416838 18.0794374 17.5008954 14.1173890  9.7611661
 [9]  5.9055055  3.1758496  1.5371112

3: observe que as categorias finais são pequenas; isso torna a distribuição qui-quadrado menos boa como uma aproximação à distribuição da estatística de teste (uma regra comum é que você queira valores esperados de pelo menos 5, embora vários trabalhos tenham mostrado que essa regra é desnecessariamente restritiva; fechar, mas a abordagem geral pode ser adaptada a uma regra mais rígida). Reduza as categorias adjacentes, para que os valores mínimos esperados fiquem pelo menos não muito abaixo de 5 (uma categoria com uma contagem esperada próxima a 1 em mais de 10 categorias não é muito ruim, duas é bastante limítrofe). Observe também que ainda não contabilizamos a probabilidade além de "10", portanto, também precisamos incorporar isso:

expec[1]=sum(expec[1:2])
expec[2:8]=expec[3:9]
expec[9]=length(y)-sum(expec[1:8])
expec=expec[1:9]
expec
sum(expec) # now adds to n

4: da mesma forma, colapsar categorias no observado:

(obs=table(y))
obs[1]=sum(obs[1:2])
obs[2:8]=obs[3:9]
obs[9]=sum(obs[10:11])
obs=obs[1:9]

$(O_i-E_i)^2/E_i$

print(cbind(obs,expec,PearsonRes=(obs-expec)/sqrt(expec),ContribToChisq=(obs-expec)^2/expec),d=4)
  obs  expec PearsonRes ContribToChisq
0   3  4.618   -0.75282      0.5667335
1   7  9.261   -0.74308      0.5521657
2  15 14.942    0.01509      0.0002276
3  19 18.079    0.21650      0.0468729
4  25 17.501    1.79258      3.2133538
5  14 14.117   -0.03124      0.0009761
6   7  9.761   -0.88377      0.7810581
7   5  5.906   -0.37262      0.1388434
8   5  5.815   -0.33791      0.1141816

$X^2 = \sum_i (E_i-O_i)^2/E_i$

(chisq = sum((obs-expec)^2/expec))
[1] 5.414413
(df = length(obs)-1-1) # lose an additional df for parameter estimate
[1] 7
(pvalue=pchisq(chisq,df))
[1] 0.3904736

Tanto o diagnóstico quanto o valor-p não mostram falta de ajuste aqui ... o que esperávamos, já que os dados que geramos na verdade eram Poisson.

Edit: aqui está um link para o blog de Rick Wicklin, que discute o enredo de Poissonness e fala sobre implementações no SAS e Matlab

http://blogs.sas.com/content/iml/2012/04/12/the-poissonness-plot-a-goodness-of-fit-diagnostic/

Edit2: Se eu entendi direito, o gráfico de Poissonness modificado da referência de 1985 seria *:

y=rpois(100,5)
n=length(y)
(x=table(y))
k=as.numeric(names(x))
x=as.vector(x)
x1 = ifelse(x==0,NA,ifelse(x>1,x-.8*x/n-.67,exp(-1)))
plot(k,log(x1)+lfactorial(k))

* Eles também ajustam a interceptação, mas eu não fiz isso aqui; isso não afeta a aparência da plotagem, mas você deve tomar cuidado se implementar alguma coisa a partir da referência (como os intervalos de confiança) se fizer isso de maneira diferente da abordagem deles.

(Para o exemplo acima, a aparência quase não muda a partir do primeiro gráfico de Poissonness.)

Execute o teste de ajuste do qui-quadrado. No caso de dados de contagem, podemos usar goodfit()incluído no pacote vcd. Observe que se o valor de p for maior que 0,05, não podemos rejeitar h0: o processo é um processo de Poisson. Ou então, não é um processo de Poisson.

# load the vcd package
library(vcd) ## loading vcd package

# generate two processes for test
set.seed(2014);y=rpois(200,5)
set.seed(2014);y=rnorm(100, 5, 0.3) # goodfit asks for non-negative values
# output the results
gf = goodfit(y,type= "poisson",method= "ML")
plot(gf,main="Count data vs Poisson distribution")
summary(gf)

# to automatically get the pvalue
gf.summary = capture.output(summary(gf))[[5]]
pvalue = unlist(strsplit(gf.summary, split = " "))
pvalue = as.numeric(pvalue[length(pvalue)]); pvalue

# to mannualy compute the pvalue
chisq = sum(  (gf$observed-gf$fitted)^2/gf$fitted )

df = length(gf$observed)-1-1
pvalue = pchisq(chisq,df)
pvalue