¿Cómo saber si un dato sigue una distribución de Poisson en R?

Soy estudiante de pregrado y tengo un proyecto para mi clase de probabilidad. Básicamente, tengo un conjunto de datos sobre los huracanes que afectaron a mi país durante una serie de años.

En mi Libro de probabilidad, (Probabilidad y estadística con R) hay un ejemplo (no completo) de cómo verificar si los datos siguen una distribución de Poisson, comienzan a tratar de probar que se siguen estos 3 criterios: (De mi libro, página 120 (criterio) página 122-123 ejemplo)

1- El número de resultados en intervalos no superpuestos es independiente. En otras palabras, el número de resultados en el intervalo de tiempo (0, t] es independiente del número de resultados en el intervalo de tiempo (t, t + h], h> 0

2- La probabilidad de dos o más resultados en un intervalo suficientemente corto es prácticamente cero. En otras palabras, siempre que h sea suficientemente pequeño, la probabilidad de obtener dos o más resultados en el intervalo (t, t + h] es insignificante en comparación con la probabilidad de obtener uno o cero resultados en el mismo intervalo de tiempo.

3- La probabilidad de exactamente un resultado en un intervalo suficientemente corto o una región pequeña es proporcional a la duración del intervalo o región. En otras palabras, la probabilidad de un resultado en un intervalo de longitud h es lambda * h.

Pero el criterio 3 se deja "como ejercicio".

A- ¿Puede alguien decirme si hay una forma más "fácil" de ver si mi conjunto de datos sigue una distribución de Poisson?

B- ¿Puede alguien explicarme los criterios 1 y 3 con algún tipo de ejemplo (si es con R, fantástico)?

¡Gracias!

Nota: Perdón por la larga publicación. Además, tengo que convertir los datos para tener una tabla como:

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

Hay un número infinito de formas para que una distribución sea ligeramente diferente de una distribución de Poisson; no puede identificar que un conjunto de datos se extrae de una distribución de Poisson. Lo que puede hacer es buscar inconsistencias con lo que debería ver con un Poisson, pero la falta de inconsistencia obvia no lo convierte en Poisson.

Sin embargo, lo que está hablando allí al verificar esos tres criterios no es verificar que los datos provienen de una distribución de Poisson por medios estadísticos (es decir, al observar los datos), sino al evaluar si el proceso por el que se generan los datos satisface condiciones de un proceso de Poisson; si todas las condiciones se mantuvieron o casi se cumplen (y eso es una consideración del proceso de generación de datos), podría tener algo de o muy cercano a un proceso de Poisson, que a su vez sería una forma de obtener datos extraídos de algo cercano a un Distribución de veneno.

Pero las condiciones no se sostienen de varias maneras ... y el más alejado de ser cierto es el número 3. No hay ninguna razón particular sobre esa base para afirmar un proceso de Poisson, aunque las violaciones pueden no ser tan malas que los datos resultantes estén lejos. de Poisson

Así que volvemos a los argumentos estadísticos que provienen del examen de los datos en sí. ¿Cómo mostrarían los datos que la distribución era Poisson, en lugar de algo así?

Como se mencionó al principio, lo que puede hacer es verificar si los datos no son obviamente inconsistentes con la distribución subyacente que es Poisson, pero eso no le dice que se extraen de un Poisson (ya puede estar seguro de que son no).

Puede hacer esta comprobación a través de pruebas de bondad de ajuste.

El chi-cuadrado que se mencionó es uno de esos, pero yo mismo no recomendaría la prueba de chi-cuadrado para esta situación **; Tiene bajo poder contra desviaciones interesantes. Si su objetivo es tener un buen poder, no lo obtendrá de esa manera (si no le importa el poder, ¿por qué lo probaría?). Su valor principal está en la simplicidad, y tiene valor pedagógico; fuera de eso, no es competitivo como prueba de bondad de ajuste.

** Agregado en una edición posterior: ahora que está claro que esto es tarea, las posibilidades de que hagas una prueba de ji cuadrado para verificar que los datos no sean inconsistentes con un Poisson aumentan bastante. Vea mi ejemplo de prueba de bondad de ajuste de chi-cuadrado, realizada debajo del primer diagrama de Poissonness

Las personas a menudo hacen estas pruebas por la razón incorrecta (por ejemplo, porque quieren decir 'por lo tanto, está bien hacer otra cosa estadística con los datos que suponen que los datos son de Poisson'). La verdadera pregunta es "¿qué tan mal podría ir eso?" ... y la bondad de las pruebas de ajuste no son de mucha ayuda con esa pregunta. A menudo, la respuesta a esa pregunta es, en el mejor de los casos, una que es independiente (/ casi independiente) del tamaño de la muestra, y en algunos casos, una con consecuencias que tienden a desaparecer con el tamaño de la muestra ... mientras que una prueba de bondad de ajuste es inútil con pequeñas muestras (donde el riesgo de violación de supuestos suele ser mayor).

