Description
HW4.1. Pairs (Xi
, Yi), i = 1, . . . , n consist of correlated standard normal random variables (mean 0, variance 1) forming a sample from a bivariate normal MVN 2(0, Σ)
distribution, with covariance matrix
Σ = ”
1 ρ
ρ 1
#
.
The density of (X, Y ) ∼ MVN 2(0, Σ) is2
f(x, y|ρ) = 1
2π
√
1 − ρ
2
exp (
−
1
2(1 − ρ
2
)
(x
2 − 2ρxy + y
2
)
)
,
with ρ as the only parameter. Take prior on ρ by assuming Jeffreys’ prior on Σ as π(Σ) =
1
|Σ|
3/2 =
1
(1−ρ
2)
3/2
, since the determinant of Σ is 1 − ρ
2
. Thus
π(ρ) = 1
(1 − ρ
2
)
3/2
1(−1 ≤ ρ ≤ 1).
(a) If (Xi
, Yi), i = 1, . . . , n are observed, write down the likelihood for ρ. Write down
the expression for the posterior, up to the proportionality constant (that is, un-normalized
posterior as the product of likelihood and prior).
(b) Since the posterior for ρ is complicated, develop a Metropolis-Hastings algorithm to
sample from the posterior. Assume that n = 100 observed pairs (Xi
, Yi) gave the following
summaries:
X
100
i=1
x
2
i = 115.9707,
X
100
i=1
y
2
i = 105.9196, and X
100
i=1
xiyi = 84.5247.
In forming a Metropolis-Hastings chain take the following proposal distribution for ρ: At
step i generate a candidate ρ
′
from the uniform U(ρi−1 − 0.1, ρi−1 + 0.1) distribution. Why
the proposal distribution cancels in the acceptance ratio expression?
(c) Simulate 51000 samples from the posterior of ρ and discard the first 1000 samples
(burn in). Plot two figures: the histogram of ρs and the realizations of the chain for the last
1000 simulations. What is the Bayes estimator of ρ based on the simulated chain?
(d) Replace the proposal distribution from (b) by the uniform U(−1, 1) (independence
proposal). Comment on the results of MCMC.
HW4.2. Exponentially distributed lifetimes have constant hazard rate equal to the
rate parameter λ. When λ is a constant hazard rate, a simple way to model heterogeneity
of hazards is to introduce a multiplicative frailty parameter µ, so that lifetimes Ti have
distribution3
Ti ∼ f(ti
|λ, µ) = λµ exp{−λµti}, ti > 0, λ, µ > 0.
The prior on (λ, µ) is
π(λ, µ) ∝ λ
c−1µ
d−1
exp{−αλ − βµ},
that is, λ and µ are apriori independent with distributions Ga(c, α) and Ga(d, β), respectively.
The hyperparameters c, d, α and β are known (elicited) and positive.
Assume that lifetimes t1, t2, . . . , tn are observed.
(a) Show that full conditionals for λ and µ are gamma,
[λ|µ, t1, . . . , tn] ∼ Ga
n + c, µXn
i=1
ti + α
!
,
and by symmetry,
[µ|λ, t1, . . . , tn] ∼ Ga
n + d, λXn
i=1
ti + β
!
.
(b) Using the result from (a) develop Gibbs Sampler algorithm that will sample 21000
pairs (λ, µ) from the posterior and burn-in the first 1000 simulations. Assume that n = 20
and that the sum of observed lifetimes is P20
i=1 ti = 522.
Assume further that the priors are specified by hyperparameters c = 3, d = 1, α = 100,
and β = 5. Start the chain with µ = 0.1.
(c) From the produced chain, plot the scatterplot of (λ, µ) as well as histograms of
individual components, λ, and µ. Estimate posterior means and variances for λ and µ. Find
95% equitailed credible sets for λ and µ.
(d) A frequentist statistician estimates the product λµ as P n
n
i=1 ti
= 20/522 = 0.0383.
What is the Bayes estimator of this product? (Hint: It is not the product of averges, it is
the average of products, so you will need to save products in the MCMC loop).
This is a toy example. In realistic applications the multiplicative frailty depends on the i, or on a subclass
of population as µj(i)
.
