Description

(a) By observing n samples (Xi,Yi),i = 1,…,n, the likelihood can be written as

By choosing prior on , we have the posterior as

(b) By choosing the given proposal distribution, the algorithm for generating samples from posterior is shown below. Note that in the following algorithm, θ denotes the target parameter that we hope to obtain from posterior, which is ρ. To avoid confusion, we use γ(θn,θ0) to denote the acceptance ratio (which is same as what Brani used ρ(θn,θ0) in the lecture video).
• Step 1. Start with arbitrary θ0.
• Step 2. At stage n + 1, generate proposal θ0 from uniform U(θn − 0.1,θn + 0.1).
• Step 3. Determine θn+1 as
– θn+1 = θ0 with probability γ(θn,θ0); – θn+1 = θn with probability 1 − γ(θn,θ0).
• Step 4. Set n = n + 1 and go to Step 2.
We now find the acceptance ratio γ(θn,θ0) as follows.
Based on the proposal distribution for θ, we know that

Similarly, we also have
.
Notice that θn − 0.1 ≤ θ0 ≤ θn + 0.1 and θ0 − 0.1 ≤ θn ≤ θ0 + 0.1 are indeed equivalent, and thus
.
Thus, the acceptance ratio can be computed as
.
By plugging in the values for and based on n = 100 observations, the
acceptance ratio can be computed as

(c) We attached the histogram of ρs and the realizations of the chain for the last 1000 simulations shown in Figure 1 and 2, respectively. The Bayes estimator can be computed from the mean of the ρs obtained. We have ρˆBayes = 0.7224.
(d) By replacing the proposal distribution, we plot the histogram of ρs and the realizations of the chain for the last 1000 simulations shown in Figure 3 and 4, respectively. The Bayes estimator of ρ based on proposal distribution U(−1,1) is computed as
ρˆBayes = 0.7227.

Figure 1: Histogram of ρs given proposal distribution U(ρi−1 − 0.1,ρi−1 + 0.1)
We notice that the Bayes estimator with proposal distribution U(−1,1) is similar to that with U(ρi−1− 0.1,ρi−1 + 0.1). The acceptance ratio with proposal distribution U(−1,1) is also lower (can be seen from the parallel segments from the behavior of the chain).
2 Gibbs Sampler and Lifetimes with Multiplicative Frailty.
(a)

Hence λ|µ,t1,..,tn respects to Gamma distribution with parameters n+c and
α. By symmetry, µ|λ,t1,..,tn respects to Gamma distribution with parameters n + d and .
(b) Gibbs sampler algorithm for sampling (λ,µ)|t1,…,tn is as follows:
• Step 1. Start with µ0 = 0.1.

Figure 2: Realizations of the chain for the last 1000 simulations given given proposal distribution U(ρi−1 − 0.1,ρi−1 + 0.1)
• Step 2. For m = 1,…,51000
– λm is a sample generated from
– µm is a sample generated from
• Step 3. Discard (λm,µm), m = 1,…,1000.
(c) Scatterplot of (λ,µ): Histograms of individual components λ and µ:
• Posterior means for λ = 0.0547 and for µ = 0.6911.
• Posterior variances for λ = 0.000358 and for µ = 0.06579.
• 95% equatailed credible sets for λ = (0.0255,0.0983) and for µ = (0.3185,1.3139).

Figure 3: Histogram of ρs given proposal distribution U(−1,1)

Figure 4: Realizations of the chain for the last 1000 simulations given given proposal distribution U(−1,1)

Figure 5: Scatter plot of (λ,µ)

Figure 6: Histograms of individual components λ and µ
A MATLAB code for problem 1
1
2 close all force
3 clear all
4
5 theta = 0; %initial values
6 thetas = [theta];
7
8 sum_x2 = 114.9707;
9 sum_y2 = 105.9196;
10 sum_xy = 82.5247;
11
12 n = 100;
13
14 tic
15 for i = 1:51000
16 %theta_prop = theta-0.1+2*0.1*rand();
17 theta_prop = -1+2*rand();
18 %————————————————————————–
19 if (theta≥-1 && theta≤1) 20 indicator_theta = 1;
21 else
22 indicator_theta = 0;
23 end
24 if (theta_prop≥-1 && theta_prop≤1)
25 indicator_theta_prop = 1;
26 else
27 indicator_theta_prop = 0;
28 end
29 %————————————————————————-30 rr = ((1-theta_prop^2)^(-(n+3)/2)* …
31 exp(-(sum_x2-2*theta_prop*sum_xy+sum_y2)/(2*(1-theta_prop^2))) * …
32 indicator_theta_prop)/((1-theta^2)^(-(n+3)/2)* …
33 exp(-(sum_x2-2*theta*sum_xy+sum_y2)/(2*(1-theta^2)))*indicator_theta);
34 %————————————————————————–
35 rho = min( rr ,1);
36 if (rand < rho)
37 theta = theta_prop;
38 end
39 thetas = [thetas theta];
40 end
41 toc
42 %Burn in 1000
43 thetas = thetas(1000:end);
44 figure(1)
45 hist(thetas, 50)
46 xlabel(‘theta’)
B R code for problem 2
1 alpha=100; beta=5
2 c=3; d=1
3 sumt=512
4 n=20
5 M=51000
6
7 mu=0.1
8 set.seed(1)
9 record = matrix(0, nrow=2, ncol=M)
10 for(m in 1:M){
11 lambda=rgamma(1, shape=(n+c), rate=(mu*sumt+alpha))
12 mu=rgamma(1, shape=(n+d), rate=(lambda*sumt+beta))
13 record[,m] = c(lambda, mu)
14 }
15 lambda=record[1,-(1:1000)]
16 mu=record[2,-(1:1000)]
17 ## scatterplot
18 par(mfrow=c(1, 1))
19 plot(lambda, mu, main=”scattor plot”, xlab =”lambda”, ylab=”mu”, cex=0.3)
20
21 ## histogram of lambda and mu
22 par(mfrow=c(1, 2))
23 hist(lambda, xlab=”lambda”, main=”histogram of lambda”, breaks=50) 24 hist(mu, xlab=”mu”, main=”histogram of mu”, breaks=50)
25
26 ## mean of lambda and mu
27 mean(lambda)
28 mean(mu)
29
30 var(lambda)
31 var(mu)
32
33 ## 95% equitailed credible set of lambda and mu
34 slambda=sort(lambda)
35 smu=sort(mu)
36 slambda[c(1250, 48750)]
37 smu[c(1250, 48750)]
38
39 ## mean of lambda*mu
40 mean(lambda*mu)