En los ejemplos que hemos visto siempre hemos necesitado una distribución de referencia que nos permita calcular los cuantiles correspondientes. Sin embargo, es posible que ese supuesto no se cumpla, o que la distribución no sea normal por ejemplo. En estos casos primero se debe de considerar una prueba de hipótesis para revisar las condiciones de nuestros datos. Una familia de problemas no paramétricos surgen a partir de esta observación.
En los temas vistos hasta ahora hemos asumido que nuestros datos siguen una distribución en particular, a menudo la distribución normal. Hemos utilizado esta suposición para realizar inferencias algunas inferencias (intervalos de confianza y pruebas de hipótesis). Pero, ¿qué pasa si nuestros datos no siguen la distribución que supusimos? ¿Cómo sabríamos si nuestra suposición es correcta?
Aquí es donde entran en juego las pruebas de bondad de ajuste. Estas pruebas nos permiten comparar nuestros datos con una distribución de referencia teórica y hacer una decisión objetiva sobre si la distribución de nuestros datos se ajusta suficientemente bien a la distribución teórica.
Si encontramos que nuestros datos no se ajustan a la distribución supuesta, entonces tendríamos que considerar el uso de métodos estadísticos no paramétricos. De esta forma se puede evitar algunos supuestos sobre nuestro problema.
En este capítulo, revisaremos toda una familia de problemas y métodos estadísticos no paramétricos diseñados específicamente para trabajar con datos que no cumplen con las suposiciones de las técnicas estadísticas paramétricas estándar.
Imagina que tienes datos categóricos. Esta es una clase de datos donde la variable en cuestión puede caer en un número finito de categorías o estados.
Los datos categóricos son una clase de datos donde la variable en cuestión puede caer en un número finito de categorías o estados.
En un escenario como este, cada categoría tiene una probabilidad asociada, representada por \(p_i\), donde \(i\) es el número de la categoría. Así, tenemos \(k\) categorías y \(p_i\) para cada \(i\) de 1 a \(k\). El total de estas probabilidades debe sumar 1, es decir, \(\displaystyle \sum_{i=1}^kp_i = 1\).
Además, se pueden tener supuestos teóricos sobre esas probabilidades para cada categoría. Podemos definir, \(p_1^0,\dotsc,p_k^0\), y nuevamente, su suma también debe ser 1.
Así, puedes plantear una prueba de hipótesis:
Ejemplo 13.1 Por ejemplo, supongamos que estás estudiando la distribución de tipos de sangre en la población. La variable puede caer en una de las cuatro categorías posibles:
| Categoría | Tipo de sangre |
|---|---|
| 1 | A |
| 2 | B |
| 3 | AB |
| 4 | O |
Si sigues con nuestro ejemplo, podrías querer probar la siguiente hipótesis sobre la distribución de los tipos de sangre:
| Categoría | Tipo de sangre | Hipótesis (\(p_i ^{0}\)) |
|---|---|---|
| 1 | A | 1/3 |
| 2 | B | 1/8 |
| 3 | AB | 1/24 |
| 4 | O | 1/2 |
En este contexto, suponga que recoges una muestra de \(n\) individuos. Deja \(N_i\) ser el número de individuos en tu muestra que caen en la categoría \(i\). Por lo tanto, \(\sum _{i=1}^kN_i = n\).
Es importante notar que el conjunto de estas cuentas de categorías, \((N_1,\dots,N_k)\), sigue una distribución multinomial con parámetros \((p_1, \dots, p_k)\), i.e., \[ (N_1,\dots,N_k) \sim \text{Multinomial}(n,p_1,\dots,p_k). \]
Observación. Una distribución multinomial tiene la siguiente forma:
\[\begin{align*} \operatorname{Pr}\left(X_{1}=x_{1}, \ldots, X_{k}=x_{k}\right) = \left\{\begin{array}{rlr} \binom {n} {x_{1}, \ldots, x_{k}} p_{1}^{x_{1}} \cdots p_{k}^{x_{k}} & \text { if } x_{1}+\cdots+x_{k}=n \\ 0 & \text { en otro caso. } \end{array}\right. \end{align*}\]
donde
\[\begin{equation*} \binom{n}{x_{1}, \ldots, x_{k}} =\frac{n !}{x_{1} ! x_{2} ! \cdots x_{k} !} \end{equation*}\]
El número esperado de elementos en la categoría \(i\) es \(n\cdot p_i^0\). Si la diferencia \(N_i -np_i^0\) es cercana a 0 para cada categoría, es un indicio de que nuestra hipótesis nula \(H_0\) podría ser cierta.
El estadístico \(\chi^2\) se define como:
\[Q = \sum_{i=1}^k\dfrac{(N_i-np_i^0)^2}{np_i^0}.\]
En 1900, Karl Pearson demostró que si el tamaño de la muestra \(n\) es grande y el número de categorías \(k\) es “relativamente” pequeño en comparación con \(n\), bajo la hipótesis nula \(H_0\), la estadística \(Q\) sigue una distribución chi-cuadrada con \(k-1\) grados de libertad:
\[Q \xrightarrow[H_0]{}\chi^2_{k-1}.\]
Para un nivel de significancia dado \(\alpha\), rechazamos la hipótesis nula si \(Q\geq c\), donde \(c\) es el valor crítico tal que:
\[\mathbb P_{H_0}[Q\geq c]\le \alpha\implies c = F^{-1}_{\chi^2_{k-1}}(1-\alpha)\]
La fórmula de Q se puede reescribir como:
\[Q = \sum_{i=1}^k\dfrac{(\text{observado}_{i} - \text{esperado}_{i})^2}{\text{esperado}_{i}}.\]
Es una medida que cuantifica cuánto varían los datos observados con respecto a lo que se esperaría si la hipótesis nula fuera cierta. Para entenderla mejor vamos a desglosar cada uno de sus componentes:
Ejemplo 13.2 Vamos a continuar con el ejemplo de los tipos de sangre en una muestra de 6004 personas de raza blanca en California. Después de recolectar los datos, se construye esta tabla:
| Grupo | Observado (\(n_i\)) | Teórico (\(p_i ^{0}\)) |
|---|---|---|
| A | 2162 | 1/3 |
| B | 738 | 1/8 |
| AB | 228 | 1/24 |
| O | 2876 | 1/2 |
Vamos a probar la hipótesis nula \(H_0: p_i = p_i^0\) para \(i=1,2,3,4\).
Los valores esperados bajo la hipótesis nula serían:
Calculamos el estadístico \(Q\):
\[\begin{multline*} Q = \dfrac{(2162-2001.3)^2}{2001.3} + \dfrac{(738-750.5)^2}{750.5} \\ - \dfrac{(228-250.2)^2}{250.2} + \dfrac{(2876-3002)^2}{3002} = 20.37. \end{multline*}\]
El valor-p asociado a esto estadístico \(Q\) es \(F_{\chi^2_{3}}(20.37) = 1.42\times 10^{-4}\).
Podemos realizar este cálculo en R con la función chisq.test:
observado <- c(2162, 738, 228, 2876)
probabilidad_hipotetica <- c(1 / 3, 1 / 8, 1 / 24, 1 / 2)
chisq.test(x = observado, p = probabilidad_hipotetica)
Chi-squared test for given probabilities
data: observado
X-squared = 20.359, df = 3, p-value = 0.000143Ejemplo 13.3 Consideremos un conjunto de 100 observaciones, \(X_i\), con \(i=1,2,\dots,100\) y \(0 < X_i < 1\). Supongamos que estas observaciones siguen una función de densidad \(f\), la cual es continua. Nos gustaría probar las siguientes hipótesis: \[H_0: f=\text{Unif}(0,1) \text{ vs } H_1: f \ne\text{Unif}(0,1). \]
Para esto, dividimos el intervalo [0,1] en 20 niveles, de tal manera que una observación \(X_j\) pertenece al nivel \(i\) si se cumple que:
\[\dfrac{i-1}{20}\leq X_j <\dfrac{i}{20}.\]
Es decir, construimos una nueva estructura de frecuencias de la forma:
| \(\text{Nivel}\) | \(1\) | \(2\) | \(\cdots\) | \(20\) |
|---|---|---|---|---|
| \(\text{Frecuencia}\) | \(N_1\) | \(N_2\) | \(\cdots\) | \(N_{20}\) |
donde \(N_i\) es el número de observaciones que están en el intervalo \(i\).
| \(i\) | \(X_i\) | \(\text{Grupo}\) |
|---|---|---|
| \(1\) | \(X_1\) | \(2\) |
| \(2\) | \(X_2\) | \(4\) |
| \(3\) | \(X_3\) | \(17\) |
| \(\vdots\) | \(\vdots\) | \(\vdots\) |
| \(100\) | \(X_{100}\) | \(20\) |
Las hipótesis anteriores son equivalentes a suponer que cada nivel tiene la misma probabilidad de ser escogido, es decir,
\[H_0: p_i = \dfrac{1}{20}, \;i=1,\dots,20.\]
El número esperado de observaciones en cada nivel bajo la hipótesis nula es \(np_i^0 = 100\cdot\dfrac 1{20} = 5,;i = 1,\dots,20\).
Definimos el estadístico de prueba como
\[Q = \sum_{i=1}^{20}\dfrac{(N_i-5)^2}{5}.\]
Rechazamos la hipótesis \(f = \text{Unif}(0,1)\) si \(Q>\chi^2_{19, 1-\alpha}\).
Este método es aplicable para cualquier tipo de distribución. Para aplicarlo, debemos seguir estos pasos:
Ejemplo 13.4 Supongamos que tenemos 23 observaciones correspondientes al tiempo de vida útil de ciertas partes mecánicas de automóviles. Trabajaremos con el logaritmo del tiempo de vida de estos dispositivos. Los datos son los siguientes:
x <- c(
17.88, 28.92, 33, 41.52, 42.12, 45.6, 48.8, 51.84, 51.96, 54.12, 55.56,
67.8, 68.44, 68.64, 68.88, 84.12, 93.12, 98.64, 105.12, 105.84, 127.92,
128.04, 173.4
)
log_x <- log(x)hist(x)
hist(log_x)
De acuerdo con estos gráficos, nos gustaría probar las siguientes hipótesis:
\[H_0: f = N(\log(50),0.25) \quad vs \quad H_1: f \neq N(\log(50),0.25)\]
Elegimos un número de grupos \(k\) tal que la probabilidad de que el log-tiempo pertenezca al i-ésimo intervalo sea mayor o igual a \(\dfrac 5{23}\approx \dfrac 14\). Es decir:
\[ p_i^0 = \mathbb P[\text{log-tiempo perteneza al }i\text{-ésimo intervalo}]\geq \dfrac 5{23}\approx \dfrac 14. \]
En este caso, podemos tomar \(k = 4\) grupos de tamaño regular. Los intervalos de estos grupos son:
Entonces solo para efectos de construir la partición
cortes <- qnorm(
p = c(0, 1 / 4, 2 / 4, 3 / 4, 1),
mean = log(50),
sd = sqrt(0.25)
)
(intervalos <- cut(c(0, 10), breaks = cortes))[1] (-Inf,3.57] (4.25, Inf]
Levels: (-Inf,3.57] (3.57,3.91] (3.91,4.25] (4.25, Inf](conteos <- cut(log_x, breaks = cortes)) [1] (-Inf,3.57] (-Inf,3.57] (-Inf,3.57] (3.57,3.91] (3.57,3.91] (3.57,3.91]
[7] (3.57,3.91] (3.91,4.25] (3.91,4.25] (3.91,4.25] (3.91,4.25] (3.91,4.25]
[13] (3.91,4.25] (3.91,4.25] (3.91,4.25] (4.25, Inf] (4.25, Inf] (4.25, Inf]
[19] (4.25, Inf] (4.25, Inf] (4.25, Inf] (4.25, Inf] (4.25, Inf]
Levels: (-Inf,3.57] (3.57,3.91] (3.91,4.25] (4.25, Inf]resumen_conteos <- summary(conteos)
names(resumen_conteos) <- c("G_1", "G_2", "G_3", "G_4")
knitr::kable(resumen_conteos, align = "c")| x | |
|---|---|
| G_1 | 3 |
| G_2 | 4 |
| G_3 | 8 |
| G_4 | 8 |
Calculando el estadístico de prueba:
\[Q = \dfrac{(3-23\cdot1/4)^2}{23\cdot 1/4} + \dfrac{(4-23\cdot1/4)^2}{23\cdot 1/4}+\dfrac{(8-23\cdot1/4)^2}{23\cdot 1/4} + \dfrac{(8-23\cdot1/4)^2}{23\cdot 1/4}= 3.609.\]
El valor-\(p\) corresponde a \(F_{\chi^2_3}(3.609) = 0.307\).
chisq.test(resumen_conteos)
Chi-squared test for given probabilities
data: resumen_conteos
X-squared = 3.6087, df = 3, p-value = 0.3069Con un nivel de significancia de 30%, no rechazamos la hipótesis de que los datos siguen una distribución normal con los parámetros especificados.
La función chisq.test si no se llama con ninguna hipótesis nula p, esta asume que \(p = 1/n\) para cada categoría. En este caso como son 4 categorías sería \(1/4\).
Una elección diferente de parámetros podría llevar a la aceptación de la hipótesis nula.
En la sección anterior, hemos realizado una prueba para verificar si nuestros datos siguen una distribución normal con media \(\log(50)\) y desviación estándar \(0.25\). Los resultados de la prueba nos indicaron que estos parámetros no son adecuados para nuestros datos.
Una pregunta natural que surge es: ¿qué parámetros serían los correctos? ¿Pertenecen los datos a una distribución normal?
Para abordar este problema, vamos a considerar una técnica que nos permita parametrizar nuestras hipótesis. En lugar de asumir que cada probabilidad \(p_i\) es fija, vamos a considerar que cada \(p_i\) es una función de algunos parámetros \(\theta = (\theta_1,\dots,\theta_s)\), es decir,
\[p_i = \pi_i(\theta),\quad i=1,\dots,k.\]
Asumimos que \(s<k-1\) y que las componentes de \(\theta\) no son redundantes, es decir, ninguna de las entradas de \(\theta\) se puede expresar como una función de las otras. Además, la suma de las probabilidades parametrizadas es igual a 1, es decir, \(\sum \pi_i(\theta) = 1\).
Nuestras hipótesis se vuelven:
La prueba estadística que usaremos es
\[Q = \sum_{i=1}^k\dfrac{[N_i-n\pi_i(\hat\theta)]^2}{n\pi_i(\hat\theta)},\]
donde \(\hat\theta\) es el estimador de máxima verosimilitud (MLE, por sus siglas en inglés) de \(\theta\), y \(N_i\) es el número de observaciones en el i-ésimo grupo.
Teorema 13.1 Bajo la hipótesis nula, cuando el número de observaciones \(n\) tiende al infinito, la distribución de \(Q\) tiende a una distribución \(\chi^2\) con \(k-1-s\) grados de libertad.
Ejemplo 13.5 Supongamos que tenemos 3 grupos y definimos un parámetro \(0<\theta<1\). El analista hace el supuesto de que las probabilidades son:
\[\begin{align*} p_1 &= \theta^2, \\ p_2 &= 2\theta(1-\theta),\\ p_3 &= (1-\theta)^2. \end{align*}\]
Notamos que la suma de estas probabilidades es 1, es decir, \(p_1+p_2+p_3=1\), ya que \[p_1+p_2+p_3 = \theta^2 + 2\theta (1-\theta +(1-\theta)^2 =[\theta+(1-\theta)]^2 = 1.\]
Aquí \(s = 1\), y el espacio de parámetros es \(\Omega = [0,1]\).
Dado que la distribución de las observaciones \((N_1,\dots,N_k)\) sigue una distribución multinomial con parámetros \(n\) y \(p_1,\dots,p_k\) bajo la hipótesis nula, podemos escribir la función de verosimilitud como
\[L(\theta|N_1,\dots,N_k) = {n \choose {N_1\cdots N_k}}(\pi_1(\theta))^{N_1}\cdots(\pi_k(\theta))^{N_k}.\]
Tomando el logaritmo obtenemos la función de log-verosimilitud:
\[\ell = \log (L) \propto N_1\ln\pi_1(\theta)+\cdots+N_k\ln\pi_k(\theta).\]
Para nuestro ejemplo, la función de log-verosimilitud se convierte en:
\[\begin{align*} \ln L(\theta) & \propto N_1\ln \theta^2 + N_2 \ln 2\theta(1-\theta) + N_3\ln (1-\theta)^2\\ & = (2N_1+N_2)\ln \theta + (2N_3+N_2)\ln(1-\theta) + N_2\ln 2. \end{align*}\]
Tomando la derivada e igualando a cero encontramos el MLE de \(\theta\):
\[\begin{align*} \dfrac{\partial \ln L(\theta)}{\partial\theta} &= \dfrac{2N_1+N_2}{\theta}-\dfrac{2N_3+N_2}{1-\theta} = 0 \\ \implies \hat\theta &= \dfrac{2N_1+N_2}{2n}. \end{align*}\]
Finalmente, podemos calcular las probabilidades \(\pi_1(\hat \theta)\),\(\pi_2(\hat \theta)\),\(\pi_3(\hat \theta)\) y el estadístico de prueba \(Q\).
Para continuar con el ejemplo supongamos que tenemos 100 observaciones, y que los datos se dividen en los siguientes grupos:
Por lo tanto, las probabilidades son:
\[\begin{align*} \pi_1(\hat\theta) &= \hat\theta^2 = \left(\dfrac{2\cdot 20+50}{2\cdot 100}\right)^2 = 0.36,\\ \pi_2(\hat\theta) &= 2\hat\theta(1-\hat\theta) = 2\left(\dfrac{2\cdot 20+50}{2\cdot 100}\right)\left(1-\dfrac{2\cdot 20+50}{2\cdot 100}\right) = 0.48,\\ \pi_3(\hat\theta) &= (1-\hat\theta)^2 = \left(1-\dfrac{2\cdot 20+50}{2\cdot 100}\right)^2 = 0.16. \end{align*}\]
Entonces la hipotesis nula del problema es \[\begin{align*} H_{0}: & p_{1}=0.36, p_{2}=0.48, p_{3}=0.16 \\ H_{1}: &p_{1} \neq 0.36, p_{2} \neq 0.48, p_{3} \neq 0.16. \end{align*}\]
Ahora si el estadístico de prueba es \[\begin{multline*} Q = \dfrac{(20-100\cdot 0.36)^2}{100\cdot 0.36} + \dfrac{(50-100\cdot 0.48)^2}{100\cdot 0.48} \\ + \dfrac{(30-100\cdot 0.16)^2}{100\cdot 0.16} = 19.44. \end{multline*}\]
Esto se debería calcular con una distribución \(\chi^2\) con 1 grado de libertad. $^2_{1,0.95} = $3.8414588.
Rechazamos la hipótesis nula si \(Q>\chi^2_{1,0.95}\).
(Q <- (20 - 100 * 0.36)^2 / (100 * 0.36) +
(50 - 100 * 0.48)^2 / (100 * 0.48) +
(30 - 100 * 0.16)^2 / (100 * 0.16))[1] 19.44444Q > qchisq(p = 0.95, df = 1)[1] TRUEEn este caso rechazamos la hipótesis nula.
El valor-\(p\) asociado a este estadístico es \(\mathbb{P}(\chi^2_{1} > 19.44)\)
pchisq(q = Q, df = 1, lower.tail = FALSE)[1] 1.035687e-05Ejemplo 13.6 (Partes de automóvil) Retomemos el ejemplo de las partes de un automóvil, y recordemos que tenemos una muestra de datos \(X_1,\dots,X_n\) que provienen de alguna distribución \(f\). Queremos probar la hipótesis de que la distribución es normal con media \(\mu\) y varianza \(\sigma^2\), ambas desconocidas, es decir, \(H_0: f = N(\mu,\sigma^2)\).
Para hacer esto, podemos construir las funciones de probabilidad parametrizadas \(\pi\), tratando de ajustar los cuantiles que habíamos definido antes con los valores teóricos de \(\mu\) y \(\sigma\). Entonces, para cada partición \(i\) con límites \(a_i\) y \(b_i\), la función de probabilidad se convierte en:
\[\begin{align*} \pi_i(\mu,\sigma^2) &= \int_{a_i}^{b_i}(2\pi\sigma^2)^{-1/2}\exp\left(-\dfrac 1{2\sigma^2}(x-\mu)^2\right)dx \\ &=\Phi\left(\dfrac{b_i-\mu}{\sigma}\right)-\Phi\left(\dfrac{a_i-\mu}{\sigma}\right), \end{align*}\]
donde \(\Phi\) es la función de distribución acumulada de una distribución normal estándar.
En el ejemplo anterior, consideramos los siguientes 4 intervalos para la partición de los datos:
\[(-\infty,3.575],(3.575,3.912],(3.912,4.249], (4.249,+\infty).\]
Ahora recuerde que en el fondo estamos buscando maximizar la log-verosimilitud
\[\ln L(\mu,\sigma^2) = N_1\ln \pi_1(\mu,\sigma^2)+\cdots+N_4\ln\pi_4(\mu,\sigma^2)\]
A continuación, se muestra cómo podemos calcular la verosimilitud y optimizarla numéricamente en R:
# Definición de los cortes de los cuantiles
cortes <- qnorm(
p = c(0, 1 / 4, 2 / 4, 3 / 4, 1),
mean = log(50),
sd = sqrt(0.25)
)
# Definición de la función de log-verosimilitud
log_versomilitud <- function(par, cortes, log_x) {
G <- length(cortes)
mu <- par[1]
sigma <- par[2]
pi <- numeric()
for (k in 1:(G - 1)) {
pi[k] <- pnorm(q = cortes[k + 1], mean = mu, sd = sigma) -
pnorm(q = cortes[k], mean = mu, sd = sigma)
}
conteos <- cut(log_x, breaks = cortes)
conteos <- summary(conteos)
l <- -sum(conteos * log(pi))
return(l)
}
# Optimización numérica de la función de log-verosimilitud
sol <- optim(
par = c(0, 1),
fn = log_versomilitud,
cortes = cortes,
log_x = log_x
)
# Imprime los parámetros optimizados
sol$par[1] 4.0926455 0.4331326Con esta información podemos calcular las probabilidades \(\pi_1(\hat \mu,\hat \sigma^2)\),\(\pi_2(\hat \mu,\hat \sigma^2)\),\(\pi_3(\hat \mu,\hat \sigma^2)\) y \(\pi_4(\hat \mu,\hat \sigma^2)\), y el estadístico de prueba \(Q\).
# Probabilidades
G <- length(cortes)
pi <- numeric()
for (k in 1:(G - 1)) {
pi[k] <- pnorm(q = cortes[k + 1], mean = sol$par[1], sd = sol$par[2]) -
pnorm(q = cortes[k], mean = sol$par[1], sd = sol$par[2])
}
(Q2 <- sum((resumen_conteos - 23 * pi)^2 / (23 * pi)))[1] 0.446496Ahora se compara con una \(\chi^2\) con \(k-s-1 = 23 - 2 - 1 =20\) grados de libertad. Es decir, nos preguntamos si rechazamos la hipótesis nula \(H_0\) de que los datos provienen de una distribución normal con media \(\mu\) y varianza \(\sigma^2\).
Q2 > qchisq(p = 0.95, df = 20)[1] FALSEEl valor-p es:
pchisq(q = Q2, df = 20, lower.tail = FALSE)[1] 1Para otra solución, considere el siguiente teorema:
Teorema 13.2 Si \(X_1,\dots, X_n\) son variables aleatorias con distribución \(F_\theta\), donde \(\theta\) es un parámetro \(p\)-dimensional y \(\hat\theta_n\) es el estimador de máxima verosimilitud (MLE) de \(\theta\) basado en \(X_1,\dots, X_n\). Tomamos una partición de \(\mathbb R\) con \(k\) intervalos disjuntos \((I_1,\dots,I_k)\), donde \(N_i\) es la cantidad de \(X_i\) que caen en el intervalo \(I_i\) y \(\pi_i(\theta)=\mathbb P_\theta[X_i\in I_i]\). Entonces, se tiene el siguiente estadístico de prueba:
\[Q' = \sum_{i=1}^k\dfrac{[N_i-n\pi_i(\hat\theta_n)]^2}{n\pi_i(\hat\theta_n)}\]
Bajo las condiciones de regularidad del MLE, si \(n\to\infty\), la función de distribución acumulada de \(Q'\) asumiendo \(H_0\) verdadero está entre \(\chi^2_{k-p-1}\) y \(\chi^2_{k-1}\).
Ejemplo 13.7 Apliquemos este teorema a nuestros datos, calculando \(\hat\mu = \bar X_n =4.1506137\) y \(\hat\sigma^2 = \dfrac{s_n^2}{n} = 0.5332049\). Tomamos la partición de los datos en los siguientes 4 intervalos:
Es decir podemos calcular lo siguiente en R
G <- length(cortes)
mu <- mean(log_x)
sigma <- sd(log_x)
pi <- numeric()
for (k in 1:(G - 1)) {
pi[k] <- pnorm(q = cortes[k + 1], mean = mu, sd = sigma) -
pnorm(q = cortes[k], mean = mu, sd = sigma)
}
pi[1] 0.1400819 0.1871877 0.2461242 0.4266062chisq.test(resumen_conteos, p = pi)Warning in chisq.test(resumen_conteos, p = pi): Chi-squared approximation may
be incorrect
Chi-squared test for given probabilities
data: resumen_conteos
X-squared = 1.3381, df = 3, p-value = 0.7201Entonces
\[Q' = \dfrac{(3-23\cdot 0.14)^2}{23\cdot 0.14} + \dfrac{(4-23\cdot 0.187)^2}{23\cdot 0.187} + \dfrac{(8-23\cdot 0.246)^2}{23\cdot 0.246} +\dfrac{(8-23\cdot 0.4266)^2}{23\cdot 0.4266} = 1.3381.\]
\(\text{valor-}p_1 = F_{\chi^2_{4-2-1}}(1.3381) = 0.7526307\).
\(\text{valor-}p_2 = F_{\chi^2_{4-1}}(1.3381) = 0.2798937\).
Rechazamos \(H_0\) (hipótesis de normalidad) si \(\alpha<0.2798\).
ggplot(data = data.frame(x = c(2.5, 6)), aes(x)) +
geom_histogram(
data = data.frame(x = log_x),
aes(x, y = after_stat(density)),
color = "white",
bins = 8
) +
stat_function(
fun = dnorm,
args = list(mean = log(50), sd = sqrt(0.25)),
aes(color = "Hipótesis Manual"),
size = 2
) +
stat_function(
fun = dnorm,
args = list(mean = sol$par[1], sd = sol$par[2]),
aes(color = "Hipótesis con Optimización"),
size = 2
) +
stat_function(
fun = dnorm,
args = list(mean = mean(log_x), sd = sd(log_x)),
aes(color = "Hipótesis con MLE"),
size = 2
) +
cowplot::theme_cowplot()
Ejemplo 13.8 Recordemos el conjunto de datos históricos sobre el número de muertes causadas por patadas de caballo en el ejército prusiano. El conjunto de datos es el siguiente:
| \(\text{Conteos}\) | \(0\) | \(1\) | \(2\) | \(3\) | \(\ge 4\) | \(\text{Total}\) |
|---|---|---|---|---|---|---|
| \(\text{Núm. de obs.}\) | \(144\) | \(91\) | \(32\) | \(11\) | \(2\) | \(280\) |
Nos preguntamos si estos datos siguen una distribución de Poisson. Podemos visualizar los datos en un gráfico de barras utilizando el siguiente código R:
df <- data.frame(
conteos = c(0, 1, 2, 3, 4),
observaciones = c(144, 91, 32, 11, 2)
)
ggplot(df, aes(x = conteos, y = observaciones)) +
geom_col() +
labs(x = "Conteos", y = "Núm. de observaciones") +
cowplot::theme_cowplot()
Formulamos la hipótesis nula como sigue: \(H_0: f = \text{Poisson}(\theta), \theta>0\).
Podemos calcular el estimador de máxima verosimilitud (MLE) de \(\theta\) como la media ponderada del número de muertes por su frecuencia, es decir,
\[\dfrac{0\cdot 144+1\cdot91+2\cdot32+3\cdot 11+2\cdot4}{280} = \dfrac{196}{280} = 0.7\]
A continuación, calculamos las probabilidades teóricas de cada número de muertes bajo la hipótesis de Poisson con \(\hat\theta = 0.7\):
Ahora, calculamos el estadístico de la prueba de bondad de ajuste \(\chi^2\) como
\[\begin{align*} Q' & = \dfrac{(144-280\cdot0.4966)^2}{280\cdot0.4966}+\dfrac{(91-280\cdot0.3476)^2}{280\cdot0.3476}+\dfrac{(32-280\cdot0.1217)^2}{280\cdot0.1217}\\ & +\dfrac{(11-280\cdot0.0283)^2}{280\cdot0.0283} +\dfrac{(2-280\cdot0.0058)^2}{280\cdot0.0058} = 1.979. \end{align*}\]
Finalmente, calculamos los valores-p para nuestro estadístico \(Q'\) bajo las dos posibles distribuciones chi-cuadrado:
Interpretamos nuestros resultados: si tomamos un nivel de significancia del 5%, no podemos rechazar la hipótesis nula de que los datos siguen una distribución de Poisson.
Finalmente, podemos visualizar nuestra prueba de bondad de ajuste en un gráfico de barras que compara las observaciones reales con las esperadas bajo la distribución de Poisson con \(\hat\theta = 0.7\).
total <- sum(df$observaciones)
ggplot(data = data.frame(x = c(0, 4)), aes(x)) +
geom_col(data = df, aes(x = conteos, y = observaciones / total)) +
stat_function(
fun = dpois,
args = list(lambda = 0.7),
aes(color = "Hipótesis con MLE"),
size = 2,
geom = "col"
) +
labs(x = "Conteos", y = "Probabilidad") +
cowplot::theme_cowplot()