Comment savoir si une donnée suit une distribution de Poisson dans R?

Je suis un étudiant de premier cycle et j'ai un projet pour ma classe de probabilité. Fondamentalement, j'ai un ensemble de données sur les ouragans qui ont touché mon pays pendant une série d'années.

Dans mon livre de probabilités, (Probabilité et statistiques avec R) il y a un exemple (non complet) de la façon de vérifier si les données suivent une distribution de Poisson, ils commencent à essayer de prouver que ces 3 critères sont suivis: (De mon livre, page 120 (critères) page 122-123 exemple)

1- Le nombre de résultats dans des intervalles sans chevauchement sont indépendants. En d'autres termes, le nombre de résultats dans l'intervalle de temps (0, t] est indépendant du nombre de résultats dans l'intervalle de temps (t, t + h], h> 0

2- La probabilité de deux résultats ou plus dans un intervalle suffisamment court est pratiquement nulle. En d'autres termes, à condition que h soit suffisamment petit, la probabilité d'obtenir deux résultats ou plus dans l'intervalle (t, t + h] est négligeable par rapport à la probabilité d'obtenir un ou zéro résultat dans le même intervalle de temps.

3- La probabilité d'exactement un résultat dans un intervalle suffisamment court ou une petite région est proportionnelle à la longueur de l'intervalle ou de la région. En d'autres termes, la probabilité d'un résultat dans un intervalle de longueur h est lambda * h.

Mais le critère 3 est laissé "comme exercice".

A- Quelqu'un peut-il me dire s'il existe un moyen plus "facile" de voir si mon jeu de données suit une distribution de Poisson?

B- Quelqu'un peut-il m'expliquer les critères 1 et 3 avec un certain type d'exemple (si c'est avec R, fantastique)?

Merci!

Remarque: Désolé pour le long post. De plus, je dois convertir les données pour avoir un tableau comme:

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

Il existe un nombre infini de façons pour une distribution d'être légèrement différente d'une distribution de Poisson; vous ne pouvez pas identifier qu'un ensemble de données est tiré d'une distribution de Poisson. Ce que vous pouvez faire, c'est rechercher une incohérence avec ce que vous devriez voir avec un Poisson, mais un manque d'incohérence évidente n'en fait pas un poisson.

