Hay tres cuestiones aquí:
La simulación probablemente no hace exactamente lo que parece.
El análisis no evaluar si la relación entre la media y la varianza es lineal.
¿Qué se puede hacer para verificar si los datos aparecen de Poisson?
Muchas de las respuestas a #3 han aparecido en otros lugares en nuestro sitio, así que me voy a centrar en las dos primeras cuestiones.
Lo que la simulación no
Las líneas de código
pos<-rpois(100,5)
...
samp<-sample(pos, 10)
dibujar $100$ iid los valores de una distribución de Poisson con una media de $5$ y, a continuación, obtener al azar $10$ de los valores (sin reemplazo). Iterando el segundo comando es una forma de análisis de remuestreo. En realidad no prueba lo que sucede cuando los datos son repetidamente extraídas de la misma subyacente distribución de Poisson; en su lugar, se estudia la distribución de la $100$ de los valores obtenidos en la primera línea, que se fija para la duración de la simulación.
Para el estudio de una distribución de Poisson, necesitamos dibujar repetidas muestras independientes:
set.seed(17) # Create reproducible results
mu <- 3 # Mean
n <- 10 # Sample size
n.trials <- 10000 # Number of trials
sim <- replicate(n.trials, {x <- rpois(n, mu); c(mean(x), var(x))})
En este punto, sim es una matriz con dos filas: una para los medios, la siguiente, de las desviaciones. Cada columna resume separado de simulación (ergo, hay $10000$ columnas en este ejemplo).
Lo que el análisis
Aquí está en una forma compacta:
xy <- apply(sim, 1, function(x) quantile(x, probs=seq(0.001, 0.999, length.out=99)))
Esto crea una matriz de los percentiles de las dos líneas de salida, con lo que el seguimiento de sus distribuciones. Una parcela es iluminadora:
plot(xy, xlab="mean", ylab="variance",col=hsv(0, 0, 0.25, .4))
![Simulation Q-Q plot]()
Como Q-Q plot, este compara sistemáticamente los percentiles de una distribución (simulación de la varianza de las muestras de $10$ Poisson valores) a los percentiles de otra distribución (simulación de la media de las muestras de $10$ Poisson valores).
La curva no es muy recta, como podemos encontrar en la regresión de los valores de y en los valores de x con un término cuadrático incluye:
z <- xy[,1]^2 # Quadratic term
fit <- lm(xy[,2] ~ xy[,1] + z) # Regression of y-values against x and x^2
#
# Overplot the fit in red.
x0 <- seq(mu - 2 * sqrt(mu/n), mu + 2 * sqrt(mu/n), length.out=99)
u <- cbind(rep(1,length(x0)), x0, x0^2) %*% coefficients(fit)
lines(x0, u, lwd=2, col=hsv(0, .8, .8, 1))
![Q-Q plot with fit]()
El término cuadrático es muy importante y tiene un gran coeficiente:
>coefficients(fit)[3]
z
0.8779665
Animo a los lectores interesados a jugar con los parámetros (mu, ny n.trials) para ver cómo las cosas (no) cambiar. Considerar también la posibilidad de cambiar rpois en la simulación a otra cosa (como rnorm).
Algo de teoría
La media de $10$ independiente se nutre de muchas distribuciones tendrá aproximadamente una distribución normal. Poisson (con una media de $5$) no es la excepción: su media debe tener una distribución normal de media de $5$ y la varianza $5/10$.
La CDF de una distribución de Poisson con considerable de la media (es decir, $3$ o más) se lo hará a sí mismo de comenzar la aproximación de una distribución normal, con la aproximación mejora para los grandes medios. En consecuencia, podemos esperar que la varianza de una muestra de $10$ independiente de Poisson dibuja comenzará a asumir la distribución de la varianza de una muestra de $10$ independiente se basa en una comparables distribución normal. (La varianza de la distribución es proporcional a la distribución chi-squared.) Podemos comprobar mediante el cálculo de la exacta percentiles de estos aproximación de distribuciones y overplotting ellos en los valores simulados:
q <- (1:100 - 0.5) / 100 # An array of percentiles for plotting
lines(qnorm(q, mean=mu, sd=sqrt(mu/n)), mu*qchisq(q, df=n-1) / (n-1), col="Blue")
![Theoretical comparison]()
Incluso para esta bastante pequeña media ($5$) y el modesto tamaño de la muestra ($10$), el ajuste es extraordinariamente buena. Por desgracia, esto demuestra que debemos esperar un aspecto similar parcela de casi cualquier moderado tamaño de la muestra para diferentes distribuciones subyacentes, lo que demuestra que esta trama nos dice poco o nada de lo que es útil acerca de si los datos provienen de una distribución de Poisson.
¿Qué podemos hacer?
"Comprobar" si una distribución tiene una forma particular puede ser mejorado mediante la indicación de las formas en que los datos pueden apartarse de esa forma.
Por ejemplo, distribuciones de Poisson son especiales ejemplos de distribuciones binomial negativa, que puede considerarse overdispersed distribuciones de Poisson. Específicamente, si obtenemos una secuencia de Poisson valores donde el subyacente de la media de la distribución de Poisson cambios en una forma que se aproxima a una Gamma de distribución, entonces el conjunto de datos recogidos es descrito por la binomial negativa.
En tal escenario, nos puede caber la binomial negativa de parámetros (utilizando máxima verosimilitud, por ejemplo) y comprobar si el "no-Poissonness" parámetro difiere significativamente de un valor correspondiente a una distribución de Poisson. En el interés de espacio, y para mantener esta respuesta razonablemente primaria, no voy a seguir esta opción aquí.
Menos formalmente, no puede simplemente aplicar una Poisson modelo lineal generalizado y mira su diagnóstico:
set.seed(17)
y <- rnbinom(100, mu=5, size=1) # NOT Poisson
x <- rpois(100, 5) # Poisson
fit.y <- glm(y ~ 1, family=poisson())
fit.x <- glm(x ~ 1, family=poisson())
plot(fit.y) # Includes a residual Q-Q plot
plot(fit.x) # Includes a residual Q-Q plot
El residual Q-Q plot para la distribución de Poisson de datos es casi lineal, mientras que la de la binomial negativa de datos tiene algunas pronunciado extremos
![Negative binomial residuals]()
![Poisson residuals]()
Como se puede ver, se necesita un moderado tamaño de la muestra para detectar estas diferencias, pero al menos que se muestran arriba.
