Respuesta :
Answer:
See answers below
Step-by-step explanation:
# prior
mu0 <- 5
k0 <- 1
s20 <- 4
nu0 <- 2
# read in data
dat <- list()
dat[1] <- read.table("school1.dat")
dat[2] <- read.table("school2.dat")
dat[3] <- read.table("school3.dat")
n <- sapply(dat, length)
ybar <- sapply(dat, mean)
s2 <- sapply(dat, var)
# posterior
kn <- k0 + n
nun <- nu0 + n
mun <- (k0 * mu0 + n * ybar)/kn
s2n <- (nu0 * s20 + (n - 1) * s2 + k0 * n * (ybar - mu0)^2/kn)/nun
# produce a posterior sample for the parameters
s.postsample <- s2.postsample <- theta.postsample <- matrix(0, 10000, 3, dimnames = list(NULL,
c("school1", "school2", "school3")))
for (i in c(1, 2, 3)) {
s2.postsample[, i] <- 1/rgamma(10000, nun[i]/2, s2n[i] * nun[i]/2)
s.postsample[, i] <- sqrt(s2.postsample[, i])
theta.postsample[, i] <- rnorm(10000, mun[i], s.postsample[, i]/sqrt(kn[i]))
}
# posterior mean and 95% CI for Thetas (the means)
colMeans(theta.postsample)
## school1 school2 school3
## 9.294 6.943 7.814
apply(theta.postsample, 2, function(x) {
quantile(x, c(0.025, 0.975))
})
## school1 school2 school3
## 2.5% 7.769 5.163 6.152
## 97.5% 10.834 8.719 9.443
# posterior mean and 95% CI for the Standard Deviations
colMeans(s.postsample)
## school1 school2 school3
## 3.904 4.397 3.757
apply(s.postsample, 2, function(x) {
quantile(x, c(0.025, 0.975))
})
## school1 school2 school3
## 2.5% 2.992 3.349 2.800
## 97.5% 5.136 5.871 5.139