Solutions to some exercises from Bayesian Data Analysis,
second edition, by Gelman, Carlin, Stern, and Rubin
4 Mar 2012
These solutions are in progress. For more information on either the solutions or the book (published by CRC), check the website, http://www.stat.columbia.edu/∼gelman/book/For each graph and some other computations, we include the code used to create it using the S
computer language. The S commands are set off from the text and appear in typewriter font.
If you find any mistakes, please notify us by e-mailing to gelman@stat.columbia.edu. Thankyou very much.
c 1996, 1997, 2000, 2001, 2003, 2004, 2006, 2007, 2009, 2010 Andrew Gelman, John Carlin, Hal
Stern, and Rich Charnigo. We also thank Jiangtao Du for help in preparing some of these solutions
and Ewan Cameron, Rob Creecy, Xin Feng, Lei Guo, Yi Lu, Pejman Mohammadi, Fei Shi, Ken
Williams, Corey Yanovsky, and Peng Yu for finding mistakes.
We have complete (or essentially complete) solutions for the following exercises:
Chapter 1: 1, 2, 3, 4, 5, 6
Chapter 2: 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 16, 18, 19, 22
Chapter 3: 1, 2, 3, 5, 9, 10
Chapter 4: 2, 3, 4
Chapter 5: 1, 2, 3, 5, 6, 7, 8, 9, 10
Chapter 6: 1, 5, 6, 7
Chapter 7: 1, 2, 7, 15
Chapter 8: 1, 2, 4, 6, 8
Chapter 11: 1, 2
Chapter 12: 6, 7
Chapter 14: 1, 3, 4, 7
Chapter 17: 1
1.1a.
p(y) ==

Pr(θ = 1)p(y|θ = 1) + Pr(θ = 2)p(y|θ = 2)
0.5N(y|1, 22 ) + 0.5N(y|2, 22 ).
y <- seq(-7,10,.02)dens <- 0.5*dnorm(y,1,2) + 0.5*dnorm(y,2,2)
plot (y, dens, ylim=c(0,1.1*max(dens)),
type="l", xlab="y", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)

-5

0

5

10

y

1.1b.
Pr(θ = 1|y = 1) ====p(θ = 1 & y = 1)p(θ = 1 & y = 1) + p(θ = 2 & y = 1)
Pr(θ = 1)p(y = 1|θ = 1)
Pr(θ = 1)p(y = 1|θ = 1) + Pr(θ = 2)p(y = 1|θ = 2)
0.5N(1|1, 22 )0.5N(1|1, 22 ) + 0.5N(1|2, 22 )0.53.
1

1.1c. As σ → ∞, the posterior density for θ approaches the prior (the data contain no information):Pr(θ = 1|y = 1) → 12 . As σ → 0, the posterior density for θ becomes concentrated at 1: Pr(θ =1|y = 1) → 1.
1.2. (1.7): For each component ui , the univariate result (1.6) states that E(ui ) = E(E(ui |v)); thus,
E(u) = E(E(u|v)), componentwise.
(1.8): For diagonal elements of var(u), the univariate result (1.8) states that var(ui ) = E(var(ui |v))+
var(E(ui |v)). For off-diagonal elements,
E[cov(ui , uj |v)] + cov[E(ui |v), E(uj |v)]===

E[E(ui uj |v) − E(ui |v)E(uj |v)] + E[E(ui |v)E(uj |v)] − E[E(ui |v)]E[E(uj |v)]
E(ui uj ) − E[E(ui |v)E(uj |v)] + E[E(ui |v)E(uj |v)] − E[E(ui |v)]E[E(uj |v)]
E(ui uj ) − E(ui )E(uj ) = cov(ui , uj ).

1.3. Note: We will use “Xx” to indicate all heterozygotes (written as “Xx or xX” in the Exercise).
Pr(child is Xx | child has brown eyes & parents have brown eyes)===

0 · (1 − p)4 + 21 · 4p(1 − p)3 + 12 · 4p2 (1 − p)2
1 · (1 − p)4 + 1 · 4p(1 − p)3 + 34 · 4p2 (1 − p)22p(1 − p) + 2p2(1 − p)2 + 4p(1 − p) + 3p2
2p.
1 + 2pTo figure out the probability that Judy is a heterozygote, use the above posterior probability asa prior probability for a new calculation that includes the additional information that her n children
are brown-eyed (with the father Xx):2p
3 n
1+2p · 4Pr(Judy is Xx | n children all have brown eyes & all previous information) = 2p
.
3 n
1
·1
+ 1+2p
1+2p · 4

Given that Judy’s children are all brown-eyed, her grandchild has blue eyes only if Judy’s childis Xx. We compute this probability, recalling that we know the child is brown-eyed and we know
Judy’s spouse is a heterozygote:
Pr(Judy’s child is Xx | all the given information)= Pr((Judy is Xx & Judy’s child is Xx) or (Judy is XX & Judy’s child is Xx) | all the given information)

 
 
2p
3 n
1
2
1
1+2p · 4
1+2p


=
+
.
2p
2p
3 n
3 n
1
1
3
2
·
·
+
+
1+2p
4
1+2p
1+2p
4
1+2p

Given that Judy’s child is Xx, the probability of the grandchild having blue eyes is 0, 1/4, or 1/2,
if Judy’s child’s spouse is XX, Xx, or xx, respectively. Given random mating, these events haveprobability (1 − p)2 , 2p(1 − p), and p2 , respectively, and soPr(Grandchild is xx | all the given information)



3 n
1
2 2p
+ 21 1+2p
1 2
1
3 1+2p · 4

=
2p(1
−
p)
+
p
2p
3 n
1
4
2
+ 1+2p
1+2p · 4

n


3
1
2 2p
+ 21 1+2p
1
3 1+2p · 4

p .
=
n
2p
3
1
2
+ 1+2p
1+2p · 42

1.4a. Use relative frequencies: Pr(A|B) = # of cases of A and B .
# of cases of B
Pr(favorite wins | point spread=8) =8= 0.67
12

Pr(favorite wins by at least 8 | point spread=8) =5= 0.4212

Pr(favorite wins by at least 8 | point spread=8 & favorite wins) =5= 0.63.
8

1.4b. Use the normal approximation for d = (score differential − point spread): d ∼ N(0, 13.862).
Note: “favorite wins” means “score differential > 0”; “favorite wins by at least 8” means “score
differential ≥ 8.”
8.5
Pr(favorite wins | point spread=8) = Φ= 0.730
13.86


0.5
= 0.514Pr(favorite wins by at least 8 | point spread=8) = Φ
13.86
0.514
Pr(favorite wins by at least 8 | point spread=8 & favorite wins) == 0.70.
0.730
Note: the values of 0.5 and 8.5 in the above calculations are corrections for the discreteness ofscores (the score differential must be an integer). The notation Φ is used for the normal cumulative
distribution function.
1.5a. There are many possible answers to this question. One possibility goes as follows. We knowthat most Congressional elections are contested by two candidates, and that each candidate typicallyreceives between 30% and 70% of the vote. For a given Congressional election, let n be the total
number of votes cast and y be the number received by the candidate from the Democratic party. If
we assume (as a first approximation, and with no specific knowledge of this election), that y/n is
uniformly distributed between 30% and 70%, then
 1
if n is even
0.4n.
Pr(election is tied|n) = Pr(y = n/2) =0
if n is odd
If we assume that n is about 200,000, with a 1/2 chance of being even, then this approximation
1
gives Pr(election is tied) ≈ 160,000.
A national election has 435 individual elections, and so the probability of at least one of them
being tied, in this analysis, is (assuming independence, since we have no specific knowledge about
the elections),Pr(at least one election is tied) = 1 − 1 −1160,000

435

≈

435
≈ 1/370.
160,000

A common mistake here is to assume an overly-precise model such as y ∼ Bin(n, 1/2). As in the
football point spreads example, it is important to estimate probabilities based on observed outcomes
rather than constructing them from a theoretical model. This is relevant even in an example such as
this one, where almost no information is available. (In this example, using a binomial model implies
that almost all elections are extremely close, which is not true in reality.)
1.5b. An empirical estimate of the probability that an election will be decided within 100 votes is49/20,597. The event that an election is tied is (y = n/2) or, equivalently, |2y −n| = 0; and the eventthat an election is decided within 100 votes is |y−(n−y)| ≤ 100 or, equivalently, |2y−n| ≤ 100. Now,3(2y−n) is a random variable that can take on integer values. Given that n is so large (at least 50,000),
and that each voter votes without knowing the outcome of the election, it seems that the distribution
1of (2y − n) should be nearly exactly uniform near 0. Then Pr(|2y − n| = 0) = 201
Pr(|2y − n| ≤ 100),
1
49
and we estimate the probability that an election is tied as 201 20,597 . As in 1.5a, the probability that49
1
any of 435 elections will be tied is then approximately 435 201
20,597 ≈ 1/190.
(We did not make use of the fact that 6 elections were decided by fewer than 10 votes, because
it seems reasonable to assume a uniform distribution over the scale of 100 votes, on which more
information is available.)
1.6. First determine the unconditional probabilities:
Pr(identical twins & twin brother) =Pr(fraternal twins & twin brother) =

1 1
·
2 300
1 1
.
Pr(fraternal twins) Pr(both boys | fraternal twins) = ·4 125

Pr(identical twins) Pr(both boys | identical twins) =The conditional probability that Elvis was an identical twin isPr(identical twins | twin brother) ===Pr(identical twins & twin brother)Pr(twin brother)
1
1
2 · 300
1
1
1
1
2 · 300 + 4 · 1255.
112.1. Prior density:p(θ) ∝ θ3 (1 − θ)3 .
Likelihood:Pr(data|θ) =
=

 
 
 
10
10
10 2
(1 − θ)10 +θ(1 − θ)9 +θ (1 − θ)80
1
2(1 − θ)10 + 10θ(1 − θ)9 + 45θ2 (1 − θ)8 .

Posterior density:p(θ|data) ∝ θ3 (1 − θ)13 + 10θ4 (1 − θ)12 + 45θ5 (1 − θ)11 .

0.0

0.2

0.4

0.6

0.8

1.0

theta <- seq(0,1,.01)
dens <- theta^3*(1-theta)^13 + 10*theta^4*(1-theta)^12 +45*theta^5*(1-theta)^11
plot (theta, dens, ylim=c(0,1.1*max(dens)),
type="l", xlab="theta", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)

theta

2.2. If we knew the coin that was chosen, then the problem would be simple: if a coin has probability π of landing heads, and N is the number of additional spins required until a head, thenE(N |π) = 1 · π + 2 · (1 − π)π + 3 · (1 − π)2 π + · · · = 1/π.

4Let T T denote the event that the first two spins are tails, and let C be the coin that was chosen.
By Bayes’ rule,
Pr(C = C1 |T T ) ==Pr(C = C1 ) Pr(T T |C = C1 )Pr(C = C1 ) Pr(T T |C = C1 ) + Pr(C = C2 ) Pr(T T |C = C2 )16
0.5(0.4)2=.
0.5(0.4)2 + 0.5(0.6)2
52

The posterior expectation of N is thenE(N |T T ) ===E[E(N |T T, C)|T T ]Pr(C = C1 |T T )E(N |C = C1 , T T ) + Pr(C = C2 |T T )E(N |C = C2 , T T )
16 1
36 1
+
= 2.24.
52 0.6 52 0.4

2.3a. E(y) = 1000( 61 ) = 166.7, and sd(y) =q1000( 61 )( 56 ) = 11.8. Normal approximation:y <- seq(120,220,.5)
dens <- dnorm (y, 1000*(1/6), sqrt(1000*(1/6)*(5/6)))
plot (y, dens, ylim=c(0,1.1*max(dens)),
type="l", xlab="y", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)

120

140

160

180

200

220

y

2.3b. From normal approximation:
5% point is 166.7 − 1.65(11.8) = 147.2

25% point is 166.7 − 0.67(11.8) = 158.8
50% point is 166.7
75% point is 166.7 + 0.67(11.8) = 174.6
95% point is 166.7 + 1.65(11.8) = 186.1

Since y is discrete, round off to the nearest integer: 147, 159, 167, 175, 186.
2.4a.
1
has mean 83.3 and sd 8.7
121
y | θ = has mean 166.7 and sd 11.8
61
y | θ = has mean 250 and sd 13.7.
4

y|θ =

The distribution for y is a mixture of the three conditional distributions:

50

100

150

200
y

250

300

y <- seq(50,300,1)
dens <- function (x, theta){
dnorm (x, 1000*theta, sqrt(1000*theta*(1-theta)))}dens.mix <- 0.25*dens(y,1/12) + 0.5*dens(y,1/6) +0.25*dens(y,1/4)
plot (y, dens.mix, ylim=c(0,1.1*max(dens.mix)),
type="l", xlab="y", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)5

