Reference no: EM131055819

TEXTBOOK - COMPUTATIONAL STATISTICS SECOND EDITION By GEOF H. GIVENS AND JENNIFER A. HOETING

Problem 1-

Consider the data set "Rat_data.csv" (Read the data using the function read.csv()) in R. Each data point represents the weight of a rat at a certain age. More specifically, weights of 30 rats (in grams) were measured at five time points of x₁ = 8, x₂ = 15, x₃ = 22, x₄ = 29, x₅ = 36(weeks). For example, weight of rat 3 in week 15 was 214 grams. Let y_ij denote the weight of rat i at age x_j, and consider the following model for the y_ij

Y_{ij ~} N(α_i + β_i(x_j - x^-), σ²) for i = 1, . . . , 30 and j = 1, . . . , 5 and x^- = 22  

(a) Obtain the maximum likelihood estimate of α_i, β_i and σ². Show your work, and compute the final estimates.

(b) Now consider a Bayesian analysis of this model, where we assume the following priors

α_{i ~} N(α, 100) with hyperprior     α _~ N(0,10000)

β_{i ~} N(β, 1)    with hyperprior       β_~ N(0,10000)

?? = (1/σ²) _~ gamma(10^-3, 10^-3)

where _~α_i, β_i and ?? are assumed to be independent. Our aim here is to use the Gibbs sampling method to generate values from the posterior distribution of (α₁, . . . , α₃₀, β₁, . . . , β₃₀, τ).

(i) State all the required conditionals. Show your derivation very briefly. 

(ii) Write a complete algorithm for this Gibbs sampling algorithm, specific to this problem and code your algorithm, using R.

(iii) Run the chain for 50,000 iterations, and burn-in the first 10,000 generated values. Plot the chain (trace) after the burn-in for α₁, β₁, α₁₅, β₁₅, τ against the iteration number, and comment on the mixing. 

(iv) Draw the histograms for the posterior distribution of α₁, β₁, α₁₅, β₁₅, τ and provide a point estimate for each parameter.

(c) Obtain a 95% credible interval for each of the parameters α₁, β₁. Interpret your intervals in the context of the problem. 

(d) In the prior, we assumed that α₁, β₁ and τ are independent. Does the posterior support this assumption? Justify your answer.

Problem 2-

Consider simulation of the bivariate normal distribution with μ = (1, 2)^T and covariance matrix

using Metropolis-Hastings algorithm with random walk generating density. Specifically, let x^∗ = x^(t)+ε, where ε _~ N₂(0, D) with D being a diagonal matrix with diagonal elements δ₁ and δ₂. Note that δ₁ and δ₂ control the spread of values generated along the first and second coordinate axis, respectively.

(a) Give the M-H algorithm, and code it in R. 

(b) Choose δ₁ = 0.001 and δ₂ = 0.001, and generate a chain, starting at x⁽⁰⁾ = (0, 0)^t, of size 10,000 using the M-H algorithm.

(i) Compute the acceptance rate, and comment on the acceptance rate. 

(ii) Burn 1000, and draw a trace plot and the density plot for the chain for each of the variables. Comment on mixing and the density plots. 

(iii) Draw 9000 pairs generated from multivariate normal with μ = (1, 2)^T and covariance matrix

Plot these pairs on a scatterplot overlaying it with the last 9000 points obtained from the M-H chain; use different color for those generated from MH and those directly generated from the target (choosing an alpha-level for color transparency is helpful). Do the M-H generated values match the target? Explain why or why not. 

(c) Choose appropriate δ₁ and δ₂ so that the M-H algorithm would generate values that meet the target requirement. Draw a similar scatter-plot as in part a(ii) for this case. 

(d) How do you expect the acceptance rate would change, if you choose δ₁ and δ₂ to be very large? Briefly explain why?

I will need the assignment along with the R codes. R codes are extremely important in this course, please have some explanation in R file. Please try to avoid using built in packages if possible.

Attachment:- rats_data.rar