1. Modelo

Sea un conjunto de datos: \[ y_1, \dots, y_n \sim \mathcal{N}(\mu, \sigma^2) \]

Priors

Usamos un prior conjugado:

Media: \[ \mu \mid \sigma^2 \sim \mathcal{N}(\mu_0, \sigma^2 / \nu_0) \]
Varianza: \[ \sigma^2 \sim \text{Inv-Gamma}(\alpha_0, \beta_0) \]

Donde elegimos hiperparámetros pequeños (por ejemplo cercanos a 2) para mantener un prior poco informativo pero propio: - \(\mu_0 = 0\) - \(\nu_0 = 2\) - \(\alpha_0 = 2\) - \(\beta_0 = 2\)

2. Posterior conjunta

La posterior conjunta es:

\[ p(\mu, \sigma^2 \mid y) \propto p(y \mid \mu, \sigma^2)\, p(\mu \mid \sigma^2)\, p(\sigma^2) \]

3. Marginales completas

3.1 Para \(\mu \mid \sigma^2, y\)

\[ \mu \mid \sigma^2, y \sim \mathcal{N}(\mu_n, \sigma^2 / \nu_n) \]

donde:

\[ \nu_n = \nu_0 + n \]

\[ \mu_n = \frac{\nu_0 \mu_0 + n \bar{y}}{\nu_n} \]

3.2 Para \(\sigma^2 \mid \mu, y\)

\[ \sigma^2 \mid \mu, y \sim \text{Inv-Gamma}(\alpha_n, \beta_n) \]

donde:

\[ \alpha_n = \alpha_0 + \frac{n}{2} \]

\[ \beta_n = \beta_0 + \frac{1}{2} \sum (y_i - \mu)^2 \]

4. Algoritmo de Gibbs

Iterar:

Muestrear \(\mu^{(t)} \sim p(\mu \mid \sigma^{2(t-1)}, y)\)
Muestrear \(\sigma^{2(t)} \sim p(\sigma^2 \mid \mu^{(t)}, y)\)

5. Implementación en R (sin paquetes)

set.seed(123)

# Datos simulados
n <- 50
y <- rnorm(n, mean = 5, sd = 2)

# Hiperparámetros
mu0 <- 0
nu0 <- 2
alpha0 <- 2
beta0 <- 2

# MCMC
B <- 1000
burn_in <- 100

mu <- numeric(B)
sigma2 <- numeric(B)

# valores iniciales
mu[1] <- mean(y)
sigma2[1] <- var(y)

for (t in 2:B) {
  
  # ---- actualizar mu ----
  nu_n <- nu0 + n
  mu_n <- (nu0 * mu0 + n * mean(y)) / nu_n
  var_mu <- sigma2[t-1] / nu_n
  
  mu[t] <- rnorm(1, mean = mu_n, sd = sqrt(var_mu))
  
  # ---- actualizar sigma^2 ----
  alpha_n <- alpha0 + n / 2
  beta_n <- beta0 + 0.5 * sum((y - mu[t])^2)
  
  # muestreo Inv-Gamma usando Gamma
  sigma2[t] <- 1 / rgamma(1, shape = alpha_n, rate = beta_n)
}

# eliminar burn-in
mu_post <- mu[(burn_in+1):B]
sigma2_post <- sigma2[(burn_in+1):B]

# resultados
mean(mu_post)
mean(sigma2_post)

# gráficos
par(mfrow=c(2,2))

plot(mu, type="l", main="Trace de mu")
hist(mu_post, main="Posterior de mu")

plot(sigma2, type="l", main="Trace de sigma2")
hist(sigma2_post, main="Posterior de sigma2")

# Priors no informativas
n0<-1
nu0<-1
t0<-100
m0<-0
s0<-100
mp<-NULL
vp<-NULL
n<-60
y<-rnorm(60,4,3)
mp[1]<-mean(y)
vp[1]<-var(y)
B<-5000
for(i in 2:B){
  mp[i]<-rnorm(1,(n0*m0/t0+n*mean(y)/vp[i-1])/(n0/t0+n/vp[i-1]))
  vp[i]<-1/rgamma(1,(n+nu0)/2,(n*var(y)+nu0*s0)/2)
}
plot(mp, type='l',col = "green")

plot(vp, type='l',col = "green")

plot(density(mp),main = 'Media posterior')
abline(v=quantile(mp,0.025),col=2)
abline(v=quantile(mp,0.5),col=5)
abline(v=quantile(mp,0.975),col=2)

plot(density(vp),main = 'Varianza posterior')
abline(v=quantile(vp,0.025),col=2)
abline(v=quantile(vp,0.5),col=5)
abline(v=quantile(vp,0.975),col=2)

# Son las muestras posteriores i.i.d.

acf(mp)

pacf(mp)

acf(vp)

pacf(vp)

TENEMOS BUENAS MUESTRAS POSTERIORES.

Modelo Normal Bayesiano con MCMC (Gibbs)

Ejemplo