Processing math: 50%
  1. Los archivos school1.dat, school2.dat, y shool3.dat contienen datos sobre la cantidad de tiempo que los estudiantes de tres colegios dedicaron a estudiar o hacer tareas durante un período de exámenes.
  1. Explore los datos gráfica y numéricamente.

  2. Analice los datos de cada una de los colegios separadamente, utilizando un modelo Normal con una distribución previa conjugada, en la que μ0=5, σ20=4, κ0=1, ν0=2, y calcule lo siguiente:

    • Medias posteriores e intervalos de credibilidad al 95% para la media θ, la desviación estándar σ, y el coeficiente de variación σμ de cada escuela.
    • La probabilidad posterior de que θi<θj<θk para las seis permutaciones {i,j,k} de {1,2,3}, donde θi es la media del del colegio i.
    • La probabilidad posterior de que ˜yi<˜yj<˜yk para las seis permutaciones {i,j,k} de {1,2,3}, donde ˜yi es una observación de la distribución predictiva posterior de la escuela i.
    • Calcule la probabilidad posterior de que θ1 sea mayor que θ2 y θ3, y la probabilidad posterior de que ˜y1 sea mayor que ˜y2 y ˜y3.
  3. Dibuje la distribución posterior conjunta de (θ,σ2) para cada escuela.

  4. Compruebe la bondad de ajuste del modelo para cada escuela utilizando como estadísticos de prueba la media y el coeficiente de variación.

  1. Considere el modelo Normal dado por yiθ,σ2iidN(θ,σ2), para i=1,,n, con distribución previa θσ2N(θ0,κ0σ2)yσ2GI(a,b), donde θ0,κ0,a,b son los hiperparámetros del modelo.
  1. Encuentre la distribución posterior de (θ,σ2).
  2. Encuentre la distribución condicional completa de θ.
  3. Encuentre la distribución marginal de θ.
  4. Encuentre la distribución marginal de σ2.
  5. Simule n=1000 observaciones i.i.d de N(5,1). Ajuste el modelo suponiendo los siguientes escenarios previos: i. estados de conocimiento previo bastante informativo acerca de los valores reales de los parámetros, ii. estado de conocimiento previo informativo acerca de θ y difuso acerca de σ2, iii. estado de conocimiento previo informativo acerca de σ2 y difuso acerca de θ, y iv. estado de conocimiento previo difuso acerca de ambos parámetros. Caracterice la distribución posterior en cada caso.
  6. Suponga que está interesado en hacer inferencia sobre η=σ/θ. Desarrolle un algoritmo de Monte Carlo para calcular la media posterior y un intervalo de credibilidad al 95% para η. Use el algoritmo para calcular estas cantidades en todos los escenarios descritos anteriormente.
  1. Considere el modelo Normal xiθ,σ2iidN(θ,σ2), para i=1,,n, donde θ es desconocido y σ2 es conocido. Además, considere una distribución previa para θ definida por medio de una mezcla finita de previas conjugadas de la forma p(θ)=K=1wϕ(θμ,τ2) donde K es un entero positivo fijo mayor o igual que 1, w1,,wK es un sistema de pesos tales que y K=1w=1 y 0w1 para =1,,K, y ϕ(θμ,τ2) denota la densidad de la distribución Normal con media μ y varianza τ2. Una distribución previa de esta forma permite especificar estados de información previos multimodales acerca de θ.
  1. Encuentre la distribución posterior de θ.
  2. Encuentre la media posterior de θ.
  3. Encuentre la distribución predictiva previa.
  4. Encuentre la distribución predictiva posterior.
  1. Considere el modelo Normal Trucado xiθ,σ2iidN(0,)(θ,σ2), para i=1,,n, donde σ2=1. Además, considere la distribución previa para θ dada por θN(μ,τ2).
  1. Encuentre la distribución posterior de θ.
  2. ¿Este modelo se puede catalogar como un modelo conjugado?

Solución

  1. Los archivos school1.dat, school2.dat, y shool3.dat contienen datos sobre la cantidad de tiempo que los estudiantes de tres colegios dedicaron a estudiar o hacer tareas durante un período de exámenes.
  1. Explore los datos gráfica y numéricamente.
# datos
school1 <- read.table("school1.dat", quote="\"", comment.char="")
school2 <- read.table("school2.dat", quote="\"", comment.char="")
school3 <- read.table("school3.dat", quote="\"", comment.char="")
x1 <- school1$V1
x2 <- school2$V1
x3 <- school3$V1
# EDA (Exploratory Data Analysis)
eda <- function(x) c(length(x), min(x), max(x), mean(x), quantile(x, c(.25,.5,.75)), sd(x), 100*sd(x)/mean(x))
tab <- round(rbind(eda(x1), eda(x2), eda(x3)), 2)
colnames(tab) <- c("n","Min","Max","Media","25%","50%","75%","DE","CV(%)")
rownames(tab) <- paste("School", 1:3, sep = " ")
knitr::kable(x = tab, digits = 2)
n Min Max Media 25% 50% 75% DE CV(%)
School 1 25 1.55 16.38 9.46 7.44 9.61 11.65 3.89 41.05
School 2 23 0.29 15.35 7.03 3.41 6.54 10.07 4.48 63.76
School 3 20 0.95 13.63 7.95 5.13 6.87 11.51 3.78 47.55
par(mfrow=c(1,1), mar=c(3,3,1,1), mgp=c(1.75,.75,0))
plot(NA,NA, xlim = c(0,4), ylim = c(0.25,17.0), xlab = "School", ylab = "Tiempo (horas)", main = "", xaxt = "n", yaxt = "n")
axis(side = 1, at = 1:3, labels = 1:3)
for (j in 1:3) {
  x <- get(x = paste0("x", j))
  boxplot(x, at = j, col = "gray90", border = "gray40", boxwex = 0.6, add = T)
  myjitter <- jitter(rep(j, length(x)), amount = 0.15)
  points(x = myjitter, y = x, pch = 20, col = rgb(0,0,0,.6)) 
}

El análisis exploratorio no revela desviaciones importantes de un comportamiento simétrico, ni observaciones atípicas.

  1. Analice los datos de cada una de los colegios separadamente, utilizando un modelo Normal con una distribución previa conjugada, en la que μ0=5, σ20=4, κ0=1, ν0=2, y calcule lo siguiente:

Se ajusta el modelo con distribución muestral yiθ,σ2iidN(θ,σ2),i=1,,n, y distribución previa conjugada θσ2N(μ0,σ2κ0)σ2GI(ν02,ν0σ202) donde μ0,κ0,ν0,σ20 son los hiperparámetros del modelo.