Cependant, ce dont vous parlez en vérifiant ces trois critères, ce n'est pas de vérifier que les données proviennent d'une distribution de Poisson par des moyens statistiques (c'est-à-dire en regardant les données), mais en évaluant si le processus par lequel les données sont générées satisfait le conditions d'un processus de Poisson; si toutes les conditions étaient remplies ou presque (et c'est une considération du processus de génération de données), vous pourriez avoir quelque chose d'un processus de Poisson ou très proche de celui-ci, ce qui serait à son tour un moyen d'obtenir des données tirées de quelque chose de proche d'un Distribution de Poisson.

Mais les conditions ne sont pas réunies de plusieurs manières ... et le plus éloigné de la vérité est le numéro 3. Il n'y a aucune raison particulière sur cette base pour affirmer un processus de Poisson, bien que les violations ne soient pas si mauvaises que les données résultantes soient loin de Poisson.

Nous revenons donc aux arguments statistiques qui proviennent de l'examen des données elles-mêmes. Comment les données montreraient-elles que la distribution était de Poisson, plutôt que quelque chose comme ça?

Comme mentionné au début, ce que vous pouvez faire est de vérifier si les données ne sont pas manifestement incompatibles avec la distribution sous-jacente étant Poisson, mais cela ne vous dit pas qu'elles sont tirées d'un Poisson (vous pouvez déjà être sûr qu'elles sont ne pas).

Vous pouvez effectuer cette vérification via des tests d'adéquation.

Le khi carré mentionné en est un, mais je ne recommanderais pas moi-même le test du khi carré pour cette situation **; il a une faible puissance contre des écarts intéressants. Si votre objectif est d'avoir un bon pouvoir, vous ne l'obtiendrez pas de cette façon (si vous ne vous souciez pas du pouvoir, pourquoi testeriez-vous?). Sa valeur principale réside dans la simplicité et il a une valeur pédagogique; en dehors de cela, ce n'est pas compétitif comme test de qualité de l'ajustement.

** Ajouté dans une édition ultérieure: Maintenant qu'il est clair qu'il s'agit de devoirs, les chances que vous vous attendiez à faire un test du chi carré pour vérifier les données ne sont pas incompatibles avec un Poisson augmentent considérablement. Voir mon exemple de test de qualité d'ajustement chi carré, effectué sous le premier tracé Poissonness

Les gens font souvent ces tests pour la mauvaise raison (par exemple parce qu'ils veulent dire "donc c'est correct de faire quelque chose de statistique avec les données qui supposent que les données sont de Poisson"). La vraie question est «à quel point cela pourrait-il mal tourner? ... et la qualité des tests d'ajustement n'aide pas vraiment à répondre à cette question. Souvent, la réponse à cette question est, au mieux, indépendante (/ presque indépendante) de la taille de l'échantillon - et dans certains cas, une avec des conséquences qui ont tendance à disparaître avec la taille de l'échantillon ... alors qu'un test de qualité de l'ajustement est inutile avec de petits échantillons (où votre risque de violation des hypothèses est souvent le plus élevé).

Si vous devez tester une distribution de Poisson, il existe quelques alternatives raisonnables. L'une consisterait à faire quelque chose qui s'apparente à un test d'Anderson-Darling, basé sur la statistique AD mais en utilisant une distribution simulée sous la valeur nulle (pour tenir compte des problèmes jumeaux d'une distribution discrète et que vous devez estimer les paramètres).

Une alternative plus simple pourrait être un test lisse de qualité de l'ajustement - il s'agit d'une collection de tests conçus pour des distributions individuelles en modélisant les données à l'aide d'une famille de polynômes orthogonaux par rapport à la fonction de probabilité dans le zéro. Les alternatives d'ordre faible (c'est-à-dire intéressantes) sont testées en testant si les coefficients des polynômes au-dessus du polynôme de base sont différents de zéro, et ceux-ci peuvent généralement traiter l'estimation des paramètres en omettant les termes d'ordre le plus bas du test. Il y a un tel test pour le Poisson. Je peux trouver une référence si vous en avez besoin.

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

Voici un exemple de ce calcul (et tracé), effectué dans 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))

Voici la statistique que j'ai suggérée pourrait être utilisée pour un test d'adéquation d'un Poisson:

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

Bien sûr, pour calculer la valeur de p, vous devez également simuler la distribution de la statistique de test sous la valeur nulle (et je n'ai pas discuté de la façon dont on pourrait traiter les nombres nuls dans la plage de valeurs). Cela devrait donner un test raisonnablement puissant. Il existe de nombreux autres tests alternatifs.

Voici un exemple de réalisation d'un tracé Poissonness sur un échantillon de taille 50 à partir d'une distribution géométrique (p = 0,3):

Comme vous le voyez, il affiche un «kink» clair, indiquant la non-linéarité

Les références pour le complot Poissonness seraient:

David C. Hoaglin (1980),
"A Poissonness Plot",
The American Statistician
Vol. 34, n ° 3 (août,), pp. 146-149

et

Hoaglin, D. et J. Tukey (1985),
"9. Checking the Shape of Discrete Distributions",
Exploring Data Tables, Trends and Shapes ,
(Hoaglin, Mosteller & Tukey eds)
John Wiley & Sons

La deuxième référence contient un ajustement du tracé pour les petits dénombrements; vous voudrez probablement l'incorporer (mais je n'ai pas la référence à portée de main).

Exemple de test de qualité d'ajustement chi carré:

En plus de réaliser la qualité d'ajustement du chi carré, la façon dont cela devrait normalement être fait dans de nombreuses classes (mais pas comme je le ferais):

1: à partir de vos données, (que je considérerai comme étant les données que j'ai générées aléatoirement en «y» ci-dessus, générez le tableau des décomptes:

(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: calculer la valeur attendue dans chaque cellule, en supposant un Poisson ajusté par 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: notez que les catégories finales sont petites; cela rend la distribution du chi carré moins bonne en tant qu'approximation de la distribution de la statistique de test (une règle courante est que vous voulez des valeurs attendues d'au moins 5, bien que de nombreux articles ont montré que cette règle soit inutilement restrictive; je vais la prendre mais l'approche générale peut être adaptée à une règle plus stricte). Réduisez les catégories adjacentes, de sorte que les valeurs minimales attendues ne soient au moins pas trop loin en dessous de 5 (une catégorie avec un compte à rebours prévu près de 1 sur plus de 10 catégories n'est pas trop mauvaise, deux est assez limite). Notez également que nous n'avons pas encore pris en compte la probabilité au-delà de "10", nous devons donc également intégrer cela:

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: de même, effondrement des catégories sur les observés:

(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

Les diagnostics et la valeur p ne montrent aucun manque d'ajustement ici ... ce à quoi nous nous attendions, car les données que nous avons générées étaient en fait de Poisson.

Edit: voici un lien vers le blog de Rick Wicklin qui discute de l'intrigue Poissonness et parle des implémentations dans SAS et Matlab

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

Edit2: Si je comprends bien, le tracé Poissonness modifié de la référence de 1985 serait *:

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

* Ils ajustent également l'interception, mais je ne l'ai pas fait ici; cela n'affecte pas l'apparence de l'intrigue, mais vous devez faire attention si vous implémentez autre chose à partir de la référence (comme les intervalles de confiance) si vous le faites différemment de leur approche.

(Pour l'exemple ci-dessus, l'apparence ne change guère par rapport à la première intrigue Poissonness.)

Effectuez le test d'adéquation du chi carré. En cas de données de comptage, nous pouvons utiliser goodfit()inclus dans le package vcd. Notez que si la valeur p est supérieure à 0,05, nous ne pouvons pas rejeter h0: le processus est un processus de Poisson. Ou bien, ce n'est pas un processus 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