Description
1. Metropolis: The Bounded Normal Mean.
Suppose that we have information
that the normal mean θ is bounded between −m and m, for some known number m. In this
case it is natural to elicit a prior on θ with the support on interval [−m, m].
A prior with interesting theoretical properties supported on [−m, m] is Bickel-Levit prior1
,
π(θ) = 1
m
cos2
πθ
2m
!
, − m ≤ θ ≤ m.
Assume that a sample [−2, −3, 4, −7, 0, 4] is observed from normal distribution
f(y|θ) ∝
√
τ exp
−
τ
2
(y − θ)
2
,
with a known precision τ = 1/4. Assume also that the prior on θ is Bickel-Levit, with m = 2.
This combination likelihood/prior does not result in an explicit posterior (in terms of
elementary functions). Construct a Metropolis algorithm that will sample from the posterior
of θ.
(a) Simulate 10,000 observations from the posterior, after discarding first 500 observations
(burn-in), and plot the histogram of the posterior.
(b) Find Bayes estimator of θ, and 95% equitailed Credible Set based on the simulated
observations.
Suggestions:
(i) Take uniform distribution on [−m, m] as a proposal distribution since it is easy to
sample from. This is an independence proposal, the proposed θ
′ does not depend on the
current value of the chain, θ.
1When m is large, this prior is an approximation of the least favorable distribution in a Bayes-Minimax
problem with class of priors limited to symmetric and unimodal distributions.
(ii) You will need to calculate Pn
i=1(yi − θ)
2
for current θ and Pn
i=1(yi − θ
′
)
2
for the
proposed θ
′
, prior to calculating the Metropolis ratio.
Gibbs Sampler and High/Low Protein Diet in Rats. Armitage and Berry (1994,
p. 111)2
report data on the weight gain of 19 female rats between 28 and 84 days after birth.
The rats were placed in randomized manner on diets with high (12 animals) and low (7
animals) protein content.
High protein Low protein
134 70
146 118
104 101
119 85
124 107
161 132
107 94
83
113
129
97
123
We want to test the hypothesis on dietary effect: Did a low protein diet result in a
significantly lower weight gain?
The classical t test against the one sided alternative will be significant at 5% significance
level, but we will not go in there. We will do the test Bayesian way using Gibbs sampler.
Assume that high-protein diet measurements y1i
, i = 1, . . . , 12 are coming from normal
distribution N (θ1, 1/τ1), where τ1 is the precision parameter,
f(y1i
|θ1, τ1) ∝ τ
1/2
1
exp
−
τ1
2
(y1i − θ1)
2
, i = 1, . . . , 12.
The low-protein diet measurements y2i
, i = 1, . . . , 7 are coming from normal distribution
N (θ2, 1/τ2),
f(y2i
|θ2, τ2) ∝ τ
1/2
2
exp
−
τ2
2
(y2i − θ2)
2
, i = 1, . . . , 7.
Assume that θ1 and θ2 have normal priors N (θ10, 1/τ10) and N (θ20, 1/τ20), respectively. Take
prior means as θ10 = θ20 = 110 (apriori no preference) and precisions as τ10 = τ20 = 1/100.
Assume that τ1 and τ2 have the gamma Ga(a1, b1) and Ga(a2, b2) priors with shapes
a1 = a2 = 0.01 and rates b1 = b2 = 4.
(a) Construct Gibbs sampler that will sample θ1, τ1, θ2, and τ2 from their posteriors.
2Armitage, P. and Berry, G. (1994). Statistical Methods in Medical Research (3rd edition). Blackwell
2
(b) Find sample differences θ1 − θ2 . Proportion of positive differences approximates the
posterior probability of hypothesis H0 : θ1 > θ2. What is this proportion if the number of
simulations is 10,000, with burn-in of 500?
(c) Using sample quantiles find the 95% equitailed credible set for θ1 − θ2. Does this set
contain 0?
Hint: No WinBUGS should be used (except maybe to check your results). Use Octave
(MATLAB), or R, or Python here. You may want to consult Handout GIBBS.pdf from the
course web repository.
3
