1 Introducción

El objetivo consiste en aproximar cualquier cantidad asociada con la distribución posterior de \(\boldsymbol{\theta}\).

Si podemos obtener muestras de \(\boldsymbol{\theta}\) procedentes de la distribución posterior, cualquier cantidad posterior de interés se puede aproximar con un grado de precisión arbitrario usando métodos de Monte Carlo (Stanislaw Ulam, John von Neumann, Nicholas Metropolis).

Los métodos de Monte Carlo son algoritmos computacionales que se basan en un muestreo aleatorio iterativo para obtener resultados numéricos. El concepto subyacente es utilizar la aleatoriedad para resolver problemas tanto estocásticos (generación de muestras de una distribución de probabilidad) como deterministicos (optimización, integración numérica).

Cualquier distribución de probabilidad (y por ende cualquier característica de esa distribución) se puede aproximar arbitrariamente bien tomando tantas muestras aleatorias de esa distribución como sea necesario dependiendo del nivel de precisión que se requiera.

2 Ejemplo

Aproximación de una distribución Gamma con diferentes niveles de precisión.

# parametros distribucion Gamma
a <- 68
b <- 45

# tamaños
m <- c(50, 250, 1250)

# simulacion con diferentes niveles de precision
set.seed(1234)
for (j in 1:length(m))
        assign(x = paste0("theta_sim_", j), value = rgamma(n = m[j], shape = a, rate = b))
