N-N分层贝叶斯模型

y <- c(28, 8, -3, 7, -1, 1, 18, 12)
sd <- c(15, 10, 16, 11, 9, 11, 10, 18)
v <- sd * sd
tau <- c(0:3000)/100
tausq <- tau * tau
ptau.y <- rep(0, 3001)
vmu <- rep(0, 3001)
muhat <- rep(0, 3001)
for (i in (1:3001)) {
  vmu[i] <- 1/sum(1/(tausq[i] + v))
  muhat[i] <- vmu[i] * sum(y/(tausq[i] + v))
  ptau.y[i] <- sqrt(vmu[i] * prod(1/(tausq[i] + v))) * prod(exp(-0.5 *
                                                                  (y - muhat[i]) * (y - muhat[i])/(tausq[i] + v)))
}
plot(tau, ptau.y, type = "l", yaxt = "n", xlab = quote(tau))

eth.tauy <- matrix(rep(0, 24008), 8, 3001)
sdth.tauy <- matrix(rep(0, 24008), 8, 3001)
for (j in (1:8)) {
  for (i in (1:3001)) {
    eth.tauy[j, i] <- (y[j]/v[j] + muhat[i]/tausq[i])/(1/v[j] + 1/tausq[i])
    sdth.tauy[j, i] <- sqrt(((1/tausq[i])/(1/v[j] + 1/tausq[i])) *
                              ((1/tausq[i])/(1/v[j] + 1/tausq[i])) * vmu[i] + 1/(1/v[j] +
                                                                                   1/tausq[i]))
  }
}
par(mfrow = c(1, 2))
taux <- matrix(rep(tau, 8), 8, byrow = T)
matplot(t(taux), t(eth.tauy), ylim = (c(-5, 30)), type = "l", xlab = "tau",
        lty = 1:8, lwd = 1, col = 1, ylab = "Estimate Treatment Effects", main = "Conditional posterior mean")
School <- c("A", "B", "C", "D", "E", "F", "G", "H")
text(x = rep(20, 8), y = t(eth.tauy)[2400, ], School)
matplot(t(taux), t(sdth.tauy), ylim = (c(0, 20)), type = "l", xlab = "tau",
        lty = 1:8, lwd = 1, col = 1, ylab = "Posterior Standard Deviations",
        main = "Conditional posterior SD")
text(x = rep(20, 8), y = t(sdth.tauy)[2400, ], School)

m <- 200
ptau.y <- ptau.y/(sum(ptau.y))
tausamp <- sample(tau, m, replace = T, prob = ptau.y)
tausamp <- sort(tausamp)
tauid <- tausamp * 100 + 1
x <- matrix(rnorm(m * 10, 0, 1), m, 10)
x[, 1] <- tausamp
x[, 2] <- muhat[tauid] + sqrt(vmu[tauid]) * x[, 2]
for (j in (1:8)) {
  thmean <- (y[j] * x[, 1] * x[, 1] + v[j] * x[, 2])/(v[j] + x[, 1] *
                                                        x[, 1])
  thsd <- sqrt(v[j] * x[, 1] * x[, 1]/(v[j] + x[, 1] * x[, 1]))
  x[, j + 2] <- thmean + thsd * x[, j + 2]
}
par(mfrow = c(1, 2))
hist(x[, 2], breaks = c(-40:50), xlab = quote(mu), yaxt = "n", main = "")
hist(x[, 3], breaks = c(-20:60), xlab = "Effect in School A", yaxt = "n",
     main = "")