Description

Bayesian Statistics and Data Analysis

Assignment 5

General information
The recommended tool in this course is R (with the IDE R-Studio). You can download R here and R-Studio here. There are many tutorials, videos and introductions to R and R-Studio online. You can nd some initial hints from RStudio Education pages.
When working with R, we recommend writing the report using R markdown and the provided R markdown template. The template includes the formatting instructions and how to include code and gures.
Instead of R markdown, you can use other software to make the PDF report, but you should use the same instructions for formatting. These instructions are also available in the PDF produced from the R markdown template.
We supply a Google Colab notebook that you can also use for the assignments. We have included the installation of all necessary R packages; hence, this can be an alternative to using your own local computer. You can nd the notebook here. You can also open the notebook in Colab here.
Report all results in a single and anonymous pdf. Note that no other formats are allowed.
The course has its own R package bsda with data and functionality to simplify coding. To install the package, just run the following (upgrade=”never” skips question about updating other packages):
1. install.packages(“remotes”)
2. remotes::install_github(“MansMeg/BSDA”, subdir = “rpackage”, upgrade=”never”)
Many of the exercises can be checked automatically using the R package markmyassignment. you can nd information on how to install and use the package here. There is no need to include markmyassignment results in the report.
You can nd common questions and answers regarding the installation and technical problems in Frequently Asked Questions (FAQ).
You can nd deadlines and information on how to turn in the assignments in Studium.
If you have any suggestions or improvements to the course material, please post in the course chat feedback channel, create an issue, or submit a pull request to the public repository here.
It is mandatory to include the following parts in all assignments (these are included already in the template):
1. Time used for reading: How long time took the reading assignment (in hours)
2. Time used for the assignment: How long time took the basic assignment (in hours)
3. Good with assignment: Write one-two sentences of what you liked with the assignment/what we should keep for next year.
4. Things to improve in the assignment: Write one-two sentences of what you think can be improved in the assignment. Can something be clari ed further? Did you get stuck on stu unrelated to the content of the assignment etc.
To pass (G) the assignment, you need 70% of the total points. To pass with distinction (VG), you need 90% of the total points. See the grading information on the point allocations for each assignment.

Information on this assignment
This assignment is related to Chapters 10 and 11.
Reading instructions: Chapter 10 and 11 in BDA3, see reading instructions.
Reporting accuracy: For posterior statistics of interest, only report digits for which the Monte Carlo standard error (MCSE) is zero. Example: If you estimate E(µ) = 1.234 with MCSE(E(µ)) = 0.01, you should report E(µ) = 1.2.
To use markmyassignment for this assignment, run the following code in R:
library(markmyassignment) assignment_path <paste(“https://github.com/MansMeg/BSDA/”,
“blob/main/assignments/tests/assignment5.yml”, sep=””)
set_assignment(assignment_path)
# To check your code/functions, just run
mark_my_assignment()
Don’t include markmyassignment results in the report.
Generalized linear model: Bioassay with Metropolis
Metropolis algorithm: Replicate the computations for the bioassay example of section 3.7 in BDA3 using the Metropolis algorithm. The Metropolis algorithm is described in BDA3 Chapter 11.2. More information on the bioassay data can be found in Section 3.7 in BDA3.
1. Implement the Metropolis algorithm as an R function for the bioassay data. Use the
Gaussian prior
, where µ and .
a) Start by implementing a function called density_ratio to compute the density ratio function, r in Eq. (11.1) in BDA3. Below is an example on how the function should work. You can test the function using markmyassignment.
library(bsda) data(“bioassay”)
density_ratio(alpha_propose = 1.89, alpha_previous = 0.374,
beta_propose = 24.76, beta_previous = 20.04, x = bioassay$x, y = bioassay$y, n = bioassay$n)
## [1] 1.305179
density_ratio(alpha_propose = 0.374, alpha_previous = 1.89,
beta_propose = 20.04, beta_previous = 24.76, x = bioassay$x, y = bioassay$y, n = bioassay$n)
## [1] 0.7661784
Hint! Compute with log-densities. Reasons are explained on page 261 of BDA3.
Remember that p1/p0 = exp(log(p1) − log(p0)). For your convenience we have provided functions that will evaluate the log-likelihood for given α and β (see bioassaylp() in the bsda package). Notice that you still need to add the prior yourself and remember the unnormalized log posterior is simply the sum of loglikelihood and log-prior. For evaluating the log of the Gaussian prior you can use the function dmvnorm from package bsda.
b) Now implement a function called Metropolis_bioassay() which implements the Metropolis algorithm using the density_ratio().
Hint! Use a simple (normal) proposal distribution. Example proposals are α∗ ∼ N(αt−1,σ = 1) and β∗ ∼ N(βt−1,σ = 5). There is no need to try to nd optimal proposal but test some di erent values for the jump scale (σ). Remember to report the one you used. E cient proposals are dicussed in BDA3 p. 295 297 (not part of the course). In real-life a pre-run could be made with an automatic adaptive control to adapt the proposal distribution.
2. Include in the report the following:
a) Describe in your own words in one paragraph the basic idea of the Metropolis algorithm (see BDA3 Section 11.2).
b) The proposal distribution (related to jumping rule) you used. Describe brie y in words how you chose the nal proposal distribution you used for the reported results.
c) The initial points of your Metropolis chains (or the explicit mechanism for generating them).
d) Report the chain length or the number of iterations for each chain. Run the simulations long enough for approximate convergence (see BDA Section 11.4).
e) Report the warm-up length (see BDA Section 11.4).
f) The number of Metropolis chains used. It is important that multiple Metropolis chains are run for evaluating convergence (see BDA Section 11.4, and lecture
5.2).
g) Plot all chains for α in a single line-plot. Overlapping the chains in this way helps in visually assessing whether chains have converged or not.
h) Do the same for β.
3. In complex scenarios, visual assessment is not su cient and Rb is a more robust indicator of convergence of the Markov chains. Use Rb for convergence analysis. You can either use Eq. (11.4) in BDA3 or the more recent version described here. You should specify which Rb you used. In R the best choice is to use function Rhat from package rstan. Remember to remove the warm-up samples before computing Rb. Report the Rb values for α and β separately. Report the values for the proposal distribution you nally used.
a ) Describe brie y in your own words the basic idea of Rb and how to to interpret the obtained Rb values.
b ) Tell whether you obtained good Rb with rst try, or whether you needed to run more iterations or how did you modify the proposal distribution.
4. Plot the draws for α and β (scatter plot) and include this plot in your report. You can compare the results to Figure 3.3b in BDA3 to verify that your code gives sensible results. Notice though that the results in Figure 3.3b are generated from posterior with a uniform prior, so even when if your algorithm works perfectly, the results will look slightly di erent (although fairly similar).