school1 <- read.table(file = "school1.dat",
quote = "/",
comment.char = "")
school2 <- read.table(file = "school2.dat",
quote = "/",
comment.char= "")
school3 <- read.table(file = "school3.dat",
quote = "/",
comment.char= "")
school1 <- school1 %>%
dplyr::select(V1) %>%
as.matrix() %>%
as.numeric()
school2 <- school2 %>%
dplyr::select(V1) %>%
as.matrix() %>%
as.numeric()
school3 <- school3 %>%
dplyr::select(V1) %>%
as.matrix() %>%
as.numeric()The kernel part derived in HomeCopy-II is non-standard or unfamiliar. So, we can’t sample from it directly. Instead we will implement Monte-Carlo simulation to approximate \((\sigma^2 | y_1, ..., y_n)\) and \((\theta | \sigma^2, y_1, ..., y_n)\).
First, we sample \(\sigma^2\) from it’s marginal posterior and then implement those values to sample \(\theta\) from its full conditional distribution.
We use the obtained Monte-Carlo samples to compute necessary means and \(95\)% confidence intervals. Since, posterior distributions are not purely theoretical quantities, we will compute (empirical) quantile-based confidence intervals from the produced samples.
Note that, the theoretical mean of the full conditional distribution, \((\theta | \sigma^2, y_1, ..., y_n)\), is NOT the posterior mean for \(\theta\). The “posterior mean” refers to the marginal posterior’s mean. But that is an unrecognized distribution, at least for \(\theta\) it is. Hence, we have to approximate it by averaging the relevant Monte-Carlo samples for each school.
# Hyper-parameters
mu0 <- 5; sigma20 <- 4; k0 <- 1; v0 <- 2
# sample means
y1bar <- school1 %>% mean()
y2bar <- school2 %>% mean()
y3bar <- school3 %>% mean()
# sample sizes
n1 <- school1 %>% length()
n2 <- school2 %>% length()
n3 <- school3 %>% length()
# posterior parameters
kn1 <- k0 + n1
kn2 <- k0 + n2
kn3 <- k0 + n3
vn1 <- v0 + n1
vn2 <- v0 + n2
vn3 <- v0 + n3
##############################################################
# sampling from 1st inverse-Gamma posterior
`sigma1^2|data` <- 1 / rgamma(
n = 1000,
shape = vn1 / 2,
rate = (vn1/2)*(
v0*(sigma20^2) +
(n1-1)*var(school1) +
(k0*n1/kn1)*(y1bar-mu0)
)
)
# sampling from 2nd inverse-Gamma posterior
`sigma2^2|data` <- 1 / rgamma(
n = 1000,
shape = vn2 / 2,
rate = (vn2/2)*(
v0*(sigma20^2) +
(n2-1)*var(school2) +
(k0*n2/kn2)*(y2bar-mu0)
)
)
# sampling from 3rd inverse-Gamma posterior
`sigma3^2|data` <- 1 / rgamma(
n = 1000,
shape = vn3 / 2,
rate = (vn3/2)*(
v0*(sigma20^2) +
(n3-1)*var(school3) +
(k0*n3/kn3)*(y3bar-mu0)
)
)
##############################################################
# theoretical means from full conditional distribution of theta
mu1 <- weighted.mean(x = c(mu0, y1bar), # school 1
w = c(k0, n1))
mu2 <- weighted.mean(x = c(mu0, y2bar), # school 2
w = c(k0, n2))
mu3 <- weighted.mean(x = c(mu0, y3bar), # school 3
w = c(k0, n3))
##############################################################
# sampling theta from 1st full conditional posterior
`theta|sigma1^2, data` <- sapply(
X = `sigma1^2|data`,
FUN = function(sig){
rnorm(n = 1, mean = mu1, sd = sig)
}
)
# sampling theta from 2nd full conditional posterior
`theta|sigma2^2, data` <- sapply(
X = `sigma2^2|data`,
FUN = function(sig){
rnorm(n = 1, mean = mu2, sd = sig)
}
)
# sampling theta from 2nd full conditional posterior
`theta|sigma3^2, data` <- sapply(
X = `sigma3^2|data`,
FUN = function(sig){
rnorm(n = 1, mean = mu3, sd = sig)
}
)
##############################################################
# posterior means approximated for theta
post.mean.theta.1 <- mean(`theta|sigma1^2, data`)
post.mean.theta.2 <- mean(`theta|sigma2^2, data`)
post.mean.theta.3 <- mean(`theta|sigma3^2, data`)
# posterior means approximated for sigma
post.mean.sigma.1 <- mean(sqrt(`sigma1^2|data`))
post.mean.sigma.2 <- mean(sqrt(`sigma2^2|data`))
post.mean.sigma.3 <- mean(sqrt(`sigma3^2|data`))
##############################################################
# quantile-based 95% confidence interval for theta
post.ci.theta.1 <- quantile(x = `theta|sigma1^2, data`,
probs = c(0.025, 0.975))
post.ci.theta.2 <- quantile(x = `theta|sigma2^2, data`,
probs = c(0.025, 0.975))
post.ci.theta.3 <- quantile(x = `theta|sigma3^2, data`,
probs = c(0.025, 0.975))
# quantile-based 95% confidence interval for sigma
post.ci.sigma.1 <- quantile(x = `sigma1^2|data`,
probs = c(0.025, 0.975)) %>% sqrt()
post.ci.sigma.2 <- quantile(x = `sigma2^2|data`,
probs = c(0.025, 0.975)) %>% sqrt()
post.ci.sigma.3 <- quantile(x = `sigma3^2|data`,
probs = c(0.025, 0.975)) %>% sqrt()theta1 <- `theta|sigma1^2, data`
theta2 <- `theta|sigma2^2, data`
theta3 <- `theta|sigma3^2, data`
(abc <- mean(theta1 < theta2 & theta2 < theta3))## [1] 0.127(acb <- mean(theta1 < theta3 & theta3 < theta2))## [1] 0.195(bac <- mean(theta2 < theta1 & theta1 < theta3))## [1] 0.177(bca <- mean(theta2 < theta3 & theta3 < theta1))## [1] 0.203(cab <- mean(theta3 < theta1 & theta1 < theta2))## [1] 0.146(cba <- mean(theta3 < theta2 & theta2 < theta1))## [1] 0.152Now, we need to simulate data from the posterior predictive distributions. The form of the posterior predictive distribution is too complicated and time-consuming to derive if not impossible. So, we’re going to employ Monte-Carlo simulation to approximate it instead.
# sampling from school1's posterior predictive distribution
y1pred <- sapply(
X = 1:length(`sigma1^2|data`),
FUN = function(i){
rnorm(n = 1,
mean = `theta|sigma1^2, data`[i],
sd = `sigma1^2|data`[i])
}
)
# sampling from school2's posterior predictive distribution
y2pred <- sapply(
X = 1:length(`sigma2^2|data`),
FUN = function(i){
rnorm(n = 1,
mean = `theta|sigma2^2, data`[i],
sd = `sigma2^2|data`[i])
}
)
# sampling from school3's posterior predictive distribution
y3pred <- sapply(
X = 1:length(`sigma3^2|data`),
FUN = function(i){
rnorm(n = 1,
mean = `theta|sigma3^2, data`[i],
sd = `sigma3^2|data`[i])
}
)
(abc <- mean(y1pred < y2pred & y2pred < y3pred))## [1] 0.151(acb <- mean(y1pred < y3pred & y3pred < y2pred))## [1] 0.178(bac <- mean(y2pred < y1pred & y1pred < y3pred))## [1] 0.183(bca <- mean(y2pred < y3pred & y3pred < y1pred))## [1] 0.187(cab <- mean(y3pred < y1pred & y1pred < y2pred))## [1] 0.161(cba <- mean(y3pred < y2pred & y2pred < y1pred))## [1] 0.14(ans1 <- mean(theta1 > theta2 & theta1 > theta3))## [1] 0.355(ans2 <- mean(y1pred > y2pred & y1pred > y3pred))## [1] 0.327What is sensitivity analysis?? Articulate.