Al finalizar este capítulo, el estudiante será capaz de:
Cuando trabajamos con un modelo estadístico con un parámetro desconocido \(\theta\), es común en el enfoque bayesiano considerar a \(\theta\) como una variable aleatoria. La distribución que describe nuestras creencias o conocimientos previos sobre \(\theta\) antes de observar cualquier dato se llama densidad previa y se denota por \(\pi(\theta)\).
Observación. La función de distribución previa \(\pi\) tiene las siguientes características.
Hagamos algunos ejemplos para entender como funciona la densidad previa.
Ejemplo 4.1 (Componentes eléctricos) Consideremos el tiempo de vida de ciertos componentes que sigue una distribución exponencial \(X_1,\dots, X_n \sim \text{Exp}(\theta)\).
Si suponemos que \(\theta\) es una variable aleatoria que sigue una distribución Gamma con parámetros \(\alpha = 1\) y \(\beta = 2\), entonces la densidad previa de \(\theta\) es:
\[\begin{equation*} \pi(\theta) = \dfrac{1}{\Gamma(\alpha)}\beta^\alpha\theta^{\alpha-1}e^{-\beta\theta} = 2e^{-2\theta}, \quad \theta > 0 \end{equation*}\]
Ejemplo 4.2 (Lanzamiento de una moneda) Supongamos que estamos interesados en averiguar la probabilidad de obtener cara al lanzar una moneda, en un experimento donde hay dos tipos de monedas: las monedas justas y las monedas alteradas. Sin embargo, hay incertidumbre sobre qué moneda se va a utilizar en la siguiente lanzada, podría ser una moneda justa o una moneda alterada que siempre cae en cara, pero se sabe que su densidad previa sigue una distribución Gamma, por lo que antes de observar cualquier lanzamiento, modelamos nuestras creencias sobre \(\theta\), asumiendo la densidad previa calculada en el ejercicio anterior, de la siguiente manera:
Esto significa que si tuviéramos 100 monedas, esperaríamos que 20 de ellas fueran alterada y 80 fueran justas.
Ejemplo 4.3 (Componentes eléctricos (cont.)) Supongamos que estamos interesados en el tiempo de vida de un componente eléctrico. Sabemos que este tiempo se puede modelar con una distribución exponencial con parámetro \(\theta\). Además, basándonos en la experiencia de un experto, creemos que \(\theta\) tiene una distribución previa Gamma. El experto nos proporciona la siguiente información:
\[ \mathbb{E}[\theta] = 0.0002, \quad \sqrt{\text{Var}(\theta)} = 0.0001 \]
Dada la distribución Gamma, sus momentos están definidos en función de sus parámetros \(\alpha\) (forma) y \(\beta\) (tasa). Específicamente, para una variable aleatoria $ $ que sigue una distribución Gamma, sus momentos están dados por:
\[\begin{align*} \mathbb{E}[\theta] &= \frac{\alpha}{\beta} = 0.0002 \\ \text{Var}(\theta) &= \frac{\alpha}{\beta^2} = 0.0001^2 \end{align*}\]
Por lo tanto despejando \(\alpha\) y \(\beta\) tenemos que:
\[\begin{align*} \beta &= \frac{\alpha}{0.0002} \\ \alpha &= 0.0001^2 \beta^2 \end{align*}\]
Sustituyendo en ambas ecuaciones, obtenemos que:
\[\begin{align*} \alpha &= 4 \\ \beta &= 20000. \end{align*}\]
La distribución previa \(\theta \sim \Gamma(4, 20000)\) luce así:
alpha_previa <- 4
beta_previa <- 20000
ggplot(data = data.frame(x = c(0, 0.001)), aes(x)) +
stat_function(
fun = dgamma,
args = list(shape = alpha_previa, rate = beta_previa),
linewidth = 2
) +
xlab("θ") +
ylab("Densidad") +
ggtitle("Distribución previa Γ(4, 20000)") +
theme_minimal()
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma as gamma_dist
alpha_previa_py = 4
beta_previa_py = 20000
x_plot = np.linspace(0, 0.001, 500)
y_plot = gamma_dist.pdf(x_plot, a=alpha_previa_py, scale=1 / beta_previa_py)
fig, ax = plt.subplots()
ax.plot(x_plot, y_plot, linewidth=2)
ax.set_xlabel("θ")
ax.set_ylabel("Densidad")
ax.set_title("Distribución previa Γ(4, 20000)")
plt.tight_layout()
plt.show()
Ejemplo 4.4 (Continuación de Ejemplo 4.3) Considerando el contexto del componente eléctrico discutido anteriormente, si tomamos una muestra \(X = (X_1, \dots, X_n)\) donde cada \(X_i\) sigue una distribución exponencial con parámetro \(\theta\) y \(X_i > 0\), la función de densidad conjunta condicionada en \(\theta\) se describe como:
\[\begin{align*} f_n(X|\theta) &= \prod_{i=1}^n \theta e^{-\theta X_i} \\ &= \theta^n e^{-\theta\sum_{i=1}^n X_i}. \end{align*}\]
Una vez que observamos los datos, actualizamos nuestras creencias sobre \(\theta\) usando el teorema de Bayes. La distribución resultante se llama densidad posterior y se denota por \(\pi(\theta|X)\).
Teorema 4.1 Dada una muestra aleatoria compuesta por las observaciones \(X_1, X_2, \dots, X_n\), podemos definir la densidad posterior de \(\theta\) como:
\[ \pi(\theta|X) = \frac{f(X_1|\theta) f(X_2|\theta) \dots f(X_n|\theta) \pi(\theta)}{g_n(X)} \]
Aquí, \(g_n(X)\) actúa como una constante de normalización que asegura que la densidad posterior integre a uno.
Demostración. Para derivar la expresión anterior, aplicamos directamente el Teorema de Bayes. La fórmula de Bayes nos dice que:
\[ \pi(\theta|X) = \frac{\pi(X|\theta) \pi(\theta)}{\pi(X)} \]
Expandiendo la densidad conjunta y reorganizando términos, obtenemos:
\[ \pi(\theta|X) = \frac{f_n(X|\theta) \pi(\theta)}{\int \pi(\theta,X) \, d\theta} \]
Esto se simplifica a:
\[ \pi(\theta|X) = \frac{f(X_1|\theta) f(X_2|\theta) \dots \times f(X_n|\theta) \times \pi(\theta)}{g_n(X)} \]
Esta expresión nos da la densidad posterior de \(\theta\) dada la muestra observada.
Ejemplo 4.5 (Componentes eléctricos (cont.)) Supongamos que tenemos una muestra \(X = (X_1, X_2, \dots, X_n)\) donde cada observación \(X_i\) sigue una distribución exponencial. La función de verosimilitud para esta muestra está dada por:
\[ f_n(X|\theta) = \theta^n e^{-\theta y}. \]
Aquí, \(y = \sum_{i=1}^n X_i\) representa la suma total de las observaciones y actúa como un estadístico para \(\theta\).
Para calcular la densidad posterior de \(\theta\), necesitamos considerar tanto la función de verosimilitud como la densidad previa. Desglosemos esto en dos partes: el numerador y el denominador de la densidad posterior \(\pi(\theta\vert X)\):
\[\begin{equation*} f_n(X|\theta)\pi(\theta) = \underbrace{\theta^n e^{-\theta y}}_{f_n(X|\theta)} \cdot \underbrace{\dfrac{20000^4}{3!}\theta^3e^{-20000\cdot\theta}}_{\pi(\theta)} = \dfrac{20000^4}{3!}\theta^{n+3}e^{-(20000+y)\theta} \end{equation*}\]
\[\begin{equation*} g_n(x) = \int_{0}^{+\infty}\theta^{n+3}e^{-(20000+y)\theta}\;d\theta = \dfrac{\Gamma(n+4)}{(20000+y)^{n+4}}. \end{equation*}\]
Con estos componentes, la densidad posterior de \(\theta\) es:
\[ \pi(\theta|X) = \frac{\theta^{n+3}e^{-(20000+y)\theta}}{\Gamma(n+4)}(20000+y)^{n+4} \]
Esta densidad posterior sigue una distribución \(\Gamma(n+4, 20000+y)\).
Supongamos que observamos 5 tiempos de fallo: 2911, 3403, 3237, 3509, y 3118. Entonces, \(y = 16478\) y \(n = 5\). Con estos datos, la distribución posterior de \(\theta\) es \(\Gamma(9, 36178)\).
La distribución posterior es sensible al tamaño de la muestra. Es decir, que una muestra grande inicial, implica un efecto menor de la previa.
x <- c(2911, 3403, 3237, 3509, 3118)
alpha_previa <- 4
beta_previa <- 20000
(y <- sum(x))[1] 16178(n <- length(x))[1] 5(alpha_posterior <- n + alpha_previa)[1] 9(beta_posterior <- beta_previa + y)[1] 36178ggplot(data = data.frame(x = c(0, 0.001)), aes(x)) +
stat_function(
fun = dgamma,
args = list(shape = alpha_previa, rate = beta_previa),
aes(color = "Previa"),
linewidth = 2
) +
stat_function(
fun = dgamma,
args = list(shape = alpha_posterior, rate = beta_posterior),
aes(color = "Posterior"),
linewidth = 2
) +
guides(color = guide_legend(title = "Densidad")) +
geom_vline(
xintercept = alpha_posterior / beta_posterior,
linetype = "dashed",
color = "#FC8D62",
linewidth = 1.2
) +
annotate(
"label",
x = alpha_posterior / beta_posterior,
y = 6000,
label = paste0("Media = ", round(alpha_posterior / beta_posterior, 6)),
hjust = -0.1,
fill = "#FC8D62"
) +
geom_vline(
xintercept = qgamma(0.5, shape = alpha_posterior, rate = beta_posterior),
linetype = "dashed",
color = "#8DA0CB",
linewidth = 1.2
) +
annotate(
"label",
x = qgamma(0.5, alpha_posterior, beta_posterior),
y = 4000,
label = paste0(
"Mediana = ",
round(qgamma(0.5, alpha_posterior, beta_posterior), 6)
),
hjust = -0.1,
fill = "#8DA0CB"
) +
geom_vline(
xintercept = (alpha_posterior - 1) / beta_posterior,
linetype = "dashed",
color = "#66C2A5",
linewidth = 1.2
) +
annotate(
"label",
x = (alpha_posterior - 1) / beta_posterior,
y = 2000,
label = paste0("Moda = ", round((alpha_posterior - 1) / beta_posterior, 6)),
hjust = -0.1,
fill = "#66C2A5"
) +
ylab("") +
cowplot::theme_cowplot()
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma as gamma_dist
x_data = np.array([2911, 3403, 3237, 3509, 3118])
alpha_previa_py = 4
beta_previa_py = 20000
y_sum = np.sum(x_data)
n = len(x_data)
alpha_posterior_py = n + alpha_previa_py
beta_posterior_py = beta_previa_py + y_sum
print(f"Posterior: Γ({alpha_posterior_py}, {beta_posterior_py})")Posterior: Γ(9, 36178)print(f"Media posterior = {alpha_posterior_py / beta_posterior_py:.6f}")Media posterior = 0.000249print(
f"Mediana posterior = {gamma_dist.ppf(0.5, a=alpha_posterior_py, scale=1 / beta_posterior_py):.6f}"
)Mediana posterior = 0.000240print(f"Moda posterior = {(alpha_posterior_py - 1) / beta_posterior_py:.6f}")Moda posterior = 0.000221x_plot = np.linspace(0, 0.001, 500)
fig, ax = plt.subplots()
ax.plot(
x_plot,
gamma_dist.pdf(x_plot, a=alpha_previa_py, scale=1 / beta_previa_py),
label="Previa Γ(4, 20000)",
linewidth=2,
)
ax.plot(
x_plot,
gamma_dist.pdf(x_plot, a=alpha_posterior_py, scale=1 / beta_posterior_py),
label=f"Posterior Γ({alpha_posterior_py}, {beta_posterior_py})",
linewidth=2,
)
ax.axvline(
alpha_posterior_py / beta_posterior_py,
color="#FC8D62",
linestyle="--",
label="Media",
)
ax.axvline(
gamma_dist.ppf(0.5, a=alpha_posterior_py, scale=1 / beta_posterior_py),
color="#8DA0CB",
linestyle="--",
label="Mediana",
)
ax.set_xlabel("θ")
ax.set_ylabel("Densidad")
ax.legend()
ax.set_title("Previa vs Posterior — Componentes eléctricos")
plt.tight_layout()
plt.show()
Observación. Le decimos Hiperparámetros a los parámetros de la previa o posterior.
De ahora en adelante vamos a entender un modelo como la unión de los siguientes componentes:
Estos dos últimos resumen el comportamiento de los datos.
Para identificar y trabajar con este modelo, lo desglosamos en las siguientes partes:
Información previa (\(\pi(\theta)\)): Esta es la información adicional o basada en la experiencia que tenemos sobre el modelo.
Verosimilitud: Los datos representan la información observada. La función de densidad \(f\) refina y mejora la información previa. La función \(f_n(X|\theta)\) se conoce como verosimilitud o función de verosimilitud y se define como:
\[ f_n(X\|\theta) = \prod\_{i=1}^n f(X_i\|\theta). \]
Densidad posterior: Es una combinación de la información previa y los datos observados, proporcionando una versión más informada de la distribución del parámetro. Se calcula como:
\[ \pi(\theta\|X) \propto f_n(x\|\theta)\pi(\theta). \]
Predicciones con datos nuevos: Si asumimos que los datos se obtienen secuencialmente, podemos calcular la distribución posterior de la siguiente manera:
\[\begin{align*} \pi(\theta|X_1) & \propto \pi(\theta) f(X_1|\theta)\\ \pi(\theta|X_1,X_2) &\propto \pi(\theta) f(X_1,X_2|\theta) \\ & = \pi(\theta) f(X_1|\theta) f(X_2|\theta) \text{ (por independencia condicional)}\\ & = \pi(\theta|X_1)f(X_2|\theta)\\ \vdots & \\ \pi(\theta|X_1,\dots,X_n) & \propto f(X_n|\theta)\pi(\theta|X_1,\dots, X_{n-1}) \end{align*}\]
Bajo la suposición de independencia condicional, no hay diferencia en la posterior si los datos se obtienen secuencialmente. Esta distribución nos ayuda a calcular probabilidades sobre datos nuevos. Por lo tanto:
\[\begin{align*} g_n(X) & = \int\_{\Omega} f(X_n\|\theta) \pi(\theta\|X_1,\dots, X_{n-1})\;d\theta\\ & = P(X_n\|X_1,\dots,X_{n-1}) \text{ (Predicción para }X_n) \end{align*}\]
Observación. En el caso de una función de verosimilitud, el argumento es \(\theta\) y los datos en \(X\) son fijos.
El proceso de inferencia bayesiana sigue siempre la misma estructura:
flowchart LR A["Conocimiento previo\n π(θ)"] --> C["Distribución posterior\n π(θ|X)"] B["Datos observados\n f_n(X|θ)"] --> C C --> D["Estimador bayesiano\n δ*(X)"]
La posterior combina lo que sabíamos antes de ver los datos (previa) con la información aportada por los datos (verosimilitud). Cuantos más datos tengamos, menor será la influencia de la previa.
Ejemplo 4.6 (Proporción de componentes defectuosos) En ejemplos anteriores, abordamos los tiempos de vida de componentes defectuosos asumiendo un valor fijo para \(\theta\). Ahora, exploraremos cómo modelar este parámetro utilizando estadística bayesiana.
Consideremos un componente que puede ser defectuoso o no. Es natural modelar su comportamiento con una variable aleatoria \(X_i\) definida como:
\[\begin{equation*} X_i = \begin{cases} 0 & \text{falló} \\ 1 & \text{no falló} \end{cases} \end{equation*}\]
Donde \(\theta\) representa la proporción de componentes defectuosos, con \(\theta \in [0,1]\). Una elección adecuada para modelar esta variable aleatoria es asumir que \(X_i \sim \mathrm{Bernoulli}(\theta)\).
Siguiendo los pasos anteriores, definamos
\[ \pi(\theta) = 1_{\{0\leq\theta\leq 1\}} \] Es importante notar que una distribución uniforme es equivalente a una distribución \(\mathrm{Beta}(1,1)\).
\[ f_n(X\|\theta) = \theta^{\sum_{i=1}^{n} X_i}(1-\theta)^{n-\sum_{i=1}^{n} X_i} = \theta^{y}(1-\theta)^{n-y} \] Donde \(y = \sum_{i=1}^{n} X_i\) es el número total de componentes que no fallaron.
\[ \pi(\theta\|X) \propto \theta^y (1-\theta)^{n-y} = \theta^{y+1-1}(1-\theta)^{n-y+1-1} \] De esta expresión, deducimos que, dada la información \(X\), \(\theta\) sigue una distribución \(\text{Beta}(y+1,n-y+1)\).
Observación. La distribución posterior anteriormente es exactamente
\[\begin{equation*} \pi(\theta|X) = \frac{\Gamma(n+1)}{\Gamma(y+1)\Gamma(n-y+1)}\theta^{y}(1-\theta)^{n-y}. \end{equation*}\]
Sin embargo el operador \(\propto\) indica que el lado izquierdo es igual al lado derecho salvo por una constante que no depende de la variable (\(\theta\)). Esto es muy útil para evitar escribir mucho en las fórmulas.
Ejemplo 4.7 (Componentes eléctricos (cont.)) En el ejemplo Ejemplo 4.5, discutimos la modelación del tiempo de vida de componentes. Ahora supongamos que se nos presenta un nuevo dato, \(X_6\). Estamos interesados en determinar la probabilidad de que este nuevo componente tenga un tiempo de vida superior a 3000, es decir, \(P(X_6>3000|X_1,X_2,X_3,X_4,X_5)\).
Para calcular esta probabilidad, necesitamos determinar la función de densidad condicional \(f(X_6|X)\). Dado que la distribución posterior de \(\theta\) está dada por:
\[\begin{align*} \pi(\theta|X) &= \frac{36178^9}{\Gamma(9)} \theta^8 e^{-36178\theta} \\ &= 2.633\times 10^{36}\theta^8 e^{-36178\theta} \end{align*}\]
Podemos expresar la función de densidad condicional como:
\[\begin{align*} f(X_6|X) &= \int_{0}^\infty f(X_6|\theta) \pi(\theta|X)\;d\theta \\ &= \int_{0}^\infty \theta e^{-\theta X_6} \pi(\theta|X)\;d\theta \\ &= 2.633\times 10^{36} \int_{0}^\infty \theta e^{-\theta X_6} \theta^8 e^{-36178\theta}\;d\theta \\ &= 2.633\times 10^{36} \int_{0}^\infty \theta^9 e^{-\theta (X_6+36178)} \;d\theta \\ &= 2.633\times 10^{36} \dfrac{9!}{(X_6+36178)^{10}} \\ &= \dfrac{9.55 \times 10^{41}}{(X_6+36178)^{10}} \end{align*}\]
Con esta función en mano, podemos calcular la probabilidad deseada:
\[ P(X_6>3000\mid X) = \int_{3000}^{\infty} \dfrac{9.55\times10^{41}}{(X_6+36178)^{10}}\; dX_6 = 0.4882 \]
Observación. La vida media del componente se define como el tiempo en el cual la probabilidad de que el componente falle es del 50%. Matemáticamente, esto se expresa como:
\[ \dfrac{1}{2} = P(X_6>u|X) \]
En este caso, lo que tendríamos es que \[\begin{align*} \dfrac{1}{2} &= \int_{u}^{\infty} \dfrac{9.55\times10^{41}}{(X_6+36178)^{10}}\; dX_6 \\ &= \dfrac{ 9.55\times10^{41}}{-9(X_6+36178)^{9}}\Big|_{u}^{\infty} \\ &= \dfrac{ 9.55\times10^{41}}{9(u+36178)^{9}} \\ \end{align*}\]
Despejando tenemos que \[\begin{align*} \dfrac{1}{2} &= \dfrac{ 9.55\times10^{41}}{9(u+36178)^{9}} \\ (u + 36178)^9 &= \dfrac{ 9.55\times10^{41}}{9\cdot 0.5} \\ u + 36178 &= \sqrt[9]{\dfrac{ 9.55\times10^{41}}{9\cdot 0.5}} \\ u &= \sqrt[9]{\dfrac{ 9.55\times10^{41}}{9\cdot 0.5}} - 36178 \\ u &= 2894 \end{align*}\]
Los siguientes bloques ilustran la actualización bayesiana con los 5 datos originales y, luego, con un 6to dato adicional (\(X_6 = 3000\)).
x <- c(2911, 3403, 3237, 3509, 3118)
alpha_previa <- 4
beta_previa <- 20000
(y <- sum(x))[1] 16178(n <- length(x))[1] 5(alpha_posterior <- n + alpha_previa)[1] 9(beta_posterior <- beta_previa + y)[1] 36178# Actualización con X6 = 3000
(alpha_previa2 <- alpha_posterior)[1] 9(beta_previa2 <- beta_posterior)[1] 36178(alpha_posterior2 <- alpha_previa2 + 1)[1] 10(beta_posterior2 <- beta_previa2 + 3000)[1] 39178ggplot(data = data.frame(x = c(0, 0.001)), aes(x)) +
stat_function(
fun = dgamma,
args = list(shape = 4, rate = 20000),
aes(color = "Previa original"),
linewidth = 2
) +
stat_function(
fun = dgamma,
args = list(shape = alpha_posterior, rate = beta_posterior),
aes(color = "Posterior (n=5)"),
linewidth = 2
) +
stat_function(
fun = dgamma,
args = list(shape = alpha_posterior2, rate = beta_posterior2),
aes(color = "Posterior (n=6)"),
linewidth = 2
) +
labs(color = "Distribución") +
cowplot::theme_cowplot()
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma as gamma_dist
x_data = np.array([2911, 3403, 3237, 3509, 3118])
alpha_previa_py = 4
beta_previa_py = 20000
y_sum = np.sum(x_data)
n = len(x_data)
alpha_post_py = n + alpha_previa_py
beta_post_py = beta_previa_py + y_sum
# Actualización con X6 = 3000
alpha_post2_py = alpha_post_py + 1
beta_post2_py = beta_post_py + 3000
x_plot = np.linspace(0, 0.001, 500)
fig, ax = plt.subplots()
ax.plot(
x_plot,
gamma_dist.pdf(x_plot, a=4, scale=1 / 20000),
label="Previa original",
linewidth=2,
)
ax.plot(
x_plot,
gamma_dist.pdf(x_plot, a=alpha_post_py, scale=1 / beta_post_py),
label=f"Posterior (n=5): Γ({alpha_post_py}, {beta_post_py})",
linewidth=2,
)
ax.plot(
x_plot,
gamma_dist.pdf(x_plot, a=alpha_post2_py, scale=1 / beta_post2_py),
label=f"Posterior (n=6): Γ({alpha_post2_py}, {beta_post2_py})",
linewidth=2,
)
ax.set_xlabel("θ")
ax.set_ylabel("Densidad")
ax.legend()
ax.set_title("Actualización secuencial de la posterior")
plt.tight_layout()
plt.show()
En estadística bayesiana, las familias conjugadas desempeñan un papel fundamental al simplificar el proceso de actualización de nuestras creencias a medida que se obtienen nuevos datos. A continuación, se presenta una introducción detallada a este concepto.
Dada una muestra \(X_1, \dots, X_n\) que es independiente e idénticamente distribuida (i.i.d.) condicionalmente dado un parámetro \(\theta\) con densidad \(f(x\mid \theta)\), consideremos,\(\psi\) como la familia de posibles densidades previas sobre \(\Omega\).
Definición 4.1 Si, independientemente de los datos observados, la densidad posterior pertenece a la misma familia que \(\psi\), entonces decimos que \(\psi\) es una familia conjugada de previas.
Ejemplo 4.8
Para el Ejemplo Ejemplo 4.3, la familia Gamma es una familia conjugada para muestras exponenciales.
En el Ejemplo Ejemplo 4.6, la familia Beta es una familia conjugada para muestras según una Binomial.
Para el caso de datos que siguen una distribución Poisson, si \(X_1,\dots,X_n\sim\mathrm{Poisson}(\lambda)\),la familia Gamma es conjugada. La función de densidad de una Poisson es:
\[\begin{equation*} P(X_i = k) = e^{-\lambda}\dfrac{\lambda^k}{k!} \end{equation*}\]
La verosimilitud, que es el producto de las densidades de los datos observados, se da por: \[\begin{equation*} f_n(X|\lambda) = \prod_{i=1}^{n}e^{-\lambda}\dfrac{\lambda^{X_i}}{X_i!} = \dfrac{e^{-n\lambda}\lambda^{y}}{\prod_{i=1}^n X_i}. \end{equation*}\] Si la previa de \(\lambda\) está definida por \(\pi(\lambda)\propto\lambda^{\alpha-1}e^{-\beta\lambda}\), entonces la posterior es: \[\begin{equation*} \pi(\lambda|X) \propto \lambda^{y+\alpha-1}e^{-(\beta+n)\lambda} \implies \lambda|X \sim \Gamma(y+\alpha,\beta+n) \end{equation*}\]
En el caso de datos que siguen una distribución normal, \(X_1,\dots,X_n\sim N(\theta,\sigma^2)\), la familia normal es conjugada si la varianza \(\sigma^2\) es conocida. Si \(\theta \sim N(\mu_0,V_0^2)\) entonces la posterior para \(\theta\) data la muestra es \(\sim N(\mu_1, V_1^2)\) donde,
\[\begin{align*} \mu\_1 &= \dfrac{\sigma^2\mu_0 + nV_0^2 \bar{X}_n}{\sigma^2 + nV_0^2} = \dfrac{\sigma^2}{\sigma^2 + nV_0^2}\mu\_0 + \dfrac{nV_0^2}{\sigma^2 + nV_0^2}\bar{X}*n\\ V*1^2 &= \frac{\sigma^2 V_0^2}{\sigma^2+n V_0^2}. \end{align*}\]
Este resultado nos muestra cómo la información previa (representada por \(\mu_0\) y \(V_0^2\)) se combina con la información de la muestra (representada por \(\bar{X}_n\)) para producir una estimación posterior de \(\theta\).
En resumen, la tabla a continuación muestra algunas de las familias conjugadas más comunes en estadística bayesiana.
| Distribución | Previa conjugada |
|---|---|
| \(\text{Bernoulli}(p)\) | \(p \sim Beta(\alpha,\beta)\) |
| \(\text{Binomial}(n,p)\) | \(p \sim Beta(\alpha,\beta)\) |
| \(\text{Binomial Negativa}(n,p)\) | \(p \sim Beta(\alpha,\beta)\) |
| \(\text{Poisson}(\lambda)\) | \(\lambda \sim Gamma(\alpha,\beta)\) |
| \(\text{Exponencial}(\theta)\) | \(\theta \sim Gamma(\alpha,\beta)\) |
| \(\text{Normal}(\mu,\sigma^2)\) | \(\mu \sim N(\mu_0,\sigma^{2}_{0}) \quad \text{ y } \quad \sigma^{2} \sim GammaInversa(\alpha,\beta)\) |
Ejemplo 4.9 Considere un conjunto de datos que sigue una distribución \(\mathrm{Poisson}(\lambda)\). Supongamos que nuestra creencia previa sobre \(\lambda\) está dada por:
\[\begin{equation*} \pi(\lambda) = \begin{cases}2e^{-2\lambda} & \lambda> 0 \\ 0 & \lambda \geq 0\end{cases} . \end{equation*}\]
Esta previa de \(\lambda\) sigue una distribución \(\Gamma(1,2)\).
Si tenemos una muestra aleatoria de tamaño \(n\), nos interesa determinar cuántas observaciones necesitamos para que la varianza de nuestra estimación de \(\lambda\) sea, a lo sumo, 0.01.
Para abordar esta pregunta, primero definimos \(y\) como la suma de nuestras observaciones, es decir,\(y=\sum_{i=1}^{n}X_i\). Utilizando el teorema de Bayes, podemos expresar nuestras distribuciones de la siguiente manera:
\[\begin{align*} \pi(\lambda) &\propto \lambda^{\alpha-1} e^{-\beta\lambda} \\ f_n(x\vert \lambda) & \propto e^{-n\lambda} \lambda^{y} \\ \pi(\lambda\vert X) & \propto \lambda^{\alpha+y-1} e^{-(\beta+n)\lambda}. \end{align*}\]
Para este caso, la distribución posterior es \(\lambda|x \sim \Gamma(y+1,n+2)\). La varianza de la distribución Gamma está dada por: \[\begin{equation*} \mathrm{Var}(\lambda) = \dfrac{\alpha}{\beta^2} = \dfrac{y+1}{(n+2)^2}. \end{equation*}\]
Sustituyendo los valores que tenemos:
\[\begin{eqnarray*} & \frac{y+1}{(n+2)^2} &\leq 0.01 \\ \implies& \frac{1}{(n+2)^2} &\leq \frac{y+1}{(n+2)^2} \leq 0.01 \\ \implies& 100 &\leq (n+2)^2 \\ \implies& n &\geq 8 \end{eqnarray*}\]
Por lo tanto, necesitamos al menos 8 observaciones para alcanzar la varianza deseada de 0.01 en nuestra estimación de \(\lambda\).
Retomemos el caso normal para entender cómo se construye la posterior en este caso.
Teorema 4.2 Si \(X_1,\dots,X_n \sim N(\theta, \sigma^2)\) con \(\sigma^2\) conocido y la previa es \(\theta \sim N(\mu_0,V_0^2)\), entonces \(\theta|X\sim N(\mu_1,V_1^2)\) donde \[\begin{equation*} \mu_1 = \dfrac{\sigma^2\mu_0 + nV_0^2 \bar{X}_n}{\sigma^2 + nV_0^2}, \quad V_1^2 = \dfrac{\sigma^2V_0^2}{\sigma^2 + nV_0^2} \end{equation*}\]
Demostración. Construyamos cada elemento por separado:
Entonces
\[ f_n(X|\theta) \propto \exp\left[-\dfrac{n}{2\sigma ^2}(\bar{X} - \theta )^2\right]. \]
Previa: \[ \pi(\theta) \propto \exp\left[-\dfrac{1}{2V_0^2}(\theta - \mu_0)^2\right]. \]
Posterior: \[ \pi(\theta|X) \propto \exp\left[-\dfrac{n}{2\sigma ^2}(\bar{X} - \theta )^2-\dfrac{1}{2V_0^2}(\theta - \mu_0)^2\right]. \]
Con \(\mu_1\) y \(V_1^2\) definidos anteriormente, se puede comprobar la siguiente identidad:
\[ -\dfrac{n}{\sigma ^2}(\bar{X} - \theta )^2-\dfrac{1}{V_0^2}(\theta - \mu_0)^2= \dfrac{1}{V_1^2}(\theta-\mu_1)^2 + \underbrace{\dfrac{n}{\sigma^2 + nV_0^2}(\bar{X}_n- \mu_0)^2}_{\text{Constante con respecto a }\theta} \]
Por lo tanto,
\[ \pi(\theta|X) \propto \exp\left[-\dfrac{1}{2V_1^2}(\theta -\mu_1)^2\right] \]
Corolario 4.1 Si definimos
\[ \mu_1 = \underbrace{\dfrac{\sigma^2}{\sigma^2 + nV_0^2}}_{W_1}\mu_0 + \underbrace{\dfrac{nV_0^2}{\sigma^2 + nV_0^2}}_{W_2} \bar{X}_n \]
Entonces,
Si \(V_0^2\) y \(\sigma^2\) son fijos, entonces \(W_1 \xrightarrow[n\to \infty]{}0\) (la importancia de la media empírica crece conforme aumenta \(n\)).
Si \(V_0^2\) y \(n\) son fijos, entonces \(W_2 \xrightarrow[\sigma^2\to \infty]{}0\) (la importancia de la media empírica decrece conforme la muestra es menos precisa).
Si \(\sigma^2\) y \(n\) son fijos, entonces \(W_2 \xrightarrow[V_0^2\to \infty]{}1\) (la importancia de la media empírica crece conforma la previa es menos precisa).
Ejemplo 4.10 Consideremos el conjutno de datos que vienen de una distribución normal \(N(\theta, \sigma^2)\) con varianza conocida \(\sigma^2=0.25\). Si asumimos que la media poblacional \(\theta\) tiene una distribución previa \(N(\mu_{0}, V_{0}^2)\) con \(\mu_0=2\) y \(V_0^2=1\). Se nos pide usar actualizar nuestr previa con los siguientes datos:
x <- c(
1.28,
0.72,
1.17,
1.31,
1.16,
1.45,
1.08,
1.22,
0.60,
1.32,
1.32,
1.47,
1.24,
1.44,
0.71,
0.51,
0.49,
1.49,
1.38,
1.33,
1.20,
0.86,
0.78,
0.57,
0.95,
1.79,
2.20,
2.27,
1.78,
1.87,
1.83,
2.94,
1.26,
1.16,
1.73,
1.45,
1.31,
1.51,
1.80,
1.47,
1.15,
1.06,
0.97,
2.01,
1.12,
1.39
)
ggplot(data.frame(x = x), aes(x = x)) +
geom_histogram(fill = "white", color = "black", bins = 15) +
xlab("x") +
ylab("frecuencia") +
cowplot::theme_cowplot()
Entonces \(\bar{X}_n = 1.329\) y tenemos que \(n = 46\)
\[\begin{align*} \mu_1 &= \frac{0.25}{0.25+46\times 1} 2 + \frac{46 \times 1}{0.25 + 46 \times 1} 1.329 = 1.333 \end{align*}\]
Comparando la previa, la verosimilitud y la posterior, tenemos que:
Ejemplo 4.11 (Determinación de n) Sean \(X_1,\dots, X_n \sim N(\theta,1)\) y \(\theta\sim N(\mu_0,4)\). Sabemos que \[V_1^2 = \dfrac{\sigma^2V_0^2}{\sigma^2 + nV_0^2}. \] Buscamos que \(V_1\leq 0.01\), entonces \[ \dfrac{4}{4n+1}\leq 0.01 \implies n\geq 99.75 \text{ (al menos 100 observaciones)}\]
Los siguientes ejemplos ilustran cómo el tamaño de la muestra y la varianza de los datos determinan si la posterior se acerca más a la previa o a la verosimilitud.
Caso 1: pocos datos (\(n=3\)) — la previa domina
set.seed(42)
x <- rnorm(n = 3, mean = 10, sd = 1)
(mu <- mean(x))[1] 10.3898(sigma <- sd(x))[1] 0.9681038(n <- length(x))[1] 3(mu_previa <- -10)[1] -10(sigma_previa <- 1)[1] 1mu_posterior <- ((sigma^2) / (sigma^2 + n * sigma_previa^2)) *
mu_previa +
((n * sigma_previa^2) / (sigma^2 + n * sigma_previa^2)) * mu
mu_posterior[1] 5.536168sigma2_posterior <- (sigma^2 * sigma_previa^2) / (sigma^2 + n * sigma_previa^2)
sigma2_posterior[1] 0.238042ggplot(data = data.frame(x = c(-15, 15)), aes(x)) +
stat_function(
fun = dnorm,
args = list(mean = mu_previa, sd = sigma_previa),
aes(color = "Previa"),
linewidth = 2
) +
stat_function(
fun = dnorm,
args = list(mean = mu_posterior, sd = sqrt(sigma2_posterior)),
aes(color = "Posterior"),
linewidth = 2
) +
stat_function(
fun = dnorm,
args = list(mean = mu, sd = sigma),
aes(color = "Verosimilitud"),
linewidth = 2
) +
cowplot::theme_cowplot()
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
rng = np.random.default_rng(42)
x_small = rng.normal(loc=10, scale=1, size=3)
mu_obs = np.mean(x_small)
sigma_obs = np.std(x_small, ddof=1)
n_obs = len(x_small)
mu_previa_n = -10.0
sigma_previa_n = 1.0
sigma2_obs = sigma_obs**2
sigma2_prev = sigma_previa_n**2
mu_post = (sigma2_obs * mu_previa_n + n_obs * sigma2_prev * mu_obs) / (
sigma2_obs + n_obs * sigma2_prev
)
sigma2_post = (sigma2_obs * sigma2_prev) / (sigma2_obs + n_obs * sigma2_prev)
print(f"Media posterior = {mu_post:.4f}, Var = {sigma2_post:.4f}")Media posterior = 5.5128, Var = 0.2246x_plot = np.linspace(-15, 15, 500)
fig, ax = plt.subplots()
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_previa_n, scale=sigma_previa_n),
label="Previa",
linewidth=2,
)
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_obs, scale=sigma_obs),
label="Verosimilitud",
linewidth=2,
)
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_post, scale=np.sqrt(sigma2_post)),
label="Posterior",
linewidth=2,
)
ax.set_xlabel("θ")
ax.set_ylabel("Densidad")
ax.legend()
ax.set_title("Normal conjugada — n=3 (previa domina)")
plt.tight_layout()
plt.show()
Caso 2: muchos datos (\(n=100\)) — la verosimilitud domina
set.seed(42)
x <- rnorm(n = 100, mean = 10, sd = 1)
(mu <- mean(x))[1] 10.03251(sigma <- sd(x))[1] 1.041357(n <- length(x))[1] 100(mu_previa <- 0)[1] 0(sigma_previa <- 1)[1] 1mu_posterior <- ((sigma^2) / (sigma^2 + n * sigma_previa^2)) *
mu_previa +
((n * sigma_previa^2) / (sigma^2 + n * sigma_previa^2)) * mu
mu_posterior[1] 9.924887sigma2_posterior <- (sigma^2 * sigma_previa^2) / (sigma^2 + n * sigma_previa^2)
sigma2_posterior[1] 0.01072791ggplot(data = data.frame(x = c(-5, 15)), aes(x)) +
stat_function(
fun = dnorm,
args = list(mean = mu_previa, sd = sigma_previa),
aes(color = "Previa"),
linewidth = 2
) +
stat_function(
fun = dnorm,
args = list(mean = mu_posterior, sd = sqrt(sigma2_posterior)),
aes(color = "Posterior"),
linewidth = 2
) +
stat_function(
fun = dnorm,
args = list(mean = mu, sd = sigma),
aes(color = "Verosimilitud"),
linewidth = 2
) +
cowplot::theme_cowplot()
rng = np.random.default_rng(42)
x_large = rng.normal(loc=10, scale=1, size=100)
mu_obs100 = np.mean(x_large)
sigma_obs100 = np.std(x_large, ddof=1)
n_obs100 = len(x_large)
mu_previa_n2 = 0.0
sigma_previa_n2 = 1.0
sigma2_obs100 = sigma_obs100**2
sigma2_prev2 = sigma_previa_n2**2
mu_post100 = (sigma2_obs100 * mu_previa_n2 + n_obs100 * sigma2_prev2 * mu_obs100) / (
sigma2_obs100 + n_obs100 * sigma2_prev2
)
sigma2_post100 = (sigma2_obs100 * sigma2_prev2) / (
sigma2_obs100 + n_obs100 * sigma2_prev2
)
print(f"Media posterior = {mu_post100:.4f}, Var = {sigma2_post100:.6f}")Media posterior = 9.8901, Var = 0.005996x_plot = np.linspace(-5, 15, 500)
fig, ax = plt.subplots()
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_previa_n2, scale=sigma_previa_n2),
label="Previa",
linewidth=2,
)
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_obs100, scale=sigma_obs100),
label="Verosimilitud",
linewidth=2,
)
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_post100, scale=np.sqrt(sigma2_post100)),
label="Posterior",
linewidth=2,
)
ax.set_xlabel("θ")
ax.set_ylabel("Densidad")
ax.legend()
ax.set_title("Normal conjugada — n=100 (verosimilitud domina)")
plt.tight_layout()
plt.show()
Caso 3: datos muy variables (\(n=1000\), \(\sigma=10\)) — la previa recupera influencia
set.seed(42)
x <- rnorm(n = 1000, mean = 10, sd = 10)
(mu <- mean(x))[1] 9.741756(sigma <- sd(x))[1] 10.02521(n <- length(x))[1] 1000(mu_previa <- 0)[1] 0(sigma_previa <- 1)[1] 1mu_posterior <- ((sigma^2) / (sigma^2 + n * sigma_previa^2)) *
mu_previa +
((n * sigma_previa^2) / (sigma^2 + n * sigma_previa^2)) * mu
mu_posterior[1] 8.852079sigma2_posterior <- (sigma^2 * sigma_previa^2) / (sigma^2 + n * sigma_previa^2)
sigma2_posterior[1] 0.09132616ggplot(data = data.frame(x = c(-5, 15)), aes(x)) +
stat_function(
fun = dnorm,
args = list(mean = mu_previa, sd = sigma_previa),
aes(color = "Previa"),
linewidth = 2
) +
stat_function(
fun = dnorm,
args = list(mean = mu_posterior, sd = sqrt(sigma2_posterior)),
aes(color = "Posterior"),
linewidth = 2
) +
stat_function(
fun = dnorm,
args = list(mean = mu, sd = sigma),
aes(color = "Verosimilitud"),
linewidth = 2
) +
cowplot::theme_cowplot()
rng = np.random.default_rng(42)
x_noisy = rng.normal(loc=10, scale=10, size=1000)
mu_obs_n = np.mean(x_noisy)
sigma_obs_n = np.std(x_noisy, ddof=1)
n_obs_n = len(x_noisy)
mu_previa_n3 = 0.0
sigma_previa_n3 = 1.0
sigma2_obs_n = sigma_obs_n**2
sigma2_prev3 = sigma_previa_n3**2
mu_post_n = (sigma2_obs_n * mu_previa_n3 + n_obs_n * sigma2_prev3 * mu_obs_n) / (
sigma2_obs_n + n_obs_n * sigma2_prev3
)
sigma2_post_n = (sigma2_obs_n * sigma2_prev3) / (sigma2_obs_n + n_obs_n * sigma2_prev3)
print(f"Media posterior = {mu_post_n:.4f}, Var = {sigma2_post_n:.6f}")Media posterior = 8.8455, Var = 0.089133x_plot = np.linspace(-5, 15, 500)
fig, ax = plt.subplots()
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_previa_n3, scale=sigma_previa_n3),
label="Previa",
linewidth=2,
)
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_obs_n, scale=sigma_obs_n),
label="Verosimilitud",
linewidth=2,
)
ax.plot(
x_plot,
norm.pdf(x_plot, loc=mu_post_n, scale=np.sqrt(sigma2_post_n)),
label="Posterior",
linewidth=2,
)
ax.set_xlabel("θ")
ax.set_ylabel("Densidad")
ax.legend()
ax.set_title("Normal conjugada — n=1000, σ=10 (alta varianza)")
plt.tight_layout()
plt.show()
Definición 4.2 Sea \(\pi\) una función positiva cuyo dominio está en \(\Omega\). Si la ingtegral de \(\pi(\theta)\) sobre \(\Omega\) es infinity, i.e. \(\int\pi(\theta)\;d\theta=\infty\); entonces decimos que \(\pi\) es una densidad impropia. Por ejemplo \(\theta \sim \text{Unif}(\mathbb{R})\), \(\lambda \sim \text{Unif}(0,\infty)\).
Observación. Una técnica para seleccionar distribuciones impropia es sustituir los hiperparámetros previos por 0.
Ejemplo 4.12 La siguiente tabla muestra el número de soldados prusianos muertos por una patada de caballo en 14 unidades de combate por 20 años, para contar 280 unidades. La siguiente tabla resume, cuántas de esas unidades tuvieron este tipo de muertes:
| Unidades | Ocurrencias |
|---|---|
| 144 | 0 |
| 91 | 1 |
| 32 | 2 |
| 11 | 3 |
| 2 | 4 |
Verosimilitud: Podemos modelar el número de muertes como una distribución Poisson: \(X_1 = 0, X_2 = 1, X_3 = 1,\dots, X_{280} = 0 \sim \text{Poisson}(\lambda)\).
Previa: Para hacer que la posterior conjugue, entonces \(\lambda \sim \Gamma(\alpha, \beta)\).
Posterior: De acuerdo a lo que hemos visto anteriormente \(\lambda|X \sim \Gamma(y+\alpha, n+\beta) = \Gamma(196 + \alpha, 280 + \beta)\).
El problema acá es que no sabemos cuál valor de \(\alpha\) o \(\beta\) deberíamso tomar. Una opción es tomar \(\alpha = \beta = 0\), lo cual nos da una distribución impropia. En este caso:
\[\begin{align*} \pi(\lambda) &= \dfrac{1}{\Gamma(\alpha)}\beta^\alpha\lambda^{\alpha-1}e^{-\beta\lambda} \\ & \propto \lambda^{\alpha-1}e^{-\lambda\beta} \\ &=\dfrac{1}{\lambda} \end{align*}\] donde \(\displaystyle\int_{0}^{\infty}\dfrac{1}{\lambda} d\lambda = \infty\).
Por teorema de Bayes, \[ \theta|X \sim \Gamma(196,280) \]
Gráficamente tenemos la versoimilitud, la previa y la posterior:
Una vez que obtenemos la distribución posterior \(\pi(\theta|X)\), surge la pregunta: ¿cómo resumir toda esa distribución en un único número que sirva como estimación de \(\theta\)? La respuesta depende del contexto del problema: no es lo mismo subestimar que sobreestimar un parámetro. Las funciones de pérdida formalizan este costo asimétrico y nos permiten seleccionar el estimador que minimiza el error esperado según las consecuencias de equivocarse.
Definición 4.3 Sean \(X_1,\dots, X_n\) datos observables cuyo modelo está indexado por \(\theta\in\Omega\). Un estimador de \(\theta\) es cualquier estadístico \(\delta(X_1,\dots, X_n)\).
Definición 4.4
Definición 4.5 Una función de pérdida es una función de dos variables: \[ L(\theta,a), \quad \theta \in\Omega\] con \(a\) un número real. Es decir, es lo que se pierde cuando el parámetro es \(\theta\) y se quiere estimar con el valor \(a\).
Asuma que \(\theta\) tiene una previa. La pérdida esperada es \[ \mathbb{E}[L(\theta,a)] = \int_{\Omega}L(\theta, a) \pi(\theta)\;d\theta\] la cual es una función de \(a\), que a su vez es función de \(X_1,\dots,X_n\). Asuma que \(a\) se selecciona el minimizar esta esperanza. A ese estimador \(a = \delta^*(X_1,\dots, X_n)\) se le llama estimador bayesiano, si ponderamos los parámetros con respecto a la posterior.
\[ \mathbb{E}[L(\theta, \delta^*)|X] = \int_{\Omega}L(\theta, a) \pi(\theta|X)\;d\theta = \min_a \mathbb{E}[L(\theta,a)|X]. \]
Si \(L(\theta,a) = (\theta-a)^2\), entonces la pérdida esperada es
\[\begin{align*} \mathbb{E}[L(\theta,a)] &= \int_{\Omega}(\theta-a)^2 \pi(\theta\mid X)\;d\theta \\ &= \int_{\Omega}\theta^2\pi(\theta\mid X)\;d\theta - 2a\int_{\Omega}\theta\pi(\theta\mid X)\;d\theta + a^2\int_{\Omega}\pi(\theta \mid X)\;d\theta \\ &= \mathbb{E}[\theta^2\mid X] - 2a\mathbb{E}[\theta\mid X] + a^2 \\ \end{align*}\]
Note que la pérdida esperada es una función cuadrática en \(a\). Para minimizarla, derivamos con respecto a \(a\) e igualamos a 0.
\[\begin{align*} \dfrac{\partial}{\partial a} \mathbb{E}[L(\theta,a)] &= -2\mathbb{E}[\theta\mid X] + 2a = 0 \\ \implies a &= \mathbb{E}[\theta\mid X] \end{align*}\]
En otras palabras el estimador bayesiano es la esperanza posterior. Formalmente:
\[\begin{equation*} \delta^*(X_1,\dots, X_n) = \mathbb{E}[\theta|X] \text{ cuando } L(\theta,a) = (\theta-a)^2. \end{equation*}\]
Ejemplo 4.13 Sea \(X_1,\dots, X_n \sim \text{Ber}(\theta)\), \(\theta \sim \text{Beta}(\alpha,\beta) \implies \theta|X \sim \text{Beta}(\alpha+y,\beta+n-y)\).
El estimador de \(\theta\) es \[\begin{align*} \delta(X_1,\dots, X_n) &= \dfrac{\alpha+y}{\alpha + \beta + n} \\ &= \overbrace{\dfrac{\alpha}{\alpha + \beta} }^{\text{Esperanza previa}}\cdot \dfrac{\alpha +\beta}{\alpha +\beta + n} + \overbrace{\dfrac{y}{n}}^{\bar{X}}\cdot \dfrac{n}{\alpha +\beta + n}. \end{align*}\]
Si \(L(\theta,a) = |\theta-a|\), entonces la pérdida esperada es
\[\begin{align*} \mathbb{E}[L(\theta,a)|X] &= \int_{-\infty}^{+\infty}|\theta-a|\pi(\theta|X)\;d\theta \\ &= \int_{a}^{+\infty}(\theta-a)\pi(\theta|X)\;d\theta + \int_{-\infty}^{a}(a-\theta)\pi(\theta|X)\;d\theta\\ &= -\int_{+\infty}^{a}(\theta-a)\pi(\theta|X)\;d\theta + \int_{-\infty}^{a}(a-\theta)\pi(\theta|X)\;d\theta\\ \end{align*}\]
La regla de Leibniz dice que si \(F(x) = \int_{-\infty}^{x}f(x,t)\;dt\), entonces
\[\begin{equation*} \dfrac{d}{dx}F(x) = f(x,x) + \int_{-\infty}^{x}\dfrac{\partial}{\partial x}f(x,t)\;dt \end{equation*}\]
Por lo tanto para este caso tenemos
\[\begin{align*} \dfrac{d}{da}\mathbb{E}[L(\theta,a)|X] &= -\left((a-a)\pi(a|X) - \int_{+\infty}^{a}\pi(\theta|X)\;d\theta\right) \\ & \qquad + (a-a) \pi(a|X) + \int_{-\infty}^{a}\pi(\theta|X)\;d\theta \\ &= \int_{-\infty}^{a} \pi(\theta|X)\;d\theta - \int_{a}^{+\infty} \pi(\theta|X)\;d\theta \\ \end{align*}\]
Igualando esta última expresión a 0, tenemos que
\[\begin{align*} \int_{-\infty}^{a} \pi(\theta|X)\;d\theta &= \int_{a}^{+\infty} \pi(\theta|X)\;d\theta \\ \implies F_{\theta\mid X}(a) &= \int_{-\infty}^{a} \pi(\theta|X)\;d\theta = \dfrac{1}{2} \\ \end{align*}\]
La mediana es el punto de \(m\) tal que \(F_{\theta\mid X}(m) = \dfrac{1}{2}\).
Observación. Bajo la función de pérdida absoluta, el estimador bayesiano es la mediana posterior.
Ejemplo 4.14 Para el caso Bernoulli. \[ \dfrac{1}{\text{Beta}(\alpha+y, \beta+n-y)}\int_{-\infty}^{X_{0.5}}\theta^{\alpha+y-1} (1-\theta)^{\beta+n-y-1}\;d\theta = \dfrac12\] Encontrar el valor que \(X_{0.5}\) que resuelve esta ecuación es díficil y se debe hacer usando métodos numéricos.
Por ejemplo, si tuvieramos que la posterior fuera una \(\mathrm{Beta}(22,18)\) entonces podríamos resolverlo en R de la siguiente forma
encuentra_mediana_beta <- function(q, alpha, beta) {
abs(pbeta(q, alpha, beta) - 1 / 2)
}
opt <- optimize(encuentra_mediana_beta, c(0, 1), alpha = 22, beta = 18)
opt$minimum[1] 0.5508349Note que la media sería \(\frac{22}{40}=0.55\).
La función de pérdida 0-1 penaliza cualquier error de estimación de forma uniforme: no importa cuán lejos esté la estimación del verdadero valor, la pérdida es siempre 1 si hay error y 0 si se acierta exactamente. Formalmente, para un parámetro continuo se define mediante un umbral \(\epsilon > 0\):
\[ L(\theta, a) = \begin{cases} 0 & \text{si } |\theta - a| < \epsilon \\ 1 & \text{si } |\theta - a| \geq \epsilon \end{cases} \]
La pérdida esperada posterior es:
\[\begin{align*} \mathbb{E}[L(\theta, a) \mid X] &= \int_\Omega L(\theta, a)\, \pi(\theta \mid X)\, d\theta \\ &= \int_{|\theta - a| \geq \epsilon} \pi(\theta \mid X)\, d\theta \\ &= 1 - \int_{a - \epsilon}^{a + \epsilon} \pi(\theta \mid X)\, d\theta \\ &= 1 - P(a - \epsilon < \theta < a + \epsilon \mid X). \end{align*}\]
Para minimizar esta expresión, debemos maximizar \(P(a - \epsilon < \theta < a + \epsilon \mid X)\), es decir, encontrar el valor \(a\) alrededor del cual la densidad posterior concentra la mayor probabilidad en un intervalo de radio \(\epsilon\).
Tomando el límite cuando \(\epsilon \to 0\), la integral se aproxima por continuidad como:
\[ \int_{a-\epsilon}^{a+\epsilon} \pi(\theta \mid X)\, d\theta \approx 2\epsilon\, \pi(a \mid X). \]
Esta cantidad se maximiza cuando \(\pi(a \mid X)\) es máxima, es decir, cuando \(a\) es el máximo de la densidad posterior.
Observación. Bajo la función de pérdida 0-1, el estimador bayesiano es la moda posterior (también llamada estimación MAP, Maximum A Posteriori):
\[ \delta^*(X_1, \dots, X_n) = \arg\max_\theta\, \pi(\theta \mid X). \]
Ejemplo 4.15 (Pérdida 0-1 en el caso Bernoulli) Sea \(X_1, \dots, X_n \sim \text{Ber}(\theta)\) y \(\theta \mid X \sim \text{Beta}(\alpha + y,\, \beta + n - y)\). La moda de una distribución \(\text{Beta}(a, b)\) con \(a > 1\) y \(b > 1\) es:
\[ \text{Moda} = \frac{a - 1}{a + b - 2}. \]
Por lo tanto, el estimador bayesiano bajo pérdida 0-1 es:
\[ \delta^*(X_1, \dots, X_n) = \frac{\alpha + y - 1}{\alpha + \beta + n - 2}. \]
Si la previa es uniforme (\(\alpha = \beta = 1\)), entonces:
\[ \delta^* = \frac{y}{n} = \bar{X}_n, \]
es decir, bajo previa uniforme y pérdida 0-1, el estimador bayesiano coincide con el estimador de máxima verosimilitud.
\(L(\theta,a) = |\theta-a|^k\), \(k\ne 1,2\), \(0<k<1\).
\(L(\theta,a) = \lambda(\theta)|\theta-a|^2\) (\(\lambda(\theta)\) penaliza la magnitud del parámetro).
\(L(\theta,a)=\begin{cases}3(\theta-a)^2 & \theta\leq a \text{ (sobreestima)}\\ (\theta-a)^2&\theta\geq a \text{ (subestima)} \end{cases}\)
Retomando el ejemplo de los componentes eléctricos, la distribución posterior fue \(\theta|X \sim \Gamma(9, 36178)\). Calculamos los estimadores bajo cada función de pérdida.
alpha <- 9
beta <- 36178
# Pérdida cuadrática: media posterior
(theta_media <- alpha / beta)[1] 0.00024877cat(
"Tiempo medio estimado (pérdida cuadrática):",
round(1 / theta_media, 1),
"unidades\n"
)Tiempo medio estimado (pérdida cuadrática): 4019.8 unidades# Pérdida absoluta: mediana posterior
encuentra_mediana_gamma <- function(q, alpha, beta) {
abs(pgamma(q, alpha, beta) - 1 / 2)
}
opt <- optimize(encuentra_mediana_gamma, c(0, 0.01), alpha = alpha, beta = beta)
theta_mediana <- opt$minimum
cat(
"Tiempo medio estimado (pérdida absoluta):",
round(1 / theta_mediana, 1),
"unidades\n"
)Tiempo medio estimado (pérdida absoluta): 100.8 unidades# Pérdida 0-1: moda posterior (MAP)
# La moda de Γ(α, β) es (α - 1) / β para α ≥ 1
(theta_moda <- (alpha - 1) / beta)[1] 0.0002211289cat(
"Tiempo medio estimado (pérdida 0-1 / MAP):",
round(1 / theta_moda, 1),
"unidades\n"
)Tiempo medio estimado (pérdida 0-1 / MAP): 4522.2 unidadesObservación. Los tres estimadores difieren porque la distribución Gamma no es simétrica. En distribuciones asimétricas siempre se cumple: moda \(\leq\) mediana \(\leq\) media (cuando la cola es hacia la derecha), o la relación inversa cuando la cola es hacia la izquierda. Cuanto más asimétrica sea la posterior, mayor será la discrepancia entre ellos.
from scipy.stats import gamma as gamma_dist
from scipy.optimize import minimize_scalar
import numpy as np
alpha_py = 9
beta_py = 36178
# Pérdida cuadrática: media posterior
theta_media_py = alpha_py / beta_py
print(
f"Media posterior (pérdida cuadrática): θ = {theta_media_py:.6f} → 1/θ = {1 / theta_media_py:.1f} unidades"
)Media posterior (pérdida cuadrática): θ = 0.000249 → 1/θ = 4019.8 unidades# Pérdida absoluta: mediana posterior
result = minimize_scalar(
lambda q: abs(gamma_dist.cdf(q, a=alpha_py, scale=1 / beta_py) - 0.5),
bounds=(0, 0.01),
method="bounded",
)
theta_mediana_py = result.x
print(
f"Mediana posterior (pérdida absoluta): θ = {theta_mediana_py:.6f} → 1/θ = {1 / theta_mediana_py:.1f} unidades"
)Mediana posterior (pérdida absoluta): θ = 0.009996 → 1/θ = 100.0 unidades# Pérdida 0-1: moda posterior (MAP)
# La moda de Γ(α, β) con α ≥ 1 es (α - 1) / β
theta_moda_py = (alpha_py - 1) / beta_py
print(
f"Moda posterior (pérdida 0-1 / MAP): θ = {theta_moda_py:.6f} → 1/θ = {1 / theta_moda_py:.1f} unidades"
)Moda posterior (pérdida 0-1 / MAP): θ = 0.000221 → 1/θ = 4522.2 unidadesEjemplo 4.16 Suponga que se tiene un examen con ítemes malos (proporción: \(\theta\)), \(\theta \in [0,1]\). Se decide usar una función de pérdida cuadrática. El tamaño de muestra son \(n=100\) ítemes, de los cuales \(y=10\) están malos.
\[ X_1,\dots,X_n\sim \text{Ber}(\theta) \]
Usemos diferentes previas para ver su efecto en los estimadores
\[ \mathbb{E}[\theta|X] = \dfrac{\alpha+y}{\alpha+\beta+n} = \dfrac{1+10}{2+100} = 0.108 \]
\[ \mathbb{E}[\theta|X] = \dfrac{1+10}{1+2+100} = \dfrac{11}{103}=0.107 \]
Las previas son muy diferentes una de la otra, pero el estimador final es muy parecido debido a la cantidad de datos. Compare esto con la media empírica \(\bar{X}_n = \dfrac{10}{100} = 0.1\).
Definición 4.6 Un estimador de \(\theta\) \(\delta(X_1,\dots, X_n)\) es consistente si \[\delta(X_1,\dots, X_n)\xrightarrow[n\to \infty]{\mathbb{P}}\theta.\]
Bajo pérdida cuadrática, \(\mathbb{E}[\theta|X] = W_1\mathbb{E}[\theta] + W_2\bar{X}_n = \delta^*\). Sabemos, por ley de grandes números, que \(\bar{X}_n \xrightarrow[n\to \infty]{\mathbb{P}}\theta\). Además, \(W_1\xrightarrow[n\to \infty]{}0\) y \(W_2\xrightarrow[n\to \infty]{}1\).
En los ejemplos que hemos analizado \[ \delta^* \xrightarrow[n\to \infty]{\mathbb{P}}\theta \]
Teorema 4.3 Bajo condiciones mencionadas anteriormente, los estimadores bayesianos son consistentes.
De forma general podemos definir un estimador de la siguiente forma
Definición 4.7 (Estimador) Si \(X_1,\dots, X_n\) es una muestra en un modelo indexado por \(\theta\), \(\theta \in \Omega\) (\(k\)-dimensiones), sea
\[ h:\Omega \to H \subset \mathbb{R}^d. \]
Defina \(\psi = h(\theta)\). Un estimador de \(\psi\) es una función o estadístico \(\delta(X_1,\dots, X_n) \in H\).
Observación. Ojo, la definición de estimador solo requiere que el parámetro \(\psi\) y el estimador \(\delta\) pertenenzca el mismo espacio. Puede ser que un \(\delta_1\) sea mejor que un \(\delta_2\).
Ejemplo 4.17 \(X_1,\dots, X_n \sim \text{Exp}(\theta)\), \(\theta|X \sim \Gamma(\alpha,\beta) = \Gamma (4,8.6)\). La característica de interés es \(\psi = \dfrac{1}\theta\), el valor esperado del tiempo de fallo.
Es estimador se calcula de la siguiente manera:
\[\begin{align*} \delta(x) = \mathbb{E}[\psi|x] & = \int_{0}^\infty \dfrac{1}\theta\pi(\theta|x)\;d\theta\\ & = \int_{0}^\infty \dfrac{1}\theta \dfrac{8.6^4}{\Gamma(4)} \theta^3e^{-8.6\theta}\;d\theta\\ &=\dfrac{8.6^4}{6} \int_{0}^\infty \theta^2 e^{-8.6\theta}\;d\theta\\ &= \dfrac{8.6^4}{6}\frac{\Gamma(3)}{8.6^3} \\ & = \dfrac{8.6^4}{6}\dfrac{2}{8.6^3} \\ & = 2.867 \text{ unidades de tiempo.} \end{align*}\]
Por otro lado, vea que \(\mathbb{E}(\theta|X) = \dfrac{4}{8.6}\). El estimador plug-in correspondería a \[\dfrac{1}{\mathbb{E}(\theta|X)} = \dfrac{8.6}{4} = 2.15.\]
Suponga que se que quiere averiguar si los estudiantes de cierto colegio duermen más de 8 horas o menos de 8 horas.
Para esto primero cargaremos el siguiente paquete,
library(LearnBayes)Suponga que se hace una encuesta a 27 estudiantes y se encuentra que 11 dicen que duermen más de 8 horas diarias y el resto no. Nuestro objetivo es encontrar inferencias sobre la proporción \(p\) de estudiantes que duermen al menos 8 horas diarias. El modelo más adecuado es
\[ f(x \vert p) \propto p^s (1-p)^f \]
donde \(s\) es la cantidad de estudiantes que duermen más de 8 horas y \(f\) los que duermen menos de 8 horas.
Una primera aproximación para la previa es usar una distribución discreta. En este caso, el investigador asigna una probabilidad a cierta cantidad de horas de sueño, según su experiencia. Así, por ejemplo:
p <- seq(0.05, 0.95, by = 0.1)
prior <- c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior <- prior / sum(prior)
plot(p, prior, type = "h", ylab = "Probabilidad Previa")
El paquete LearnBayes tiene la función pdisc que estima la distribución posterior para una previa discreta binomial. Recuerde que el valor 11 representa la cantidad de estudiantes con más de 8 horas de sueño y 16 lo que no duermen esa cantidad.
data <- c(11, 16)
post <- pdisc(p, prior, data)
round(cbind(p, prior, post), 2) p prior post
[1,] 0.05 0.03 0.00
[2,] 0.15 0.18 0.00
[3,] 0.25 0.28 0.13
[4,] 0.35 0.25 0.48
[5,] 0.45 0.16 0.33
[6,] 0.55 0.07 0.06
[7,] 0.65 0.02 0.00
[8,] 0.75 0.00 0.00
[9,] 0.85 0.00 0.00
[10,] 0.95 0.00 0.00Y podemos ver la diferencia entre la previa (negro) y la posterior (roja),
plot(p, post, type = "h", col = "red")
lines(p + 0.01, prior, type = "h")
¿Qué se puede deducir de estos resultados?
Ejercicio 4.1 Suponga que se tiene la base de datos studentdata. Realice los cálculos anteriores con esos datos,
data("studentdata")
horas_sueno <- studentdata$WakeUp - studentdata$ToSleep
horas_sueno <- na.omit(horas_sueno)
summary(horas_sueno) Min. 1st Qu. Median Mean 3rd Qu. Max.
2.500 6.500 7.500 7.385 8.500 12.500 hist(horas_sueno, main = "")
Ahora supongamos que se tiene quiere ajustar una previa continua a este modelo. Para esto usaremos una distribución Beta con parámetros \(\alpha\) y \(\beta\), de la forma
\[ pi(p\vert \alpha, \beta) \propto p^{1-\alpha} (1-p)^{1-\beta}. \]
El ajuste de los parámetro de la Beta depende mucho de la información previa que se tenga del modelo. Una forma fácil de estimarlo es a través de cuantiles con los cuales se puede reescribir estos parámetros. Para una explicación detallada revisar https://stats.stackexchange.com/a/237849
En particular, suponga que se cree que el \(50\%\) de las observaciones la proporción será menor que 0.3 y que el \(90\%\) será menor que 0.5.
Para esto ajustaremos los siguientes parámetros
quantile2 <- list(p = .9, x = .5)
quantile1 <- list(p = .5, x = .3)
(ab <- beta.select(quantile1, quantile2))[1] 3.26 7.19a <- ab[1]
b <- ab[2]
s <- 11
f <- 16En este caso se obtendra la distribución posterior Beta con paramétros \(\alpha + s\) y \(\beta + f\),
curve(
dbeta(x, a + s, b + f),
from = 0,
to = 1,
xlab = "p",
ylab = "Densidad",
lty = 1,
lwd = 4
)
curve(dbeta(x, s + 1, f + 1), add = TRUE, lty = 2, lwd = 4)
curve(dbeta(x, a, b), add = TRUE, lty = 3, lwd = 4)
legend(
.7,
4,
c("Previa", "Verosimilitud", "Posterior"),
lty = c(3, 2, 1),
lwd = c(3, 3, 3)
)
(p_hat <- 11 / 27)[1] 0.4074074(p_posterior <- (a + s) / (a + s + b + f))[1] 0.3807744En este capítulo se estudió el enfoque bayesiano para la estimación de parámetros. Los conceptos centrales son:
Flujo bayesiano: \[ \underbrace{\pi(\theta|X)}_{\text{posterior}} \propto \underbrace{\pi(\theta)}_{\text{previa}} \times \underbrace{f_n(X|\theta)}_{\text{verosimilitud}} \]
Familias conjugadas más importantes:
| Modelo de datos | Previa conjugada | Posterior |
|---|---|---|
| \(\text{Exponencial}(\theta)\) | \(\Gamma(\alpha, \beta)\) | \(\Gamma(\alpha+n,\, \beta+\sum X_i)\) |
| \(\text{Bernoulli}(\theta)\) | \(\text{Beta}(\alpha,\beta)\) | \(\text{Beta}(\alpha+y,\, \beta+n-y)\) |
| \(\text{Poisson}(\lambda)\) | \(\Gamma(\alpha,\beta)\) | \(\Gamma(\alpha+y,\, \beta+n)\) |
| \(\text{Normal}(\theta,\sigma^2)\), \(\sigma^2\) conocido | \(N(\mu_0, V_0^2)\) | \(N(\mu_1, V_1^2)\) |
Estimadores bayesianos según función de pérdida:
| Función de pérdida | Estimador óptimo |
|---|---|
| \(L(\theta,a) = (\theta-a)^2\) (cuadrática) | \(\delta^* = \mathbb{E}[\theta\mid X]\) (media posterior) |
| \(L(\theta,a) =\vert \theta-a \vert\) (absoluta) | \(\delta^* = \text{mediana de } \pi(\theta\mid X)\) |
| \(L(\theta,a) = \mathbf{1}_{\vert \theta-a\vert \geq\epsilon}\) (0-1) | \(\delta^* = \arg\max_\theta\,\pi(\theta\mid X)\) (moda posterior / MAP) |
Propiedades asintóticas: