BSDA-Homework 4 Solved

Bioassay model
In this exercise, you will use a dose-response relation model that is used in Section 3.7 of the course book. The used likelihood is the same, but instead of uniform priors, we will use a bivariate normal distribution as the joint prior distribution of the parameters α and β.

a)     In the prior distribution for (α,β), the marginal distributions are α ∼ N(0,22) and β ∼ N(10,102), and the correlation between them is corr(α,β) = 0.6. Report the mean (vector of two values) and covariance (two by two matrix) of the bivariate normal distribution.

Hint! The mean and covariance of the bivariate normal distribution are a length 2 vector and a 2 × 2 matrix. The elements of the covariance matrix can be computed using the relation of correlation and covariance.

b)     You are given 4000 independent draws from the posterior distribution of the model. Load the draws with data("bioassay_posterior"). Report the mean as well as 5 % and 95 % quantiles separately for both α and β. Report also the Monte Carlo standard errors (MCSEs) for the mean and quantile estimates. Report as many digits for the mean and quantiles as the MCSEs allow. In other words, leave out digits where MCSE is nonzero (Example: if posterior mean is 2.345678 and MCSE is 0.0012345, report two digits after the decimal sign, taking into account the usual rounding rule, so you would report 2.35. Further digits do not contain useful information due to the Monte Carlo uncertainty.). Explain in words what does Monte Carlo standard error mean and how you decided the number of digits to show.

Note! The answer is graded as correct only if the number of digits reported is correct! The number of signi cant digits can be di erent for the mean and quantile estimates. In some other cases, the number of digits reported can be less than MCSE allows for practical reasons.

Hint! Quantiles can be computed with the quantile function. With S draws,

 the MCSE for E[θ] is pVar[θ]/S. MCSE for the quantile estimates can be computed with the mcse_quantile function from the bsda package.

Importance sampling
Now we discard our posterior draws and switch to importance sampling.

c) Implement a function for computing the log importance ratios (log importance weights) when the importance sampling target distribution is the posterior distribution, and the proposal distribution is the prior distribution from a). Below is a test example, the functions can also be tested with markmyassignment. Explain in words why it’s better to compute log ratios instead of ratios.

Note! The values below are only a test case. In this c) part, you only need to report the source code of your function, as it will be needed in later parts.

Hints! Use the function rmvnorm from the bsda package for sampling. Nonlog importance ratios are given by equation (10.3) in the course book. The fact that our proposal distribution is the same as the prior distribution makes this task easier. The logarithm of the likelihood can be computed with the bioassaylp function from the bsda package. The data required for the likelihood can be loaded with data("bioassay").

alpha <- c(1.896, -3.6, 0.374, 0.964, -3.123, -1.581) beta <- c(24.76, 20.04, 6.15, 18.65, 8.16, 17.4) round(log_importance_weights(alpha, beta),2)

## [1] -8.95 -23.47 -6.02 -8.13 -16.61 -14.57

d)     Implement a function for computing normalized importance ratios from the unnormalized log ratios in c). In other words, exponentiate the log ratios and scale them such that they sum to one. Explain in words what is the e ect of exponentiating and scaling so that sum is one. Below is a test example, the functions can also be tested with markmyassignment.

Note! The values below are only a test case. In this d) part, you only need to report the source code of your function, as it will be needed in later parts.

alpha

## [1] 1.896 -3.600 0.374 0.964 -3.123 -1.581

beta

## [1] 24.76 20.04 6.15 18.65 8.16 17.40

round(normalized_importance_weights(alpha = alpha, beta = beta),3)

## [1] 0.045 0.000 0.852 0.103 0.000 0.000

e)     Sample 4000 draws of α and β from the prior distribution from a). Compute and plot a histogram of the 4000 normalized importance ratios. Use the functions you implemented in c) and d).

f)      Using the importance ratios, compute the importance sampling e ective sample size Se and report it.

Note! The values below are only a test case, you need to use 4000 draws for alpha and beta in the nal report.

alpha

## [1] 1.896 -3.600 0.374 0.964 -3.123 -1.581

beta

## [1] 24.76 20.04 6.15 18.65 8.16 17.40

round(S_eff(alpha = alpha, beta = beta), 3)

## [1] 1.354

Hint! Equation (10.4) in the course book.

Note! BDA3 1st (2013) and 2nd (2014) printing have an error for w˜(θs) used in the e ective sample size equation (10.4). The normalized weights equation should not have the multiplier S (the normalized weights should sum to one). Errata for the book can be found here: http://www.stat.columbia.edu/ ~gelman/book/errata_bda3.txt. The later printings and slides have the correct equation.

g)     Explain in your own words what the importance sampling e ective sample size represents. Also explain how the e ective sample size is seen in the histogram of the weights that you plotted in e).

h)     Implement a function for computing the posterior mean using importance sampling, and compute the mean using your 4000 draws. Explain in your own words the computation for importance sampling. Below is an example how the function would work with the example values for alpha and beta above. Report the means for alpha and beta, and also the Monte Carlo standard errors (MCSEs) for the mean estimates. Report the number of digits for the means based on the MCSEs.

Note! The values below are only a test case, you need to use 4000 draws for alpha and beta in the nal report.

 Hint! Use the same equation for the MCSE of E[θ] as earlier (pVar[θ]/S), but now replace S with Se . To compute Var[θ] with importance sampling, use the identity Var[θ] = E[θ2] − E[θ]2.

alpha

## [1] 1.896 -3.600 0.374 0.964 -3.123 -1.581 beta

## [1] 24.76 20.04 6.15 18.65 8.16 17.40 round(posterior_mean(alpha = alpha, beta = beta),3)

## [1] 0.503 8.275