# hiperparametros
mu0 <- 5  
k0  <- 1
s20 <- 4 
nu0 <- 2
# n. de simulaciones
S <- 1000  
# distribucion posterior
THETA <- matrix(NA, S, 3)
SIG2  <- matrix(NA, S, 3)
set.seed(1)
for (j in 1:3) {
  # datos
  x    <- get(x = paste0("x", j))
  n    <- length(x)
  # estadisticos
  ybar <- mean(x)
  s2   <- var(x)
  # parametros posterior
  kn  <- k0 + n
  mun <- (k0*mu0 + n*ybar)/kn
  nun <- nu0 + n
  s2n <- (nu0*s20 +(n-1)*s2 + k0*n*(ybar-mu0)^2/kn)/nun
  # simulacion
  SIG2 [,j] <- 1/rgamma(n = S, shape = nun/2, rate = nun*s2n/2)
  THETA[,j] <- rnorm(n = S, mean = mun, sd = sqrt(SIG2[,j]/kn))
}
for (j in 1:3) {
  tab <- rbind(
  c(mean(THETA[,j]), quantile(THETA[,j], c(0.025,0.975))),
  c(mean(sqrt(SIG2)[,j]), quantile(sqrt(SIG2[,j]), c(0.025,0.975))),
  c(mean(sqrt(SIG2[,j])/THETA[,j]), quantile(sqrt(SIG2[,j])/THETA[,j], c(0.025,0.975))))
  colnames(tab) <- c("Media","Q2.5%","Q97.5%")
  rownames(tab) <- c("theta","sigma","cv")
  assign(x = paste0("tab",j), value = tab)
}

School 1:

knitr::kable(x = tab1, digits = 2, align = "c")
Media Q2.5% Q97.5%
theta 9.28 7.64 10.84
sigma 3.93 3.02 5.12
cv 0.43 0.31 0.59

School 2:

knitr::kable(x = tab2, digits = 2, align = "c")
Media Q2.5% Q97.5%
theta 6.91 5.09 8.64
sigma 4.34 3.35 5.84
cv 0.64 0.45 0.99

School 3:

knitr::kable(x = tab3, digits = 2, align = "c")
Media Q2.5% Q97.5%
theta 7.83 6.32 9.47
sigma 3.73 2.83 5.15
cv 0.48 0.34 0.71

La situación más probable es θ2<θ3<θ1.

