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
)
df <- as.data.frame(x)16 Pruebas de Kolmogorov-Smirnov
En los capítulos anteriores, implementamos la prueba \(\chi^2\) para determinar si un conjunto de datos se ajustaba a una distribución continua. Sin embargo, existe una prueba más adecuada para este propósito: la prueba Kolmogorov-Smirnov.
Supongamos que tenemos una serie de observaciones \(X_1,\dots,X_n\), que siguen una distribución continua \(F\). También asumimos que los valores observados \(x_1,\dots,x_n\) son todos diferentes.
Definición 16.1 Dados los valores observados \(x_1,\dots,x_n\) de la muestra aleatoria \(X_1,\dots, X_n\), para cada \(x\) definimos \(F_n(x)\) como la proporción de valores observados en la muestra que son menores o iguales a \(x\). Es decir, si hay \(k\) valores observados menores o iguales a \(x\), \[ F_n(x) = \frac{k}{n}. \]
La función \(F_n(x)\) se conoce como la función de distribución de la muestra (empírica). \(F_n\) es una distribución escalonada con saltos de magnitud \(\dfrac 1n\) en los valores \(x_1,\dots,x_n\). Se puede expresar como \[F_n(x) = \begin{cases} 0 & \text{si } x<X_{(1)}\\ \displaystyle \dfrac 1n \sum_{i=1}^n 1_{\{x_i\leq x\}} & \text{si } X_{(1)}\leq x<X_{(n)}\\ 1 & \text{si} X_{(n)}\geq x \end{cases} \]
Dado que \(\{x_i\}_{i=1}^n\) son independientes, \(\{1_{{x_i\leq x}}\}_{i=1}^n\) también son independientes. Por lo tanto, por la ley de los grandes números
\[ F_n(x) = \dfrac 1n \sum_{i=1}^n 1_{\{x_i\leq x\}} \xrightarrow[n\to\infty]{\mathbb P} \mathbb E[1_{X_i\leq x}] = F(x) \]
por lo que \(F_n(x)\) es consistente.
Ejemplo 16.1 Suponga que se tiene el siguiente conjunto de datos:
Para este ejemplo tenemos que \(F_n(x)\) se ve de la forma:
ggplot() +
stat_ecdf(
data = df,
mapping = aes(x, color = "F_n(x)"),
linewidth = 2
) +
cowplot::theme_cowplot()
En los capítulos anteriores, determinamos que los parámetros \(\mu =3.912\) y \(\sigma ^{2} = 0.25\) se ajustaban bien a los log-valores de nuestros datos. Ahora, podemos verificar este ajuste con la siguiente gráfica:
ggplot() +
stat_ecdf(
data = df,
mapping = aes(log(x), color = "F_n(x)"),
linewidth = 2
) +
stat_function(
data = data.frame(x = c(3, 5)),
fun = pnorm,
aes(color = "F(x)"),
linewidth = 2,
args = list(mean = 3.912, sd = sqrt(0.25))
) +
cowplot::theme_cowplot()
La prueba de Kolmogorov-Smirnov se basa en el siguiente teorema, conocido como el Lema de Glivenko-Cantelli.
Lema 16.1 (Glivenko-Cantelli) Sea \(F_n(x)\) la distribución empírica de una muestra \(X_1, \dots , X_n\) provenientes de la distribución \(F\). Defina \[ D_n = \sup_{-\infty<x<\infty}|F_n(x)-F(x)| \] Entonces \[D_n\xrightarrow[]{\mathbb P} 0\]
Esto implica que si una distribución empírica \(F_n(x)\) se obtiene realmente de la distribución teórica \(F(x)\), entonces la diferencia entre estas dos convergerá en probabilidad a 0 cuando \(n \to \infty\).
La demostración de este teorema está más allá del alcance de este curso, por lo que no se incluirá en esta sección.
16.1 Prueba de Kolmogorov-Smirnov para una muestra
La prueba de Kolmogorov-Smirnov se utiliza para responder a la pregunta de si una distribución empírica de datos, \(F\), es igual a una distribución teórica, \(F^\star\). Se establecen las siguientes hipótesis:
- Hipótesis nula, \(H_0\): \(F = F^{\star}\)
- Hipótesis alternativa, \(H_1\): \(F\neq F^{\star}\)
El estadístico de prueba se define como \[D_n^{\star} = \sup_{-\infty<x<\infty}|F_n(x)-F^(x)|\].
Importante
Si la hipótesis nula es verdadera, entonces el valor de \(D_n^{\star}\) no depende de \(F^{\star}\).
En términos prácticos, si la distribución empírica \(F_n(x)\) está cerca de la distribución teórica \(F^{\star}(x)\), entonces el valor de \(D_n^{\star}\) será cercano a cero. Por lo tanto, rechazamos la hipótesis nula, \(H_0\), si \(n^{\frac{1}{2}}D_n^{\star}\geq c\) para algún valor crítico \(c\).
Este valor crítico \(c\) se obtiene a partir de la distribución de Kolmogorov-Smirnov. De acuerdo con el teorema de Kolmogorov-Smirnov (1930), si la hipótesis nula es verdadera, para \(t>0\),
Teorema 16.1 (Kolmogorov-Smirnov (1930)) Si \(H_0\) es cierto, para \(t>0\), \[\lim_{n\to \infty} \mathbb P(n^{1/2}D_n^{\star}\leq t) = 1-2\sum_{i=1}^\infty (-1)^{i-1}e^{-2i^2t^2} = H(t).\]
Rechazamos \(H_0\) si \(n^{1/2}D_n^{\star}\geq c\), \(n\) grande. Para un nivel de significancia \(\alpha_0\), \(c = H^{-1}(1-\alpha_0)\), donde \(H\) denota el valor de la parte derecha de la ecuación anterior.
La función \(H(t)\) es algo complicada de estimar, y sus cuantiles lo son aún más. Estos normalmente son definidos a través de métodos núméricos que están fuera del alcance del este curso. La siguiente tabla muestra el conjunto de valores estimados para cada \(t\)
La función de distribución se ve de la siguiente forma:
H <- function(t) {
i <- 1:1000
1 - 2 * sum((-1)^(i - 1) * exp(-2 * i^2 * t^2))
}
t <- seq(0.1, 2, length.out = 1000)
data.frame(t, H = Vectorize(H)(t)) %>%
ggplot(aes(t, H)) +
geom_line(linewidth = 2) +
cowplot::theme_cowplot()
Los valores más comunes de cuantiles para las pruebas son
| \(\alpha\) | \(H^{-1}(1-\alpha)\) |
|---|---|
| 0.01 | 1.63 |
| 0.05 | 1.36 |
| 0.1 | 1.22 |
Solo como ilustración, la funcion ks.test tiene programada la estimación de \(\mathbb P (D_n \leq t)\). En este gráfico se presentan algunos ejemplos con \(n\) igual a 10, 50 y 100.
pkolmogorov1x <- function(x, n) {
## Probability function for the one-sided one-sample Kolmogorov
## statistics, based on the formula of Birnbaum & Tingey (1951).
if (x <= 0) {
return(0)
}
if (x >= 1) {
return(1)
}
j <- seq.int(from = 0, to = floor(n * (1 - x)))
1 - x * sum(exp(lchoose(n, j)
+ (n - j) * log(1 - x - j / n)
+ (j - 1) * log(x + j / n)))
}
t <- seq(0.001, 0.999, length.out = 1000)
df <- rbind(
data.frame(t, ks = Vectorize(pkolmogorov1x)(t, 10), n = 10),
data.frame(t, ks = Vectorize(pkolmogorov1x)(t, 50), n = 50),
data.frame(t, ks = Vectorize(pkolmogorov1x)(t, 100), n = 100)
)ggplot(df, aes(x = t, y = ks, color = as.factor(n))) +
geom_line(linewidth = 2) +
labs(color = "n") +
cowplot::theme_cowplot()
Ejemplo 16.2 En el caso de las partes mecánicas quisiéramos saber si los log-valores siguen o no una distribución normal de parámetros \(N (\hat{\mu } , s^{2})\).
Para calcular manualmente el valor de \(D_{n}^{\star }\) se deben usar los siguientes dos estimadores:
\[\begin{align*} D_n^{+} &= \sup_{1\leq i\leq n}\left\{ \dfrac{i}{n}-F(X_{(i)})\right\} \\ D_n^{-} &= \sup_{1\leq i\leq n}\left\{ F(X_{(i)})-\dfrac{i-1}{n}\right\} \end{align*}\]
Es decir, se debe tomar en cuenta si la distribución empírica está por encima o por debajo de la distribución teórica.
Una vez con esto se puede calcular el valor de \(D_{n}^{\star }\) como
\[ D_{n}^{\star }=\max \left\{ D_{n}^{+},D_{n}^{-}\right\} \]
df <- data.frame(logx = sort(log(x)))
F_t <- pnorm(sort(log(x)), mean = mean(log(x)), sd = sd(log(x)))
F_n <- seq(1, length(F_t)) / length(F_t)
df <- cbind(df, F_n, F_t)
df <- df %>%
mutate(
D_n_plus = abs(F_n - F_t),
D_n_minus = abs(F_t - lag(F_n)),
D_n_star = pmax(D_n_plus, D_n_minus)
)
knitr::kable(df,
digits = 3,
col.names = c(
"$\\log(x)$", "$F_n(\\log(x))$",
"$F(\\log(x))$", "$D_n^+$", "$D_n^-$", "$D_n^\\star$"
)
)| \(\log(x)\) | \(F_n(\log(x))\) | \(F(\log(x))\) | \(D_n^+\) | \(D_n^-\) | \(D_n^\star\) |
|---|---|---|---|---|---|
| 2.884 | 0.043 | 0.009 | 0.035 | NA | NA |
| 3.365 | 0.087 | 0.070 | 0.017 | 0.027 | 0.027 |
| 3.497 | 0.130 | 0.110 | 0.020 | 0.023 | 0.023 |
| 3.726 | 0.174 | 0.213 | 0.039 | 0.083 | 0.083 |
| 3.741 | 0.217 | 0.221 | 0.004 | 0.047 | 0.047 |
| 3.820 | 0.261 | 0.268 | 0.007 | 0.050 | 0.050 |
| 3.888 | 0.304 | 0.311 | 0.007 | 0.050 | 0.050 |
| 3.948 | 0.348 | 0.352 | 0.004 | 0.048 | 0.048 |
| 3.950 | 0.391 | 0.354 | 0.038 | 0.006 | 0.038 |
| 3.991 | 0.435 | 0.382 | 0.052 | 0.009 | 0.052 |
| 4.017 | 0.478 | 0.401 | 0.077 | 0.033 | 0.077 |
| 4.217 | 0.522 | 0.549 | 0.027 | 0.071 | 0.071 |
| 4.226 | 0.565 | 0.556 | 0.009 | 0.034 | 0.034 |
| 4.229 | 0.609 | 0.558 | 0.050 | 0.007 | 0.050 |
| 4.232 | 0.652 | 0.561 | 0.091 | 0.048 | 0.091 |
| 4.432 | 0.696 | 0.701 | 0.006 | 0.049 | 0.049 |
| 4.534 | 0.739 | 0.764 | 0.025 | 0.068 | 0.068 |
| 4.591 | 0.783 | 0.796 | 0.013 | 0.057 | 0.057 |
| 4.655 | 0.826 | 0.828 | 0.002 | 0.045 | 0.045 |
| 4.662 | 0.870 | 0.831 | 0.038 | 0.005 | 0.038 |
| 4.851 | 0.913 | 0.906 | 0.007 | 0.036 | 0.036 |
| 4.852 | 0.957 | 0.906 | 0.051 | 0.007 | 0.051 |
| 5.156 | 1.000 | 0.970 | 0.030 | 0.014 | 0.030 |
ggplot(df) +
geom_step(aes(logx, F_n, color = "F_n"), linewidth = 2) +
geom_line(aes(logx, F_t, color = "F"), linewidth = 2) +
geom_segment(
aes(x = logx, xend = logx, y = F_n, yend = F_t),
color = "purple",
linewidth = 2
) +
geom_segment(
aes(x = logx, xend = logx, y = F_t, yend = lag(F_n)),
color = "gold",
linewidth = 2
) +
cowplot::theme_cowplot()
(D_n_star <- max(df$D_n_star, na.rm = TRUE))[1] 0.09124598El valor se estimaría usando la función vista anteriormente
n <- length(x)
sqrt(n) * D_n_star[1] 0.4376003A un nivel de significancia del 5% se rechaza la hipótesis nula si
sqrt(n) * D_n_star > 1.36[1] FALSE1 - H(sqrt(n) * D_n_star)[1] 0.9908784library(KSgeneral)
cont_ks_c_cdf(D_n_star, n)[1] 0.9814592Otra forma de hacer es usando la función ks.test
ks.test(
x = log(x),
y = "pnorm",
mean = mean(log(x)),
sd = sd(log(x))
)
Exact one-sample Kolmogorov-Smirnov test
data: log(x)
D = 0.091246, p-value = 0.9815
alternative hypothesis: two-sidedcont_ks_test(log(x), "pnorm", mean(log(x)), sd(log(x)))
One-sample Kolmogorov-Smirnov test
data: log(x)
D = 0.091246, p-value = 0.9815
alternative hypothesis: two-sidedObservación. Los valores de los valores-p son diferentes ya que la función ks.test usa varias aproximaciones para calcularlos.
Note que esta localización es muy importante ya que si se quisiera comparar con una distribución \(N(0,1)\) el resultado es diferente.
ks.test(
x = log(x),
y = "pnorm",
mean = 0,
sd = 1
)
Exact one-sample Kolmogorov-Smirnov test
data: log(x)
D = 0.99803, p-value = 4.441e-16
alternative hypothesis: two-sided16.2 Prueba de 2 muestras
Suponga que se tiene \(X_1,\dots,X_m\sim N(\mu_1,\sigma^2)\) y \(Y_1,\dots,Y_n\sim N(\mu_2,\sigma^2)\) y se desea saber si ambas muestras tienen la misma distribución.
Una opción es probar que \[ H_0:\mu_1 = \mu_2 \text{ vs } H_1: \mu_{1} \neq \mu_{2} \]
Uno de los supuestos fuertes para este tipo de pruebas es la normalidad.
Otra opción es decir que
\[ H_0: F_1 = F_2 \text{ vs } H_1: F_{1} \neq F_{2} \]
donde \(F_1\) es la distribución de \(X\) y \(F_2\) la de \(Y\).
Igual en este caso estamos probando dos distribuciones normales, pero ¿Es posible del todo quitar el supuesto de normalidad?
Es decir, para \(X_1,\dots,X_m\sim F\) y \(Y_1,\dots,Y_m\sim G\) continuas, sin valores en común, probar \(H_0: F(x) = G(x)\), \(x \in \mathbb R\).
Considere \[ D_{mn}^{\star} = \sup_{-\infty<x<\infty}|F_m(x)-G_n(x)| \]
Se tiene por el teorema de Glivenko-Cantelli que \(D_{mn}^{\star}\xrightarrow[]{\mathbb P} 0\), \(m,n\to\infty\) cuando \(H_0\) es verdadera
Para el caso de dos muestras se puede probar que si \(H(t)\) es la distribución límite en el caso de una muestra y \(t>0\), entonces se cumple que
\[\lim_{m,n\to \infty} \mathbb P \left( \left( \dfrac{mn}{m+n}\right)^{\frac 12} D_{mn}^{\star}\leq t\right) = H(t)\]
En este caso se rechaza la hipótesis nula si \(\left(\dfrac{mn}{m+n}\right)^{\frac 12}D_{mn}^{\star} \geq H^{-1}(1-\alpha_0)\).
Ejemplo 16.3 Suponga que se tienen dos grupos de personas a las cuales a unas se les dio un tratamiento para la presión arterial y al otro se le dio un placebo.
Al se midieron las diferencias entre las presiones arteriales al cabo de 12 semanas de tratamiento.
Los resultados fueron estos
Medicina <- c(7, -4, 18, 17, -3, -5, 1, 10, 11, -2)
Placebo <- c(-1, 12, -1, -3, 3, -5, 5, 2, -11, -1, -3)La pregunta es si ambos conjuntos de datos vienen de la misma distribución.
Podemos repetir el mismo procedimiento anterior, pero modificado para dos muestras. Primero necesitamos estimar el \(D_{nm}^{\star}\) usando estos dos estimadores:
\[\begin{align*} D_{nm}^{+} &= \sup_{1\leq i\leq n}\left\{ F_{m}(X_{(i)})-G_{n}(X_{(i)})\right\} \\ D_{nm}^{-} &= \sup_{1\leq i\leq n}\left\{ G_{n}(X_{(i)})-F_{m}(X_{(i-1)})\right\} \end{align*}\]
Estimamos \(D_{nm}^{\star}\) como
\[ D_{nm}^{\star }=\max \left\{ D_{nm}^{+},D_{nm}^{-}\right\} \]
Medicina <- sort(Medicina)
Placebo <- sort(Placebo)
df <- data.frame(Medicina,
F_Medicina = ecdf(Medicina)(Medicina),
F_Placebo = ecdf(Placebo)(Medicina)
)
df <- df %>%
mutate(
D_mn_plus = abs(F_Medicina - F_Placebo),
D_mn_minus = abs(F_Placebo - lag(F_Medicina)),
D_mn_star = pmax(D_mn_plus, D_mn_minus)
)
knitr::kable(df,
digits = 3,
col.names = c(
"Medicina",
"$F_\\text{Medicina}$",
"$F_\\text{Placebo}$", "$D_{mn}^+$", "$D_{mn}^-$", "$D_{mn}^\\star$"
)
)| Medicina | \(F_\text{Medicina}\) | \(F_\text{Placebo}\) | \(D_{mn}^+\) | \(D_{mn}^-\) | \(D_{mn}^\star\) |
|---|---|---|---|---|---|
| -5 | 0.1 | 0.182 | 0.082 | NA | NA |
| -4 | 0.2 | 0.182 | 0.018 | 0.082 | 0.082 |
| -3 | 0.3 | 0.364 | 0.064 | 0.164 | 0.164 |
| -2 | 0.4 | 0.364 | 0.036 | 0.064 | 0.064 |
| 1 | 0.5 | 0.636 | 0.136 | 0.236 | 0.236 |
| 7 | 0.6 | 0.909 | 0.309 | 0.409 | 0.409 |
| 10 | 0.7 | 0.909 | 0.209 | 0.309 | 0.309 |
| 11 | 0.8 | 0.909 | 0.109 | 0.209 | 0.209 |
| 17 | 0.9 | 1.000 | 0.100 | 0.200 | 0.200 |
| 18 | 1.0 | 1.000 | 0.000 | 0.100 | 0.100 |
ggplot(df) +
geom_step(aes(Medicina, F_Medicina, color = "F_Medicina"), linewidth = 2) +
geom_step(aes(Medicina, F_Placebo, color = "F_Placebo"), linewidth = 2) +
geom_segment(
aes(x = Medicina, xend = Medicina, y = F_Medicina, yend = F_Placebo),
color = "purple",
linewidth = 2
) +
geom_segment(
aes(x = Medicina, xend = Medicina, y = F_Placebo, yend = lag(F_Medicina)),
color = "gold",
linewidth = 2
) +
cowplot::theme_cowplot()
(D_mn_star <- max(df$D_mn_star, na.rm = TRUE))[1] 0.4090909n <- length(Medicina)
m <- length(Placebo)
sqrt(n * m / (n + m)) * D_mn_star[1] 0.93628171 - H(sqrt(n * m / (n + m)) * D_mn_star)[1] 0.3446214La prueba se realiza con el comando ks.test:
ks.test(Medicina, Placebo)
Exact two-sample Kolmogorov-Smirnov test
data: Medicina and Placebo
D = 0.40909, p-value = 0.2422
alternative hypothesis: two-sidedEn este caso rechazamos la hipótesis nula de que ambas distribuciones son iguales con un nivel de \(\alpha \geq 0.346\).
df <- rbind(
data.frame(x = Medicina, Tratamiento = "Medicina"),
data.frame(x = Placebo, Tratamiento = "Placebo")
)
ggplot(df) +
stat_ecdf(aes(x, color = Tratamiento), linewidth = 2) +
cowplot::theme_cowplot()