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.Explore los datos gráfica y numéricamente.
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:
Dibuje la distribución posterior conjunta de (θ,σ2) para cada escuela.
Compruebe la bondad de ajuste del modelo para cada escuela utilizando como estadísticos de prueba la media y el coeficiente de variación.
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.# 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.
Se ajusta el modelo con distribución muestral yi∣θ,σ2iid∼N(θ,σ2),i=1,…,n, y distribución previa conjugada θ∣σ2∼N(μ0,σ2κ0)σ2∼GI(ν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
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))
}
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)
}
}
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)}
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.
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:
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")
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.
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*}
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}
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*}
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.