# permutaciones
per <- gtools::permutations(n = 3, r = 3)
# probabilidades
tab <- NULL
for (i in 1:nrow(per)) {
  tab <- rbind(tab, mean((THETA[,per[i,1]] < THETA[,per[i,2]]) & (THETA[,per[i,2]] < THETA[,per[i,3]])))
}
tab <- cbind(per, tab)
colnames(tab) <- c("i","j","k","P(theta i < theta j < theta k")
knitr::kable(x = tab, digits = 3, align = "c")
i j k P(theta i < theta j < theta k
1 2 3 0.008
1 3 2 0.004
2 1 3 0.087
2 3 1 0.679
3 1 2 0.013
3 2 1 0.209

La situación más probable es nuevamente ˜y2<˜y3<˜y1, aunque las probabilidades no son tan disimiles con los parámetros de los modelos.

# distribucion predictiva posterior
YP <- matrix(NA, S, 3)
set.seed(1) 
for (j in 1:3) {
  YP[,j] <- rnorm(n = S, mean = THETA[,j], sd = sqrt(SIG2[,j]))
}
# probabilidades
tab <- NULL
for (i in 1:nrow(per)) {
  tab <- rbind(tab, mean((YP[,per[i,1]] < YP[,per[i,2]]) & (YP[,per[i,2]] < YP[,per[i,3]])))
}
tab <- cbind(per, tab)
colnames(tab) <- c("i","j","k","P(y* i < y* j < y* k)")
knitr::kable(x = tab, digits = 3, align = "c")
i j k P(y* i < y* j < y* k)
1 2 3 0.108
1 3 2 0.114
2 1 3 0.186
2 3 1 0.274
3 1 2 0.133
3 2 1 0.185

La probabilidad posterior de que θ1 sea mayor que θ2 y θ3 es 0.888.

mean((THETA[,1] > THETA[,2]) & (THETA[,1] > THETA[,3]))
## [1] 0.888

La probabilidad posterior de que ˜y1 sea mayor que ˜y2 y ˜y3 es 0.459.

mean((YP[,1] > YP[,2]) & (YP[,1] > YP[,3]))
## [1] 0.459
  1. Dibuje la distribución posterior conjunta de (θ,σ2) para cada escuela.
par(mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(1.75,.75,0))
for (j in 1:3) {
  plot(THETA[,j], SIG2[,j], xlab = expression(theta), ylab = expression(sigma^2), main = paste0("School ", j), pch = 16, cex = 0.8, col = j, xlim = c(3,12), ylim = c(5,45))
}

  1. Compruebe la bondad de ajuste del modelo para cada escuela utilizando como estadísticos de prueba la media y el coeficiente de variación.
par(mfrow=c(3,2), mar=c(3,3,1,1), mgp=c(1.75,.75,0))
for (j in 1:3) {
  # datos
  x <- get(x = paste0("x", j))
  n <- length(x)
  # estadisticos observados
  EP0 <- c(mean(x), sd(x)/mean(x))
  # distribucion predictiva
  EP <- NULL
  set.seed(1)
  for (b in 1:S) {
    xx <- rnorm(n = n, mean = THETA[b,j], sd = sqrt(SIG2[b,j]))
    EP <- rbind(EP, c(mean(xx), sd(xx)/mean(xx)))
  }
  # graficos
  display <- c("Media", "CV")
  for (i in 1:2) {
    ppp <- round(mean(EP[,i] > EP0[i]), 2)
    hist(x = EP[,i], freq = F, col = "gray90", border = "gray90", xlab = display[i], ylab = "Densidad", main = paste0("School ",j,", ", display[i],", ppp = ", ppp))
    abline(v = EP0[i], col = 2)
    abline(v = quantile(EP[,i], c(.025,.975)), col = 1, lty = 2)
  }
}

  1. Considere el modelo Normal dado por yiθ,σ2iidN(θ,σ2), para i=1,,n, con distribución previa θσ2N(θ0,κ0σ2)yσ2GI(a,b), donde θ0,κ0,a,b son los hiperparámetros del modelo.
  1. Encuentre la distribución posterior de (θ,σ2).
  2. Encuentre la distribución condicional completa de θ.
  3. Encuentre la distribución marginal de θ.
  4. Encuentre la distribución marginal de σ2.

Tenemos que la distribución muestral es \begin{align*} f(\boldsymbol{x}|\,\boldsymbol{\theta}) &= \prod_{i=1}^n f(x_i|\,\boldsymbol{\theta})\\ &=\prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma}\,\exp{\left\{ -\frac{1}{2\sigma^2}(x_i - \theta)^2 \right\}}\\ &=\frac{1}{(2\pi)^{n/2}\sigma^n}\,\exp{\left\{ -\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i - \theta)^2 \right\}}\\ &=\frac{1}{(2\pi)^{n/2}\sigma^n}\,\exp{\left\{ -\frac{1}{2\sigma^2}\sum_{i=1}^n (x_i^2 - 2\theta x_i + \theta^2)\right\}}\\ &=\frac{1}{(2\pi)^{n/2}\sigma^n}\,\exp{\left\{ -\frac{1}{2\sigma^2}\left[\sum_{i=1}^n x_i^2 - 2n\bar{x}\theta+ n\theta^2\right] \right\}}\\ &=\frac{1}{(2\pi)^{n/2}\sigma^n}\,\exp{\left\{ -\frac{1}{2}\left[\frac{\sum x_i^2}{\sigma^2} - 2\theta\frac{n\bar{x}}{\sigma^2} + \frac{n\theta^2}{\sigma^2}\right] \right\}}\\ &=\frac{1}{(2\pi)^{n/2}\sigma^n}\,\exp{\left\{ -\frac{1}{2} \left[\theta^2\frac{n}{\sigma^2} - 2\theta\frac{n\bar{x}}{\sigma^2} \right] \right\}}\, \exp{\left\{ -\frac{1}{\sigma^2} \frac{\sum x_i^2}{2} \right\}} \end{align*} y la distribución previa es \begin{align*} \pi(\boldsymbol{\theta})&=\pi_{\theta|\sigma^2}(\theta|\sigma^2)\,\pi_{\sigma^2}(\sigma^2)\\ &= \frac{1}{\sqrt{2\pi}\,\sqrt{\kappa_0}\,\sigma}\,\exp{\left\{ -\frac{1}{2\kappa_0\sigma^2}(\theta-\theta_0)^2 \right\}} \cdot \frac{b^2}{\Gamma(a)}\,(\sigma^2)^{-(a+1)} \exp{ \left\{ -\frac{b}{\sigma^2} \right\} } \\ &= \frac{b^2}{\sqrt{2\pi}\,\sqrt{\kappa_0}\,\sigma^{2a+3}\,\Gamma(a)}\, \exp{\left\{ -\frac{1}{2}\left[\theta^2\frac{1}{\kappa_0\sigma^2} - 2\theta\frac{\theta_0}{\kappa_0\sigma^2}+\frac{\theta_0^2}{\kappa_0\sigma^2}\right] \right\}}\, \exp{ \left\{ -\frac{b}{\sigma^2} \right\} } \\ &= \frac{b^2}{\sqrt{2\pi}\,\sqrt{\kappa_0}\,\sigma^{2a+3}\,\Gamma(a)}\, \exp{\left\{ -\frac{1}{2}\left[\theta^2\frac{1}{\kappa_0\sigma^2} - 2\theta\frac{\theta_0}{\kappa_0\sigma^2}\right] \right\}}\, \exp{ \left\{ -\frac{1}{\sigma^2 }\left[b + \frac{\theta_0^2}{2\kappa_0} \right] \right\} } \end{align*} y por lo tanto, la distribución posterior es \begin{align*} \pi(\boldsymbol{\theta}|\,\boldsymbol{x})&\propto \frac{1}{\sigma^n}\,\exp{\left\{ -\frac{1}{2} \left[\theta^2\frac{n}{\sigma^2} - 2\theta\frac{n\bar{x}}{\sigma^2} \right] \right\}}\, \exp{\left\{ -\frac{1}{\sigma^2} \frac{\sum x_i^2}{2} \right\}}\\ &\hspace{2cm}\cdot\frac{1}{\sigma^{2a+3}}\, \exp{\left\{ -\frac{1}{2}\left[\theta^2\frac{1}{\kappa_0\sigma^2} - 2\theta\frac{\theta_0}{\kappa_0\sigma^2}\right] \right\}}\, \exp{ \left\{ -\frac{1}{\sigma^2 }\left[b + \frac{\theta_0^2}{2\kappa_0} \right] \right\} } \\ &\propto \frac{1}{\sigma^{2a+n+3}}\,\exp{\left\{ -\frac{1}{2} \left[\theta^2\left(\frac{n}{\sigma^2} + \frac{1}{\kappa_0\sigma^2}\right) - 2\theta\left(\frac{n\bar{x}}{\sigma^2} + \frac{\theta_0}{\kappa_0\sigma^2}\right) \right] \right\}}\\ &\hspace{2cm}\cdot\exp{\left\{ -\frac{1}{\sigma^2} \left[\frac{\sum x_i^2}{2} + b + \frac{\theta_0^2}{2\kappa_0} \right] \right\}}\\ &\propto \frac{1}{\sigma^{2a+n+3}}\,\exp{\left\{ -\frac{1}{2} \left[\frac{\theta^2 }{\tau^2} - \frac{2\theta\mu}{\tau^2} \right] \right\}}\, \exp{\left\{ -\frac{B}{\sigma^2} \right\}} \end{align*} donde \begin{align*} &\frac{1}{\tau^2}\equiv \frac{n}{\sigma^2} + \frac{1}{\kappa_0\sigma^2}= \frac{n\kappa_0+1}{\kappa_0\sigma^2}& &\Rightarrow\quad \boxed{ \tau^2= \frac{\kappa_0\sigma^2}{n\kappa_0+1} }\\ &\frac{\mu}{\tau^2}\equiv \frac{n\bar{x}}{\sigma^2} + \frac{\theta_0}{\kappa_0\sigma^2}= \frac{n\bar{x}\kappa_0+\theta_0}{\kappa_0\sigma^2}& &\Rightarrow\quad \boxed{ \mu=\frac{n\bar{x}\kappa_0+\theta_0}{n\kappa_0+1} }\\ &B\equiv \frac{\sum x_i^2}{2} + b + \frac{\theta_0^2}{2\kappa_0} & &\Rightarrow\quad \boxed{ B= \frac{\kappa_0\sum x_i^2 + 2b\kappa_0 + \theta_0^2}{2\kappa_0} } \end{align*} Así, completando el cuadrado y reorganizando términos, se sigue que \begin{align*} \pi(\boldsymbol{\theta}|\,\boldsymbol{x}) &\propto \frac{1}{\sigma^{2a+n+3}}\,\exp{\left\{ -\frac{1}{2\tau^2} \left(\theta- \mu \right)^2 \right\}}\, \exp{\left\{ -\frac{B}{\sigma^2} + \frac{\mu^2}{2\tau^2} \right\}}\\ &\propto\frac{1}{\tau}\,\exp{\left\{ -\frac{1}{2\tau^2} \left(\theta- \mu \right)^2 \right\}}\,\frac{\tau}{\sigma^{2a+n+3}}\, \exp{\left\{ -\frac{1}{\sigma^2}\left[B - \frac{\mu^2}{2}\frac{n\kappa_0+1}{\kappa_0}\right] \right\}}\\ &\propto\frac{1}{\tau}\,\exp{\left\{ -\frac{1}{2\tau^2} \left(\theta- \mu \right)^2 \right\}}\,\frac{\sqrt{\kappa_0}\sigma}{\sqrt{ n\kappa_0+1 }}\,\frac{1}{\sigma^{2a+n+3}}\, \exp{\left\{ -\frac{1}{\sigma^2}\left[B - \frac{\mu^2}{2}\frac{n\kappa_0+1}{\kappa_0}\right] \right\}}\\ &\propto\,\exp{\left\{ -\frac{1}{2\tau^2} \left(\theta- \mu \right)^2 \right\}}\, \left( \sigma^2 \right)^{-(\alpha+1)}\, \exp{\left\{ -\frac{\beta}{\sigma^2} \right\}} \end{align*} donde \begin{align*} &\quad\boxed{\alpha\equiv a+\frac{n}2 }\\ \beta\equiv B - \frac{\mu^2}{2}\frac{n\kappa_0+1}{\kappa_0}\quad\Rightarrow&\quad \boxed{ \beta=B - \frac{(n\bar{x}\kappa_0 + \theta_0)^2}{2\kappa_0(n\kappa_0 + 1)} } \end{align*} Normalizando, se tiene que \boxed{ \pi(\boldsymbol{\theta}|\,\boldsymbol{x}) =\frac{1}{\sqrt{2\pi}\tau}\,\exp{\left\{ -\frac{1}{2\tau^2} \left(\theta- \mu \right)^2 \right\}}\, \frac{\beta^\alpha}{\Gamma(\alpha)}\,\left( \sigma^2 \right)^{-(\alpha+1)}\, \exp{\left\{ -\frac{\beta}{\sigma^2} \right\}} } lo que significa que p(\boldsymbol{\theta}|\,\boldsymbol{x})=\pi_1(\theta|\,\sigma^2, \boldsymbol{x})\,\pi_2(\sigma^2|\,\boldsymbol{x}) con \pi_1(\theta|\,\sigma^2, \boldsymbol{x})=N(\mu,\tau^2) y \pi_2(\sigma^2|\,\boldsymbol{x})=IG(\alpha,\beta).

