Description

(a) As we know that X|θ ∼ Poi(θ) (i.e., ), and the likelihood can be found as
.
Since we have

we know that the posterior follows . Thus the Bayes estimator for θ can be computed as
.
Based on the four realizations 1,2,0, and 1, we know that the Bayes estimator for θ is
.
MLE for . This shows that θBayes > θMLE.
(b) We use the following code to find the 95% equitailed credible set
> gaminv(0.025, 4+0.5, 1/4) [1] 0.3375
> gaminv(0.975, 4+0.5, 1/4)
[1] 2.3778
The 95% equitailed credible set is computed as [0.3375, 2.3778].
(c) We find the 95% HPD as [0.23781, 2.17410]. The sketch of HPD and the distribution of posterior is shown in Figure 1. The MATLAB code for finding HPD is attached in appendix.
(d) We compute the mode of posterior, which is also MAP estimator of θ, as
.
(e) We use the following code to compute the probability for H1.

Figure 1: Sketch of HPD
> gamcdf(1, 4+0.5, 1/4)
[1] 0.4659
This means that the probability for H0 is 1 − 0.4659 = 0.5341. Hence, hypothesis H0 is more favored.
2 Lady Guessing Coin Flips.
(a) According to Bayes formula,

and p1 = 1 − p0 = 0.9621. The Bayes factor is

Since K < 1, the result is negative (supports H1)
(b) Since the Bayes factor K is less than 1, the result is negative (supports H1). The experiment is convincing that the lady possesses ESP.
A MATLAB for finding HPD in problem 1 (c)
1 %%
2 % search for k
3 format long 4 a = 4.5; b = 4;
5 for k=0.1:0.000001:0.8
6 ff=@(x) 1/gamma(a) * x.^(a-1) .* b.^a .* exp(- b * x) – k;
7 a1=fzero(ff, 0.5);
8 a2=fzero(ff, 2);
9 c=gamcdf(a2, a, 1/b) – gamcdf(a1, a, 1/b);
10 if (abs(c-0.95)<0.00001)
11 break;
12 end
13 end
14 k
15
16 %%
17 % find HPD
18 format long
19 k=0.111507;
20 a=4.5; b=4;
21 ff=@(x) 1/gamma(a) * x.^(a-1) .* b.^a .* exp(- b * x) – k
22 a1=fzero(ff, 0.5) % 0.23781
23 a2=fzero(ff, 2) % 2.17410
24 gamcdf(a2, a, 1/b) – gamcdf(a1, a, 1/b) % 0.95
25 format short
26 lengthhpd = a2 – a1 % 1.9363
27
28 xx=0.01:0.001:3.5 ;
29 plot(xx, f(xx, a, b),’k-‘,’linewidth’,2)
30 hold on
31 plot(xx, k*ones(size(xx)),’r-‘)
32 plot(0.23781, f(0.23781, a, b), ‘o’)
33 plot(0.23781, 0, ‘ro’)
34 plot(2.17410, f(2.17410, a, b), ‘o’)
35 plot(2.17410, 0, ‘ro’)
36 plot([0.23781 2.17410],[0 0], ‘r-‘,’linewidth’,8)
37 hold off