Modelo nomral multivairado bayesiano
# Datos: puntajes de comprensión de lectura antes y después de la instrucción
Y <- matrix(
c(
59, 43, 34, 32, 42, 38, 55, 67, 64, 45, 49, 72, 34, 70, 34, 50, 41, 52, 60, 34, 28, 35,
77, 39, 46, 26, 38, 43, 68, 86, 77, 60, 50, 59, 38, 48, 55, 58, 54, 60, 75, 47, 48, 33
),
nrow = 22, ncol = 2,
dimnames = list(NULL, c("pre_test", "post_test"))
)
# Dimensiones
(n <- nrow(Y))
## [1] 22
(p <- ncol(Y))
## [1] 2
# Inspeccionar datos
summary(Y)
## pre_test post_test
## Min. :28.00 Min. :26.00
## 1st Qu.:34.25 1st Qu.:43.75
## Median :44.00 Median :52.00
## Mean :47.18 Mean :53.86
## 3rd Qu.:58.00 3rd Qu.:60.00
## Max. :72.00 Max. :86.00
df<-data.frame(y1=Y[,1],y2=Y[,2])
plot(Y[,1],Y[,2])
# EstadÃsticos suficientes
yb <- c(colMeans(Y))
round(yb, 1)
## pre_test post_test
## 47.2 53.9
SS <- cov(Y)
round(SS, 1)
## pre_test post_test
## pre_test 182.2 148.4
## post_test 148.4 243.6
# Previa
mu0 <- c(50,50)
L0 <- matrix(data = c(278,139,139,278), nrow = 2, ncol = 2)
nu0 <- 4
S0 <- matrix(data = c(278,139,139,278), nrow = 2, ncol = 2)
# Inicialización
# Solve ya que la previa conjugada es Wishart inversa
Sigma <- solve(rWishart::rWishart(n = 1, df = nu0, Sigma = S0)[,,1])
# Número de muestras
B <- 10000
# Almacenamiento
THETA <- NULL
SIGMA <- NULL
YS <- NULL
LL <- numeric(B)
# Cálculo de cantidades constantes previas l procedimiento iterativo
iL0 <- solve(L0)
Lm0 <- iL0 %*% mu0
nun <- nu0 + n
SSn <- S0 + (n - 1) * SS
# Muestreador de Gibbs
set.seed(1234)
for (b in 1:B) {
# Actualización de theta
iSigma <- solve(Sigma)
Ln <- solve(iL0 + n * iSigma)
theta <- c(mvtnorm::rmvnorm(n = 1, mean = Ln %*% (Lm0 + n * (iSigma %*% yb)), sigma = Ln))
# Actualización de Sigma
Sigma <- MCMCpack::riwish(v = nun, S = SSn + n * ((yb - theta) %*% t(yb - theta)))
# Muestra de la predictiva posterior
YS <- rbind(YS, mvtnorm::rmvnorm(n = 1, mean = theta, sigma = Sigma))
# Cálculo de log-verosimilitud
LL[b] <- sum(apply(X = Y, MARGIN = 1, FUN = function(y) mvtnorm::dmvnorm(y, mean = theta, sigma = Sigma, log = T)))
# Almacenamiento de resultados
THETA <- rbind(THETA, theta)
SIGMA <- rbind(SIGMA, c(Sigma))
# Mostrar progreso cada 10%
if (b %% floor(B / 10) == 0)
cat(100 * round(b / B, 1), "% completado...\n", sep = "")
}
## 10% completado...
## 20% completado...
## 30% completado...
## 40% completado...
## 50% completado...
## 60% completado...
## 70% completado...
## 80% completado...
## 90% completado...
## 100% completado...
# Etiquetado de columnas para claridad
colnames(THETA) <- c("theta1", "theta2")
colnames(SIGMA) <- c("sigma1^2", "sigma21", "sigma12", "sigma2^2")
# Tamaño efectivo de la muestra
neff <- coda::effectiveSize(cbind(THETA, SIGMA))
round(neff)
## theta1 theta2 sigma1^2 sigma21 sigma12 sigma2^2
## 10000 10000 9059 9013 9013 8951
# Error estándar de Monte Carlo
EMC <- apply(X = cbind(THETA, SIGMA), MARGIN = 2, FUN = sd)/sqrt(neff)
round(EMC, 4)
## theta1 theta2 sigma1^2 sigma21 sigma12 sigma2^2
## 0.0284 0.0322 0.6120 0.6007 0.6007 0.8182
# Coeficiente de variación de Monte Carlo (%)
CVMC <- 100*EMC/abs(colMeans(cbind(THETA, SIGMA)))
round(CVMC, 4)
## theta1 theta2 sigma1^2 sigma21 sigma12 sigma2^2
## 0.0603 0.0601 0.3295 0.4076 0.4076 0.3340
plot(THETA[,1],main="theta 1",type="l")
plot(density(THETA[,1]),main="theta 1")
qqnorm(THETA[,1])
acf(THETA[,1])
plot(THETA[,2],main="theta 2",type="l")
plot(density(THETA[,2]),main="theta 2")
qqnorm(THETA[,2])
acf(THETA[,2])
plot(SIGMA[,1],main="sigma1^2",type="l")
plot(density(SIGMA[,1]),main="SIGMA$sigma1^2")
qqnorm(sqrt(SIGMA[,1])*sample(c(-1,1),10000,replace = TRUE))
acf(sqrt(SIGMA[,1])*sample(c(-1,1),10000,replace = TRUE))
colnames(SIGMA)
## [1] "sigma1^2" "sigma21" "sigma12" "sigma2^2"