Teniendo en cuenta esta distribución posterior, se tiene directamente que \begin{align*} &p(\sigma^2|\, \boldsymbol{x})= \int p(\boldsymbol{\theta}|\,\boldsymbol{x})\,d\theta=\int \pi_1(\theta|\,\sigma^2, \boldsymbol{x})\,\pi_2(\sigma^2|\,\boldsymbol{x})\,d\theta=\pi_2(\sigma^2|\,\boldsymbol{x})\\ &\Rightarrow\quad \boxed{p(\sigma^2|\, \boldsymbol{x})=\pi_2(\sigma^2|\,\boldsymbol{x})}\quad\therefore\quad\boxed{(\sigma^2|\, \boldsymbol{x}) \sim IG(\alpha,\beta)} \end{align*} Además, \begin{align*} &p(\theta|\,\sigma^2, \boldsymbol{x})= \frac{p(\theta,\sigma^2|\,\boldsymbol{x})}{p(\sigma^2|\, \boldsymbol{x})}=\frac{\pi_1(\theta|\,\sigma^2, \boldsymbol{x})\,\pi_2(\sigma^2|\,\boldsymbol{x})}{\pi_2(\sigma^2|\,\boldsymbol{x})}=\pi_1(\theta|\,\sigma^2, \boldsymbol{x})\\ &\Rightarrow\quad \boxed{p(\theta|\,\sigma^2, \boldsymbol{x})=\pi_1(\theta|\,\sigma^2, \boldsymbol{x})}\quad\therefore\quad\boxed{(\theta|\,\sigma^2, \boldsymbol{x}) \sim N\left(\mu,\tau^2\right)} \end{align*} Integrando la distribución posterior respecto a \sigma^2, se sigue que \begin{align*} p(\theta|\,\boldsymbol{x})&=\int p (\boldsymbol{\theta}|\,\boldsymbol{x})\,d\sigma^2=\int \pi_1(\theta|\,\sigma^2, \boldsymbol{x})\,\pi_2(\sigma^2|\,\boldsymbol{x}) \,d\sigma^2\\ &= \int \frac{1}{\sqrt{2\pi}\tau}\,\exp{\left\{ -\frac{1}{2\tau^2} \left(\theta- \mu \right)^2 \right\}}\, \frac{\beta^\alpha}{\Gamma(\alpha)}\,\left( \sigma^2 \right)^{-(\alpha+1)}\, \exp{\left\{ -\frac{\beta}{\sigma^2} \right\}}\,d\sigma^2\\ &= \frac{\beta^\alpha}{\sqrt{2\pi}\,\Gamma(\alpha)}\int \frac{\left( \sigma^2 \right)^{-(\alpha+1)}}{\tau}\, \exp{\left\{ -\frac{1}{2\tau^2} \left(\theta- \mu \right)^2 \right\}}\, \exp{\left\{ -\frac{\beta}{\sigma^2} \right\}}\,d\sigma^2\\ &= \frac{c_1}{\sqrt{c_2}} \int \frac{\left( \sigma^2 \right)^{-(\alpha+1)}}{\sigma}\, \exp{\left\{ -\frac{1}{2\,c_2\,\sigma^2} \left(\theta- \mu \right)^2 \right\}}\, \exp{\left\{ -\frac{\beta}{\sigma^2} \right\}}\,d\sigma^2 \end{align*} donde \begin{align*} c_1\equiv\frac{\beta^\alpha}{\sqrt{2\pi}\,\Gamma(\alpha)}\quad\text{and}\quad &c_2\equiv\frac{\kappa_0}{n\kappa_0+1} \end{align*} Así, reconociendo términos, se sigue que \begin{align*} p(\theta|\,\boldsymbol{x}) &= \frac{c_1}{\sqrt{c_2}} \int \sigma^{-(2\alpha+3)}\, \exp{\left\{ -\frac{1}{\sigma^2}\left[ \frac{\left(\theta- \mu \right)^2}{2\,c_2} + \beta\right] \right\}}\,d\sigma^2\\ &= \frac{c_1}{\sqrt{c_2}}\, \frac{\Gamma(\lambda)}{\gamma^\lambda} \, \int \frac{\gamma^\lambda}{\Gamma(\lambda)}\, \left(\sigma^2\right)^{-(\lambda+1)}\, \exp{\left\{ -\frac{\gamma}{\sigma^2}\right\}}\,d\sigma^2\\ &= \frac{c_1}{\sqrt{c_2}}\, \frac{\Gamma(\lambda)}{\gamma^\lambda} \end{align*} con \boxed{\gamma\equiv\frac{\left(\theta- \mu \right)^2}{2\,c_2} + \beta} \quad\text{and}\quad \boxed{ \lambda\equiv \alpha+\frac12} Entonces, reconociendo términos, se sigue que \begin{align*} p(\theta|\,\boldsymbol{x}) &= \frac{1}{\sqrt{2\pi c_2}}\,\,\frac{\Gamma(\alpha+1/2)}{\Gamma(\alpha)}\,\, \frac{\beta^\alpha}{\gamma^{\alpha+1/2}}\\ &= \frac{1}{\sqrt{2\pi c_2}}\,\,\frac{\Gamma(\alpha+1/2)}{\Gamma(\alpha)}\,\frac{1}{\beta^{1/2}} \,\left(\frac{\gamma}{\beta}\right)^{-(\alpha+1/2)}\\ &= \frac{1}{\sqrt{\pi\, (2\beta c_2)}}\,\,\frac{\Gamma(\alpha+1/2)}{\Gamma(\alpha)}\, \,\left(\frac{(\theta-\mu)^2}{2c_2\beta} + 1\right)^{-(\alpha+1/2)}\\ &= \frac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2)}\,\frac{1}{\sqrt{\pi\,\nu\psi^2}} \,\left(\frac{1}{\nu}\frac{\,(\theta-\mu)^2}{\psi^2} + 1\right)^{-(\nu+1)/2} \end{align*} con \boxed{ \psi^2\equiv2\beta c_2 }\quad\text{and}\quad \boxed{\nu\equiv2\alpha} lo que significa que (\theta|\,\boldsymbol{x}) sigue una distribución t Student con parámetro de localización \mu, parámetro de escala \psi y grados de libertad \nu: \therefore\quad\boxed{(\theta|\,\boldsymbol{x}) \sim t(\mu, \psi, \nu)}

  1. Simule n=1000 observaciones i.i.d de \textsf{N}(5,1). Ajuste el modelo suponiendo los siguientes escenarios previos: i. estados de conocimiento previo bastante informativo acerca de los valores reales de los parámetros, ii. estado de conocimiento previo informativo acerca de \theta y difuso acerca de \sigma^2, iii. estado de conocimiento previo informativo acerca de \sigma^2 y difuso acerca de \theta, y iv. estado de conocimiento previo difuso acerca de ambos parámetros. Caracterice la distribución posterior en cada caso.

  2. Suponga que está interesado en hacer inferencia sobre \eta=\sigma/\theta. Desarrolle un algoritmo de Monte Carlo para calcular la media posterior y un intervalo de credibilidad al 95% para \eta. Use el algoritmo para calcular estas cantidades en todos los escenarios descritos anteriormente.

