Assignment 8

Md Shafiqul Islam

8.1
Consider the hierarchical normal model
yij | θj , σ 2 ∼ N (θj , σ 2 ), θj | µ, τ 2 ∼ N (µ, τ 2 ),
with hyperparameters µ, τ 2 , σ 2 .

(a)
We compare Var(yij | θj , σ 2 ) and Var(yij | µ, τ 2 ).
Conditioning on θj and σ 2 gives a simple normal model, so
Var(yij | θj , σ 2 ) = σ 2 .
If instead we condition only on µ, τ 2 , the group mean θj is still random. Using the law
of total variance,
Var(yij | µ, τ 2 ) = E Var(yij | θj , σ 2 , µ, τ 2 ) | µ, τ 2 + Var E[yij | θj , σ 2 , µ, τ 2 ] | µ, τ 2 .
 

The first term is E[σ 2 | µ, τ 2 ] = σ 2 , and the second term is Var(θj | µ, τ 2 ) = τ 2 . Hence
Var(yij | µ, τ 2 ) = σ 2 + τ 2 ,
which is larger than σ 2 because it combines within-group (σ 2 ) and between-group (τ 2 ) vari-
ation.

(b)
Now consider two observations from the same group, yi1 j and yi2 j with i1 ̸= i2 .
Given θj , σ 2 , the data within group j are independent:
Cov(yi1 j , yi2 j | θj , σ 2 ) = 0.
If we condition only on µ, τ 2 , the group mean θj is random, and the two observations
share the same random mean. Using the law of total covariance,
Cov(yi1 j , yi2 j | µ, τ 2 ) = E Cov(yi1 j , yi2 j | θj , σ 2 , µ, τ 2 ) | µ, τ 2 +Cov E[yi1 j | θj , σ 2 , µ, τ 2 ], E[yi2 j | θj , σ 2 , µ, τ 2 ]
 

The first term is zero (conditional independence given θj ), and the second term is Cov(θj , θj |
µ, τ 2 ) = τ 2 , so
Cov(yi1 j , yi2 j | µ, τ 2 ) = τ 2 > 0.
Thus observations in the same group are positively correlated when we condition only on
the hyperparameters.

1
(c)
Pnj
Let ȳ·j = n−1
j i=1 yij be the sample mean in group j.

Conditional on θj , σ 2 . Since the yij are independent N (θj , σ 2 ) given θj , σ 2 ,

σ2
Var(yij | θj , σ 2 ) = σ 2 , Var(ȳ·j | θj , σ 2 ) = ,
nj

and, as in part (b),

Cov(yi1 j , yi2 j | θj , σ 2 ) = 0 (i1 ̸= i2 ).

Conditional on µ, τ 2 . Using the law of total variance and the computations from part
(a),
Var(yij | µ, τ 2 ) = σ 2 + τ 2 .
For the group mean,

Var(ȳ·j | µ, τ 2 ) = E Var(ȳ·j | θj , σ 2 , µ, τ 2 ) | µ, τ 2 + Var E[ȳ·j | θj , σ 2 , µ, τ 2 ] | µ, τ 2

 

 2 
σ
=E 2
| µ, τ + Var(θj | µ, τ 2 )
nj
2
σ
= + τ 2.
nj

The covariance of two distinct observations in the same group is, as in part (b),

Cov(yi1 j , yi2 j | µ, τ 2 ) = τ 2 .

(d)
We are asked to examine the conditional distribution of µ given all other quantities and the
data. Starting from Bayes’ rule,

p(µ | θ1 , . . . , θm , σ 2 , τ 2 , y1 , . . . , ym ) ∝ p(y1 , . . . , ym | θ1 , . . . , θm , σ 2 ) p(θ1 , . . . , θm | µ, τ 2 ) p(µ) p(σ 2 ) p(τ 2 ).

In this expression, the likelihood term p(y1 , . . . , ym | θ1 , . . . , θm , σ 2 ) and the prior on σ 2

do not involve µ. Treating these as constants with respect to µ, we can absorb them into
the normalizing constant. Hence

p(µ | θ1 , . . . , θm , σ 2 , τ 2 , y1 , . . . , ym ) = p(µ | θ1 , . . . , θm , τ 2 ).

That is, once the group means θ1 , . . . , θm and the between-group variance τ 2 are known,
the data y1 , . . . , ym and the within-group variance σ 2 provide no additional information
about µ.

2
8.3
Eight files ([Link] through [Link]) provide weekly homework hours for students
in eight schools. We fit a hierarchical normal model:

yij | θj , σ 2 ∼ N (θj , σ 2 ), θj | µ, τ 2 ∼ N (µ, τ 2 ),

with priors matching Hoff:

µ ∼ N (µ0 , γ02 ), σ 2 ∼ IG(ν0 /2, ν0 σ02 /2), τ 2 ∼ IG(η0 /2, η0 τ02 /2).

µ0 = 7, γ02 = 5, τ02 = 10, η0 = 2, σ02 = 10, ν0 = 2.

We run Gibbs sampling to obtain posteriors of {θ1 , . . . , θ8 , µ, σ 2 , τ 2 }.

Data Loading
import numpy as np
import pandas as pd
import glob, os

# Load school data

