Solutions 9
Bayesian methods
1. The odds on the culprit being of Chinese origin, prior to the genotype result, are
3:1. Genotype result gives a likelihood ratio
Pr(genotype|Chinese) 1
= 0.0031/0.048 ≈
Pr(genotype|Indian) 15
Bayes theorem gives posterior odds on the culprit being Chinese as approximately
3/15 = 1:5 (5 to 1 against). Approximately, posterior probability that culprit is
Chinese is 1/6, Indian 5/6. Perhaps police should transfer attention to the Indian
community.
2. Odds that mother is a carrier are 1:1 (prior to information on status of her sons).
Each son inherits single X chromosome from his mother. If the mother is not a
carrier, the son is unaffected with probability 1. If the mother is a carrier, each son
has probability 1/2 of inheriting the ‘bad’ gene. Given that both sons are
unaffected, the likelihood ratio for mother a carrier/not a carrier is (1/2)2 . Using
Bayes theorem to combine this with the prior odds, the posterior odds that the
woman is a carrier is 1/4 (4 to 1 against), corresponding to a probability of 1/5.
If subsequently the lady has a third son who is unaffected, we can use the posterior
distribution derived above as a prior for the new situation. The new data provides a
likelihood ratio of 1/2, combined with prior odds of 1/4 gives posterior odds of 1:8
(probability that mother is carrier = 1/9).
3. (a) Given 10 heads, 10 tails, posterior is a beta distribution with parameters
10 + 2 = 12, 10 + 2 = 12.
(b) Posterior distribution of θ/(1 − θ) is F (24, 24). From Cambridge Table 12 (c),
97.5% point is 2.269, 2.5% point is 1/2.269 = 0.441, so (0.441, 2.269) is a 95%
Bayesian confidence interval for θ/(1 − θ). Interval for θ has end-points at
0.441/1.441 = 0.3060, 2.269/3.269 = 0.6941. Cambridge Table 29 is based on a
different prior so enter the Tables at n = 22, r = 11 to get (0.3059, 0.6941). In
this case, the equi-tailed interval and interval of highest posterior density coincide,
because posterior distribution is symmetric.
(c) To plot prior and posterior distributions,
lims ¡- qbeta(c(0.025,0.975), 12, 12)
curve(dbeta(x,12,12), from = lims[1], to = lims[2])
curve(dbeta(x,2,2), add=TRUE, lty=”dashed”)
4. Here Y + a = n − Y + b = 12, so, approximately, posterior distribution of logit(θ) is
normal with zero √
mean and variance 1/12 + 1/12 = 1/6. A 95% CI for the logit is
therefore ±1.96/ 6 = ±0.8. To transform this into a CI for θ, apply the inverse
transformation exp(x)/(1 + exp(x)) to ±0.8 to obtain (0.31, 0.69).
5. The ML estimate can be obtained as in week 8:
minuslogL ¡- function(x) – -125*log(2+x) - 38*log(1-x) -34*log(x) ˝
library(stats4)
mle(minuslogL, start = list(x = 0.5))
This gives θ̂ = 0.627 with s.e. 0.051.
To sample from posterior, distinguish between two types of animal in the first
group. One occurs with probability θ/4 and the other with probability 1/2. Let n0
be the number of animals of the first type. We treat this as an extra unknown
parameter. The conditional distributions required for Gibbs sampling are those of
n0 , given θ, and of θ, given n0 .
(a) Choose starting value for θ.
(b) Generate n0 as binomial with index n1 , parameter θ/(2 + θ).
(c) Generate θ as beta with parameters n0 + n3 + 1, n2 + 1.
(d) Repeat steps (b) and (c) N times, where N is a large number.
n1 ¡- 125; n2 ¡- 38; n3 ¡- 34
N ¡- 1500
gibbs ¡- numeric(N)
theta ¡- 0.5 # step (a)
for (i in 1:N) –
n0 ¡- rbinom(1,n1,theta/(2+theta)) # step (b)
theta ¡- rbeta(1, n0+n3+1, n2+1) # step (c)
gibbs[i] ¡- theta ˝
Plot of the Gibbs sample as a time series showed satisfactory mixing from the start,
so I did not reject any initial values as ‘burn-in’. Sample mean was 0.622 and
standard deviation 0.049, in close agreement with ML. Since Gibbs sampling
involves simulation, your results will differ by small random amounts.
6. Given the data (bull is black, offspring is black), the aa genotype is excluded for
both animals, and the ‘forward’ probabilities Pr(offspring = y |bull = x) are simply
Offspring
AA Aa
AA 1/2 1/2
Bull Aa 1/3 2/3
For ‘backward’ probabilities Pr(bull = y |offspring = x), use same reasoning as for
pedigree in lecture notes: Pr(bull = y |offspring = x) is proportional to
Pr(bull = y |bull’s parents Aa × Aa)×Pr(offspring = x|offspring’s parents y × Aa)
Working this through for x = AA, Aa and y =AA, Aa and normalising probabilities,
we find that the same probabilities apply backwards as well as forwards:
Bull
AA Aa
AA 1/2 1/2
Offspring Aa 1/3 2/3
Let θ be the probability that the (parent) bull is a carrier. The table of probabilities
shows that the probability that the offspring is a carrier is
(1/2)(1 − θ) + (2/3)θ = (1/2) + (1/6)θ
Applying this probability to the second table we find that the probability that the
bull is a carrier is
(2/3)(1/2 + θ/6) + (1/2)(1/2 − θ/6) = (7/12) + (1/36)θ
At equilibrium we must have θ = (7/12) + (1/36)θ, and solving this equation gives
θ = 3/5, the result we obtained in week 2. To simulate,
[Link] ¡- matrix( c(1/2,1/2,1/3,2/3), 2, 2, byrow=TRUE)
N ¡- 1000
gibbs ¡- numeric(N)
bull ¡- 1 # starting value
for (i in 1:N) –
offspring ¡- [Link]( 2, 1, prob = [Link][bull,] )
bull ¡- [Link]( 2, 1, prob = [Link][offspring,] )
gibbs[i] ¡- bull ˝
table(gibbs)
There should be approximately 400 ones (AA) and 600 twos (Aa). Of course this
problem can be easily solved directly with Bayes theorem but it is instructive to see
Gibbs sampling at work in a simple problem.