Se simulan n = 1000 observaciones iid a partir de N(5, 1) y se ajusta el modelo a estos datos, bajo los siguientes escenarios:

  1. Previas informativas sobre \theta y \sigma^2 alrededor de los valores verdaderos de los parámetros.
  2. Previa informativa sobre \theta y previa difusa sobre \sigma^2.
  3. Previa informativa sobre \sigma^2 y previa difusa sobre \theta.
  4. Previas difusas sobre \theta y \sigma^2.

Para lograr la especificación anterior los hiperparámetros para \theta se eligen de forma que:

Similarmente para \sigma^2. La siguiente tabla muesta los valores de los hiperparámetros de las previas en cada caso:

Caso \theta_0 \kappa_0 a b
1 5.0 0.5 26.0 25.0
2 5.0 0.5 10.0 0.01
3 5.0 2.0 26.0 25.0
4 5.0 2.0 10.0 0.01

Estas previas se muestran en el siguiente gráfico:

# datos
n <- 1000
theta <- 5
sigma <- 1
set.seed(123456789)
data <- rnorm(n, theta, sigma)
# grafico
par(mfrow=c(1,2), mar=c(3,3,2,1), mgp=c(1.75,.75,0))
# previa sigma^2
l <- 0.01; s <- 2
a <- 25; b <- a+1; fun <- function(x) pscl::densigamma(x, a, b)
curve(fun, l, s, col="blue", las=1, lwd=1, n=1000, cex.axis=0.75, main=expression( paste( 'p(' , sigma^2 , ')' ) ), ylab="Densidad", xlab=expression(sigma^2))
a <- 10; b <- 0.1; fun <- function(x) pscl::densigamma(x, a, b)
curve(fun, l, s, col="red", lwd=1, n=1000, add=TRUE)
abline(v=sigma, lty = 2, col="black", lwd=1)
legend("toprigh", c("Informativa","Difusa"), col = c("blue", "red"), lty=1, lwd=2, cex=0.7, bty="n")
# previa theta | sigma^2
l <- -1; s <- 11
t0 <- 5; k0 <- 0.5; fun <- function(x) dnorm( x, t0, k0*sigma^2 )
curve(fun, l, s, col="blue", las=1, lwd=1, n=1000, cex.axis=0.75, main=expression( paste( 'p(' ,theta,'|' ,sigma^2 , ')' ) ), ylab="Densidad", xlab=expression(theta) )
t0 <- 5; k0 <- 2; fun <- function(x) dnorm( x, t0, k0*sigma^2 )
curve(fun, l, s, col="red", lwd=1, n=1000, add=TRUE)
abline(v=theta, lty = 2, col="black", lwd=1)
legend("toprigh", c("Informativa","Difusa"), col = c("blue", "red"),lty=1, lwd=2, cex=0.7, bty="n")

Se lleva a cabo una simulación de Monte Carlo con S=100000 iteraciones siguiendo las expresiones de las distribuciones (\sigma^2|\,\boldsymbol{x}) y (\theta|\,\sigma^2,\boldsymbol{x}). La siguiente rutina lleva a cabo el ajuste del modelo y el resumen posterior de los parámetros correspondientes.

fit <- function(n, theta, sigma, data, t0, k0, a, b, S) {
    ### posterior ###
    alpha <- a + n/2
    sum.x2 <- sum( data^2 ) # sum.x2 := sum_{i=1}^n x_i^2
    xb <- mean(data)        # xb := sample mean of x
    B <- ( k0*sum.x2 + 2*b*k0 + t0^2 ) / ( 2*k0 )
    beta  <- B - ( n*xb*k0 + t0 )^2 / ( 2*k0*(n*k0 + 1) )
    mu <- ( n*xb*k0 + t0 ) / ( n*k0 + 1 )
    k <- sqrt( k0 / ( n*k0 +1 ) )
    ### Monte Carlo ###
    suppressMessages(suppressWarnings(library(pscl)))
    sig.mc <- rigamma(S, alpha, beta)
    the.mc <- rep(NA, S)
    eta.mc <- rep(NA, S)
    for(i in 1:S)   {
        tau <-  k*sqrt( sig.mc[i] )
        the.mc[i] <- rnorm( n = 1, mean = mu, sd = tau )
        eta.mc[i] <- sqrt( sig.mc[i] ) / the.mc[i]
    }
    ###resumen posterio ###
    ### sigma^2 ###
    sig.me <- mean( sig.mc )
    sig.sd <- sd( sig.mc )
    sig.LI <- quantile( sig.mc, probs=c(0.025) )
    sig.LS <- quantile( sig.mc, probs=c(0.975) )
    ### theta ###
    the.me <- mean( the.mc )
    the.sd <- sd( the.mc )
    the.LI <- quantile( the.mc, probs=c(0.025) )
    the.LS <- quantile( the.mc, probs=c(0.975) )
    ### eta ###
    eta.me <- mean( eta.mc )
    eta.sd <- sd( eta.mc )
    eta.LI <- quantile( eta.mc, probs=c(0.025) )
    eta.LS <- quantile( eta.mc, probs=c(0.975) )
    ### resumen ###
    summa <- as.matrix( c( alpha, beta, mu,
        sig.me, sig.sd, sig.LI, sig.LS,
        the.me, the.sd, the.LI, the.LS,
        eta.me, eta.sd, eta.LI, eta.LS
        ) )
    rownames(summa) <- c( 'alpha', 'beta', 'mu',
        'sigma Media', 'sigma DE', 'sigma Q2.5%', 'sigma Q97.5%',
        'theta Media', 'theta DE', 'theta Q2.5%', 'theta Q97.5%',
        'eta   Media', 'eta   DE', 'eta   Q2.5%', 'eta   Q97.5%'
    )
    list( sig.mc, the.mc, eta.mc, summa )
}

La siguiente tabla muestra un resumen posterior de las cantidades de interés \theta, \sigma, y \eta=\sigma/\theta.