data = []
for i in range(1,9):
arr = [Link](f"school{i}.dat")
[Link](arr)

ns = [Link]([len(d) for d in data])

ys = [Link]([[Link](d) for d in data])
vars_ = [Link]([[Link](d, ddof=1) for d in data])

Summary statistics are:

ȳ = (9.464, 7.033, 7.953, 6.232, 10.766, 6.205, 6.133, 7.381)

n = (25, 23, 20, 24, 24, 22, 22, 20).

(a)
We sample from the full conditionals:

θj | · ∼ N (Eθj , vθj ),
!
2 2
P P
ν σ + (y − θ )
P
ν 0 + n j 0 0 j i ij j
σ 2 | · ∼ IG , ,
2 2

3
 µ | · ∼ N (Eµ , vµ ),
!
2 2
P
η 0 + m η0 τ 0 + j (θ j − µ)
τ 2 | · ∼ IG , .
2 2

Python Code for Gibbs Sampler

[Link](1)

# Construct combined dataset Y (school_id, value)

Y_list = []
for i in range(1,9):
vals = [Link](f"school{i}.dat")
ids = np.full_like(vals, i, dtype=float)
Y_list.append(np.column_stack([ids, vals]))
Y = [Link](Y_list)

m = 8
n = [Link]([len(Y[Y[:,0]==j]) for j in range(1,9)])
ybar = [Link]([Y[Y[:,0]==j,1].mean() for j in range(1,9)])
sv = [Link]([Y[Y[:,0]==j,1].var(ddof=1) for j in range(1,9)])

# Prior parameters
mu0 = 7.0
gamma0_sq = 5.0
tau0_sq = 10.0
eta0 = 2.0
sigma0_sq = 10.0
nu0 = 2.0

# Initial values
theta = [Link]()
sigma_sq = [Link]()
mu = [Link]()
tau_sq = [Link](ddof=1)

# MCMC
S = 5000
THETA = [Link]((S,m))
MST = [Link]((S,3)) # columns: mu, sigma_sq, tau_sq
N_total = int([Link]())

for s in range(S):
# _j updates
for j in range(m):
v_theta = 1.0/(n[j]/sigma_sq + 1.0/tau_sq)
E_theta = v_theta*(n[j]*ybar[j]/sigma_sq + mu/tau_sq)
theta[j] = [Link](E_theta, [Link](v_theta))

# ² update
nu_n = nu0 + N_total
ss = nu0*sigma0_sq
for j in range(m):
yj = Y[Y[:,0]==(j+1),1]

4
 ss += [Link]((yj - theta[j])**2)
sigma_sq = 1.0/[Link](nu_n/2, 2.0/ss)

# update
v_mu = 1.0/(m/tau_sq + 1.0/gamma0_sq)
E_mu = v_mu*(m*[Link]()/tau_sq + mu0/gamma0_sq)
mu = [Link](E_mu, [Link](v_mu))

# ² update
etam = eta0 + m
ss_tau = eta0*tau0_sq + [Link]((theta-mu)**2)
tau_sq = 1.0/[Link](etam/2, 2.0/ss_tau)

THETA[s,:] = theta
MST[s,:] = [mu, sigma_sq, tau_sq]

Trace Plots

Figure 1: Trace plots for µ, σ 2 , and τ 2 .

Effective Sample Sizes

Estimated ESS:

ESS(µ) = 4164, ESS(σ 2 ) = 4789, ESS(τ 2 ) = 3995.

5
 All exceed 1000.

(b)
From the MCMC draws:

post_means = [Link](axis=0)
cis = [Link](MST, [0.025,0.975], axis=0)
print(post_means)
print(cis)

Results:

E[µ | Y ] = 7.57, 95% CI = [5.96, 9.09],

E[σ 2 | Y ] = 14.38, 95% CI = [11.68, 17.73],
E[τ 2 | Y ] = 5.48, 95% CI = [1.94, 14.56].
Posterior uncertainty for all parameters is substantially reduced relative to the priors.

(c)
S = [Link][0]

sigma_sq_pri = [Link](shape=nu0/2, scale=2*sigma0_sq/nu0, size=S)

tau_sq_pri = [Link](shape=eta0/2, scale=2*tau0_sq/eta0, size=S)

R_pri = tau_sq_pri/(sigma_sq_pri + tau_sq_pri)

R_pos = MST[:,2]/(MST[:,1] + MST[:,2])

Prior and posterior means:

E[R]prior = 0.50, E[R]post = 0.26.

6
 Figure 2: Prior and posterior densities of R.

Posterior mass shifts toward smaller values: only about 26% of total variance is between
schools.

(d)
p_theta7_lt_theta6 = [Link](THETA[:,6] < THETA[:,5])
p_theta7_smallest = [Link](THETA[:,6] <= [Link](axis=1))

Results:

P (θ7 < θ6 | Y ) = 0.524, P (θ7 = min θj | Y ) = 0.322.

j

(e)
theta_bar = [Link](axis=0)

7
 Figure 3: Posterior expectations of θj vs sample means ȳj .

Schools with extreme sample means are pulled toward the overall mean (partial pooling).
Overall sample mean:

ȳoverall = 7.691, E[µ | Y ] = 7.571.