En la parte anterior consideramos un número fijo de resultados numéricos de experimentos aleatorios, por ejemplo, cuando \(X\) el resultado de una tirada de dado. En este caso, un modelo de probabilidad para \(X\) asigna una probabilidad dada a cada posible resultado, por ejemplo
\[P(X=1) = 1/6\]
e igualmente \(P(X=2)=\cdots = P(X=6) = 1/6\). En muchos casos, la cantidad \(X\) que nos interesa puede tomar valores numéricos arbitrarios, y en este esquema no está claro cómo asignaríamos probabilidades.
Supongamos que giramos una ruleta con una flecha indicadora, y el resultado del experimento es el ángulo en grados final de la flecha. ¿Cómo podríamos poner, por ejemplo \(P(X= 92.7)\)? ¿Qué pasa si podemos medir con resultados con varios decimales de exactitud?
En estos casos, en lugar de considerar eventos de la forma \(X=a\), podemos considerar eventos resultantes de la forma \(a\leq X \leq b\), es decir, buscamos asignar probabilidades a eventos de la forma
\[P(X \in [a,b]).\]
Por ejemplo, si \(Y\) es la estatura adulta de una persona que acaba de nacer, podríamos preguntarnos cómo asignar probabilidades a eventos como
\[P(Y\in [150,170]) = P(150\leq Y \leq 170),\]
y quizá también otros como
\[P(Y < 180).\]
Los modelos más simple para una medición continua \(X\) son los modelos uniformes.
Para nuestra ruleta, por ejemplo, \(X\) puede tomar valores en el intervalo \([0, 360)\). Si la ruleta es justa, entonces la probabilidad de que la flecha caiga en cualquier sector \([a,b]\) debe ser igual. Una manera de lograr esto usar como probabilidad la proporción de la longitud de \([a,b]\) con respecto al total de \([0, 360)\):
\[P(X\in [a,b]) = \frac{b-a}{360}.\]
Supongamos que una variable aleatoria puede tomar valores en el intervalo \([L,U]\). La variable aleatoria es uniforme en \([L,U]\) cuando
\[P(X \in [a,b]) = \frac{b-a}{L-U}\]
Ahora supongamos que nuestra ruleta no está del todo balanceada. Por ejemplo, podría ser estuviera colgada en una pared, y al girar la flecha es un poco más probable que la flecha apunte hacia el piso en lugar de hacia el cielo.
En este caso, si la dirección hacia arriba es 90 grados y hacia abajo es 270 grados, quisiéramos por ejemplo que
\[P(260 < X <280) > P(80 < X < 100)\]
Y nótese que debe ser posible asignar probabilidades a cualquier sector de la ruleta con el nuevo modelo que propongamos. ¿Cómo podríamos modificar nuestra asignación de probabilidades?
Una de las maneras más fáciles es pensando que nuestra probabilidad se obtiene integrando una funcion constante:
\[P(X\in [a,b]) = \int_a^b \frac{1}{360} \,dx = \frac{b-a}{360}\]
De forma que si \(f(x)= 1/360\) para valores \(0 \leq x < 360\), nuestra probabilidad se escribe como la integral
\[P(X\in [a,b]) = \int_a^b f(x) \,dx \]
Para generalizar la idea es la siguiente:
\[P(X\in [0, 360]) = \int_0^{360} f(x)\,dx = 1\]
Podríamos utilizar por ejemplo:
f_dens <- function(x){
x_rad <- 2 * pi * x / 360
(1/360) + 0.0002 * cos(x_rad - 3 * pi / 2)
}
graf_1_tbl <- tibble(angulo = seq(0, 360, 1), tipo = "uniforme",
f = 1 / 360)
graf_2_tbl <- tibble(angulo = seq(0, 360, 1), tipo = "colgada") %>%
mutate(f = f_dens(angulo))
graf_tbl <- bind_rows(graf_1_tbl, graf_2_tbl)
ggplot(graf_tbl, aes(x = angulo, y = f, colour = tipo)) +
geom_line() +
ylim(c(0, 0.003)) + facet_wrap(~tipo, nrow = 1)
El cálculo se hace ahora con área bajo la curva. Para calcular la probabilidad
\[P(X\in [50, 130]),\]
integramos la función \(f\) correspondiente, que corresponde a calcular área bajo la curva:
ggplot(graf_tbl, aes(x = angulo, y = f, colour = tipo)) +
geom_line() +
ylim(c(0, 0.003)) + facet_wrap(~tipo, nrow = 1) +
geom_area(aes(x = ifelse(angulo > 50 & angulo < 130, angulo, 0)),
fill="salmon", alpha = 0.5)
Y ahora vemos que para la versión perturbada, más de la probabilidad se concentra alrededor de 270 grados que alrededor de 90. Por las propiedades de la integral, todas las propiedades usuales de probabilidad se cumplen.
Cuando trabajamos con mediciones de tipo continuo, es más conveniente definir asignaciones de probabilidad utilizando funciones de densidad de probabilidad:
Una función \(f(x)\) no negativa cuya integral es igual a 1 es una función de densidad de probabilidad. Las probabilidades asociadas se calculan integrando:
\[P(X\in [a,b]) = \int_a^b f(x)\,dx\]
En este caso decimos que \(f\) es la función de densidad de probabilidad asociada a la variable aleatoria \(X\). A este tipo de variables aleatorias les llamamos continuas.Supongamos que tenemos una variable aleatoria que tiene mediana 2, y puede tomar valores entre 0 y 4. Podríamos definir una densidad como sigue: Si \(x\) está entre 0 y 2, entonces
\[f(x) = \frac{x}{4}\]
y si \(x\) está entre 2 y 4, entonces
\[f(x) = 1 - \frac{x}{4}\]
dens_triangular <- function(x){
(x > 0) * (x < 4) * ifelse(x < 2, x/4, 1 - x/4)
}
triangular_tbl <- tibble(x = seq(-1, 5, 0.001)) %>%
mutate(f = dens_triangular(x))
ggplot(triangular_tbl, aes(x = x, y = f)) +
geom_line()
Supongamos que una variable \(X\) tiene mediana 2 y rango de 0 a 4, con densidad triangular. ¿Cuál es la probabilidad \(P(X>1)\)?
Solución: Por reglas usuales de probabilidad, \(P(X>1) = P(1<X<2) + P(X\geq2)\). Tenemos que \(P(X\geq 2) = 0.5\). Ahora usamos la fórmula de la densidad triangular para obtener
\[P(1<X<2) = \int_{1}^{2} f(x)\,dx = \int_1^2 \frac{x}{4}\,dx = \left [\frac{x^2}{8}\right ]_1^2 = 1/2 - 1/8 = 3/8 = 0.375\]
de modo que
\[P(X>1) = 0.375 +0.500 = 0.875\]
En general, podemos dar una fórmula para una densidad triangular en el intervalo \([A,B]\) con mediana en \((A + B)/2\). ¿Cómo sería la fórmula?
Antes vimos la definición de cuantiles para datos numéricos. Podemos definirlos también para variables aleatorias numéricas:
Sea \(p\in (0,1)\). El cuantil-\(p\) de la variable \(X\) con función de densidad \(f(x)\) es el valor \(x(p)\) tal que
\[\int_{-\infty}^{x(p)} f(x)\,dx = p\]
Observación: nótese que usamos como límite inferior \(-\infty\) para indicar que integramos \(f\) sobre toda la densidad que esté a la izquierda de \(x(p)\).
Supongamos que \(X\) tiene la densidad triangular mostrada arriba. Calcula el cuartil inferior y superior (es decir, los cuantiles 0.25 y 0.75). Para el cuartil superior, por ejemplo, buscamos al \(x(0.75)\) de la siguiente gráfica:
source("R/triangular.R")
ggplot(triangular_tbl, aes(x = x, y = f)) +
geom_line() +
geom_area(aes(x = ifelse(x > 0 & x < qtri(0.75, 0, 4), x, 0)),
fill="salmon", alpha = 0.5) +
ylim(c(0, 0.7)) +
annotate("text", x = qtri(0.75, 0, 4), y = 0.03, label = "x(0.75)") +
annotate("point", x = qtri(0.75, 0, 4), y = 0.0) 
Comenzaremos por el cuartil inferior. Buscamos una \(x(0.25)\) tal que
\[\int_0^{x(0.25)} f(x)\,dx = 0.25\]
Sabemos que \(x(0.25)< 2\), pues la integral hasta 2 es 0.5, así que
\[\int_0^{x(0.25)} f(x)\,dx = \int_0^{x(0.25)} x/4 \,dx = \left [ x^2/8\right]_0^{x(0.25)} = (x(0.25))^2/8\]
Si queremos que este valor sea igual a 0.25, entonces despejando obtenemos
\[x(0.25) = \sqrt{0.25(8)} = \sqrt{2}\approx 1.4142\]
Ahora podríamos calcular la otra integral, pero por simetría podemos concluir que
\[x(0.75) = 2 + (2 - 1.4142) \approx 2.5858\]
y concluimos que los cuartiles inferiores y superiores son aproximadamente 1.41 y 2.59
Calcula la mediana, y los percentiles 0.10 y 0.90 de una variable uniforme en \([0, 10]\).
Los cuantiles que vimos en la parte de descriptivos para datos numéricos se llaman usualmente cuantiles empíricos. Estos cuantiles podemos compararlos con cuantiles teóricos para ver qué tan similares son, y si el modelo teórico describe adecuadamente los datos.
Simularemos 500 datos uniformes en \([0, 10]\):
x_sim_u <- runif(500, 0, 10)Podríamos calcular algunos cuantiles empíricos:
quantile(x_sim_u, c(0.10, 0.50, 0.90)) 10% 50% 90%
0.9108561 4.9737151 9.1764574 Por el ejercicio anterior sabemos cuáles son los cuantiles teóricos correspondientes a una uniforme en \([0,10]\). Podemos calcularlos también como sigue:
qunif(c(0.10, 0.5, 0.90), 0, 10)[1] 1 5 9Y vemos que son muy similares los cuantiles empíricos y teóricos. Una mejor manera de considerar esta similitud es graficando todos los cuantiles empíricos y comparándolos con los teóricos:
ggplot(tibble(x = x_sim_u), aes(sample = x)) +
geom_abline(colour = "red") +
geom_qq(distribution = stats::qunif, dparams = list(min = 0, max = 10)) +
xlab("Cuantiles Teóricos U(0,10)") + ylab("Cuantiles de datos")
Y vemos que la forma de las dos distribuciones es muy similar: los cuantiles empíricos son muy similares a los teóricos. Existen algunas fluctuaciones debidas al muestreo aleatorio.
Repetimos para la distribución triangular. Los cuantiles que calculamos arriba son:
qtri(c(0.25, 0.75), a = 0, b = 4)[1] 1.414214 2.585786x_sim_tri <- rtri(500, 0, 4)
ggplot(tibble(x = x_sim_tri), aes(sample = x)) +
geom_abline(colour = "red") +
geom_qq(distribution = qtri, dparams = list(a = 0, b = 4)) +
xlab("Cuantiles Teóricos triangular(0,4)") + ylab("Cuantiles de datos")
Nótese que otra vez, los cuantiles teóricos se alinean bien con los teóricos.
La distribución normal es una que aparece naturalmente en al teoría de probabilidad.
Consideremos que tiramos 40 dados justos de 6 lados, y consideramos su promedio \(\bar{X}\) como resultado de nuestro experimento aleatorio. ¿Cómo se ve la distribución de probabilidades de esta variable \(\bar{X}\)?
Comenzamos haciendo simulacion:
simular_bolsa <- function(num_dados = 40){
dados <- sample(1:6, num_dados, replace = TRUE)
media <- mean(dados)
media
}
simular_bolsa()[1] 3.475Veamos cómo se ven los resultados si hacemos este experimento un número grande de veces:
set.seed(23)
sims_dados <- map_dbl(1:10000, ~ simular_bolsa())
ggplot(tibble(resultado = sims_dados),
aes(x = resultado)) +
geom_histogram(binwidth = 0.1)
Y notamos una forma de “campana” interesante. Esto se explica porque típicamente tendremos tiros bajos y altos, de modo que muchos resultados de este experimento se concentran alrededor de la media de un dado (3.5). Además, existen fluctuaciones aleatorios, y a veces tenemos un poco más de tiros altos o de tiros bajos, de forma que existe dispersión alrededor de 3.5.
Sin embargo, estas desviaciones de 3.5 no pueden ser muy grandes: por ejemplo, para tener un promedio de 1, todas las tiradas de los 40 dados tendrían que dar 1, y eso es muy poco probable. Igualmente, para que el promedio fuera cercano a 6, la gran mayoría de los 40 dados deberían de dar 6, lo cual otra vez es muy poco probable. Esto explica al menos la forma general de la forma de las colas derecha e izquierda de esta distribución.
Los dados podrían ser diferentes (por ejemplo, un poco cargados a 1 o 6, o más cargados a valores centrales), y los argumentos de arriba también se cumplirían. Lo que es sorprendente es que, independientemente de cómo sean las particularidades de los dados, la forma analítica de esta distribución que acabamos de mostrar es la misma.
Esta forma está descrita por la densidad normal estándar:
\[f(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}\]
cuya gráfica presentamos a continuación:
tibble(x = seq(-3.5, 3.5, 0.01)) %>%
mutate(f = (1/sqrt(2*pi)) * exp(-x^2 / 2)) %>%
ggplot(aes(x = x, y = f)) + geom_line()
A una variable \(Z\) que tiene esta densidad le llamamos una variable con distribución normal estándar. Comparemos cuantiles en nuestro ejemplo:
ggplot(tibble(resultado = sims_dados),
aes(sample = resultado)) +
geom_qq(distribution = stats::qnorm) +
xlab("Normal estándar teórica") +
ylab("Promedio de 40 dados")
Y notamos que los cuantiles no corresponden, pero el espaciamiento entre los cuantiles de los datos y los teóricos de la normal estándar es el mismo. Quiere decir que estas dos distribuciones tienen la misma forma, aunque estén escaladas y centradas en distintos valores.
Probemos con promedios de 20 observaciones triangulares en \((0,1)\) por ejemplo. El resultado es el mismo:
set.seed(23)
sims_tri <- map_dbl(1:10000, ~ mean(rtri(20, 0 ,1)))
ggplot(tibble(resultado = sims_tri),
aes(x = resultado)) +
geom_histogram(binwidth = 0.01)
ggplot(tibble(resultado = sims_tri),
aes(sample = resultado)) +
geom_qq(distribution = stats::qnorm) +
xlab("Normal estándar teórica") +
ylab("Promedio de 20 triangulares (0,1)")
Otra vez, la forma general es la misma, aún cuando los datos están centrados y escalados de manera distinta.
Como expicamos, la densidad normal estándar está dada por
\[f(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}},\]
cuya gráfica es como sigue:
normal_std_graf <- tibble(x = seq(-3.5, 3.5, 0.01)) %>%
mutate(f = dnorm(x, 0, 1))
ggplot(normal_std_graf, aes(x = x, y = f)) +
geom_line()
Las probabilidades asociadas a una normal estándar se calculan integrando esta curva (que tiene que hacerse de forma numérica). Por ejemplo, para calcular
\[P(Z < 1.5),\]
podemos usar
pnorm(1.5)[1] 0.9331928Que es el área bajo la curva mostrada abajo:
normal_std_graf <- tibble(x = seq(-3.5, 3.5, 0.005)) %>%
mutate(f = dnorm(x, 0, 1))
ggplot(normal_std_graf, aes(x = x, y = f)) +
geom_line() +
geom_area(aes(x = ifelse(x > -3.5 & x < 1.5, x, 0)),
fill="salmon", alpha = 0.5) +
ylim(c(0,0.4)) +
scale_x_continuous(breaks = seq(-3.5, 3.5, 0.5))
Esta es la forma de la densidad estándar. Podemos centrar esta campana en otro valor \(\mu\) y aumentar la dispersión por un factor \(\sigma\). Si \(Z\) es una variable normal estándar, la variable
\[X = \mu + \sigma Z\]
es una variable normal con parámetros \((\mu, \sigma)\), o de manera más compacta, decimos que \(X\) es \(N(\mu, \sigma)\). La distribución normal estándar es \(N(0,1)\).
Por ejemplo, si escogemos \(\mu=5\) y \(\sigma = 0.5\), obtenemos:
normal_graf <- tibble(x = seq(3, 7, 0.01)) %>%
mutate(f = dnorm(x, 5, 0.5))
ggplot(normal_graf, aes(x = x, y = f)) +
geom_line()
Podemos mostrar juntas estas dos distribuciones:
densidad_normal <- tibble(x = seq(3, 7, 0.1)) %>%
mutate(densidad = dnorm(x, 5, 0.5))
densidad_normal_estandar <- tibble(x = seq(-3, 3, 0.1)) %>%
mutate(densidad = dnorm(x))
g_2 <- ggplot(densidad_normal_estandar, aes(x = x, y = densidad)) + geom_line()
g_3 <- g_2 + xlim(c(-3, 7)) + ylim(c(0, 1))
g_1 <- ggplot(densidad_normal, aes(x = x, y = densidad)) + geom_line() + xlim(c(-3, 7)) + ylim(c(0, 1))
g_3 + g_1
Se puede demostrar que:
Distribución normal
Con un poco de cálculo podemos ver qué tan fuertemente se concentra la densidad alrededor de la media para una distribución normal. La regla es la siguiente:
grafica_concentracion <- function(mu, sigma, z){
x <- seq(mu - 3.1 * sigma, mu + 3.1 * sigma, 0.01)
valor <- dnorm(x, mu, sigma)
datos <- tibble(x = x, `f(x)`=valor)
texto <- round(100*(pnorm(z) - pnorm(-z)), 1)
texto <- paste0(texto, "%")
ggplot(datos, aes(x = x, y = `f(x)`)) +
geom_area(data = filter(datos, x < mu + z*sigma, x > mu - z*sigma),
fill = "salmon") +
geom_line() +
annotate("text", x = mu, y = 0.1, label = texto) +
scale_x_continuous(breaks = mu + sigma*seq(-3, 3, 1)) +
theme_minimal() + ylab("")
}
g_68 <- grafica_concentracion(10, 2, 1)
g_95 <- grafica_concentracion(10, 2, 2)
g_997 <- grafica_concentracion(10, 2, 3)
paneles <- g_68 + g_95 + g_997
paneles + plot_annotation(title = "Concentración alrededor de la media (normal)")
Nota: esto aplica para cualquier densidad normal, independientemente de los parámetros.
Obsérvese que esto nos da una interpretación natural de la desviación estándar de una distribución normal en términos de percentiles de los datos, y la manera usual con la que entendemos la desviación estándar de la distribución normal.
Una de las razones por las que el modelo normal es tan importante es el siguiente resultado:
Si \(X_1,X_2,\ldots, X_n\) son variables aleatorias independientes con media \(\mu\) y desviación estándar \(\sigma\) con una densidad \(f(x)\):
Esto nos permite, por ejemplo, considerar nuestro primera técnica de estimación por intervalos:
\[\bar{x} = \frac{1}{n}(x_1+\cdots + x_n) = \frac{1}{n}\sum_i x_i,\]
\[P(\mu - 2\sigma \leq \bar{x} \leq\mu + 2\sigma)\approx 0.95\]
Despejando \(\mu\) obtenemos
\[P(\bar{x} - 2\sigma \leq \mu \leq\bar{x} + 2\sigma)\approx 0.95\]
Finalmente, no conocemos \(\sigma\), pero la estimamos con
\[\hat{\sigma}^2 = \frac{1}{n}((x_1 - \bar{x})^2 + \dots + (x_n - \bar{x})^2) = \frac{1}{n}\sum_i(x_i - \bar{x})^2\]
Y aproximamos sustituyendo nuestra estimación:
\[P(\bar{x} - 2\hat{\sigma } \leq \mu \leq\bar{x} + 2\hat{\sigma})\approx 0.95\]
Esto nos da un intervalo (llamado el intervalo de Wald) con 95% de confianza para la media poblacional:
\[[\bar{x} - 2\hat{\sigma } , \bar{x} + 2\hat{\sigma}]\]
Notas:
Supongamos que queremos queremos modelar la estatura de cantantes tenores y sopranos en una población, y tomamos una muestra aleatoria. Nuestro primer supuesto será que los datos de ambas poblaciones se distribuyen de manera aproximadamente normal, con medias distintas pero la misma dispersión.
Nuestro modelo es entonces:
Si tomamos muestras aleatorias independientes de estos dos grupos, y denotaremos las estaturas registradas como:
\[y_1, y_2, \ldots, y_n\]
y la tesitura del cantante como
\[s_1,s_2,\ldots, s_n\]
donde \(s_i=1\) indica soprano y \(s_i=0\) indica tenor. Escribimos nuestra verosimilitud, que en este caso consiste en multipicar las funciones de densidad normal correspondientes. Por ejemplo, si tenemos \(n=3\), \(y_1 = 160, s_1=0\) y \(y_2 = 150, y_3 = 162, s_2=1, s_3 = 1,\) y \(f(x;\mu,\sigma)\) es la función de densidad normal, podríamos escribir la verosimilitud como
\[L(\mu_1, \mu_2, \sigma; y) = f\left(160; \mu_t,\sigma\right)f\left(150; \mu_s,\sigma\right)f\left(162; \mu_s, \sigma\right)\]
En general, podemos escribir la log verosimilitud como (aunque con algo de trabajo puede simplificarse a una versión numéricamente más conveniente):
crear_log_verosim <- function(y, s){
log_verosim <- function(pars){
mu_s <- pars[1]
mu_t <- pars[2]
sigma <- pars[3]
res <- sum(dnorm(y, mean = mu_s * s + mu_t * (1 - s), sigma, log = TRUE))
res
}
log_verosim
}Nuestra muestra es:
set.seed(32413)
cantantes <- lattice::singer |>
mutate(estatura_cm = round(2.54 * height)) |>
filter(str_detect(voice.part, "Tenor|Soprano")) |>
group_by(voice.part) |>
slice_sample(n = 20) |>
ungroup() |>
mutate(tesitura = ifelse(str_detect(voice.part, "Tenor"), "tenor", "soprano")) |>
select(tesitura, estatura_cm)
cantantes |> ggplot(aes(x = tesitura, y = estatura_cm)) + geom_boxplot()
Estimamos las medias y desviación estándar con máxima verosimilitud:
y <- cantantes$estatura_cm
s <- (cantantes$tesitura == "soprano") |> as.numeric()
log_vero <- crear_log_verosim(y, s)
res_mv <- optim(par = c(150, 150, 30), log_vero,
control = list(fnscale = -1, maxit = 2000))
est_mv <- res_mv$par
est_mv |> round(2)[1] 162.33 176.32 6.08Y construimos intervalos de confianza utilizando el bootstrap paramétrico. Sin embargo, deberíamos checar nuestro modelo antes de proseguir. Si el ajuste del modelo es malo, es posible que nuestros intervalos no tengan una cobertura cercana a la nominal.
Una estrategia general para checar ajuste de modelos es generar datos a partir del modelo estimado. Así, estos datos generados por los modelos pueden ser contrastados con los datos reales para encontrar desviaciones importantes que puedan afectar la inferencia.
La simulación debe incluir la incertidumbre que tenemos acerca de los parámetros del modelo, de forma que podemos usar el bootstrap paramétrico para incorporar esta incertidumbre.
En nuestro caso, tenemos las siguentes funciones
generar_muestra <- function(s, mu_s, mu_t, sigma){
y <- rnorm(length(s), mu_s * s + mu_t * (1 - s), sigma)
tibble(s = s, estatura_cm = y)
}
generar_muestra_boot <- function(s, est_mv){
# muestrear del modelo ajustado
muestra_boot <- generar_muestra(s, est_mv[1], est_mv[2], est_mv[3])
# reestimar parámetros
y <- muestra_boot$estatura_cm
s <- muestra_boot$s
log_boot_vero <- crear_log_verosim(y, s)
res_boot_mv <- optim(par = c(150, 150, 30), log_boot_vero,
control = list(fnscale = -1, maxit = 1000))
est_mv <- res_boot_mv$par
generar_muestra(s, est_mv[1], est_mv[2], est_mv[3])
}Y las utilizamos para generar muestras. Nótese que redondeamos para usar la misma precisión de los datos originales:
muestras_tbl <- map_df(1:19, function(rep) {
generar_muestra_boot(s, est_mv) |>
mutate(estatura_cm = 2.54 * round(estatura_cm/2.54)) |>
mutate(rep = rep)
})muestras_tbl <- bind_rows(muestras_tbl, cantantes |>
mutate(rep = 20, s = ifelse(tesitura == "soprano",1, 0)))
ggplot(muestras_tbl, aes(sample = estatura_cm, colour = factor(s))) +
geom_qq() + facet_wrap(~ rep)
Notamos que no hay desjustes grandes de nuestro modelo. En primer lugar, el supuesto de normalidad parece consistente con los datos a grandes rasgos. Nuestra hipótesis de igualdad de varianzas parece sostenerse razonablemente bien, aunque es posible que la dispersión en los tenores sea un poco mayor que en las sopranos.
Prueba con un modelo que no suponga la misma varianza para los dos grupos, y discute ventajas y desventajas.
Finalmente, construimos intervalos del 95% para cada uno de los parámetros:
parametros_boot <- function(s, est_mv){
# muestrear del modelo ajustado
muestra_boot <- generar_muestra(s, est_mv[1], est_mv[2], est_mv[3])
# reestimar parámetros
y <- muestra_boot$estatura_cm
s <- muestra_boot$s
log_boot_vero <- crear_log_verosim(y, s)
res_boot_mv <- optim(par = c(150, 150, 30), log_boot_vero,
control = list(fnscale = -1, maxit = 1000))
est_boot_mv <- res_boot_mv$par
tibble(mu_s = est_boot_mv[1], mu_t = est_boot_mv[2], sigma = est_boot_mv[3])
}
remuestreo_tbl <- map_df(1:1000, ~ parametros_boot(s, est_mv))remuestreo_tbl |> pivot_longer(mu_s:sigma, names_to = "parametro", values_to = "valor") |>
group_by(parametro) |>
summarise(inf_5 = quantile(valor, 0.05), inf_95 = quantile(valor, 0.95)) |>
mutate(across(where(is.numeric), ~ round(.x, 1)))# A tibble: 3 × 3
parametro inf_5 inf_95
<chr> <dbl> <dbl>
1 mu_s 161. 164.
2 mu_t 175. 178.
3 sigma 5.2 6.7En las secciones anteriores, definimos modelos de probabilidad para datos observados (de tipo discreto y continuo).
Estos modelos dependen generalmente de parámetros desconocidos poblacionales (como probabilidad de ser seropositivo, o medias y dispersión de la población de cantantes).
Las cantidades de interés ahora son esos parámetros del modelo de probabilidad, y queremos hacer inferencia sobre ellos.
Introdujimos el método de máxima verosimilitud como una manera general de encontrar estimadores de esos parámetros dada una muestra, y vimos cómo construir intervalos de confianza usando el bootstrap paramétrico.
Mostramos un ejemplo simple de cómo checar los supuestos del modelo para saber si podemos confiar en la inferencia.
Nota: en lugar de utilizar el bootstrap paramétrico, es posible utilizar resultados asintóticos (como el teorema central del límite) para hacer inferencia.