Si debe probar una distribución de Poisson, hay algunas alternativas razonables. Una sería hacer algo similar a una prueba de Anderson-Darling, basada en la estadística AD pero usando una distribución simulada bajo nulo (para tener en cuenta los problemas gemelos de una distribución discreta y que debe estimar los parámetros).

Una alternativa más simple podría ser una prueba de suavidad para la bondad de ajuste: se trata de una colección de pruebas diseñadas para distribuciones individuales modelando los datos utilizando una familia de polinomios que son ortogonales con respecto a la función de probabilidad en el nulo. Las alternativas de bajo orden (es decir, interesantes) se prueban probando si los coeficientes de los polinomios por encima de la base son diferentes de cero, y estos generalmente pueden ocuparse de la estimación de parámetros omitiendo los términos de orden más bajo de la prueba. Hay tal prueba para el Poisson. Puedo desenterrar una referencia si la necesitas.

También puede usar la correlación (o, para ser más parecido a una prueba de Shapiro-Francia, quizás ) en un gráfico de Poissonness, por ejemplo, un gráfico de vs (ver Hoaglin, 1980) - como estadística de prueba.log ( x k ) + log ( k ! ) $n(1-r^2)$$\log(x_k)+\log(k!)$$k$

Aquí hay un ejemplo de ese cálculo (y gráfico), hecho en 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))

Aquí está la estadística que sugerí que podría usarse para una prueba de bondad de ajuste de un Poisson:

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

Por supuesto, para calcular el valor p, también necesitaría simular la distribución del estadístico de prueba bajo el valor nulo (y no he discutido cómo se podría tratar con los recuentos cero dentro del rango de valores). Esto debería producir una prueba razonablemente poderosa. Existen numerosas otras pruebas alternativas.

Aquí hay un ejemplo de cómo hacer un diagrama de Poissonness en una muestra de tamaño 50 a partir de una distribución geométrica (p = .3):

Como puede ver, muestra un 'pliegue' claro, que indica no linealidad

Las referencias para el diagrama de Poissonness serían:

David C. Hoaglin (1980),
"Una trama de Poissonness",
The American Statistician
Vol. 34, núm. 3 (agosto, pág. 146-149).

y

Hoaglin, D. y J. Tukey (1985),
"9. Comprobación de la forma de las distribuciones discretas",
Exploración de tablas de datos, tendencias y formas ,
(Hoaglin, Mosteller y Tukey eds)
John Wiley & Sons

La segunda referencia contiene un ajuste a la trama para recuentos pequeños; probablemente quieras incorporarlo (pero no tengo la referencia a mano).

Ejemplo de hacer una prueba de bondad de ajuste de chi-cuadrado:

Además de realizar la bondad de ajuste de chi-cuadrado, la forma en que normalmente se esperaría que se hiciera en muchas clases (aunque no de la forma en que lo haría):

1: comenzando con sus datos, (que tomaré como los datos que generé aleatoriamente en 'y' arriba, genere la tabla de conteos:

(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 el valor esperado en cada celda, suponiendo un 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: tenga en cuenta que las categorías finales son pequeñas; esto hace que la distribución de chi-cuadrado sea menos buena como una aproximación a la distribución del estadístico de prueba (una regla común es que desea valores esperados de al menos 5, aunque numerosos documentos han demostrado que esa regla es innecesariamente restrictiva; lo tomaré cerca, pero el enfoque general se puede adaptar a una regla más estricta). Contraiga las categorías adyacentes, de modo que los valores mínimos esperados no estén al menos muy por debajo de 5 (una categoría con una cuenta regresiva esperada cerca de 1 de más de 10 categorías no es tan mala, dos son bastante dudosas). También tenga en cuenta que aún no hemos tenido en cuenta la probabilidad más allá de "10", por lo que también debemos incorporar eso:

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 manera similar, colapsar las categorías observadas:

(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]

5: Poner en una tabla, (opcionalmente) junto con la contribución a chi-cuadrado y el residuo de Pearson (la raíz cuadrada firmada de la contribución), estos pueden ser útiles al intentar ver dónde no encaja muy bien:$(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 el diagnóstico como el valor p no muestran falta de ajuste aquí ... lo que esperaríamos, ya que los datos que generamos en realidad eran Poisson.

Editar: aquí hay un enlace al blog de Rick Wicklin que analiza la trama de Poissonness y habla sobre implementaciones en SAS y Matlab

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

Edit2: si tengo razón, la gráfica de Poissonness modificada de la referencia de 1985 sería *:

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

* En realidad también ajustan la intercepción, pero no lo he hecho aquí; no afecta la apariencia de la trama, pero debe tener cuidado si implementa algo más de la referencia (como los intervalos de confianza) si lo hace de manera diferente a su enfoque.

(Para el ejemplo anterior, la apariencia apenas cambia desde la primera gráfica de Poissonness).

Realice la prueba de bondad de ajuste chi-cuadrado. En el caso de los datos de conteo, podemos usarlos goodfit()incluidos en el paquete vcd. Tenga en cuenta que si el valor p es mayor que 0.05, no podemos rechazar h0: el proceso es un proceso de Poisson. O de lo contrario, no es un proceso 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