### resumen ###
c1 <- fit(n, theta, sigma, data, t0=5, k0=0.5, a=26, b=25,  S=100000)
c2 <- fit(n, theta, sigma, data, t0=5, k0=0.5, a=10, b=0.1, S=100000)
c3 <- fit(n, theta, sigma, data, t0=5, k0=2.0, a=26, b=25,  S=100000)
c4 <- fit(n, theta, sigma, data, t0=5, k0=2.0, a=10, b=0.1, S=100000)
tab <- cbind(c1[[4]] , c2[[4]] , c3[[4]] , c4[[4]])
colnames(tab) <- paste("Caso ", 1:4)
knitr::kable(x = tab, digits = 3, align = "c")
Caso 1 Caso 2 Caso 3 Caso 4
alpha 526.000 510.000 526.000 510.000
beta 538.513 513.613 538.513 513.613
mu 5.012 5.012 5.012 5.012
sigma Media 1.026 1.009 1.026 1.009
sigma DE 0.045 0.045 0.045 0.045
sigma Q2.5% 0.942 0.925 0.942 0.925
sigma Q97.5% 1.117 1.101 1.117 1.101
theta Media 5.012 5.012 5.012 5.012
theta DE 0.032 0.032 0.032 0.032
theta Q2.5% 4.950 4.950 4.949 4.950
theta Q97.5% 5.075 5.074 5.075 5.075
eta Media 0.202 0.200 0.202 0.200
eta DE 0.005 0.005 0.005 0.005
eta Q2.5% 0.193 0.192 0.193 0.192
eta Q97.5% 0.211 0.210 0.211 0.210

A continuación se presenta una gráfico de la distribución posterior de las cantidades de interés \theta, \sigma, y \eta=\sigma/\theta. Dado que el tamaño de la muestra es grande, se observa que el impacto de la información previa sobre la distribución posterior es débil en todos los casos.

par(mfrow=c(1,3), mar=c(3,3,2,1), mgp=c(1.75,.75,0))
### sigma^2 ###
sig.mc1 <- c1[[1]]; sig.mc2 <- c2[[1]]; sig.mc3 <- c3[[1]]; sig.mc4 <- c4[[1]] 
l <- min(sig.mc1, sig.mc2, sig.mc3, sig.mc4)
s <- max(sig.mc1, sig.mc2, sig.mc3, sig.mc4)
plot(density(sig.mc1), las=1, cex.axis=0.75, 
    xlim=c(l, s),
    xlab=expression(sigma^2), sub=" ", ylab='Densidad',
    main=expression( paste( 'p(' , sigma^2 , '|x)' ) ))
    lines(density(sig.mc2), col='blue')
    lines(density(sig.mc3), col='red')
    lines(density(sig.mc4), col='green')
abline(v=sigma, lty = 2, col='grey', lwd=1)
legend("toprigh", legend = paste("Caso", 1:4), 
    col = c('black',"blue", "red", 'green'), 
    lty=1, lwd=2, cex=0.7, bty="n")
### theta ###
the.mc1 <- c1[[2]]; the.mc2 <- c2[[2]]; the.mc3 <- c3[[2]]; the.mc4 <- c4[[2]] 
l <- 4.85
s <- 5.15
plot(density(the.mc1), las=1, cex.axis=0.75, 
    xlim=c(l, s),
    xlab=expression(theta), sub=" ", ylab='Densidad',
    main=expression( paste( 'p(',theta , '|',sigma^2 , ',x)' ) ))
    lines(density(the.mc2), col='blue')
    lines(density(the.mc3), col='red')
    lines(density(the.mc4), col='green')
abline(v=theta, lty = 2, col='grey', lwd=1)
legend("toprigh", legend = paste("Caso", 1:4), 
    col = c('black',"blue", "red", 'green'), 
    lty=1, lwd=2, cex=0.7, bty="n")
### eta ###
eta.mc1 <- c1[[3]]; eta.mc2 <- c2[[3]]; eta.mc3 <- c3[[3]]; eta.mc4 <- c4[[3]] 
l <- 0.17
s <- 0.23
plot(density(eta.mc1), las=1, cex.axis=0.75, 
    xlim=c(l, s),
    xlab=expression(eta), sub=" ", ylab='Densidad',
    main=expression( paste( 'p(',eta , '|',theta,',',sigma^2 , ',x)' ) ))
    lines(density(eta.mc2), col='blue')
    lines(density(eta.mc3), col='red')
    lines(density(eta.mc4), col='green')
abline(v=sigma/theta, lty = 2, col='grey', lwd=1)
legend("toprigh", legend = paste("Caso", 1:4), 
    col = c('black',"blue", "red", 'green'), 
    lty=1, lwd=2, cex=0.7, bty="n")

  1. Considere el modelo Normal x_i\mid\theta,\sigma^2 \stackrel{\text{iid}}{\sim} \textsf{N}(\theta,\sigma^2), para i=1,\ldots,n, donde \theta es desconocido y \sigma^2 es conocido. Además, considere una distribución previa para \theta definida por medio de una mezcla finita de previas conjugadas de la forma p(\theta) = \sum_{\ell=1}^K w_\ell\,\phi(\theta\mid\mu_\ell,\tau^2)\, donde K es un entero positivo fijo mayor o igual que 1, w_1,\ldots,w_K es un sistema de pesos tales que y \sum_{\ell=1}^K w_\ell = 1 y 0\leq w_\ell\leq 1 para \ell=1,\ldots,K, y \phi(\theta\mid\mu,\tau^2) denota la densidad de la distribución Normal con media \mu y varianza \tau^2. Una distribución previa de esta forma permite especificar estados de información previos multimodales acerca de \theta.
  1. Encuentre la distribución posterior de \theta.

