Bayesian MAS3301 Statistics Problems 3 and Solutions
Semester 2 2008-9
Problems 3
1. In a small survey, a random sample of 50 people from a large population is selected. Each person is asked a question to which the answer is either Yes or No. Let the proportion in the population who would answer Yes be . Our prior distribution for is a beta(1.5, 1.5) distribution. In the survey, 37 people answer Yes. (a) Find the prior mean and prior standard deviation of . (b) Find the prior probability that < 0.6. (c) Find the likelihood. (d) Find the posterior distribution of . (e) Find the posterior mean and posterior standard deviation of . (f) Plot a graph showing the prior and posterior probability density functions of on the same axes. (g) Find the posterior probability that < 0.6. Notes: The probability density function of a beta(a, b) distribution is f (x) = kxa1 (1 x)b1 where k is a constant. If X beta(a, b) then the mean of X is E(X) = and the variance of X is var(X) = a a+b
ab . (a + b + 1)(a + b)2
If X beta(a, b) then you can use a command such as the following in R to nd Pr(X < c). pbeta(c,a,b) To plot the prior and posterior probability densities you may use R commands such as the following. theta<-seq(0.01,0.99,0.01) prior<-dbeta(theta,a,b) posterior<-dbeta(theta,c,d) plot(theta,posterior,xlab=expression(theta),ylab="Density",type="l") lines(theta,prior,lty=2)
1
2. The populations, ni , and the number of cases, xi , of a disease in a year in each of six districts are given in the table below. Population n 120342 235967 243745 197452 276935 157222 Cases x 2 5 3 5 3 1
We suppose that the number Xi in a district with population ni is a Poisson random variable with mean ni /100000. The number in each district is independent of the numbers in other districts, given the value of . Our prior distribution for is a gamma distribution with mean 3.0 and standard deviation 2.0. (a) Find the parameters of the prior distribution. (b) Find the prior probability that < 2.0. (c) Find the likelihood. (d) Find the posterior distribution of . (e) Find the posterior mean and posterior standard deviation of . (f) Plot a graph showing the prior and posterior probability density functions of on the same axes. (g) Find the posterior probability that < 2.0. Notes: The probability density function of a gamma(a, b) distribution is f (x) = kxa1 exp(bx) where k is a constant. If X gamma(a, b) then the mean of X is E(X) = a/b and the variance of X is var(X) = a/(b2 ). If X gamma(a, b) then you can use a command such as the following in R to nd Pr(X < c). pgamma(c,a,b) To plot the prior and posterior probability densities you may use R commands such as the following. lambda<-seq(0.00,5.00,0.01) prior<-dgamma(lambda,a,b) posterior<-dgamma(lambda,c,d) plot(lambda,posterior,xlab=expression(lambda),ylab="Density",type="l") lines(lambda,prior,lty=2) 3. Geologists note the type of rock at xed vertical intervals of six inches up a quarry face. At this quarry there are four types of rock. The following model is adopted. The conditional probability that the next rock type is j given that the present type is i and 4 given whatever has gone before is pij . Clearly j=1 pij = 1 for all i. The following table gives the observed (upwards) transition frequencies.
2
From rock
1 2 3 4
1 56 15 20 6
To rock 2 3 13 24 93 22 25 153 35 11
4 4 35 11 44
Our prior distribution for the transition probabilities is as follows. For each i we have a uniform distribution over the space of possible values of pi1 , . . . , pi4 . The prior distribution of pi1 , . . . , pi4 is independent of that for pk1 , . . . , pk4 for i = k. Find the matrix of posterior expectations of the transition probabilities. Note that the integral of xn1 xn2 xn3 xn4 over the region such that xj > 0 for j = 1, . . . , 4 and 1 2 3 4 4 xj = 1, where n1 , . . . , n4 are positive is j=1
1
x n1 1
0 0
1x1
xn2 2
0
1x1 x2
xn3 (1 x1 x2 x3 )n4 .dx3 .dx2 .dx1 3
=
(n1 + 1)(n2 + 1)(n3 + 1)(n4 + 1) (n1 + n2 + n3 + n4 + 4)
4. A biologist is interested in the proportion, , of badgers in a particular area which carry the infection responsible for bovine tuberculosis. The biologists prior distribution for is a beta(1, 19) distribution. (a) i. Find the biologists prior mean and prior standard deviation for . ii. Find the cumulative distribution function of the biologists prior distribution and hence nd values 1 , 2 such that, in the biologists prior distribution, Pr( < 1 ) = Pr( > 2 ) = 0.05.
(b) The biologist captures twenty badgers and tests them for the infection. Assume that, given , the number, X, of these carrying the infection has a binomial(20, ) distribution. The observed number carrying the infection is x = 2. i. ii. iii. iv. Find the likelihood function. Find the biologists posterior distribution for . Find the biologists posterior mean and posterior standard deviation for . Use R to plot a graph showing the biologists prior and posterior probability density functions for .
5. A factory produces large numbers of packets of nuts. As part of the quality control process, samples of the packets are taken and weighed to check whether they are underweight. Let the true proportion of packets which are underweight be and assume that, given , the packets are independent and each has probability of being underweight. A beta(1, 9) prior distribution for is used. (a) The procedure consists of selecting packets until either an underweight packet is found, in which case we stop and note the number X of packets examined, or m = 10 packets are examined and none is underweight, in which we case we stop and note this fact. i. Find the posterior distribution for when X = 7 is observed. ii. Find the posterior distribution for when no underweight packets are found out of m = 10. (b) Now consider varying the value of m. Use R to nd the posterior probability that < 0.02 when no underweight packets are found out of i. m = 10, ii. m = 20, 3
iii. m = 30. 6. The numbers of patients arriving at a minor injuries clinic in 10 half-hour intervals are recorded. It is supposed that, given the value of a parameter , the number Xj arriving in interval j has a Poisson distribution Xj Poisson() and Xj is independent of Xk for j = k. The prior distribution for is a gamma(a, b) distribution. The prior mean is 10 and the prior standard deviation is 5. (a) i. Find the values of a and b. ii. Let W 2 . Find values w1 , w2 such that Pr(W < w1 ) = Pr(W > w2 ) = 0.025. 2a Hence nd values l1 , l2 such that, in the prior distribution, Pr( < l1 ) = Pr( > l2 ) = 0.025. iii. Using R (or otherwise) nd a 95% prior highest probability density interval for . iv. Compare these two intervals.
(b) The data are as follows. 9 12 16 12 16 11 18 13 12 19
i. Find the posterior distribution of . ii. Using R (or otherwise) nd a 95% posterior highest probability density interval for . 7. The numbers of sales of a particular item from an Internet retail site in each of 20 weeks are recorded. Assume that, given the value of a parameter , these numbers are independent observations from the Poisson() distribution. Our prior distribution for is a gamma(a, b) distribution. (a) Our prior mean and standard deviation for are 16 and 8 respectively. Find the values of a and b. (b) The observed numbers of sales are as follows. 14 19 14 21 22 33 15 13 16 19 27 22 27 21 16 25 14 23 22 17 Find the posterior distribution of . (c) Using R or otherwise, plot a graph showing both the prior and posterior probability density functions of . (d) Using R or otherwise, nd a 95% posterior hpd interval for . (Note: The R function hpdgamma is available from the Module Web Page). 8. In a medical experiment, patients with a chronic condition are asked to say which of two treatments, A, B, they prefer. (You may assume for the purpose of this question that every patient will express a preference one way or the other). Let the population proportion who prefer A be . We observe a sample of n patients. Given , the n responses are independent and the probability that a particular patient prefers A is . Our prior distribution for is a beta(a, a) distribution with a standard deviation of 0.25. (a) Find the value of a. (b) We observe n = 30 patients of whom 21 prefer treatment A. Find the posterior distribution of . (c) Find the posterior mean and standard deviation of . (d) Using R or otherwise, plot a graph showing both the prior and posterior probability density functions of . (e) Using R or otherwise, nd a symmetric 95% posterior probability interval for . (Hint: The R command qbeta(0.025,a,b) will give the 2.5% point of a beta(a, b) distribution). 4
9. The survival times, in months, of patients diagnosed with a severe form of a terminal illness are thought to be well by modelled an exponential() distribution. We observe the survival times of n such patients. Our prior distribution for is a gamma(a, b) distribution. (a) Prior beliefs are expressed in terms of the median lifetime, m. Find an expression for m in terms of . (b) In the prior distribution, the lower 5% point for m is 6.0 and the upper 5% point is 46.2. Find the corresponding lower and upper 5% points for . Let these be k1 , k2 respectively. (c) Let k2 /k1 = r. Find, to the nearest integer, the value of such that, in a 2 distribution, the 95% point divided by the 5% point is r and hence deduce the value of a. (d) Using your value of a and one of the percentage points for , nd the value of b. (e) We observe n = 25 patients and the sum of the lifetimes is 502. Find the posterior distribution of . (f) Using the relationship of the gamma distribution to the 2 distribution, or otherwise, nd a symmetric 95% posterior interval for . Note: The R command qchisq(0.025,nu) will give the lower 2.5% point of a 2 distribution on nu degrees of freedom.
Homework 3
Solutions to Questions 7, 8, 9 of Problems 3 are to be submitted in the Homework Letterbox no later than 4.00pm on Monday March 9th.
5
Solutions
1. (a) In the prior a = 1.5 and b = 1.5. So the mean is 1.5 a = = 0.5. a+b 3.0 The variance is (a + so the standard deviation is ab b)2 (a + b + 1) = 1.5 1.5 1 = 24 3 16
1 = 0.25. 4
(b) Using R the prior probability that < 0.6 is 0.62647. > pbeta(0.6,1.5,1.5) [1] 0.62647 (c) The likelihood is 50 37 37 (1 )13 .
(d) The prior density is proportional to 1.51 (1 )1.51 The likelihood is proportional to 37 (1 )13 Hence the posterior density is proportional to 38.51 (1 )14.51 The posterior distribution is beta(38.5, 14.5). (e) In the posterior a = 38.5 and b = 14.5. So the mean is 38.5 a = = 0.7264. a+b 53.0 The variance is ab 38.5 14.5 = = 3.6803 103 (a + b)2 (a + b + 1) 532 54 so the standard deviation is 0.06067. (f) See Figure 1. > > > > > theta<-seq(0.01,0.99,0.01) prior<-dbeta(theta,1.5,1.5) posterior<-dbeta(theta,38.5,14.5) plot(theta,posterior,xlab=expression(theta),ylab="Density",type="l") lines(theta,prior,lty=2)
(g) Using R the posterior probability that < 0.6 is 0.02490528. > pbeta(0.6,38.5,14.5) [1] 0.02490528 2. (a) The mean is a/b = 3 and the variance is a/b2 = 4. So 9 a2 /b2 = = a, 4 a/b2 giving a = 2.25 and b= 2.25 = 0.75. 3
6
Density
0 0.0
1
2
3
4
5
6
0.2
0.4
0.6
0.8
1.0
Figure 1: Prior (dashes) and posterior (solid) pdfs for Question 1. (b) Using R the prior probability that < 2.0 is 0.3672305. > pgamma(2,2.25,0.75) [1] 0.3672305 (c) The likelihood is ei xi i xi ! i=1
n
=
e
i
xi i (ni /100000)xi xi !
xi !
= en/100000 S where n = ni = 1231663 and S = xi = 19. This is proportional to e12.31663 19 .
(d) The prior density is proportional to 2.251 e0.75 19 12.316663 The likelihood is proportional to e Hence the posterior density is proportional to 21.251 e13.06663 The posterior distribution is gamma(21.25, 13.06663). (e) In the posterior a = 21.25 and b = 13.06663. So the mean is a 21.25 = = 1.6262. b 13.06663 The standard deviation is a 21.25 = = 0.3528. b 13.06663 7
Density
0.0 0
0.2
0.4
0.6
0.8
1.0
2
4
6
8
Figure 2: Prior (dashes) and posterior (solid) pdfs for Question 2. (f) See Figure 2. > > > > > lambda<-seq(0.05,8.0,0.05) prior<-dgamma(lambda,2.25,0.75) posterior<-dgamma(lambda,21.25,13.06663) plot(lambda,posterior,xlab=expression(lambda),ylab="Density",type="l") lines(lambda,prior,lty=2)
(g) Using R the posterior probability that < 2.0 is 0.8551274. > pgamma(2,21.25,13.06663) [1] 0.8551274 3. Since the prior distribution is uniform the prior density is a constant. Therefore the posterior density is proportional to the likelihood. the likelihood is
4 4
L=
i=1 j=1
pijij
n
where nij is the observed number of transitions from rock i to rock j. The posterior density is therefore
4
fi (pi1 , pi2 , pi3 , pi4 )
i=1
(1)
where fi (pi1 , pi2 , pi3 , pi4 ) = k1i pni1 pni2 pni3 pni4 i1 i2 i3 i4 is the posterior density of pi1 , pi2 , pi3 , pi4 . 8
(1)
Since fi (pi1 , pi2 , pi3 , pi4 ) dpi1 dpi2 dpi3 = 1,
R (1)
we must have
1 k1i
=
R
pni1 pni2 pni3 pni4 dpi1 dpi2 dpi3 i1 i2 i3 i4 (ni1 + 1)(ni2 + 1)(ni3 + 1)(ni4 + 1) (Ni + 4)
=
where Ni = ni1 + ni2 + ni3 + ni4 and the integrals are taken over the region R in which (pi1 , pi2 , pi3 , pi4 ) must lie and pi4 = 1 pi1 pi2 pi3 . Now, the posterior mean of pi1 , for example, is
E(1) (pi1 )
=
R
pi1 fi (pi1 , pi2 , pi3 , pi4 ) dpi1 dpi2 dpi3 ki1 pni1 +1 pni2 pni3 pni4 dpi1 dpi2 dpi3 i2 i3 i4 i1
R
(1)
= = k1i k2i where k2i = So E(1) (pi1 ) = = = In general
(ni1 + 2)(ni2 + 1)(ni3 + 1)(ni4 + 1) . (Ni + 5)
(ni1 + 2)(ni2 + 1)(ni3 + 1)(ni4 + 1) (Ni + 4) (ni1 + 1)(ni2 + 1)(ni3 + 1)(ni4 + 1) (Ni + 5) (Ni + 4) (ni1 + 1)(ni1 + 1) (ni1 + 1) (Ni + 4)(Ni + 4) ni1 + 1 Ni + 4
E(1) (pij ) = The table of posterior means is as follows.
nij + 1 . Ni + 4
From rock
1 2 3 4
1 0.5644 0.0947 0.0986 0.0700
To rock 2 3 0.1386 0.2475 0.5562 0.1361 0.1221 0.7230 0.3600 0.1200
4 0.0495 0.2130 0.0563 0.4500
4. 5. 6.
9
7. (a) Prior mean:
a = 16, b a = 64. b2
Prior variance: Hence a = 4 and b = 0.25.
(1 mark) (b) From the data s = ...