---
title: "Assignment 5"
author: anonymous # <-- hand in anonymously
format:
  html:
    toc: true
    code-tools: true
    code-line-numbers: true
    number-sections: true
    mainfont: Georgia, serif
    page-layout: article
  pdf:
    geometry:
    - left=1cm,top=1cm,bottom=1cm,right=7cm
    number-sections: true
    code-annotations: none
editor: source
---


# General information

:::: {.content-hidden when-format="pdf"}
::: {.callout-warning collapse=false}

## Setup


*This block will only be visible in your HTML output, but will be hidden when rendering to PDF with quarto for the submission.*
**Make sure that this does not get displayed in the PDF!**




This is the template for [assignment 5](assignment5.html). You can download the [qmd-file](https://avehtari.github.io/BDA_course_Aalto/assignments/template5.qmd) or copy the code from this rendered document after clicking on `</> Code` in the top right corner.

**Please replace the instructions in this template by your own text, explaining what you are doing in each exercise.**

The following will set-up [`markmyassignment`](https://github.com/MansMeg/markmyassignment) to check your functions at the end of the notebook:

```{r}
if(!require(markmyassignment)){
    install.packages("markmyassignment")
    library(markmyassignment)
}
assignment_path = paste("https://github.com/avehtari/BDA_course_Aalto/",
"blob/master/assignments/tests/assignment5.yml", sep="")
set_assignment(assignment_path)
```

The following installs and loads the `aaltobda` package:
```{r}
if(!require(aaltobda)){
    install.packages("aaltobda", repos = c("https://avehtari.github.io/BDA_course_Aalto/", getOption("repos")))
    library(aaltobda)
}
```
The following installs and loads the [`latex2exp` package](https://github.com/stefano-meschiari/latex2exp), which allows us to use LaTeX in plots:
```{r}
if(!require(latex2exp)){
    install.packages("latex2exp")
    library(latex2exp)
}
```
The following installs and loads the [`posterior` package](https://github.com/stan-dev/posterior) which imports the `rhat_basic()` function:
```{r}
if(!require(posterior)){
    install.packages("posterior")
    library(posterior)
}
```
The following installs and loads the [`ggplot2` package](https://ggplot2.tidyverse.org/) and the [`bayesplot` package](https://mc-stan.org/bayesplot/index.html)
```{r}
if(!require(ggplot2)){
    install.packages("ggplot2")
    library(ggplot2)
}
if(!require(bayesplot)){
    install.packages("bayesplot")
    library(bayesplot)
}
```

:::
::::


# Generalized linear model: Bioassay model with Metropolis algorithm

## (a)


Write your answers/code here!
```{r}
# Useful functions: runif, rnorm
# bioassaylp, dmvnorm (from aaltobda)

data("bioassay")
# Start by implementing a function called `density_ratio` to
# compute the density ratio function, $r$ in Eq. (11.1) in BDA3:
density_ratio <- function(alpha_propose, alpha_previous, beta_propose, beta_previous, x, y, n){
    # Do computation here, and return as below.
    # Below are the correct return values for two different calls of this function:

    # 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.305179

    # 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
    0.7661784

}
# Then implement a function called `metropolis_bioassay()` which
# implements the Metropolis algorithm using the `density_ratio()`:
metropolis_bioassay <- function(alpha_initial, beta_initial, alpha_sigma, beta_sigma, no_draws, x, y, n){
    # Do computation here, and return as below.
    # Below are "wrong" values (unlikely to actually occur)
    # in the "correct" format (such that they work with the plotting functions further down).
    data.frame(
        alpha=c(alpha_initial, alpha_initial+alpha_sigma, alpha_initial-alpha_sigma),
        beta=c(beta_initial, beta_initial+beta_sigma, beta_initial-beta_sigma)
    )

}
df = metropolis_bioassay(0, 0, 1, 1, 1000, bioassay$x, bioassay$y, bioassay$n)
```



## (b)


Write your answers/code here!

Have a look at [`bayesplot` trace plot examples](http://mc-stan.org/bayesplot/reference/MCMC-traces.html#examples) and tune your plot if wanted/needed. Don't forget to include a title/caption/description.

**The below example plot only includes a single chain, but your report should include a plot with multiple chains overlayed!**
```{r}
# Useful functions: mcmc_trace (from bayesplot)
mcmc_trace(df, pars=c("alpha", "beta"))
```



## (c)


Write your answers/code here!
```{r}
# Useful functions: rhat_basic (from posterior)
```



## (c)


Write your answers/code here!

Have a look at [`bayesplot` scatter plot examples](https://mc-stan.org/bayesplot/reference/MCMC-scatterplots.html#examples) and tune your plot if wanted/needed. Don't forget to include a title/caption/description.
```{r}
# Useful functions: mcmc_scatter (from bayesplot)
mcmc_scatter(df, pars=c("alpha", "beta"))
```



:::: {.content-hidden when-format="pdf"}
::: {.callout-warning collapse=false}

## markmyassignment


*This block will only be visible in your HTML output, but will be hidden when rendering to PDF with quarto for the submission.*
**Make sure that this does not get displayed in the PDF!**



The following will check the functions for which `markmyassignment` has been set up:

```{r}
mark_my_assignment()
```

:::
::::