Se tiene que \begin{align*} \pi(\theta|\boldsymbol{x}) &\propto f(\boldsymbol{x}|\theta)\pi(\theta)\\ &\propto \left[\prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma}\,\exp{\left\{ -\frac{1}{2\sigma^2}(x_i-\theta)^2 \right\} }\right]\cdot\\ &\hspace{4cm}\left[ \sum_{l=1}^K w_l\, \frac{1}{\sqrt{2\pi}\tau} \exp{\left\{ -\frac{1}{2\tau^2}(\theta-\mu_l)^2 \right\} }\right]\\ &\propto \left[\exp{\left\{ -\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\theta)^2 \right\} }\right]\cdot \left[ \sum_{l=1}^K w_l\, \exp{\left\{ -\frac{1}{2\tau^2}(\theta-\mu_l)^2 \right\} }\right]\\ &\propto \sum_{l=1}^K w_l\, \exp{\left\{ -\frac{1}{2}\left(\frac{1}{\sigma^2}\sum_{i=1}^n (x_i-\theta)^2 +\frac{1}{\tau^2}(\theta-\mu_l)^2\right) \right\} }\\ &\propto \sum_{l=1}^K w_l\, \exp{\left\{ -\frac{1}{2}\left(\frac{\sum x_i^2 - 2n\bar{x}\theta+ n\theta^2}{\sigma^2} +\frac{\theta^2 -2\mu_l\theta+\mu_l^2}{\tau^2}\right) \right\} }\\ &\propto \sum_{l=1}^K w_l\, \exp{\left\{ -\frac{1}{2}\left[ \theta^2\left(\frac{n}{\sigma^2}+\frac{1}{\tau^2}\right) - 2\theta\left(\frac{n\bar{x}}{\sigma^2} + \frac{\mu_l}{\tau^2}\right) + \left(\frac{\sum x_i^2}{\sigma^2} +\frac{ \mu_l^2}{\tau^2} \right) \right] \right\} }\\ &\propto \sum_{l=1}^K w_l\, \exp{\left\{ -\frac{1}{2}\left[ \frac{\theta^2}{\upsilon^2} - 2\theta\frac{\xi_l}{\upsilon^2} + \psi_l \right] \right\} } \end{align*} donde (\upsilon^2)^{-1}=\frac{n}{\sigma^2}+\frac{1}{\tau^2}=\frac{n\tau^2+\sigma^2}{\sigma^2\tau^2}\quad\Rightarrow\quad \boxed{\upsilon^2=\frac{\sigma^2\tau^2}{n\tau^2+\sigma^2}} \frac{\xi_l}{\upsilon^2}=\frac{n\bar{x}}{\sigma^2} + \frac{\mu_l}{\tau^2}=\frac{n\bar{x}\tau^2+\sigma^2\mu_l}{\sigma^2\tau^2}\quad\Rightarrow\quad \boxed{\xi_l=\frac{n\bar{x}\tau^2+\sigma^2\mu_l}{n\tau^2+\sigma^2}} \boxed{\psi_l=\frac{ \mu_l^2}{\tau^2}} Entonces, \begin{align*} \pi(\theta|\boldsymbol{x}) &\propto \sum_{l=1}^K w_l\, \exp{\left\{ -\frac{1}{2}\left[ \frac{\theta^2}{\upsilon^2} - 2\theta\frac{\xi_l}{\upsilon^2} + \frac{\xi_l^2}{\upsilon^2} \right] \right\} }\, \exp{\left\{-\frac12\left[\psi_l -\frac{\xi_l^2}{\upsilon^2}\right]\right\}}\\ &\propto \sum_{l=1}^K w_l^*\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} } \end{align*} donde \boxed{w_l^*=w_l\,\exp{\left\{-\frac12\left[\psi_l -\frac{\xi_l^2}{\upsilon^2}\right]\right\}}} Así, \begin{align*} \pi(\theta|\boldsymbol{x}) &= \frac{\sum_{l=1}^K w_l^*\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} }}{\int \sum_{l=1}^K w_l^*\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} } d\theta}\\ &= \frac{\sum_{l=1}^K w_l^*\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} }} {\sqrt{2\pi}\upsilon\sum_{l=1}^K w_l^*\, \int \frac{1}{\sqrt{2\pi}\upsilon}\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} } d\theta}\\ &= \frac{1}{\sqrt{2\pi}\upsilon\sum_{l=1}^K w_l^*}\sum_{l=1}^K w_l^*\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} } \\ &= \sum_{l=1}^K \omega_l\, \frac{1}{\sqrt{2\pi}\upsilon}\exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} } \\ &= \sum_{l=1}^K \omega_l\, \phi(\theta|\xi_l,\upsilon^2) \\ %&\\ &\Rightarrow \boxed{\pi(\theta|\boldsymbol{x})=\sum_{l=1}^K \omega_l\, \phi(\theta|\xi_l,\upsilon^2)} \end{align*} donde \boxed{\omega_l=\frac{w_l^*}{\sum_{j=1}^K w_j^*}} y \phi(\theta|\xi_l,\upsilon^2) denota la función de densidad de una variable aleatoria con distribución Normal con media \xi_l y varianza \upsilon^2.

  1. Encuentre la media posterior de \theta.

La media posterior de \theta está dada por \begin{align*} \mathbb{E}\left(\theta|\boldsymbol{x}\right)&=\int \theta\,\pi(\theta|\boldsymbol{x}) d\theta\\ &=\int \theta\,\left[ \sum_{l=1}^K \omega_l\, \phi(\theta|\xi_l,\upsilon^2)\right] d\theta\\ &=\sum_{l=1}^K \omega_l \left[ \int \theta\, \phi(\theta|\xi_l,\upsilon^2) d\theta\right]\\ &=\sum_{l=1}^K \omega_l \xi_l\\ %&\\ &\Rightarrow\boxed{\mathbb{E}\left(\theta|\boldsymbol{x}\right)=\sum_{l=1}^K \omega_l \xi_l} \end{align*}

  1. Encuentre la distribución predictiva previa.

La distribución predicrtiva previa asociada con este modelo es \begin{align*} p(\boldsymbol{x}) &= \int f(\boldsymbol{x}|\theta)\,\pi(\theta) d\theta\\ &=\frac{1}{(2\pi)^{n/2}\sigma^n}\frac{1}{\sqrt{2\pi}\tau}\int \sum_{l=1}^K w_l^*\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} } d\theta\\ &=\frac{1}{(2\pi)^{(n+1)/2}\sigma^n\tau}\sqrt{2\pi}\upsilon\sum_{l=1}^K w_l^*\, \int \frac{1}{\sqrt{2\pi}\upsilon} \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} } d\theta\\ &=\frac{\upsilon}{(2\pi)^{n/2}\sigma^n\tau}\sum_{l=1}^K w_l^*\\ &=\frac{\upsilon}{(2\pi)^{n/2}\sigma^n\tau}\sum_{l=1}^K w_l\, \exp{\left\{-\frac12\left[\frac{\tau^2\sum x_i^2 + \sigma^2 \mu_l^2}{\sigma^2\tau^2} -\frac{\xi_l^2}{\upsilon^2}\right]\right\}}\\ %&\\ &\Rightarrow\boxed{ p(\boldsymbol{x})=(2\pi)^{-n/2}\frac{\upsilon}{\sigma^n\tau}\exp{\left\{-\frac12\left[\frac{\tau^2\sum x_i^2 + \sigma^2 \mu_l^2}{\sigma^2\tau^2} -\frac{\xi_l^2}{\upsilon^2}\right]\right\}}} \end{align*} donde \upsilon^2=\frac{\sigma^2\tau^2}{n\tau^2+\sigma^2} \,\,\, \text{and} \,\,\, \xi_l=\frac{n\bar{x}\tau^2+\sigma^2\mu_l}{n\tau^2+\sigma^2}

  1. Encuentre la distribución predictiva posterior.

