The files school1.dat, schoo12. dat and schoo13.dat contain data on the amount of time students from three high schools spent on studying or homework during an exam period. Analyze data from each of these schools separately, using the normal model with a conjugate prior -5, σ1-4,K0-1Po-2) and compute or distribution, in which f40 approximate the following: a) posterior means and 95% confidence intervals for the mean θ and standard deviation σ from each school b) the posterior probability that θǐ 〈 θj 〈 θk for all six permutations i, j,k of 1,2,3); c) the posterior probability that Yi 〈 Yk for all six permutations i, j, k} of {1,2,3], where , is a sample from the posterior predictive distribution of school i, d) Compute the posterior probability that is bigger than both θ2 and , and the posterior probability that Y1 is bigger than both Yo and 3

School1.dat:

2.11
9.75
13.88
11.3
8.93
15.66
16.38
4.54
8.86
11.94
12.47
11.11
11.65
14.53
9.61
7.38
3.34
9.06
9.45
5.98
7.44
8.5
1.55
11.45
9.73

School2.dat:

0.29
1.13
6.52
11.72
6.54
5.63
14.59
11.74
9.12
9.43
10.64
12.28
9.5
0.63
15.35
5.31
8.49
3.04
3.77
6.22
2.14
6.58
1.11

School3.dat:

4.33
7.77
4.15
5.64
7.69
5.04
10.01
13.43
13.63
9.9
5.72
5.16
4.33
12.9
11.27
6.05
0.95
6.02
12.22
12.85

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

take it another