Esta respuesta presenta dos soluciones: Las correcciones de Sheppard y una estimación de máxima verosimilitud. Ambos coinciden en una estimación de la desviación típica: $7.70$ para la primera y $7.69$ para el segundo (cuando se ajusta para que sea comparable al estimador "insesgado" habitual).
Correcciones de Sheppard
Las "correcciones de Sheppard" son fórmulas que ajustan los momentos calculados a partir de datos binned (como éstos) donde
-
se supone que los datos se rigen por una distribución soportada en un intervalo finito $[a,b]$
-
ese intervalo se divide secuencialmente en intervalos iguales de anchura común $h$ que es relativamente pequeño (ningún contenedor contiene una gran proporción de todos los datos)
-
la distribución tiene una función de densidad continua.
Se derivan de la fórmula de la suma de Euler-Maclaurin, que aproxima las integrales en términos de combinaciones lineales de valores del integrando en puntos regularmente espaciados, y por lo tanto de aplicación general (y no sólo a las distribuciones Normales).
Aunque en sentido estricto una distribución Normal es no apoyado en un intervalo finito, a una aproximación extremadamente cercana lo es. Esencialmente, toda su probabilidad está contenida dentro de siete desviaciones estándar de la media. Por lo tanto, las correcciones de Sheppard son aplicables a los datos que se supone proceden de una distribución Normal.
Las dos primeras correcciones de Sheppard son
-
Utilice la media de los datos divididos en cubos para la media de los datos (es decir, no es necesaria ninguna corrección para la media).
-
Resta $h^2/12$ a partir de la varianza de los datos agrupados para obtener la varianza (aproximada) de los datos.
¿Dónde está $h^2/12$ ¿De dónde viene? Es igual a la varianza de una variante uniforme distribuida en un intervalo de longitud $h$ . Intuitivamente, entonces, la corrección de Sheppard para el segundo momento sugiere que Al dividir los datos por intervalos -sustituyéndolos por el punto medio de cada intervalo- parece que se añade un valor distribuido de forma aproximadamente uniforme que oscila entre los siguientes valores $-h/2$ y $h/2$ por lo que infla la varianza en $h^2/12$ .
Hagamos los cálculos. Utilizo R para ilustrarlas, empezando por especificar los recuentos y los intervalos:
counts <- c(1,2,3,4,1)
bin.lower <- c(40, 45, 50, 55, 70)
bin.upper <- c(45, 50, 55, 60, 75)
La fórmula adecuada para realizar los recuentos es la siguiente replicando las anchuras de las casillas por las cantidades dadas por los recuentos; es decir, los datos divididos en casillas equivalen a
42.5, 47.5, 47.5, 52.5, 52.5, 57.5, 57.5, 57.5, 57.5, 72.5
Sin embargo, su número, media y varianza pueden calcularse directamente sin tener que expandir los datos de esta forma: cuando un contenedor tiene punto medio $x$ y un recuento de $k$ entonces su contribución a la suma de cuadrados es $kx^2$ . Esto conduce a la segundo de las fórmulas de Wikipedia citadas en la pregunta.
bin.mid <- (bin.upper + bin.lower)/2
n <- sum(counts)
mu <- sum(bin.mid * counts) / n
sigma2 <- (sum(bin.mid^2 * counts) - n * mu^2) / (n-1)
La media ( mu ) es $1195/22 \approx 54.32$ (que no necesita corrección) y la varianza ( sigma2 ) es $675/11 \approx 61.36$ . (Su raíz cuadrada es $7.83$ como se indica en la pregunta). Dado que la anchura común de los contenedores es $h=5$ restamos $h^2/12 = 25/12 \approx 2.08$ de la varianza y sacar su raíz cuadrada, obteniendo $\sqrt{675/11 - 5^2/12} \approx 7.70$ para la desviación típica.
Estimaciones de máxima verosimilitud
Un método alternativo consiste en aplicar una estimación de máxima verosimilitud. Cuando la distribución subyacente supuesta tiene una función de distribución $F_\theta$ (en función de los parámetros $\theta$ a estimar) y el bin $(x_0, x_1]$ contiene $k$ de un conjunto de valores independientes e idénticamente distribuidos de $F_\theta$ entonces la contribución (aditiva) a la logverosimilitud de este bin es
$$\log \prod_{i=1}^k \left(F_\theta(x_1) - F_\theta(x_0)\right) = k\log\left(F_\theta(x_1) - F_\theta(x_0)\right)$$
(véase MLE/verosimilitud del intervalo con distribución lognormal ).
La suma de todos los intervalos da la probabilidad logarítmica $\Lambda(\theta)$ para el conjunto de datos. Como de costumbre, encontramos una estimación $\hat\theta$ que minimiza $-\Lambda(\theta)$ . Esto requiere una optimización numérica y eso se agiliza suministrando buenos valores de partida para $\theta$ . Los siguientes R hace el trabajo para una distribución Normal:
sigma <- sqrt(sigma2) # Crude starting estimate for the SD
likelihood.log <- function(theta, counts, bin.lower, bin.upper) {
mu <- theta[1]; sigma <- theta[2]
-sum(sapply(1:length(counts), function(i) {
counts[i] *
log(pnorm(bin.upper[i], mu, sigma) - pnorm(bin.lower[i], mu, sigma))
}))
}
coefficients <- optim(c(mu, sigma), function(theta)
likelihood.log(theta, counts, bin.lower, bin.upper))$par
Los coeficientes resultantes son $(\hat\mu, \hat\sigma) = (54.32, 7.33)$ .
Recuerde, sin embargo, que para distribuciones normales la estimación de máxima verosimilitud de $\sigma$ (cuando los datos se dan exactamente y no se dividen) es el población de los datos, no la estimación "corregida de sesgo" más convencional en la que la varianza se multiplica por $n/(n-1)$ . Corrijamos entonces (por comparación) la MLE de $\sigma$ , encontrando $\sqrt{n/(n-1)} \hat\sigma = \sqrt{11/10}\times 7.33 = 7.69$ . Esto se compara favorablemente con el resultado de la corrección de Sheppard, que fue $7.70$ .
Verificación de los supuestos
Para visualizar estos resultados podemos trazar la densidad Normal ajustada sobre un histograma:
hist(unlist(mapply(function(x,y) rep(x,y), bin.mid, counts)),
breaks = breaks, xlab="Values", main="Data and Normal Fit")
curve(dnorm(x, coefficients[1], coefficients[2]),
from=min(bin.lower), to=max(bin.upper),
add=TRUE, col="Blue", lwd=2)
![Figure]()
Para algunos, esto podría no parecer una buena opción. Sin embargo, como el conjunto de datos es pequeño (sólo $11$ ), pueden producirse desviaciones sorprendentemente grandes entre la distribución de las observaciones y la verdadera distribución subyacente.
Comprobemos más formalmente el supuesto (realizado por el MLE) de que los datos se rigen por una distribución Normal. Se puede obtener una prueba aproximada de bondad de ajuste a partir de un $\chi^2$ prueba: los parámetros estimados indican la cantidad esperada de datos en cada recipiente; la $\chi^2$ compara los recuentos observados con los esperados. He aquí una prueba en R :
breaks <- sort(unique(c(bin.lower, bin.upper)))
fit <- mapply(function(l, u) exp(-likelihood.log(coefficients, 1, l, u)),
c(-Inf, breaks), c(breaks, Inf))
observed <- sapply(breaks[-length(breaks)], function(x) sum((counts)[bin.lower <= x])) -
sapply(breaks[-1], function(x) sum((counts)[bin.upper < x]))
chisq.test(c(0, observed, 0), p=fit, simulate.p.value=TRUE)
El resultado es
Chi-squared test for given probabilities with simulated p-value (based on 2000 replicates)
data: c(0, observed, 0)
X-squared = 7.9581, df = NA, p-value = 0.2449
El programa ha realizado una prueba de permutación (que es necesaria porque la estadística de la prueba no sigue exactamente una distribución chi-cuadrado: véase mi análisis en Cómo entender los grados de libertad ). Su valor p de $0.245$ que no es pequeño, muestra muy pocos indicios de desviación de la normalidad: tenemos motivos para confiar en los resultados de máxima verosimilitud.