Para encontrar la distribución predictiva posterior asociada con este mode, se sigue el procedimiento anterior para obtener que \begin{align*} g(x^*|\boldsymbol{x}) &= \int f(x^*|\theta)\, \pi(\theta|\boldsymbol{x})\, d\theta\\ &= \frac{1}{2\pi\sigma\upsilon} \int \left[\,\exp{\left\{ -\frac{1}{2\sigma^2}(x^*-\theta)^2 \right\} }\right]\cdot \left[ \sum_{l=1}^K \omega_l\, \exp{\left\{ -\frac{1}{2\upsilon^2}(\theta-\xi_l)^2 \right\} }\right] \,d\theta\\ &= \frac{1}{2\pi\sigma\upsilon} \sum_{l=1}^K \omega_l\,\int \exp{\left\{ -\frac{1}{2}\left(\frac{1}{\sigma^2} (x^*-\theta)^2 +\frac{1}{\upsilon^2}(\theta-\xi_l)^2\right) \right\} }d\theta\\ &= \frac{1}{2\pi\sigma\upsilon} \sum_{l=1}^K \omega_l\, \int \exp{\left\{ -\frac{1}{2}\left(\frac{x^{*2} - 2x^*\theta+ n\theta^2}{\sigma^2} +\frac{\theta^2 -2\xi_l\theta+\xi_l^2}{\upsilon^2}\right) \right\} }d\theta\\ &= \frac{1}{2\pi\sigma\upsilon} \sum_{l=1}^K \omega_l\, \int \exp{\left\{ -\frac{1}{2}\left[ \theta^2\left(\frac{n}{\sigma^2}+\frac{1}{\upsilon^2}\right) - 2\theta\left(\frac{x^*}{\sigma^2} + \frac{\xi_l}{\upsilon^2}\right) + \left(\frac{x^{*2}}{\sigma^2} +\frac{ \xi_l^2}{\upsilon^2} \right) \right] \right\} }d\theta\\ &= \frac{1}{2\pi\sigma\upsilon} \sum_{l=1}^K \omega_l\,\int \exp{\left\{ -\frac{1}{2}\left[ \frac{\theta^2}{s^2} - 2\theta\frac{m_l}{s^2} + c_l \right] \right\} } d\theta \end{align*} donde (s^2)^{-1}=\frac{n}{\sigma^2}+\frac{1}{\upsilon^2}=\frac{n\upsilon^2+\sigma^2}{\sigma^2\upsilon^2}\quad\Rightarrow\quad \boxed{s^2=\frac{\sigma^2\upsilon^2}{n\upsilon^2+\sigma^2}} \frac{m_l}{s^2}=\frac{x^*}{\sigma^2} + \frac{\xi_l}{\upsilon^2}=\frac{x^*\upsilon^2+\sigma^2\xi_l}{\sigma^2\upsilon^2}\quad\Rightarrow\quad \boxed{m_l=\frac{x^*\upsilon^2+\sigma^2\xi_l}{n\upsilon^2+\sigma^2}} \boxed{c_l=\frac{ \xi_l^2}{\upsilon^2}} Entonces, \begin{align*} g(x^*|\boldsymbol{x}) &= \frac{1}{2\pi\sigma\upsilon} \sum_{l=1}^K \omega_l\, \int \exp{\left\{ -\frac{1}{2}\left[ \frac{\theta^2}{s^2} - 2\theta\frac{m_l}{s^2} + \frac{m_l^2}{s^2} \right] \right\} }\, \exp{\left\{-\frac12\left[c_l -\frac{m_l^2}{s^2}\right]\right\}} d\theta\\ &= \frac{1}{2\pi\sigma\upsilon} \sqrt{2\pi}s \sum_{l=1}^K \omega_l\,\exp{\left\{-\frac12\left[c_l -\frac{m_l^2}{s^2}\right]\right\}}\, \int \frac{1}{\sqrt{2\pi}s} \, \exp{\left\{ -\frac{1}{2s^2}(\theta-m_l)^2 \right\} }d\theta\\ &= \frac{s}{\sqrt{2\pi}\sigma\upsilon} \sum_{l=1}^K \omega_l\,\exp{\left\{-\frac12\left[c_l -\frac{m_l^2}{s^2}\right]\right\}}\\ %&\\ &\Rightarrow \boxed{g(x^*|\boldsymbol{x})=\frac{s}{\sqrt{2\pi}\sigma\upsilon} \sum_{l=1}^K \omega_l\,\exp{\left\{-\frac12\left[\frac{\upsilon^2 x^{*2} + \sigma^2 \xi_l^2}{\sigma^2\upsilon^2} -\frac{m_l^2}{s^2}\right]\right\}}} \end{align*}

  1. Considere el modelo Normal Trucado x_i\mid\theta,\sigma^2 \stackrel{\text{iid}}{\sim} \textsf{N}_{(0,\infty)}(\theta,\sigma^2), para i=1,\ldots,n, donde \sigma^2=1. Además, considere la distribución previa para \theta dada por \theta\sim\textsf{N}_{(0,\infty)}(\mu,\tau^2).
  1. Encuentre la distribución posterior de \theta.
  2. ¿Este modelo se puede catalogar como un modelo conjugado?

La distribución posterior es proporcional a \begin{align*} \pi(\theta|x) &\propto f(\boldsymbol{x}|\theta)\,\pi(\theta)\\ &\propto \left[ \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\Phi(\theta)}\,\exp{\left\{ -\frac12(x_i-\theta)^2\right\}}I_{[0,\infty)}(x_i) \right] \cdot \frac{1}{\sqrt{2\pi}\tau}\,\exp{\left\{ -\frac{1}{2\tau^2}(\theta-\mu)^2\right\}}\\ &\propto \frac{1}{[\Phi(\theta)]^n}\, \exp{\left\{ -\frac12\left[\sum_{i=1}^n (x_i-\theta)^2 + \frac{(\theta-\mu)^2}{\tau^2}\right] \right\}}\\ &\propto \frac{1}{[\Phi(\theta)]^n}\, \exp{\left\{ -\frac12\left[\sum_{i=1}^nx_i^2 -2\theta n\bar{x}+ \theta^2 + \frac{\theta^2-2\theta\mu+\mu^2}{\tau^2}\right] \right\}}\\ &\propto \frac{1}{[\Phi(\theta)]^n}\, \exp{\left\{ -\frac12\left[ \theta^2\left(1+\frac{1}{\tau^2}\right) -2\theta\left(n\bar{x}+ \frac{\mu}{\tau^2}\right)\right] \right\}}\, \exp{\left\{ -\frac12\left[\sum_{i=1}^nx_i^2 + \frac{\mu^2}{\tau^2}\right] \right\}}\\ &\propto \frac{1}{[\Phi(\theta)]^n}\, \exp{\left\{ -\frac12 \left[ \frac{\theta^2}{\upsilon^2} -\frac{2\theta\xi}{\upsilon^2} + \frac{\xi^2}{\upsilon^2}\right] \right\}} \exp{\left\{ \frac12 \frac{\xi^2}{\upsilon^2} \right\}}\\ &\propto \frac{1}{[\Phi(\theta)]^n}\, \exp{\left\{ -\frac{1}{2\upsilon^2} (\theta-\xi)^2 \right\}}\\ %&\\ &\Rightarrow\boxed{\frac{1}{[\Phi(\theta)]^n}\, \exp{\left\{ -\frac{1}{2\upsilon^2} (\theta-\xi)^2 \right\}}} \end{align*} donde (\upsilon^2)^{-1}= 1+\frac{1}{\tau^2} =\frac{\tau^2+1}{\tau^2}\quad\Rightarrow\quad \boxed{\upsilon^2=\frac{\tau^2}{\tau^2+1}} \frac{\xi}{\upsilon^2}=n\bar{x}+ \frac{\mu}{\tau^2}= \frac{n\bar{x}\tau^2+\mu}{\tau^2}\quad\Rightarrow\quad \boxed{\xi=\frac{n\bar{x}\tau^2+\mu}{\tau^2+1}} Esta no es una previa conjugada porque se tiene un término extra, a saber (1/[\Phi(\theta)]^n), en el núcleo de la distribución posterior que lo hace diferente del núcleo de una distribución Nomrl, lo que significa que la distribución posterior no es miembro de la familia de distribuciones Normales, y por lo tanto la previa no es conjugada.