2.4b. Because the three humps of the distribution have very little overlap, 14 of the distribution ofy is in the first hump, 12 is in the second hump, and 41 is in the third hump.
1)): 83.3 − (0.84)8.7 = 75.9,
The 5% point of p(y) is the 20% point of the first hump (p(y | θ = 12
round to 76. (−0.84 is the 20% point of the standard normal distribution.)
The 25% point of p(y) is between the first and second humps (approximately 120, from thegraph).
The 50% point of p(y) is at the middle of the second hump: 166.7, round to 167.
The 75% point of p(y) is between the second and third humps (approximately 205 or 210, fromthe graph).
The 95% point of p(y) is the 80% point of the first hump: 250 + (0.84)13.7 = 261.5, round to
262.
2.5a.
Pr(y = k) =Z 1

Pr(y = k | θ)dθZ 1 n k=
θ (1 − θ)n−k dθk0 n Γ(k + 1)Γ(n − k + 1)=
Γ(n + 2)
k
1
=
.
n+1
0

(1)
(2)
(3)

R1
To go from (1) to (2), use the identity 0 θα−1 (1 − θ)β−1 dθ = Γ(α)Γ(β)Γ(α+β) ; that is, the beta density has
an integral of 1. To go from (2) to (3), use the fact that Γ(x) = (x − 1)!.
α+yα
2.5b. Posterior mean is α+β+n. To show that it lies between α+β
and ny , we will write it asα+yyα
α+β+n = λ α+β + (1 − λ) n , and show that λ ∈ (0, 1). To do this, solve for λ:


y
y
α
α+y
=
+λ
−
α+β+n
n
α+β
n


α+y
y
y
α
−
= λ
−
α+β+n nα+βnnα − αy − βynα − αy − βy= λ(α + β + n)n(α + β)n
α+β
λ =
,
α+β+n

which is always between 0 and 1. So the posterior mean is a weighted average of the prior mean andthe data.
1
2.5c. Uniform prior distribution: α = β = 1. Prior variance is (α+β)2αβ(α+β+1) = 12 .

(1 + y)(1 + n − y)(2 + n)2 (3 + n)




1+y
1+n−y
1
=
.
2+n
2+n
3+n

Posterior variance =(4)

The first two factors in (4) are two numbers that sum to 1, so their product is at most 41 . And, since1n > 1, the third factor is less than 31 . So the product of all three factors is less than 12.
6

2.5d. There is an infinity of possible correct solutions to this exercise. For large n, the posteriorvariance is definitely lower, so if this is going to happen, it will be for small n. Try n = 1 and y = 1(1 success in 1 try). Playing around with low values of α and β, we find: if α = 1, β = 5, then priorvariance is 0.0198, and posterior variance is 0.0255.
2.7a. The binomial can be put in the form of an exponential family with (using the notation ofθthelast part of Section 2.4) f (y) = ny , g(θ) = (1−θ)n , u(y) = y and natural parameter φ(θ) = log( 1−θ).
A uniform prior density on φ(θ), p(φ) ∝ 1 on the entire real line, can be transformed to give theprior density for θ = eφ /(1 + eφ ):
 φ 


e
d
θ
q(θ) = p∝ θ−1 (1 − θ)−1 .
log
1 + eφdθ1−θ
2.7b. If y = 0 then p(θ|y) ∝ θ−1 (1 − θ)n−1 which has infinite integral over any interval near θ = 0.
Similarly for y = n at θ = 1.
2.8a.
θ|y ∼ Nn1
402 180 + 202 150
, 1
n
n
1
402 + 202
402 + 202

 1

2.8b.
ỹ|y ∼ N 1

n1
402 180 + 202 150
, 1
1
n
n
+
+
402
202
402
202

+ 202



2.8c.
95% posterior interval for θ | y = 150, n = 10:

95% posterior interval for ỹ | y = 150, n = 10:

150.7 ± 1.96(6.25) = [138, 163]

150.7 ± 1.96(20.95) = [110, 192]

2.8d.
95% posterior interval for θ | y = 150, n = 100:95% posterior interval for ỹ | y = 150, n = 100:

[146, 154]
[111, 189]

2.9a. From (A.3) on p. 582:
α+β

=

α =
β

0.0

0.2

0.4

0.6
theta

0.8

1.0

=

E(θ)(1 − E(θ))− 1 = 1.67
var(θ)(α + β)E(θ) = 1(α + β)(1 − E(θ)) = 0.67
theta <- seq(0,1,.001)
dens <- dbeta(theta,1,.67)
plot (theta, dens, xlim=c(0,1), ylim=c(0,3),
type="l", xlab="theta", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)
lines (c(1,1),c(0,3),col=0)lines (c(1,1),c(0,3),lty=3)

The density blows up at θ = 1 but has a finite integral.
72.9b. n = 1000, y = 650. Posterior distribution is p(θ|y) = Beta(α+650, β+350) = Beta(651, 350.67).
The data dominate the prior distribution. E(θ|y) = 0.6499, sd(θ|y) = 0.015.

0.60

0.62

0.64
0.66
theta

0.68

theta <- seq(0,1,.001)
dens <- dbeta(theta,651,350.67)
cond <- dens/max(dens) > 0.001
plot (theta[cond], dens[cond],
type="l", xlab="theta", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)0.70

2.10a.
p(data|N ) =p(N |data) ∝1
N

if N ≥ 203
otherwise0p(N )p(data|N )1(0.01)(0.99)N −1 for N ≥ 203
N
1(0.99)N for N ≥ 203.
N=∝
2.10b.

1(0.99)N .
N
P
We need to compute the normalizing constant, c.
N p(N |data) = 1, sop(N |data) = c

∞
X
1
1
=
(0.99)N .
c
N
N =203

P202
P∞This sum can be computed analytically (as N =0 N1 (0.99)N − N =0 N1 (0.99)N ), but it is easier todo the computation numerically on the computer (the numerical method is also more general andcan be applied even if the prior distribution does not have a simple, analytically-summable form).
Approximation on the computer:
Error in the approximation:

1000
X

N =203∞
X

N =10011(0.99)NN=

0.046581(0.99)N
N

<

1
1001

=
=
So 1c = 0.04658 and c = 21.47 (to a good approximation).
E(N |data) ==∞X

N =203∞Xc

N p(N |data)(0.99)N

N =203=21.47(0.99)203
1 − 0.99
8∞X(0.99)N

N =10011001

1 (0.99)
1001 1 − 0.99
4.3 × 10−6 (very minor).

=sd(N |data) =
≈
=

279.1
v
u ∞
u X
1
t
(N − 279.1)2 c (0.99)NN
N =203v
u 1000u X
1t(N − 279.1)2 21.47 (0.99)NN
N =20379.6.

2.10c. Many possible solutions here (see Jeffreys, 1961; Lee, 1989; and Jaynes, 1996). One ideathat does not work is the improper discrete uniform prior density on N
P:∞p(N ) ∝ 1. This density
leads to an improper posterior density: p(N ) ∝ N1 , for N ≥ 203. ( N =203 (1/N ) =P ∞.) Theprior density p(N ) ∝ 1/N is improper, but leads to a proper prior density, because N 1/N 2 isconvergent.
Note also that:• If more than one data point is available (that is, if more than one cable car number is observed),
then the posterior distribution is proper under all the above prior densities.
• With only one data point, perhaps it would not make much sense in practice to use a noninformative prior distribution here.

normalized density
0.4
0.8
1.2
0.0

2.11a.

dens <- function (y, th){dens0 <- NULL
for (i in 1:length(th))
dens0 <- c(dens0, prod (dcauchy (y, th[i], 1)))
dens0}y <- c(-2, -1, 0, 1.5, 2.5)
step <- .01
theta <- seq(step/2, 1-step/2, step)dens.unnorm <- dens(y,theta)
dens.norm <- dens.unnorm/(step*sum(dens.unnorm))plot (theta, dens.norm, ylim=c(0,1.1*max(dens.norm)),
type="l", xlab="theta", ylab="normalized density",
xaxs="i", yaxs="i", cex=2)

0.2

0.4

0.6

0.8

1.2

theta

0.8

1.0

[a common
mistake]
0.2

0.4

0.6

0.8

Note: a common error here is to forget toscale the y-axis from zero, thus yielding aplot as shown to the left. This is incorrect because it misleadingly implies thatthe density goes to zero at θ = 1. Whenplotting densities, the y-axis must extend
to zero!

thetas <- sample (theta, 1000, step*dens.norm,
replace=TRUE)
hist (thetas, xlab="theta", yaxt="n",breaks=seq(0,1,.05), cex=2)

2.11b.

0.0

0.2

0.4

0.6

0.8

1.0

theta

The histogram is jagged because there are only 1000 simulation draws.

9y6 <- rcauchy (length(thetas), thetas, 1)
hist (y6, xlab="new observation", yaxt="n",nclass=100, cex=2)

2.11c.

-200

0

200
400
new observation

600

Draws from a Cauchy distribution (or, in this case, a mixture of Cauchy distributions) do not fit
well onto a histogram. Compare to Figure 4.2 from the book.
2.12. The Poisson density function is p(y|θ) = θy e−θ /y!, and so J(θ) = E(−d2 log p(y|θ)/dθ2 |θ) =E(y/θ2 ) = 1/θ. This is an (improper) gamma density with a = 1/2 and b = 0.
2.13a. Let yi = number of fatal accidents in year i, for i = 1, . . . , 10, and θ = expected number ofaccidents in a year. The model for the data is yi |θ ∼ Poisson(θ).
Use the conjugate family of distributions for convenience. If the prior distribution for θ is
Gamma(α, β), then the posterior distribution is Gamma(α + 10y, β + 10).
Assume a noninformative prior distribution: (α, β) = (0, 0)—this should be ok since we haveenough information here: n = 10. Then the posterior distribution is θ|y ∼ Gamma(238, 10). Let ỹbe the number of fatal accidents in 1986. Given θ, the predictive distribution for ỹ is Poisson(θ).
Here are two methods of obtaining a 95% posterior interval for ỹ:• Simulation. Draw θ from p(θ|y) and ỹ from p(ỹ|θ):Computed 95% interval is [14, 35].

theta <- rgamma(1000,238)/10y1986 <- rpois(1000,theta)
print (sort(y1986)[c(25,976)])• Normalapproximation. From gamma distribution, E(θ|y) = √238/10 = 23.8, sd(θ|y) =√238/10 = 1.54. From Poisson distribution, E(ỹ|θ) = θ, sd(ỹ|θ) = θ.
From (1.6) and (1.7), the mean and variance of the posterior predictive distribution for ỹ are:E(ỹ|y) = E(E(ỹ|θ, y)|y)= E(θ|y) = 23.8
var(ỹ|y) = E(var(ỹ|θ, y)|y) + var(E(ỹ|θ, y)|y)= E(θ|y) + var(θ|y)= 26.2 = 5.122.
Normal approximation to p(ỹ|y) gives a 95% interval for ỹ of [23.8 ± 1.96(5.12)] = [13.8, 33.8].
But ỹ must be an integer, so interval containing at least 95% becomes [13, 34].
2.13b. Estimated numbers of passenger miles in each year: for 1976, (734/0.19)(100 million miles) =3.863 × 1011 miles; for 1977, (516/0.12)(100 million miles) = 4.300 × 1011 miles; and so forth:

10

Year
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985

Estimated number of passenger miles
3.863 × 1011
4.300 × 1011
5.027 × 1011
5.481 × 1011
5.814 × 1011
6.033 × 1011
5.877 × 1011
6.223 × 1011
7.433 × 1011
7.106 × 1011

Let xi =number of passenger miles flown in year i and θ=expected accident rate per passenger
mile. The model for the data is yi |xi , θ ∼ Poisson(xi θ).
Again use Gamma(0, 0) prior distribution for θ. Then the posterior distribution for θ isy|θ ∼ Gamma(10y, 10x) = Gamma(238, 5.716 × 1012 ).
Given θ, the predictive distribution for ỹ is Poisson(x̃θ) = Poisson(8 × 1011 θ).
Here are two methods of obtaining a 95% posterior interval for ỹ:• Simulation. Draw θ from p(θ|y) and ỹ from p(ỹ|x̃, θ):theta <- rgamma(1000,238)/5.716e12y1986 <- rpois(1000,theta*8e11)
print (sort(y1986)[c(25,976)])

Computed 95% interval is [22, 47].

(Note: your answer may differ slightly due to simulation variability. Given the wide range ofthe posterior interval, there is no practical reason to do a huge number of simulations in orderto estimate the endpoints of the interval more precisely.)• Normal approximation.
From gamma distribution, E(θ|y) = 238/(5.716 × 1012 ) = 4.132 ×√−11
× 1012 ) = 0.270 × 10−11 , From Poisson distribution, E(ỹ|θ) =10 , sd(θ|y) = 238/(5.716p
11(8 × 10 )θ, sd(ỹ|x̃, θ) = (8 × 1011 )θ.

From (1.6) and (1.7), the mean and variance of the posterior predictive distribution for ỹ are:E(ỹ|y) ==var(ỹ|y) ====E(E(ỹ|θ, y)|y)E((8 × 1011 θ|y) = (8 × 1011 )(4.164 × 10−11 ) = 33.3
E(var(ỹ|θ, y)|y) + var(E(ỹ|θ, y)|y)E((8 × 1011 )θ|y) + var((8 × 1011 )θ|y)(8 × 1011 )(4.164 × 10−11 ) + (8 × 1011 )2 (0.270 × 10−11 )2
38.0 = 6.22 .

Normal approximation to p(ỹ|y) gives a 95% interval for ỹ of [33.3 ± 1.96(6.2)] = [21.1, 45.5].
But ỹ must be an integer, so interval containing at least 95% becomes [21, 46].
The 95% interval for (b) is higher than for (a) because the model in (a) has constant expected
accident rate per year, whereas the model in (b) has expected accident rate that is proportional topassenger miles, which are increasing over time.
2.13c. Repeat analysis from (a), replacing 238 by 6919, the total number of deaths in the data.
From 1000 simulation draws: 95% posterior interval is [638, 750] deaths.

11

2.13d. Repeat analysis from (b), replacing 238 by 6919. From 1000 simulation draws: 95% posterior interval is [900, 1035] deaths.
2.13e. Just based on general knowledge, without specific reference to the data in Table 2.2, thePoisson model seems more reasonable with rate proportional to passenger miles (models (b) and(d)), because if more miles are flown in a year, you would expect more accidents. On the other
hand, if airline safety is improving, we would expect the accident rate to decline over time, so maybe
the expected number of accidents would remain roughly constant (as the number of passenger miles
flown is gradually increasing). The right thing to do here is to put a time trend in the model, as in
Exercise 3.12.
Also just based on general knowledge, we would expect accidents to be independent. But passenger deaths are not independent—they occur in clusters when airplanes crash. So the Poisson model
should be more reasonable for accidents (models (a) and (b)) than for total deaths. What if you
were interested in total deaths? Then it might be reasonable to set up a compound model: a Poisson
distribution for accidents, and then another distribution for deaths given accidents.
In addition, these models can be checked by comparing to data (see Exercise 6.2).
2.16a.
p(y) =

Zp(y|θ)p(θ)dθZ 1 
Γ(α + β) α−1n y=θ (1 − θ)n−yθ
(1 − θ)β−1 dθyΓ(α)Γ(β)
0Z
Γ(α + β) 1 y+α−1
Γ(n + 1)θ
(1 − θ)n−y+β−1 dθ=Γ(y + 1)Γ(n − y + 1) Γ(α)Γ(β) 0Γ(α + β) Γ(y + α)Γ(n − y + β)Γ(n + 1),=
Γ(y + 1)Γ(n − y + 1) Γ(α)Γ(β)Γ(n + α + β)which is the beta-binomial density (see pp. 576–577).

2.16b. We are interested only in seeing whether p(y) varies with y, so we need only look at
the factors in p(y) that depend on y; thus, we investigate under what circumstances the quantityΓ(a + y)Γ(b + n − y)/Γ(y + 1)Γ(n − y + 1) changes with y.
The preceding expression evaluates to 1 if a = b = 1; hence p(y) is constant in y if a = b = 1.
On the other hand, suppose that p(y) is constant in y. Then, in particular, p(0) = p(n) andp(0) = p(1). The first equality implies that Γ(a)Γ(b + n)/Γ(1)Γ(n + 1) = Γ(a + n)Γ(b)/Γ(n + 1)Γ(1),so that Γ(a)Γ(b + n) = Γ(a + n)Γ(b).
Recalling the formula Γ(t) = (t − 1)Γ(t − 1), we must have Γ(a)Γ(b)(b + n − 1)...(b + 1)(b) =Γ(a)Γ(b)(a + n − 1)...(a + 1)(a), which implies that (b + n − 1)...(b + 1)(b) = (a + n − 1)...(a + 1)(a).
Since each term in each product is positive, it follows that a = b.
Continuing, p(0) = p(1) implies Γ(a)Γ(b + n)/Γ(1)Γ(n + 1) = Γ(a + 1)Γ(b + n − 1)/Γ(2)Γ(n);again using Γ(t) = (t − 1)Γ(t − 1), we have b + n − 1 = na. Since a = b, this reduces to a + n − 1 = na,
whose only solution is a = 1.
Therefore a = b = 1 is a necessary as well as sufficient condition for p(y) to be constant in y.
2.18. Uniform prior distribution on θ means that the posterior distribution is θ|y ∼ N(17.653, 0.122/9).
Then
ỹ|y

∼ N(17.653, 0.122/9 + 0.122)= N(17.653, 0.1262).

99% tolerance interval is [17.653 ± 2.58(0.126)] = [17.33, 17.98].
12

√
√ √
2.19a. Let u = σ 2 , so that σ = u. Then p(σ 2 ) = p(u) = p( u)d u/du = p(σ)(1/2)u−1/2 =(1/2)p(σ)/σ, which is proportional to 1/σ 2 if p(σ) ∝ σ −1 .
2.19b. Proof by contradiction: The posterior density p(σ|data) is proportional to σ −1−n exp(−c/σ 2 ),which may be written as (σ 2 )−1/2−n/2 exp(−c/σ 2 ). The posterior density p(σ 2 |data) is proportional
to (σ 2 )−1−n/2 exp(−c/σ 2 ). In the preceding expressions, we have defined c = nv/2, and we assume
this quantity√ √to be positive.
Let ( a, b) be the 95% interval of highest density for p(σ|data). (It can be shown using calculusthat the density function is unimodal: there is only one value of σ for which dp(σ|data)/dσ = 0.
Therefore the 95% region of highest density is a single interval.) Thena−1/2−n/2 exp(−c/a) = b−1/2−n/2 exp(−c/b).
Equivalently, (−1/2 − n/2) log a − c/a = (−1/2 − n/2) log b − c/b.
Now, if (a, b) were the 95% interval of highest density for p(σ 2 |data), thena−1−n/2 exp(−c/a) = b−1−n/2 exp(−c/b).
That is, (−1 − n/2) log a − c/a = (−1 − n/2) log b − c/b. Combining the two equations, 1/2 log a =1/2 log b, so that a = b, in which case [a, b] cannot be a 95% interval, and we have a contradiction.
2.22a.
p(θ | y ≥ 100) ∝∝p(θ | y ≥ 100) =p(y ≥ 100 | θ)p(θ)

exp(−100θ)θα−1 exp(−βθ)
Gamma(θ | α, β + 100)ααand (β+100)
The posterior mean and variance of θ are β+100
2 , respectively.

2.22b.
p(θ | y = 100) ∝ p(y = 100 | θ)p(θ)∝ θ exp(−100θ)θα−1 exp(−βθ)

p(θ | y = 100) = Gamma(θ | α + 1, β + 100)α+1α+1and (β+100)
The posterior mean and variance of θ are β+100
2 , respectively.

2.22c. Identity (2.8) says on average, the variance of θ decreases given more information: in thiscase,E(var(θ|y) | y ≥ 100) ≤ var(θ | y ≥ 100).
(5)
Plugging in y = 100 to get var(θ|y = 100) is not the same as averaging over the distribution ofy | y ≥ 100 on the left side of (5).
3.1a Label the prior distribution p(θ) as Dirichlet(a1 , . . . , an ). Then the posterior distributionis p(θ|y) = Dirichlet(y1 + a1 , . . . , yn + an ). From the properties of the Dirichlet distribution (seeAppendix A), the marginal posterior distribution of (θ1 , θ2 , 1 − θ1 − θ2 ) is also Dirichlet:p(θ1 , θ2 |y) ∝ θ1y1 +a1 −1 θ2y2 +a2 −1 (1−θ1 −θ2 )yrest +arest −1 , where yrest = y3 +. . .+yJ , arest = a3 +. . .+aJ .
(This can be proved using mathematical induction, by first integrating out θn , then θn−1 , and soforth, until only θ1 and θ2 remain.)13

1
Now do a change of variables to (α, β) = ( θ1θ+θ, θ1 + θ2 ). The Jacobian of this transformation2
is |1/β|, so the transformed density isp(α, β|y) ∝=∝β(αβ)y1 +a1 −1 ((1 − α)β)y2 +a2 −1 (1 − β)yrest +arest −1

αy1 +a1 −1 (1 − α)y2 +a2 −1 β y1 +y2 +a1 +a2 −1 (1 − β)yrest +arest −1
Beta(α|y1 + a1 , y2 + a2 )Beta(β|y1 + y2 + a1 + a2 , yrest + arest ).

Since the posterior density divides into separate factors for α and β, they are independent, and, as
shown above, α|y ∼ Beta(y1 + a1 , y2 + a2 ).
3.1b. The Beta(y1 + a1 , y2 + a2 ) posterior distribution can also be derived from a Beta(a1 , a2 ) priordistribution and a binomial observation y1 with sample size y1 + y2 .
3.2. Assume independent uniform prior distributions on the multinomial parameters. Then theposterior distributions are independent multinomial:(π1 , π2 , π3 )|y∼ Dirichlet(295, 308, 39)(π1∗ , π2∗ , π3∗ )|y1, α2 =and α1 = π1π+π2π1∗
π1∗ +π2∗∼ Dirichlet(289, 333, 20),

. From the properties of the Dirichlet distribution (see Exercise 3.1),
α1 |y∼ Beta(295, 308)

α2 |y∼ Beta(289, 333).

The histogram of 1000 draws from the posterior density of α2 − α1 is attached. Based on thishistogram, the posterior probability that there was a shift toward Bush is 19%.

-0.10

-0.05
0.0
alpha_2 - alpha_1

alpha.1 <- rbeta (2000, 295, 308)
alpha.2 <- rbeta (2000, 289, 333)
dif <- alpha.2 - alpha.1
hist (dif, xlab="alpha_2 - alpha_1", yaxt="n",breaks=seq(-.12,.08,.01), cex=2)
print (mean(dif>0))0.05

Essentially the same answer may be obtained using normal approximations for the distributionsof α2 and α1 with means and standard deviations computed from the relevant beta distributions.
3.3a. Data distribution is p(y|µc , µt , σc , σt ) =distribution is

Q32

2
i=1 N(yc i |µc , σc )Q36

2i=1 N(yt i |µt , σt ).

Posteriorp(µc , µt , log σc , log σt |y) = p(µc , µt , log σc , log σt )p(y|µc , µt , log σc , log σt )
=

32
Y

i=1

|

N(yc i |µc , σc2 ){zp(µc ,log σc |y)

36
Y

i=1

}|

N(yt i |µt , σt2 ) .
{zp(µt ,log σt |y)}The posterior density factors, so (µc , σc ) are independent of (µt , σt ) in the posterior distribution.
So, under this model, we can analyze the two experiments separately. Using the results from Section
3.2, the marginal posterior distributions for µc and µt are:
µc |y
µt |y

∼
∼

t31 (1.013, 0.242/32)| {z }
0.042422

t35 (1.173, 0.20 /36).
| {z }
0.0333214

3.3b The histogram of 1000 draws from the posterior density of µt − µc is attached. Based on thishistogram, a 95% posterior interval for the average treatment effect is [0.05, 0.27].

-0.1

0.0

0.1
0.2
mu_t - mu_c

0.3

0.4

mu.c <- 1.013 + (0.24/sqrt(32))*rt(1000,31)
mu.t <- 1.173 + (0.20/sqrt(36))*rt(1000,35)
dif <- mu.t - mu.c
hist (dif, xlab="mu_t - mu_c", yaxt="n",breaks=seq(-.1,.4,.02), cex=2)
print (sort(dif)[c(25,976)])

3.5a. If the observations are treated as exact measurements then the theory of Section 3.2 applies.
Then p(σ 2 |y) ∼ Inv-χ2 (n − 1, s2 ) and p(µ|σ 2 , y) ∼ N (ȳ, σ 2 /n) with n = 5, ȳ = 10.4, s2 = 1.3.
3.5b. The posterior distribution assuming the non-informative prior distribution p(µ, σ 2 ) ∝ 1/σ 2
is



n  
1 Y
yi − 0.5 − µyi + 0.5 − µp(µ, σ 2 |y) ∝ 2
−Φ
,
Φ
σ i=1
σ
σ
where Φ is the standard normal cumulative distribution function.

3.5c. We computed contour plots on a grid on the (µ, log σ) scale (levels are 0.0001, 0.001, 0.01,
0.05–0.95). Draws from the posterior density in 3.5a can be obtained directly as described in Section
3.2. Draws from the posterior density in 3.5b were obtained by sampling from the grid approximation(again using the (µ, log σ) scale).
mean
sd 2.5% 25% 50% 75% 97.5%
Ignoring roundingµ10.4 0.74
8.9 10.1 10.4 10.8
11.8
σ
1.4 0.74
0.68 0.98
1.2
1.6
3.3
Accounting for roundingµ10.5 0.69
9.2 10.1 10.4 10.8
11.8
σ
1.4 0.73
0.61 0.90
1.2
1.6
3.3
The only notable difference is in the lower quantiles of σ.

15

post.a <- function(mu,sd,y){ldens <- 0
for (i in 1:length(y)) ldens <- ldens +log(dnorm(y[i],mu,sd))
ldens}
post.b <- function(mu,sd,y){ldens <- 0
for (i in 1:length(y)) ldens <- ldens +log(pnorm(y[i]+0.5,mu,sd) - pnorm(y[i]-0.5,mu,sd))
ldens}
summ <- function(x){c(mean(x),sqrt(var(x)),
quantile(x, c(.025,.25,.5,.75,.975)))}-2-1

log sigma
0
1
2

3

4

Posterior density, ignoring rounding

5

10
mu

15

-2

-1

log sigma
0
1
2

3

4

Posterior density, accounting for rounding

5

10
mu

15

nsim <- 2000
y <- c(10,10,12,11,9)n <- length(y)ybar <- mean(y)s2 <- sum((y-mean(y))^2)/(n-1)
mugrid <- seq(3,18,length=200)
logsdgrid <- seq(-2,4,length=200)
contours <- c(.0001,.001,.01,seq(.05,.95,.05))
logdens <- outer (mugrid, exp(logsdgrid), post.a, y)dens <- exp(logdens - max(logdens))
contour (mugrid, logsdgrid, dens, levels=contours,
xlab="mu", ylab="log sigma", labex=0, cex=2)
mtext ("Posterior density, ignoring rounding", 3)
sd <- sqrt((n-1)*s2/rchisq(nsim,4))mu <- rnorm(nsim,ybar,sd/sqrt(n))
print (rbind (summ(mu),summ(sd)))
logdens <- outer (mugrid, exp(logsdgrid), post.b, y)dens <- exp(logdens - max(logdens))
contour (mugrid, logsdgrid, dens, levels=contours,
xlab="mu", ylab="log sigma", labex=0, cex=2)
mtext ("Posterior density, accounting for rounding",
cex=2, 3)dens.mu <- apply(dens,1,sum)muindex <- sample (1:length(mugrid), nsim, replace=T,
prob=dens.mu)
mu <- mugrid[muindex]
sd <- rep (NA,nsim)
for (i in (1:nsim)) sd[i] <- exp (sample(logsdgrid, 1, prob=dens[muindex[i],]))
print (rbind (summ(mu),summ(sd)))

3.5d. Given µ and σ, the conditional distribution of zi is N(µ, σ 2 ), truncated to fall in the range(yi − 0.5, yi + 0.5). So, for each of our posterior draws of (µ, σ), we draw a vector (z1 , . . . , z5 ) fromindependent truncated normal distributions using the inverse-cdf method.
The posterior mean of(z1 − z2 )2 is 0.16.

z <- matrix (NA, nsim, length(y))
for (i in 1:length(y)){
lower <- pnorm (y[i]-.5, mu, sd)
upper <- pnorm (y[i]+.5, mu, sd)z[,i] <- qnorm (lower + runif(nsim)*(upper-lower), mu, sd)}
mean ((z[,1]-z[,2])^2)

16

3.9.
p(µ, σ 2 |y) ∝ p(y|µ, σ 2 )p(µ, σ 2 )(n − 1)s2 + n(µ − ȳ)2ν0 σ02 + κ0 (µ − µ0 )2−1 2 −(ν0 /2+1)∝ (σ 2 )−n/2 exp −
σ
(σ
)
exp
−
2σ 2
2σ 2


nκ0 (ȳ−µ0 )2µ0 κ0 +nȳ 2
2
2
+
(n
+
κ
)(µ
−
)
ν
σ
+
(n
−
1)s
+
0
0
0
n+κ0
n+κ0

∝ σ −1 (σ 2 )−((ν0 +n)/2+1) exp −2σ 2
And so,2µ, σ |y ∼ N-Inv-χ
where
σn2 =

2




µ0 κ0 + nȳ
σn22,; n + ν0 , σn ,n + κ0n + κ0(ȳ−µ0 )ν0 σ02 + (n − 1)s2 + nκ0n+κ
02

n + ν0.

3.10. From equation (3.4) on page 75, p(σj2 |y) ∝ (σj2 )−n/2−1/2 exp(−(n − 1)s2 /2σj2 ) for each j.
Thus p(1/σj2 |y) ∝ (σj2 )2 (1/σj2 )n/2+1/2 exp(−(n − 1)s2 /2σj2 ) = (1/σj2 )n/2−3/2 exp(−(n − 1)s2 /2σj2 ).
Since the proportionality is in σj2 , it is also correct to say that the latter object is proportional to((n − 1)s2 /σj2 )n/2−3/2 exp(−(n − 1)s2 /2σj2 ) = ((n − 1)s2 /σj2 )(n−1)/2−1 exp(−(n − 1)s2 /2σj2 ), which
implies that (n − 1)s2 /σj2 has a χ2n−1 distribution.
The independence assumptions imply that (n1 − 1)s21 /σ12 and (n2 − 1)s22 /σ22 are independent χ2with n1 − 1 and n2 − 1 degrees of freedom, respectively, and so we invoke the well-known result thatthe quotient of two independent χ2 random variables, each divided by their degrees of freedom, hass21 /σ12
the F distribution. In particular, s2 /(σ2 has the Fn1 −1,n2 −1 distribution. A trivial rearrangement2)2
yields the desired result.
4.2.
p(yi |θ)li = log p(yi |θ)
d2 li
dα2
d2 li
dαdβ
d2 li
dβ 2

∝ (logit−1 (α + βxi ))yi (1 − logit−1 (α + βxi ))ni −yi= constant + yi log(logit−1 (α + βxi )) + (ni − yi ) log(1 − logit−1 (α + βxi ))ni exp(α + βxi )= −(1 + exp(α + βxi ))2
ni xi exp(α + βxi )= −(1 + exp(α + βxi ))2
ni x2i exp(α + βxi )= −(1 + exp(α + βxi ))2

The prior density on (α, β) is uniform, so log p(θ|y) = constant +I(θ̂) =

P4

ni exp(α+βxi )
i=1 (1+exp(α+βxi ))2P4 ni xi exp(α+βxi)
i=1 (1+exp(α+βxi ))2P4P4

i=1 li , andni xi exp(α+βxi )
i=1 (1+exp(α+βxi ))2P4 ni x2i exp(α+βxi)
i=1 (1+exp(α+βxi ))2!

,(α,β)=(α̂,β̂)a c, the normal approximation variancesc bbaof α and β are the diagonal elements of I −1 : ab−c
2 , and ab−c2 , respectively.
where θ̂ = (α̂, β̂) is the posterior mode. Denoting I as17

4.3. Let θ = LD50; we also introduce ν = β, in anticipation of a change of coordinates. Formula(4.1) on page 101 suggests that the (asymptotic) posterior median and mode should coincide, and the(asymptotic) posterior standard variance should be the inverse of observed information, evaluatedat the posterior mode.
With some effort (or using the “generalized linear model” function in a statistics package), it
is possible to obtain decent numerical estimates for the posterior mode and standard deviation
associated with this data set; we will take these values as proxies for their asymptotic analogues.
Here is one way to proceed:R
R
Observe that p(θ|y) = p(θ, ν|y)dν = p(α, β|y)|ν|dν, where α = −θν and β = ν in the last
integral, and |ν| is the Jacobian associated with the change in coordinates. (An expression forp(α, β|y) is given as (3.21) on page 90.
Fixing θ equal to a value aroundR −0.12, which appears to be near the posterior mode (see Figure
3.5 on page 93), we may compute p(α, β|y)|ν|dν numerically (up to a proportionality constant that
does not depend on θ). The region of integration is infinite, but the integrand decays rapidly afterν passes 50, so that it suffices to estimate the integral assuming that ν runs from 0 to, say, 70.
This procedure can be repeated for several values of θ near −0.12. The values may be compareddirectly to find the posterior mode for θ. To three decimal places, we obtain −0.114.
One can fix a small value of h, such as h = 0.002, and estimate the d2 /dθ2 log p(θ|y), evaluated at θ
equal to the posterior mode, by the expression [log p(−0.114+h|y)−2 log p(−0.114|y)+log p(−0.114−h|y)]/h2 (see (12.2) on page 314).
The negative of of the preceding quantity, taken to the −0.5 power, is our estimate for theposterior standard deviation. We get 0.096, which appears quite reasonable when compared to
Figure 3.5, but which seems smaller than what we might guess from Figure 4.2(b). As remarked
on page 105, however, the actual posterior distribution of LD50 has less dispersion than what issuggested by Figure 4.2.
4.4. In the limit as n → ∞, the posterior variance approaches zero—that is, the posterior distribution becomes concentrated near a single point. Any one-to-one continuous transformation on thereal numbers is locally linear in the neighborhood of that point.
5.1a. The following results are based on a different set of simulations than contained in Section
5.5. Two noteworthy differences: the support of τ for simulation purposes was increased to [0, 40],and the number of draws from the posterior distribution is 1000. Based on the 1000 posterior
simulations, we obtain
Pr(better than school)
School Pr(best)
A
B
C
D
E
F
G
H
A
0.25
−
0.64 0.67 0.66 0.73 0.69 0.53 0.61
B
0.10
0.36
−
0.55 0.53 0.62 0.61 0.37 0.49
C
0.10
0.33 0.45−
0.46 0.58 0.53 0.36 0.45
D
0.09
0.34 0.47 0.54−
0.61 0.58 0.37 0.47
E
0.05
0.27 0.38 0.42 0.39−
0.48 0.28 0.38
F
0.08
0.31 0.39 0.47 0.42 0.52−
0.31 0.40
G
0.21
0.47 0.63 0.64 0.63 0.72 0.69
−
0.60
H
0.12
0.39 0.51 0.55 0.53 0.62 0.60 0.40−
5.1b. In the model with τ set to ∞, the school effects θj are independent inqtheir posterior

distribution with θj |y ∼ N (yj , σj2 ). It follows that Pr(θi > θj |y) = Φ((yi − yj )/ σi2 + σj2 ). The
probability that θi is the largest of the school effects can be expressed as a single integral given

18

below.
Pr(θi is the largest) =Z ∞Y−∞ j6=iΦθi − θjσjφ(θi |yi , σi )dθi

This integral can be evaluated numerically (results given below) or estimated by simulating school
effects from their independent normal posterior distributions. The results are provided below.
Pr(better than school)
School Pr(best)
A
B
C
D
E
F
G
H
A
0.556
−
0.87 0.92 0.88 0.95 0.93 0.72 0.76
B
0.034
0.13
−
0.71 0.53 0.73 0.68 0.24 0.42
C
0.028
0.08 0.29−
0.31 0.46 0.43 0.14 0.27
D
0.034
0.12 0.47 0.69−
0.70 0.65 0.23 0.40
E
0.004
0.05 0.27 0.54 0.30−
0.47 0.09 0.26
F
0.013
0.07 0.32 0.57 0.35 0.53−
0.13 0.29
G
0.170
0.28 0.76 0.86 0.77 0.91 0.87
−
0.61
H
0.162
0.24 0.58 0.73 0.60 0.74 0.71 0.39−
5.1c. The model with τ set to ∞ has more extreme probabilities. For example, in the first column,
the probability that School A is the best increases from 0.25 to 0.56. It is also true in the pairwise
comparisons. For example, the probability that School A’s program is better than School E is under
the full hierarchical model is 0.73, whereas it is 0.95 under the τ = ∞ model. The more conservative
answer under the full hierarchical model reflects the evidence in the data that the coaching programs
appear fairly similar in effectiveness. Also, noteworthy is that the preferred school in a pair canchange, so that School E is better than School C when τ is set to ∞, whereas School C is better
than School E when averaging over the posterior distribution of τ . This effect occurs only because
the standard errors, σj , differ.
5.1d. If τ = 0 then all of the school effects are the same. Thus no school is better or worse thanany other.
5.2a. Yes, they are exchangeable. The joint distribution is −1 X Y
2J
J
Y
2J

N(θp(j) |1, 1)
N(θp(j) | − 1, 1) ,p(θ1 , . . . , θ2J ) =
J
p
j=1

(6)

j=J+1

where the sum is over all permutations p of (1, . . . , 2J). The density (6) is obviously invariant topermutations of the indexes (1, . . . , 2J).
5.2b. Pick any i, j. The covariance of θi , θj is negative. You can see this because if θi is large, thenit probably comes from the N(1, 1) distribution, which means that it is more likely than not that θjcomes from the N(−1, 1) distribution (because we know that exactly half of the 2J parameters comefrom each of the two distributions), which means that θj will probably be negative. Conversely, if
θi is negative, then θj is most likely positive.
Then, by Exercise 5.3, p(θ1 , . . . , θ2J ) cannot be written as a mixture of iid components.
The above argument can be made formal and rigorous by defining φ1 , . . . , φ2J , where half ofthe φj ’s are 1 and half are −1, and then setting θj |φj ∼ N(φj , 1). It’s easy to show first thatcov(φi , φj ) < 0, and then that cov(θi , θj ) < 0 also.

195.2c. In the limit as J → ∞, the negative correlation between θi and θj approaches zero, andthe joint distribution approaches iid. To put it another way, as J → ∞, the distinction disappears
between (1) independently assigning each θj to one of two groups, and (2) picking exactly half ofthe θj ’s for each group.
5.3. Let µ(φ) = E(θj |φ). From (1.8) on page 24 (also see Exercise 1.2),
cov(θi , θj )= E(cov(θi , θj |φ)) + cov(E(θi |φ), E(θj |φ))= 0 + cov(µ(φ), µ(φ))= var(µ(φ))≥ 0.
5.5a. We want to find E(y) and var(y), where y|θ ∼ Poisson(θ) and θ ∼ Gamma(α, β). From (2.7),
E(y) = E(E(y|θ)); from properties of the Poisson and Gamma distributions, we have E(y|θ) = θ sothat E(E(y|θ)) = E(θ) = α/β.
Similarly, by formula (2.8),
var(y) = E(var(y|θ)) + var(E(y|θ)) = E(θ) + var(θ) = α/β + α/β 2 = α/β 2 (β + 1).
5.5b. This part is a little bit trickier because we have to think about what to condition on in
applying formulas (2.7) and (2.8). Some reflection (and, perhaps, a glance at formula (3.3) on page75) will lead to the choice of conditioning on σ 2 .
Throughout these computations, n, s, and ȳ are essentially treated like constants.
From (2.7) and (3.3),√√√√E( n(µ−ȳ)/s|y) = E(E( n(µ−ȳ)/s|σ, y)|y) = E(( n/s)E(µ−ȳ|σ, y)|y) = E(( n/s)·0|y) = E(0|y) = 0.
Obviously we must have n > 1 for s to be defined. But in order for the expectation to exist, wemust have n > 2. You can deduce this from inspecting the formula for p(µ|y) at the top of page 69:
the exponent must be less than negative one for the quantity to be integrable, which is true if andonly if n > 2.
Similarly, we can compute from (2.8), (3.3), and (3.5) that√√√var( n(µ − ȳ)/s|y) = var(E( n(µ − ȳ)/s|σ, y)|y) + E(var( n(µ − ȳ)/s|σ, y)|y)= var(0|y) + E((n/s2 )var(µ|σ, y)|y)= E((n/s2 )σ 2 /n|y)n−1.
= E(σ 2 |y)/s2 =n−3Note that n > 3 is required here.
5.6. Let pm (θ|y) denote the posterior density of θ corresponding to the prior density pm (θ). Thatis, for each m,
Ppm (θ|y) = pm (θ)p(y|θ)/pm (y), where pm (y) is the prior predictive density.
If
p(θ)
=
(θ), then the posterior density of θ is proportional tom λm pmPPm λm pm (y)pm (θ|y): this is a mixture of the posterior densities pm (θ|y) withm λm pm (θ)p(y|θ) =weights proportional to λm pm (y). Since each pm (θ) is conjugate for the model for y given θ, thepreceding computation demonstrates that the class of finite mixture prior densities is also conjugate.
Consider an example: p1 (θ) ∼ N(1, 0.52 ), p2 (θ) ∼ N(−1, 0.52 ), and suppose that λ1 = 0.9 amd
λ2 = 0.1. (The exact choices of λ1 and λ2 are not important. What is important is how andwhy the posterior mixture weights are different than the prior weights.) We know that p1 (θ|y) ∼20

N(1.5/14, 1/14) and that p2 (θ|y) ∼ N(−6.5/14, 1/14); see, e.g., formulas (2.9) and (2.10). We also
know from the last paragraph that
P p(θ|y) will be a weighted sum of these conditional posterior
densities with weights λm pm (y)/ k λk pk (y) for m = 1, 2.
You can compute p1 (y) and p2 (y), using convenient properties of normal distributions: p1 (y) =N(−0.25|1, 0.52 + 1/10) = 0.072, and p2 (y) = N(−0.25| − 1, 0.52 + 1/10) = 0.302.
0.9·0.072
So the weights for p1 (θ|y) and p2 (θ|y) are not 0.9 and 0.1 but are, rather, 0.9·0.072+0.1·0.302= 0.680.1·0.302
and 0.9·0.072+0.1·0.302 = 0.32.
theta <- seq(-3,3,.01)
prior <- c (0.9, 0.1)
dens <- prior[1]*dnorm(theta,1,0.5) +prior[2]*dnorm(theta,-1,0.5)
plot (theta, dens, ylim=c(0,1.1*max(dens)),
type="l", xlab="theta", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)
mtext ("prior density", cex=2, 3)prior density

-3

-2

-1

0
theta

1

2

3

marg <- dnorm(-.25,c(1,-1),sqrt(c(0.5,0.5)^2+1/10))
posterior <- prior*marg/sum(prior*marg)

posterior density

−3

−2

−1

0
theta

1

2

dens <- posterior[1]*dnorm(theta,1.5/14,sqrt(1/14)) +posterior[2]*dnorm(theta,-6.5/14,sqrt(1/14))
plot (theta, dens, ylim=c(0,1.1*max(dens)),
type="l", xlab="theta", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)
mtext ("posterior density", cex=2, 3)

3

5.7a. Consider the limit (α + β) → ∞ with α/β fixed at any nonzero value. The likelihood (seeequation (5.8)) isp(y|α, β) ∝≈=JY
[α · · · (α + yj − 1)][β · · · (β + nj − yj − 1)](α + β) · · · (α + β + nj − 1)

j=1

J
Y
αyj β nj −yj

j=1(α + β)nj

J 
Y

j=1

(7)

α
α+β

yj 

β
α+β

nj −yj,which is a constant (if we are considering y, n, and α/β to be fixed), so the prior density determines
whether the posterior density has a finite integral in this limit. A uniform prior density on log(α+ β)
has an infinite integral in this limit, and so the posterior density does also in this case.
5.7b. The Jacobian of the transformation isβ(α+β)2− 21 (α + β)−3/2

α− (α+β)21− 2 (α + β)−3/2= constant · (α + β)−5/2 .

5.7c. Note error in statement of problem in the first and second printings: “improper” should be“proper” and “yj 6= 0 for at least one experiment j and yj 6= nj for at least one experiment j”
should be “0 < yj < nj for at least one experiment j.”
There are 4 limits to consider:21

1. α → 0 with α + β fixed
2. β → 0 with α + β fixed
3. α + β → 0 with α/β fixed
4. α + β → ∞ with α/β fixed
As in Exercise 5.7a, we work with expression (5.8). We have to show that the integral is finite ineach of these limits.
Let J0 be the number of experiments with yj > 0, J1 be the number of experiments with yj < nj ,and J01 be the number of experiments with 0 < yj < nj .
1. For α → 0, proving convergence is trivial. All the factors in the likelihood (7) go to constants
except αJ0 , so the likelihood goes to 0 (if J0 > 0) or a constant (if J0 = 0). The priordistribution is a constant as α goes to 0. So whether J0 = 0 or J0 > 0, the integral of theposterior density in this limit is finite.
2. For β → 0, same proof.
α
3. First transform to ( α+β, α + β), so we only have to worry about the limit in one dimension.
Multiplying the prior distribution (5.9) by the Jacobian yieldsαp, α + β ∝ (α + β)−3/2
α+βFor α + β → 0, the likelihood looks likep(y|α, β) ∝ (α + β)−J01 ,αignoring all factors such as Γ(α + β + nj ) and α+βthat are constant in this limit. So theposterior density in this parameterization is


α
(8)
, α + β y ∝ (α + β)−3/2 (α + β)J01 .
p
α+βThe function xc has a finite integral as x → 0 if c > −1, so (8) has a finite integral if J01 > 21 .

Note: the statistical reasoning here is that if J01 = 0 (so that all the data are of “0 successes
out of nj ” or “nj successes out of nj ”), then it is still possible that α = β = 0 (corresponding
to all the θj ’s being 0 or 1), and we have to worry about the infinite integral of the improper
prior distribution in the limit of (α + β) → 0. If J01 > 0, the likelihood gives the informationthat (α + β) = 0 is not possible. In any case, values of α and β near 0 do not make muchsense in the context of the problem (modeling rat tumor rates), and it might make sense tojust constrain α ≥ 1 and β ≥ 1.
4. For α + β → ∞, the likelihood is constant, and so we just need to show that the prior density
has a finite integral (see the solution to Exercise 5.7a). As above,αp, α + β ∝ (α + β)−3/2 ,α+βwhich indeed has a finite integral as (α + β) → ∞.
5.8. We first note that, since p(µ|τ, y) and p(θ|µ, τ, y) have proper distributions, the joint posterior
density p(θ, µ, τ |y) is proper if and only if the marginal posterior density p(τ |y) from (5.21) is
proper—that is, has a finite integral for τ from 0 to ∞.
22

5.8a. Everything multiplying p(τ ) in (5.21) approaches a nonzero constant limit as τ tends to zero;
call that limit C(y). Thus the behavior of the posterior density near τ = 0 is determined by theprior density. The function p(τ ) ∝ 1/τ is not integrable for any small interval including τ = 0, andso it leads to a nonintegrable posterior density. (We can make this more formal: for any δ > 0 we
can identify an interval including zero on which p(τ |y) ≥ p(τ )(C − δ).)5.8b. The argument from 5.8a shows that if p(τ ) ∝ 1 then the posterior density is integrable nearzero. We need to examine the behavior as τ → ∞ and find an upper bound that is integrable.
Theterm is clearly less than or equal to 1. We can rewrite the remaining terms as
PJ exponential
Q( j=1 [ k6=j (σk2 + τ 2 )])−1/2 . For τ > 1 we make this quantity bigger by dropping all of the σ 2 toyield (Jτ 2(J−1) )−1/2 . An upper bound on p(τ |y) for τ large is p(τ )J −1/2 /τ J−1 . When p(τ ) ∝ 1,
this upper bound is integrable if J > 2, and so p(τ |y) is integrable if J > 2.
5.8c. There are several reasonable options here. One approach is to abandon the hierarchical modeland just fit the two schools with independent noninformative prior distributions (as in Exercise 3.3
but with the variances known).
Another approach would be to continue with the uniform prior distribution on τ and try toanalytically work through the difficulties with the resulting improper posterior distribution. In thiscase, the posterior distribution has all its mass near τ = ∞, meaning that no shrinkage will occur
in the analysis (see (5.17) and (5.20)). This analysis is in fact identical to the first approach ofanalyzing the schools independently.
If the analysis based on noninformative prior distributions is not precise enough, it might beworthwhile assigning an informative prior distribution to τ (or perhaps to (µ, τ )) based on outside
knowledge such as analyses of earlier coaching experiments in other schools. However, one would haveto be careful here: with data on only two schools, inferences would be sensitive to prior assumptions
that would be hard to check from the data at hand.
5.9a. p(θ, µ, τ |y) is proportional to p(θ, µ, τ )p(y|θ, µ, τ ), where we may note that the latter density
is really just p(y|θ). Since p(θ, µ, τ ) = p(θ|µ, τ )p(µ, τ ),p(θ, µ, τ |y) ∝ p(µ, τ )

J 
Y

j=1

Y
J
 yj

1
θj−1 (1 − θj )−1 τ −1 exp(− (logit(θj ) − µ)2 /τ 2 )θj (1 − θj )nj −yj .
2
j=1(The factor θj−1 (1 − θj )−1 is d(logit(θj ))/dθj ; we need this since we want to have p(θ|µ, τ ) rather
than p(logit(θ)|µ, τ ).)5.9b. Even though we can look at each of the J integrals individually—the integrand separates
into independent factors—there is no obvious analytic technique which permits evaluation of these
integrals. In particular, we cannot recognize a multiple of a familiar density function inside theintegrals. One might try to simplify matters with a substitution like uj = logit(θj ) or vj = θj /(1−θj ),but neither substitution turns out to be helpful.
5.9c. In order for expression (5.5) to be useful, we would have to know p(θ|µ, τ, y). Knowing it up
to proportionality in θ is insufficient because our goal is to use this density to find another densitythat depends on µ and τ . In the rat tumor example (page 128, equation (5.7)), the conjugacy ofthe beta distribution allowed us to write down the relevant “constant” of proportionality, which ofcourse appears in equation (5.8).

23

5.10 Following the hint we apply (2.7) and (2.8).

1
1y + τ12 µy + τ12 µ̂
σj2 j
σj2 j


E(θj |τ, y) = E[E(θj |µ, τ, y) | τ, y] = E
τ, y =
1
1
+ τ12
+ τ12
σ2
σ2
j

var(θj |τ, y) = E[var(θj |µ, τ, y) | τ, y] + var[E(θj |µ, τ, y) | τ, y] =

j

1
1
+ τ12
σj2



+ 1

σj2

1
τ2

+ τ12

2

 Vµ ,

where expressions for µ̂ = E(µ|τ, y) and Vµ = var(µ|τ, y) are given in (5.20) of Section 5.4.

6.1a. Under the model that assumes identical effects in all eight schools (τ = 0), the posterior
distribution for the common θ is N(7.9, 4.22); see page 139–140. To simulate the posterior predictive
distribution of y rep we first simulate θ ∼ N(7.9, 4.22 ) and then simulate yjrep ∼ N(θ, σj2 ) for j =1, . . . , 8. The observed order statistics would be compared to the distribution of the order statistics
obtained from the posterior predictive simulations. The expected order statistics provided in the
problem (from a N(8, 132 ) reference distribution) are approximately the mean order statistics wewould find in the replications. (The mean in the replications would differ because the school standard
errors vary and because the value 13 does not incorporate the variability in the common θ.)
Just by comparing the observed order statistics to the means provided, it is clear that theidentical-schools model does fit this aspect of the data. Based on 1000 draws from the posterior
predictive distribution we find that observed order statistics correspond to upper-tail probabilities
ranging from 0.20 to 0.62, which further supports the conclusion.
6.1b. It is easy to imagine that we might be interested in choosing some schools for further study
or research on coaching programs. Under the identical-effects model, all schools would be equally
good for such follow-up studies. However, most observers would wish to continue experimenting in
School A rather than School C.
6.1c. The decision to reject the identical-effects model was based on the notion that the posterior
distribution it implies (all θj ’s the same) does not make substantive sense.
6.5.
1. A “hypothesis of no difference” probably refers to a hypothesis that θ = 0 or θ1 = θ2 in somemodel. In this case, we would just summarize our uncertainty about the parameters using theposterior distribution, rather than “testing” the hypothesis θ = 0 or θ1 = θ2 (see, for example,
Exercises 3.2–3.4). The issue is not whether θ is exactly zero, but rather what we know aboutθ.
2. One can also think of a model itself as a hypothesis, in which case the model checks describedin this chapter are hypothesis tests. In this case, the issue is not whether a model is exactly
true, but rather what are the areas where the model (a) does not fit the data, and (b) is
sensitive to questionable assumptions. The outcome of model checking is not “rejection” or“acceptance,” but rather an exploration of the weak points of the model. With more data, wewill be able to find more weak points. It is up to the modeler to decide whether an area inwhich the model does not fit the data is important or not (see page 366 for an example).

24

6.6a. Under the new protocol,p(y, n|θ) =

n
Y

i=1

θyi (1 − θ)1−yi 1Pn−1 (1−y )=12 1Pn
i=1

i

i=1

(1−yi )=13

For the given data, this is just θ7 (1 − θ)13 , which is the likelihood for these data under the “stopafter 20 trials” rule. So if the prior distribution for θ is unchanged, then the posterior distribution
is unchanged.

6.6b.

0

5

10
T (y-rep)

15

20

test <- NULL
for (i in 1:1000){
theta <- rbeta (1,8,14)
y.rep <- rbinom (1,1,theta)
while (sum(y.rep==0) < 13)y.rep <- c(y.rep, rbinom(1,1,theta))
n.rep <- length(y.rep)
test <- c(test,sum (y.rep[2:n.rep] != y.rep[1:(n.rep-1)]))}
hist (test, xlab="T (y-rep)", yaxt="n",breaks=seq(-.5, max(test)+.5), cex=2)

The posterior predictive distribution differs from that of Figure 6.4 in two notable ways. First,
there is more spread, which makes sense since n is now a random variable. Second, T (y rep) ismuch more likely to be even, especially for low values of T , which is a peculiarity of this particularsequential data collection rule.
6.7a. For very large A, the posterior distribution for θ is approximately N(5.1, 1/100). The posterior predictive distribution is y rep |y ∼ N(5.1, (1/100)11t + I), where 11t is the 100 × 100 squarematrix all of whose elements are 1. We compute the distribution of T (y rep) by simulation. The
posterior predictive p-value is estimated as the proportion of simulations of T (y rep) that exceed
T (y) = 8.1; this is 14%.

6.5

7.0

7.5

8.0
T (y-rep)

8.5

9.0

9.5

test <- NULL
for (i in 1:1000){
theta <- rnorm (1,5.1,1/sqrt(100))y.rep <- rnorm (100,theta,1)test <- c(test,max(abs(y.rep)))}postscript ("fig6.6a.ps", horizontal=TRUE)
par (mar=c(5,5,4,1)+.1)
hist (test, xlab="T (y-rep)", yaxt="n",nclass=20, cex=2)
lines (rep(8.1,2), c(0,1000))
print (mean(test>8.1))

6.7b. For very large A, the prior predictive distribution for T (y rep) = maxi yirep will extend from0 out to the neighborhood of A. Thus, if A is large, the prior predictive p-value for T (y) = 8.1 will
be approximately 1.0.
6.7c. Since the prior distribution of θ is essentially noninformative (for these data), the data almost
completely determine the posterior distribution. If the data fit the model reasonably well, then it
makes sense that the posterior predictive distribution is consistent with the data. On the other
hand, we would not expect a diffuse prior distribution to be consistent with any observed data set.
In particular, any finite dataset will be unusually close to 0 compared to an improper distributionthat extends out to ∞.
25

7.1a. “Randomization implies exchangeability” is not necessarily true because if you have a randomized design but covariates x are also available, and you think the covariates might be related toy, it makes sense to model (x, y) together. For example, if the cow experiment in Section 7.5 hadreally been conducted using a single complete randomization, it would still be a good idea to model
milk production conditional on age, lactation, and initial weight. Setting up a conditional model
will yield more precise inferences for the distribution of y.
“Randomization implies exchangeability” has a kernel of truth because, under randomization,
it is acceptable to use an exchangeable model: no additional variables are needed for an adequate
summary of the data.
7.1b. “Randomization is required for exchangeability” is not necessarily true because if you have noother information on a set of units but the measurements y, then you must model them exchangeably.
For example, no randomization was used in many of the examples used in this book (the footballgames, the speed of light measurements, the selection of the 8 schools, etc.)—but with no otherinformation available, exchangeable models are the only possibility.
“Randomization is required for exchangeability” has a kernel of truth because if data are collected
nonrandomly based on additional information, and there is information available to distinguish theunits that affected the data collection procedure (as in the cow experiment), then it is not appropriate
to use an exchangeable model on the data y alone.
7.1c. “Randomization implies ignorability” is not necessarily true for two reasons:
1. Even under randomization, data collection can be affected by other factors. For example, what
if treatments are assigned randomly, but then, in addition, all values of y > 200 are censored.
Then the censoring should be in the model; the data collection cannot be ignored.
2. If θ (the parameters for the data model) and φ (the parameters for the data collection model)
are not distinct, then the design is not ignorable: even if the data are collected randomly, theinclusion pattern provides information about φ, and thus θ (see variation 2 in Section 7.8).
“Randomization implies ignorability” has a kernel of truth because if the above two issues do
not arise, then a randomized design is ignorable—because under most methods of randomization,p(I|φ, y) does not depend on ymis , so the data are missing at random. Missing at random anddistinct parameters imply ignorability (see page 204).
7.1d. “Randomization is required for ignorability” is not necessarily true because many deterministic data collection designs have the property thet inclusion does not depend on ymis (and thus, sincea deterministic design has no parameters, the design is ignorable). For example: designs in which all
the data are used (for example, the football games, the speed of light experiment), designs chosen inadvance (for example, laying out two agricultural treatments in an AABBAABBAABB pattern),
designs that have no clear rule but are based entirely on observed covariates (for example, the cow
experiment), and designs based on observed outcomes (for example, sequential designs, which are
ignorable conditional on the order of the units and the observed outcomes).
“Randomization is required for ignorability” has a kernel of truth because if a nonrandomized
design is used, with data collection depending on the (possibly unobserved) value of y or on an
unobserved variable, then the pattern of data collected is, in general, informative, and the datacollection is nonignorable. For example, consider a medical experiment in which treatment A is more
likely to be given to patients that “look sick,” but this measure of “looking sick” is not quantified
so that it is not available for the analysis. This is a censoring-like problem and is nonignorable.

267.1e. “Ignorability implies exchangeability” is not necessarily true for the same reason as theanswer to Exercise 7.1a. Even if you are allowed to use an exchangeable model on y, you might still
do better by modeling (x, y), even if x is not used in the data collection.
“Ignorability implies exchangeability” has a kernel of truth because, if the data collection is
ignorable given yobs , so that it does not depend on other variables x, then it is not necessary toinclude x in the analysis: one can apply an exchangeable model to y.
7.1f. “Ignorability is required for exchangeability” is not necessarily true: consider censored datacollection (for example, variation 3 in Section 7.8). The design is nonignorable, and the inclusion
pattern must be included in the model, but the model is still exchangeable in the yobs i ’s.
“Ignorability is required for exchangeability” has a kernel of truth because if the data collection
is nonignorable given yobs , and additional information x is available that is informative with respect
to the data collection, then x should be included in the model because p(θ|x, yobs , I) 6= p(θ|yobs ).
This is an exchangeable model on (x, y), not y.
7.2. Many possibilities: for example, in the analysis of the rat tumor experiments in Chapter 5, weare assuming that the sample sizes, nj , provide no information about the parameters θj —this is likethe assumption of distinct parameters required for ignorability. But one could imagine situations inwhich this is not true—for example, if more data are collected in labs where higher tumor rates are
expected. This would be relevant for the parameter of interest, y71 , because n71 has the relatively
low value of 14. One way the analysis could reflect this possibility is to model θj conditional on nj .
To do this would be awkward using the beta distribution on the θj ’s; it might be easier to set up alogistic regression model—basically, a hierarchical version of the model in Section 3.7. This would
take a bit of computational and modeling effort—before doing so, it would make sense to plot theestimated θj ’s vs. the nj ’s to see if there seems to be any relation between the variables.
7.7a. We can write the posterior distribution of µ, σ as,
σ 2 |y obsµ|σ 2 , y obs∼∼

Inv-χ2 (n − 1, s2 )N(ȳ obs , σ 2 /n),where ȳ obs and s2 are the sample mean and variance of the n data points yiobs .
The predictive distribution of the mean of the missing data is,ȳ mis |µ, σ 2 , y obs ∼ N(µ, σ 2 /(N − n)).
Averaging over µ yields,ȳ mis |σ 2 , y obs ∼ N(ȳ obs , (1/n + 1/(N − n))σ 2 ).
Since σ 2 has an inverse-χ2 posterior distribution, averaging over it yields a t distribution (as on p.
76):ȳ mis |y obs ∼ tn−1 (ȳ obs , (1/n + 1/(N − n))s2 ).

The complete-data mean is ȳ = n/N ȳ obs + (N − n)/N ȳ mis , which is a linear transformation ofȳ(trating n, N , and ȳ obs as constants, since they are all observed) and so it has a t posterior
distribution also:
misȳ|y obs∼==

tn−1 (ȳ obs , s2 (1/n + 1/(N − n))(N − n)2 /N 2 )tn−1 (ȳ obs , s2 (N − n)/(nN ))tn−1 (ȳ obs , s2 (1/n − 1/N )).
27

7.7b. The t distribution with given location and scale parameters tends to the corresponding
normal distribution as the degrees of freedom go to infinity:
Since (1 + z/w)w → exp(z) as w → ∞, the expression (1/σ)(1 + (1/ν)((θ − µ)/σ)2 )−(ν+1)/2
approaches (1/σ) exp(− 12 ((θ − µ)/σ)2 ) as ν → ∞.
7.15a. The assignment for the current unit depends on previous treatment assignments. Thus, thetime sequence of the assignments must be recorded. This can be expressed as a covariate as ti foreach observation i, where ti = 1, 2, . . .. A time covariate giving the actual times (not just the timeorder) of the observations would be acceptable also (that is, the design would be ignorable given
this covariate) and, in practice, might make it easier to model the outcomes, given t.
7.15b. If the outcome is y = (y1 , . . . , yn ), we would model it as (y|x, θ), along with a prior distribution, p(θ|x), where θ describes all the parameters in the model and x are the covariates, which must
include the treatment assignment and the time variable (and can also include any other covariates
of interest, such as sex, age, medical history, etc.).
7.15c. Under this design, inferences are sensitive to dependence of y on t (conditional on θ andthe other variables in x). If such dependence is large, it needs to be modeled, or the inferences
will not be appropriate. Under this design, t plays a role similar to a blocking variable or the
pre-treatment covariates in the milk production example in Section 7.5. For example, suppose thatE(y|x, θ) has a linear trend in t but that this dependence is not modeled (that is, suppose that amodel is fit ignoring t). Then posterior means of treatment effects will tend to be reasonable, butposterior standard deviations will be too large, because this design yields treatment assignments
that, compared to complete randomization, tend to be more balanced for t.
7.15d. Under the proposed design and alternative (ii), you must include the times t as one of thecovariates in your model; under alternative (i), it is not necessary to include t. This is an advantage
of alternative (i), since it reduces the data collection requirements. On the other hand, it could
be considered an advantage of the proposed design and alternative (ii), since it forces you to use
this information, which may very well lead to more relevant inferences (if, for example, treatment
efficacy increases with t).
Suppose that, under all three designs, you will model the data using the same set of covariates
x (that include treatment assignment and t). Of the three designs, alternative (ii) is most balanced
for t, alternative (i) is least balanced (on average), and the proposed design falls in between. So,
assuming smooth changes in the dependence of y on t (conditional on x and θ), we would expect
data collected under alternative (ii) to yield the most precise inferences and under alternative (i)the least precise. If for some reason there was an alternating trend (for example, if the patients are
assigned to two doctors, and the doctors tend to alternate), then alternative (ii) could become theleast precise.
Another potential difficulty in alternative (ii) is posterior predictive model checking (or, in classical terms, hypothesis testing). Since this alternative has only two possible treatment vectors, it isbasically impossible to check the fit of the model under replications of the treatment under the same
patients. (It would be, however, possible to test the fit of the model as it would apply to additional
patients.)
8.1. We assume for simplicity that the posterior distribution is continuous and that needed moments exist. The proofs can be easily adapted to discontinuous distributions.

28

8.1a. The calculus below shows that L(a|y) has zero derivative at a = E(θ|y). Also the second
derivative is positive so this is a minimizing value.
ZZdd
E(L(a|y)) =(θ − a)2 p(θ|y)dθ = −2 (θ − a)p(θ|y)dθ = −2(E(θ|y) − a) = 0 if a = E(θ|y).
dada
8.1b. We can apply the argument from 4.5c with ko = k1 = 1.
Ra
8.1c. The calculus below shows that L(a|y) has zero derivative at any a for which −∞ p(θ|y)dθ =k0 /(k0 + k1 ). Once again the positive second derivative indicates that it is a minimizing value.
Z a

Z ∞
d
d
E(L(a|y)) =k1 (a − θ)p(θ|y)dθ +k0 (θ − a)p(θ|y)dθ
da
da
−∞
a
Z a
Z ∞
= k1
p(θ|y)dθ − k0p(θ|y)dθ−∞a
Z a= (k1 + k0 )p(θ|y)dθ − k0 .
−∞8.2. Denote the posterior mean by m(y) = E(θ|y) and consider m(y) as an estimator of θ.
Unbiasedness implies that E(m(y)|θ) = θ. In this case, the marginal expectation of θm(y) is
E(θm(y)) = E[E(θm(y)|θ)] = E[θ2 ]. But we can also write E(θm(y)) = E[E(θm(y)|y)] = E[m(y)2 ].
It follows that E[(m(y) − θ)2 ] = 0. This can only hold in degenerate problems for which m(y) = θ
with probability 1.
8.4. Our goal is to show that, for sufficiently large sigma, the “Bayes estimate” (the posterior meanof θ based on the prior density p(θ) = 1 in [0, 1]) has lower mean squared error than the maximum
likelihood estimate, for any value of θ ∈ [0, 1].
The maximum likelihood estimate, restricted to the interval [0, 1], takes value 0 with probabilityΦ(−c/σ) and takes value 1 with probability 1 − Φ((1 − c)/σ); these are just the probabilities (given σand θ) that y is less than 0 or greater than 1, respectively. For very large σ, these probabilities both
approach Φ(0) = 1/2. Thus, for very large σ, the mean squared error of the maximum likelihood
estimate is approximately 1/2[(1 − θ)2 + θ2 ] = 1/2 − θ + θ2 . (We use the notation Φ for the unit
normal cumulative distribution function.)On the other hand, p(θ|y) is proportional to the density function N(θ|y, σ 2 ) for θ ∈ [0, 1]. The
R1
R1
posterior mean of θ, E(θ|y), is 0 θN (θ|y, σ 2 )dθ/ 0 p(θ|y, σ)dθ. For very large σ, N(θ|y, σ 2 ) is
approximately constant over small ranges of θ (e.g., over [0, 1]). So the Bayes estimate is close toR1
0 θdθ = 1/2. (This works because, since the true value of θ is assumed to lie in [0, 1], the observed
y will almost certainly lie within a few standard deviations of 1/2.) Hence, for large σ, the meansquared error of the posterior mean is about (θ − 1/2)2 = 1/4 − θ + θ2 .
The difference in (asymptotic, as σ → ∞) mean squared errors is independent of the true value of
θ. Also, for large σ the maximum likelihood estimate generally chooses 0 or 1, each with probability
almost 1/2, whereas the Bayes estimate chooses 1/2 with probability almost 1.
8.6. This is a tricky problem.
For this to work out, given the form of θ̂, the prior distribution should be a mixture of a spikeat θ = 0 and a flat prior distribution for θ 6= 0. It’s easy to get confused with degenerate andnoninformative prior distributions, so let’s write it asp(θ) = λN(θ|0, τ12 ) + (1 − λ)N(θ|0, τ22 ),29

work through the algebra, and then take the limit τ1 → 0 (so that θ̂ = 0 is a possibility) and τ2 → ∞(so that θ̂ = y is a possibility).
We’ll have to figure out the appropriate limit for λ by examining the posterior distribution:p(θ|y)∝ p(θ)N(y|θ, σ 2 /n)∝ λN(θ|0, τ12 )N(y|θ, σ 2 /n) + (1 − λ)N(θ|0, τ22 )N(y|θ, σ 2 /n)
!
n 2
1
2
2
σ2 y
∝ λN(y|0, τ1 + σ /n) N θ 1+,+ σn2 τ12 + σn2
τ12
1
!
n 2
y
1
2
(1 − λ)N(y|0, τ22 + σ 2 /n) N θ 1 σ n , 1+ σ2 τ 2 + σn2τ222This is a mixture of two normal densities. In the limit τ1 → 0 and τ2 → ∞, this isp(θ|y) = λN(y|0, σ 2 /n)N(θ|0, τ12 ) + (1 − λ)N(y|0, τ22 )N(θ|y, σ 2 /n).
The estimate θ̂ cannot be the posterior mean (that would not make sense, since, as defined, θ̂is a discontinuous function of y). Since the two normal densities have much different variances, itwould also not make sense to use a posterior mode estimate. A more reasonable estimate (if a point
estimate must be used) is the mode of the hump of the posterior distribution that has greater mass.
That is,Set θ̂ = 0 if:

λN(y|0, σ 2 /n) >


1 n
1
exp − 2 y 2
λp
>
2σ
2π/nσ

(1 − λ)N(y|0, τ22 )1,(1 − λ) √
2πτ2

(9)

and set θ̂ = y otherwise.
√
For condition (9) to be equivalent to “if y < 1.96σ/ n,” as specified, we must have1
0.146λ p
2π/nσ
λ
1−λ

= (1 − λ) √
=

1
2πτ2

σ
√
.
0.146 nτ2

Since we are considering the limit τ2 → ∞, this means that λ → 0. That would be acceptable (akind of improper prior distribution), but the more serious problem here is that the limiting value forλτ2 depends on n, and thus the prior distribution for θ depends on n. A prior distribution cannot
depend on the data, so there is no prior distribution for which the given estimate θ̂ is a reasonable
posterior summary.
8.8. There are many possible answers here; for example:(a) With a weak prior distribution, why bother with Bayesian inference? Because it is a direct
way to get inferences for quantities of interest. For instance, in the examples of Sections 3.5
and 3.7, the Bayesian inference is straightforward. Why not use Bayesian inference in theseexamples? Is there a method that gives better results?
(b) With a strong prior distribution, why collect new data? This question misses the point, because
any strong prior distribution must come from lots of previous data (as in the example of thesex ratio of births in Chapter 2) or a theoretical model, or some combination of the two. If
enough previous data have been collected to form a precise inference, then new data may be30collected in order to see whether the model generalizes to other situations. For example, inthe placenta previa example of Section 2.5, we cannot really use the data from all previous
births in Germany to construct an extremely strong prior distribution—we have to be open to
the possibility that the new data will have a different rate of female births.
(c) With something in between, does Bayesianism “duck the issue”? The “issue” here is: whatshould the relative weights be for the prior and data? Bayesian inference sets these weights
for different cases based on the availability of data, as is explained for the cancer rate example
in Section 2.8.
The issue is further clarified in Chapter 5. In each of the three examples in Chapter 5, there
are multiple experiments, each run under slightly different conditions and with a relatively
small sample size. Using an informative—but not extremely strong—prior distribution allows
us to get better inferences for the result of each individual experiment. The Bayesian method
does not “duck the issue” because the relative weights given to the prior distribution and thedata are determined by the data themselves.
11.1a. Suppose that θ is drawn from the density proportional to g(θ), and U is a random Uniform(0, 1)
draw. Then we can express the cumulative distribution function of draws accepted by rejection sampling as∗

Pr(θ ≤ θ |θ is accepted)

=

=

=

p(θ|y)
Pr(θ ≤ θ∗ and U ≤ Mg(θ))

p(θ|y))
Pr(U ≤ Mg(θ)
R θ∗ R p(θ|y)/(Mg(θ))g(θ)dudθ−∞ 0
R ∞ R p(θ|y)/(Mg(θ))g(θ)dudθ
−∞ 0
∗
R
θ
1
M −∞ p(θ|y)dθ,1
M

which is the cumulative density function for p(θ|y).
A similar argument works in higher dimensions.
11.1b. The above proof requires that (i) g(θ) > 0 wherever p(θ) > 0 (otherwise the final integral
will not cover the support of p) and that (ii) p(θ)/M g(θ) ≤ 1 always (otherwise the integral over u
will have the range of [0, 1] rather than [0, p/M g], and the transition to the last line of the equation
will not be valid). If p/g is unbounded, than one or both of (i) and (ii) will be violated.
11.2. As described at the top of p. 291: consider any two points θa and θb , at iteration t labeled
so that p(θb |y)Jt (θa |θb ) ≥ p(θa |y)Jt (θb |θa ). To show that the posterior distribution is a stationary
distribution, suppose that θt−1 is a draw from the posterior distribution. Then the unconditional
probability density of a transition from θa to θb isp(θt−1 = θa , θt = θb ) = p(θa |y)Jt (θb |θa ),where the acceptance probability is 1 because of our labeling of a and b so that the ratio of importanceratios is at least 1. The unconditional probability density of a transition from θb to θa isp(θa |y)/Jt (θa |θb )tt−1,p(θ = θa , θ= θb ) = p(θb |y)Jt (θa |θb )p(θb |y)/Jt (θb |θa )which is the same as the unconditional probability in the other direction. Since their joint distribution is symmetric, θt and θt−1 have the same marginal distributions, and so p(θ|y) is a stationary
distribution of the Markov chain.
31

As with the Metropolis algorithm, the stationary distribution is unique if the Markov chain is
irreducible, aperiodic, and not transient.
12.6. We find the conditional posterior distribution of γ = log σ by transforming the inverse-χ2.
22γp(γ|θ, µ, τ, y) ∝ 2e2γ (e−2γ )n/2+1 e−nσ̂ /2e .
Maximizing the log of this density (which is equivalent to maximizing the density) yields γ̂ = log σ̂.
The argument for log τ is similar.
12.7. The joint posterior density is given by the formula on page 300. If θj = µ for all j, then theexpression simplifies top(θ, µ, log σ, log τ |y) ∝
=

nj
J Y
Y
√
N(yij |θj , σ 2 )τ ( 2πτ )−J
j=1 i=1−J/2 −(J−1)

(2π)

τ

nj
J Y
Y

j=1 i=1

N(yij |θj , σ 2 ).

This is just τ −(J−1) multiplied by a function of the other parameters and the data. In this case,J = 4, so as τ → 0 with θj = µ for all j, the joint posterior density is proportional to τ −3 , whichapproaches infinity.
14.1a. We fit a linear regression with a uniform prior density on (β, log σ). We include coefficientsfor the basement indicator and for three county indicators, and no constant term. (If data on morecounties were available, or if we were interested in forecasting for additional counties, then we would
fit a hierarchical regression at the county level, as discussed in Chapter 15.)
Quantity of interest
geometric mean for Blue Earth (no basement), exp(β2 )geometric mean for Blue Earth County (basement), exp(β1 + β2 )geometric mean for Clay County (no basement), exp(β3 )geometric mean for Clay County (basement), exp(β1 + β3 )geometric mean for Goodhue County (no basement), exp(β4 )geometric mean for Goodhue County (basement), exp(β1 + β4 )factor for basement vs. no basement, exp(β1 )geometric sd of predictions, exp(σ)

Posterior quantiles
25% 50% 75%
4.1
5.0
6.5
6.1
7.1
8.2
3.8
4.7
5.8
5.6
6.5
7.6
3.9
4.9
6.2
5.8
6.8
7.9
1.1
1.4
1.7
2.1
2.2
2.4

All estimates in the table come from the S code below: applying the normal regression model tothe logarithms of the data, drawing 1000 simulations of (β, σ), and then computing quantiles of thequantities of interest.
Nontechnical summary: The first six lines of the above table gives stimates and uncertainties
for the average county radon level (separating into houses without and with basements). These
range from about 4 to 8 pCi/L. Within a county, the houses with basements have, on average, radonlevels between 10% and 70% higher than houses without basements. One can expect the radonmeasurement in any given house to vary by about a factor of 2.2 from the mean of houses of thesame “basement status” in the same county.
make.indicators <- function (x){ux <- unique(x)32mat1 <- matrix (x, nrow=length(x), ncol=length(ux))
mat2 <- matrix (ux, nrow=length(x), ncol=length(ux), byrow=T)(mat1==mat2)*1}
y.1 <- c (5.0, 13.0, 7.2, 6.8, 12.8, 5.8, 9.5, 6.0, 3.8, 14.3, 1.8, 6.9, 4.7, 9.5)
y.2 <- c (0.9, 12.9, 2.6, 3.5, 26.6, 1.5, 13.0, 8.8, 19.5, 2.5, 9.0, 13.1, 3.6, 6.9)
y.3 <- c (14.3, 6.9, 7.6, 9.8, 2.6, 43.5, 4.9, 3.5, 4.8, 5.6, 3.5, 3.9, 6.7)
basement.1 <- c(1,1,1,1,1,0,1,1,1,0,1,1,1,1)
basement.2 <- c(0,1,1,0,1,1,1,1,1,0,1,1,1,0)
basement.3 <- c(1,0,1,0,1,1,1,1,1,1,1,1,1)
counties <- rep(1:3,c(length(y.1),length(y.2),length(y.3)))y <- c(y.1,y.2,y.3)
x <- cbind (c(basement.1,basement.2,basement.3), make.indicators(counties))
ls.out <- lsfit (x, log(y), intercept=F)
lsd <- ls.diag (ls.out)
nsim <- 10000n <- nrow(x)k <- ncol(x)sigma <- rep (NA, nsim)
beta <- array (NA, c(nsim, k))for (i in 1:nsim){
sigma[i] <- lsd$std.dev*sqrt((n-k)/rchisq(1,n-k))
beta[i,] <- ls.out$coef + (sigma[i]/lsd$std.dev)*lsd$std.err*t(chol(lsd$corr))%*%rnorm(k)}output <- exp (cbind (beta[,2], beta[,1]+beta[,2], beta[,3],beta[,1] + beta[,3], beta[,4], beta[,1] + beta[,4], beta[,1], sigma))for (i in 1:ncol(output)) print (round(quantile(output[,i],c(.25,.5,.75)),1))

14.1b. Let ỹ be the radon measurement in the new house, b̃ be the variable that equals 1 if thehouse has a basement, and c̃ be the county where the house is located. We are given c̃. Then we
can draw from the posterior predictive distribution of ỹ in four steps:
1. Draw (β, σ 2 ) from the posterior distribution as in Exercise 8.1a.
2. Draw b̃ from the posterior distribution p(b|c̃).
3. Draw log ỹ from the conditional posterior distribution p(log y|b̃, c̃) from the regression modelin Exercise 8.1a (N(β3 , σ 2 ) if b = 0 or N(β1 + β3 , σ 2 ) if b = 1).
4. Compute ỹ = exp(log ỹ).
To perform step 2 above, we need a model for b|c. Let θc be the proportion of houses in Blue
Earth County that have basements. For simplicity, we shall assume a binomial model (equivalent toassuming simple random sampling of houses within each county, and a large population of houses
relative to the sample size), so that bi |ci ∼ Bernoulli, (θc ), for each house i in this county. We also,
for simplicity, assume a uniform prior distribution on θc . The posterior distribution for θc is thenBeta(3, 13). So we draw b̃ by first drawing θc ∼ Beta(3, 13), then drawing b̃ ∼ Bernoulli(θc ).

0

10
20
30
radon measurement (new house)40

theta <- rbeta (nsim, 3, 13)b <- rbinom (nsim, 1, theta)
logy.rep <- rnorm (nsim, beta[,3] + b*beta[,1], sigma)
y.rep <- exp(logy.rep)
print (round(quantile(y.rep,c(.025,.25,.5,.75,.975)),1))
hist (y.rep[y.rep<40], yaxt="n", breaks=0:40,
xlab="radon measurement (new house)", cex=2)

The 2.5%, 25%, 50%, 75%, and 97.5% posterior quantiles for ỹ are 0.8, 2.8, 5.1, 9.0, and 33.3. To
make the histogram readable, we omit the approximately 2% of simulated values that are greaterthan 40.
33

14.3. Assuming p(β|σ) = p(β) to be uniform, we have p(β|σ, y) proportional to p(y|β, σ), that is,
to exp(− 21 (y − Xβ)t (y − Xβ)/σ 2 ). This posterior density is an exponential of a quadratic form in
β and hence is normal.
To show: the mean is β̂ = (X t X)−1 X t y and the variance is (X t X)−1 σ 2 .
It suffices to demonstrate that the following two quantities are equal (up to additive terms notdepending on β, which correspond to multiplicative constants in the posterior density function):(y − Xβ)t (y − Xβ) and (β − β̂)t X t X(β − β̂).
Expanding the quadratics, (y − Xβ)t (y − Xβ) = −β t X t y − y t Xβ + β t X t Xβ + constant, and(β − β̂)t X t X(β − β̂) = β t X t Xβ − y t X(X t X)−1 (X t X)β − β t (X t X)(X t X)−1 X t y + constant =β t X t Xβ − y t Xβ − β t X t y + constant, which establishes (14.4) and (14.5).
14.4. p(σ 2 |y) = p(β, σ 2 |y)/p(β|σ 2 , y). Assuming a prior density p(β, σ 2 ) = p(β)p(σ 2 ) = σ −2 andfixing β at the value β̂ on the right side of this equation, we get something proportional toσ −(n+2) exp(− 12 (y − X β̂)t (y − X β̂)/σ 2 ).
(det((X t X)−1 σ 2 ))−1/2 exp(0)(10)

Note that det((X t X)−1 σ 2 ) = σ 2k det((X t X)−1 ), since X t X has dimension k, and this expressionis σ 2k times a constant. Thus (10) is proportional in σ to σ −(n−k+2) exp(− 21 σ −2 (y − X β̂)t (y −X β̂)/σ 2 ).
Comparing this result to the scaled inverse-χ2 density in Appendix A reveals that our posterior
density is, in fact, inverse-χ2(n − k, s2 ).
14.7. First write the joint distribution of everything, p(β, log σ, y, ỹ) = N(y|β, σ 2 I)N(ỹ|β, σ 2 I).
Now consider y and σ to be constants (that is, condition on them). The joint density is then aquadratic function in (β, ỹ). So (β, ỹ) jointly have a normal distribution, which means that themarginal distribution of ỹ is also normal (see Appendix A).
17.1a. Conditional on the parameters µ1 , τ1 , λ1 , µ2 , τ2 , λ2 (which are assumed known for this problem), the posterior distribution isp(θ|y) ∝8 X
2Y

j=1 m=1

2
λm N(θj |µm , τm)N(yj |θj , σj2 ).

We use the fact thatN(θ|µ, τ 2 )N(y|θ, σ 2 ) = N(y|µ, σ 2 + τ 2 )N θ

1
11
τ 2 µ + σ2 y
.
,
1
1
1
1
τ 2 + σ2τ 2 + σ2(This is based on the identity p(θ)p(y|θ) = p(y, θ) = p(y)p(θ|y) for the normal model.) We can thenwrite the posterior distribution for the mixture model as


1
1
2
8 X
2 µm + σ 2 y j
Y
τ
1
m
j
,
, 1
λ∗mj N θjp(θ|y) ∝1
11
τ 2 + σ2τ 2 + σ2
j=1 m=1mwhere

λ∗mj =jmj

2
λm N(yj |µm , σj2 + τm)
for m = 1, 2.
22
λ1 N(yj |µ1 , σj + τ1 ) + λ2 N(yj |µ2 , σj2 + τ22 )For the given parameter values and the data in Table 5.2, the posterior distributions are independentand are as follows:34θ1 |y ∼ 0.70 N(8.8, 8.32) + 0.30 N(24.9, 12.82)
θ2 |y ∼ 0.90 N(3.9, 7.12) + 0.10 N(8.9, 9.42)
θ3 |y ∼ 0.94 N(−0.8, 8.52) + 0.06 N(2.5, 13.72)
θ4 |y ∼ 0.91 N(3.1, 7.42) + 0.09 N(8.1, 10.12)
θ5 |y ∼ 0.97 N(−0.3, 6.82) + 0.03 N(1.3, 8.82)
θ6 |y ∼ 0.95 N(0.3, 7.52) + 0.05 N(3.1, 10.42)
θ7 |y ∼ 0.71 N(8.7, 7.22) + 0.29 N(17.6, 9.62)
θ8 |y ∼ 0.88 N(3.0, 8.72) + 0.12 N(13.1, 14.42)
17.1b. The effect of the two-component mixture model is to shrink toward µ1 when yj is low; toshrink toward µ2 when yj is high; and to look like a mixture for intermediate values of yj .

post density of theta_8 when y_8 = 0050
theta_8

100

post density of theta_8 when y_8 = 25

050
theta_8

100

lambda0 <- c(.9,.1)
mu0 <- c(0,15)
tau0 <- c(10,25)
y <- c(28.39, 7.94, -2.75, 6.82, -0.64, 0.63, 18.01, 12.16)
sigma <- c(14.9, 10.2, 16.3, 11.0, 9.4, 11.4, 10.4, 17.6)
J <- length(y)M <- length(lambda0)
lambda <- array (NA, c(J,M))
mu <- array (NA, c(J,M))
tau <- array (NA, c(J,M))
for (m in 1:M){
lambda[,m] <- lambda0[m]*dnorm(y,mu0[m],sqrt(sigma^2+tau0[m]))tau[,m] <- sqrt(1/(1/tau0[m]^2 + 1/sigma^2))mu[,m] <- (mu0[m]/tau0[m]^2 + y/sigma^2)*tau[,m]^2}for (j in 1:J)
lambda[j,] <- lambda[j,]/sum(lambda[j,])print (cbind(lambda[,1],mu[,1],tau[,1],lambda[,2],mu[,2],tau[,2]))post density of theta_8 when y_8 = 50

050
theta_8

100

post density of theta_8 when y_8 = 100

050
theta_8

100

possible <- c(0,25,50,100)for (i in 1:length(possible)){y[8] <- possible[i]for (m in 1:M){
lambda[,m] <- lambda0[m]*dnorm(y,mu0[m],sqrt(sigma^2+tau0[m]))tau[,m] <- sqrt(1/(1/tau0[m]^2 + 1/sigma^2))mu[,m] <- (mu0[m]/tau0[m]^2 + y/sigma^2)*tau[,m]^2}for (j in 1:J)
lambda[j,] <- lambda[j,]/sum(lambda[j,])theta <- seq(-40,120,.1)
dens <- lambda[8,1]*dnorm(theta,mu[8,1],tau[8,1]) +lambda[8,2]*dnorm(theta,mu[8,2],tau[8,2])
plot (theta, dens, ylim=c(0,1.1*max(dens)),
type="l", xlab="theta_8", ylab="", xaxs="i",yaxs="i", yaxt="n", bty="n", cex=2)
mtext (paste("Post density of theta_8 when","y_8 =", possible[i]), cex=2, 3)}35