# grafico
par(mfrow = c(2,3), mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
for (j in 1:length(m)) {
        hist(x = get(paste0("theta_sim_", j)), prob = T, xlim = c(.75, 2.25), ylim = c(0, 2.7), 
             xlab = expression(theta), ylab = "Densidad", main = paste0("B = ", m[j]), 
             col = "gray90", border = "gray90")
        curve(expr = dgamma(x, shape = a, rate = b), col = "blue", add = T)
}
for (j in 1:length(m)) {
        plot(x = density(get(paste0("theta_sim_", j))), xlim = c(.75, 2.25), ylim = c(0, 2.7), 
             xlab = expression(theta), ylab = "Densidad", main = paste0("B = ", m[j]),
             col = "gray75", lwd = 3)
        curve(expr = dgamma(x, shape = a, rate = b), col = "blue", add = T)
}

3 Implementación

Sea \(\theta\) un parámetro de interés y \(\boldsymbol{y} = (y_1,\dots,y_n)\) un conjunto de observaciones. Supongamos que es posible obtener una muestra aleatoria de \(B\) valores de \(\theta\) asociados con la distribución posterior \(p(\theta\mid \boldsymbol{y})\), esto es, \[ \theta^{(1)},\ldots,\theta^{(B)}\stackrel{\text{iid}}{\sim} p(\theta\mid \boldsymbol{y})\,. \]

La distribución empírica de \(\theta^{(1)},\ldots,\theta^{(B)}\) aproxima la distribución posterior \(p(\theta\mid \boldsymbol{y})\) y tal aproximación puede ser tan precisa como se quiera incrementando el valor de \(B\). Esta aproximación se conoce como la aproximación de Monte Carlo de \(p(\theta\mid \boldsymbol{y})\).

(Ley débil de los grandes números). Sea \(X_1,\ldots,X_n\) una secuencia de variables aleatorias independientes e identicamente distribuidas con media \(\mu\) y varianza finita \(\sigma^2\). Entonces, el promedio muestral \(\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i\) converge en probabilidad a \(\mu\) cuando \(n\rightarrow\infty\), i.e., para todo \(\epsilon > 0\) se tiene que \[ \lim\limits_{n\to\infty}\textsf{Pr}(|\bar{X} - \mu| > \epsilon) = 0\,. \]

La ley débil de los grandes números garantiza que, si \(\theta^{(1)},\ldots,\theta^{(B)}\stackrel{\text{iid}}{\sim} p(\theta\mid \boldsymbol{y})\), entonces \[ \frac{1}{B}\sum_{b=1}^{B} g\left(\theta^{(b)}\right)\longrightarrow \int_\Theta g(\theta)\,p(\theta\mid \boldsymbol{y})\,\text{d}\theta = \textsf{E}(g(\theta)\mid \boldsymbol{y})\,\,\text{cuando $B\rightarrow\infty$}, \] donde \(g(\theta)\) es una función arbitraria de \(\theta\).

Por lo tanto:

4 Inferencia sobre una función arbitraria de \(\theta\)

Los métodos de Monte Carlo permiten hacer fácilmente inferencia posterior sobre cualquier función arbitraria de \(\theta\), digamos \(\gamma = g(\theta)\):

La secuencia \(\gamma^{(1)},\ldots,\gamma^{(B)}\) constituye un conjunto de valores independientes de \(p(\gamma\mid \boldsymbol{y})\) con los cuales se puede hacer inferencia posterior sobre cualquier aspecto de \(\gamma\).

5 Bondad de ajuste

Los métodos de Monte Carlo también permiten examinar detalladamente la distribución predictiva posterior \(p(y^*\mid\boldsymbol{y})\), lo que hace posible chequear la bondad de ajuste interna del modelo por medio estadísticos calculados a partir de la distribución predictiva posterior:

Si \(t_0\) es un valor típico de la distribución de \(t^{(1)},\ldots,t^{(B)}\), entonces se dice que el modelo captura adecuadamente la característica de interés que representa el estadístico de prueba.

Se recomienda evaluar todos aquellos aspectos de la población que sea de interés caracterizar por medio del modelo.

6 Ejemplo

La Encuesta Social General de 1998 recopiló datos sobre el nivel educativo y el número de hijos de una muestra de mujeres que tenían 40 años de edad. En este ejemplo, se comparan las mujeres con títulos universitarios con aquellas que no los tienen, en función del número de hijos.

6.1 Aproximaciones con diferentes niveles de precisión

# modelo Gamma-Poisson
# datos mujeres con educacion superior: s = 66, n = 44
s <- 66
n <- 44

# hiperparametros
a <- 2
b <- 1

# tamaños
m <- c(5, 500, 50000)

# simulacion con diferentes niveles de precision
set.seed(1234)
for (j in 1:length(m))
        assign(x = paste0("theta_sim_", j), value = rgamma(n = m[j], shape = a+s, rate = b+n))

# media posterior exacta
(a+s)/(b+n) 
## [1] 1.511111
# media posterior aproximada
tab <- c(mean(theta_sim_1), mean(theta_sim_2), mean(theta_sim_3))
names(tab) <- paste0("B = ", m)
print(tab)
##     B = 5   B = 500 B = 50000 
##  1.475078  1.516255  1.511432
# intervalo de credibilidad al 95% exacto
qgamma(p = c(.025,.975), shape = a+s, rate = b+n)
## [1] 1.173437 1.890836
# intervalo de credibilidad al 95% aproximado
tab <- NULL
for (j in 1:length(m))
        tab <- rbind(tab, quantile(x = get(paste0("theta_sim_", j)), probs = c(.025,.975)))
colnames(tab) <- c("2.5%", "97.5%")
rownames(tab) <- paste0("B = ", m)
print(tab)
##               2.5%    97.5%
## B = 5     1.295904 1.591065
## B = 500   1.195910 1.898362
## B = 50000 1.174680 1.890198

6.2 Inferencia sobre \(\gamma = \theta_1 - \theta_2\)

# datos mujeres sin educacion superior
s1 <- 217
n1 <- 111

# datos muejeres con educacion superior
s2 <- 66
n2 <- 44   

# hiperparametros
a <- 2
b <- 1

# numero de muestras de MC 
B <- 50000

# simulacion
set.seed(1234)
theta1_mc <- rgamma(n = B, shape = a+s1, rate = b+n1)
theta2_mc <- rgamma(n = B, shape = a+s2, rate = b+n2)
gamma_mc  <- theta1_mc - theta2_mc

# estimacion puntual de gamma
mean(gamma_mc)
## [1] 0.4445117
# intervalo de credibilidad al 95% para gamma
quantile(x = gamma_mc, probs = c(0.025, 0.975))
##         2.5%        97.5% 
## -0.007322471  0.872844678
# grafico
par(mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
plot(density(gamma_mc), main = "", xlim = c(-.7,1.5), xlab = expression(gamma), 
     ylab = expression(paste("p(",gamma," | y)", sep="")), col = 4, lwd = 2)
abline(v = quantile(x = gamma_mc, probs = c(0.025, 0.975)), lty = 2, lwd = 2, col = 2)
abline(v = mean(gamma_mc), lty = 2, lwd = 2, col = 3)
legend("topright", legend = c("Posterior", "IC 95%", "Media"), 
       col = c(4, 2, 3), lty = 1, lwd = 2, bty = "n")

6.3 Distribución predictiva posterior

# distribucion predictiva posterior
set.seed(1234)
y1_mc <- rpois(n = B, lambda = theta1_mc)
y2_mc <- rpois(n = B, lambda = theta2_mc)

# media predictiva posterior exacta (binomal negativa)
c((a+s1)/(b+n1), (a+s2)/(b+n2))
## [1] 1.955357 1.511111
# media predictiva posteriori aproximada
c(mean(y1_mc), mean(y2_mc))
## [1] 1.95806 1.51564
# probabilidades de la distribucion predictiva
tab <- rbind(dnbinom(x = 0:7, size = a+s1, mu = (a+s1)/(b+n1)), table(y1_mc)[1:8]/B)
colnames(tab) <- 0:7
rownames(tab) <- c("Directo", "Simulacion")
round(tab, 3)
##                0     1     2     3     4     5     6     7
## Directo    0.143 0.277 0.269 0.176 0.086 0.034 0.011 0.003
## Simulacion 0.142 0.275 0.270 0.177 0.086 0.034 0.011 0.003
tab <- rbind(dnbinom(x = 0:7, size = a+s2, mu = (a+s2)/(b+n2)), table(y2_mc)[1:8]/B)
colnames(tab) <- 0:7
rownames(tab) <- c("Directo", "Simulacion")
round(tab, 3)
##                0     1     2     3     4     5     6     7
## Directo    0.224 0.332 0.249 0.126 0.049 0.015 0.004 0.001
## Simulacion 0.224 0.329 0.251 0.126 0.049 0.015 0.004 0.001
# probabilidad posterior de que y*_1 > y*_2
mean(y1_mc > y2_mc)
## [1] 0.48374
# distribucion predictiva posterior de d = y*_1 - y*_2
par(mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
plot(table(y1_mc - y2_mc)/B, type = "h", cex.axis = 0.9, xlab = "d", 
     ylab = "p(d | y )", col = 4, lwd = 3)

6.4 Chequeo del modelo

# estadistico observado
t_obs <- s1/n1
print(t_obs)
## [1] 1.954955
# distribucion predictiva posterior
t_mc <- NULL
set.seed(1234)
for(i in 1:B) {
        y_rep   <- rpois(n = n1, lambda = theta1_mc[i])
        t_mc[i] <- mean(y_rep)
}
# grafico
par(mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
hist(x = t_mc, freq = F, col = "gray90", border = "gray90", xlab = "t", ylab = "p(t | y)", main = "")
lines(density(t_mc), col = 4, lwd = 2)
abline(v = t_obs, col = 1, lwd = 2, lty = 1)
abline(v = quantile(x = t_mc, probs = c(0.025, 0.975)), lty = 2, lwd = 2, col = 2)
legend("topright", legend = c("Posterior", "IC 95%", "t obs"), col = c(4, 2, 1), 
       lty = 1, lwd = 2, bty = "n")

# ppp (valor p predictivo posterior)
mean(t_mc > t_obs)
## [1] 0.48274