---
title: "Assignment 3"
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 3](assignment3.html). You can download the qmd-files ([full](https://avehtari.github.io/BDA_course_Aalto/assignments/template3.qmd), [simple](https://avehtari.github.io/BDA_course_Aalto/assignments/simple_template3.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/assignment3.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)
}
```

:::
::::



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

## Showcase: Setting up advanced packages (`posterior` and `ggdist`)


*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 block showcases advanced tools, which you will be allowed and expected to use after this assignment.*
**For now, you should solve the assignment without the tools showcased herein.**




The following installs and loads the [`posterior` package](https://mc-stan.org/posterior/index.html),
which allows us to use its [`rvar` Random Variable Datatype
](https://mc-stan.org/posterior/articles/rvar.html):
```{r}
if(!require(posterior)){
    install.packages("posterior")
    library(posterior)
}
```

The following installs and loads the [`ggdist` package](https://mjskay.github.io/ggdist/)
for advanced plotting functions:
```{r}
if(!require(ggplot2)){
    install.packages("ggplot2")
    library(ggplot2)
}
ggplot2::theme_set(theme_minimal(base_size = 14))
if(!require(ggdist)){
    install.packages("ggdist")
    library(ggdist)
}
```



*This block showcases advanced tools, which you will be allowed and expected to use after this assignment.*
**For now, you should solve the assignment without the tools showcased herein.**



:::
::::


# Inference for normal mean and deviation (3 points)

Loading the library and the data.
``` {r}
data("windshieldy1")
# The data are now stored in the variable `windshieldy1`.
# The below displays the data:
windshieldy1
```
The below data is **only for the tests**, you need to change to the
full data `windshieldy1` when reporting your results.
``` {r}
windshieldy_test <- c(13.357, 14.928, 14.896, 14.820)
```



## (a)


Write your answers here!



## (b)


Write your answers and code here!

**Keep the below name and format for the functions to work with `markmyassignment`:**
```{r}
# Useful functions: mean(), length(), sqrt(), sum()
# and qtnew(), dtnew() (from aaltobda)

mu_point_est <- function(data) {
    # Do computation here, and return as below.
    # This is the correct return value for the test data provided above.
    14.5

}
mu_interval <- function(data, prob = 0.95) {
    # Do computation here, and return as below.
    # This is the correct return value for the test data provided above.
    c(13.3, 15.7)

}
```

You can plot the density as below if you implement `mu_pdf` to compute the PDF of the posterior $p(\mu|y)$ of the average hardness $\mu$.

```{r}
#| label: fig-2b-density
#| fig-cap: PDF of the posterior $p(\mu|y)$ of the average hardness $\mu$
mu_pdf <- function(data, x){
    # Compute necessary parameters here.
    # These are the correct parameters for `windshieldy_test`
    # with the provided uninformative prior.
    df = 3
    location = 14.5
    scale = 0.3817557
    # Use the computed parameters as below to compute the PDF:

    dtnew(x, df, location, scale)
}

x_interval = mu_interval(windshieldy1, .999)
lower_x = x_interval[1]
upper_x = x_interval[2]
x = seq(lower_x, upper_x, length.out=1000)
plot(
    x, mu_pdf(windshieldy1, x), type="l",
    xlab=TeX(r'(average hardness $\mu$)'),
    ylab=TeX(r'(PDF of the posterior $p(\mu|y)$)')
)
```



## (c)


Write your answers and code here!

**Keep the below name and format for the functions to work with `markmyassignment`:**
```{r}
# Useful functions: mean(), length(), sqrt(), sum()
# and qtnew(), dtnew() (from aaltobda)

mu_pred_point_est <- function(data) {
    # Do computation here, and return as below.
    # This is the correct return value for the test data provided above.
    14.5

}
mu_pred_interval <- function(data, prob = 0.95) {
    # Do computation here, and return as below.
    # This is the correct return value for the test data provided above.
    c(11.8, 17.2)

}
```

You can plot the density as below if you implement `mu_pred_pdf` to compute the PDF of the posterior predictive $p(\tilde{y}|y)$ of a new hardness observation $\tilde{y}$.

```{r}
#| label: fig-2c-density
#| fig-cap: PDF of the posterior predictive $p(\tilde{y}|y)$ of a new hardness observation $\tilde{y}$
mu_pred_pdf <- function(data, x){
    # Compute necessary parameters here.
    # These are the correct parameters for `windshieldy_test`
    # with the provided uninformative prior.
    df = 3
    location = 14.5
    scale = 0.8536316
    # Use the computed parameters as below to compute the PDF:

    dtnew(x, df, location, scale)
}

x_interval = mu_pred_interval(windshieldy1, .999)
lower_x = x_interval[1]
upper_x = x_interval[2]
x = seq(lower_x, upper_x, length.out=1000)
plot(
    x, mu_pred_pdf(windshieldy1, x), type="l",
    xlab=TeX(r'(new hardness observation $\tilde{y}$)'),
    ylab=TeX(r'(PDF of the posterior predictive $p(\tilde{y}|y)$)')
)
```


# Inference for the difference between proportions (3 points)

## (a)


Write your answers here!



## (b)


Write your answers and code here!

The below data is **only for the tests**:
``` {r}
set.seed(4711)
ndraws = 1000
p0 = rbeta(ndraws, 5, 95)
p1 = rbeta(ndraws, 10, 90)

```

**Keep the below name and format for the functions to work with `markmyassignment`:**

```{r}
# Useful function: mean(), quantile()

posterior_odds_ratio_point_est <- function(p0, p1) {
    # Do computation here, and return as below.
    # This is the correct return value for the test data provided above.
    2.650172

}
posterior_odds_ratio_interval <- function(p0, p1, prob = 0.95) {
    # Do computation here, and return as below.
    # This is the correct return value for the test data provided above.
    c(0.6796942,7.3015964)

}
```



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

## Showcase: advanced tools (`posterior`'s `rvar`, `ggdist`'s `stat_dotsinterval`)


*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 block showcases advanced tools, which you will be allowed and expected to use after this assignment.*
**For now, you should solve the assignment without the tools showcased herein.**




The `posterior` package's random variable datatype `rvar` is a
["sample-based representation of random variables"](https://mc-stan.org/posterior/articles/rvar.html#:~:text=sample%2Dbased%20representation%20of%20random%20variables)
which makes handling of random samples (of draws) such as the ones contained in the above variables `p0` and `p1` easier.
[By default, it prints as the mean and standard deviation of the draws](https://mc-stan.org/posterior/articles/rvar.html#:~:text=The%20default%20display%20of%20an%20rvar%20shows%20the%20mean%20and%20standard%20deviation%20of%20each%20element%20of%20the%20array.), **such that `rvar(p0)` prints as `r rvar(p0)` and `rvar(p1)` prints as `r rvar(p1)`**.

The datatype is ["designed to [...] be able to be used inside `data.frame()`s and `tibble()`s, and to be used with distribution visualizations in the ggdist package."](https://mc-stan.org/posterior/articles/rvar.html#:~:text=designed%20to%20interoperate%20with%20vectorized%20distributions%20in%20the%20distributional%20package%2C%20to%20be%20able%20to%20be%20used%20inside%20data.frame()s%20and%20tibble()s%2C%20and%20to%20be%20used%20with%20distribution%20visualizations%20in%20the%20ggdist%20package.)
The code below sets up an [R `data.frame()`](https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/data.frame) with the draws in `p0` and `p1` wrapped in an `rvar`, and uses that data frame to visualize the draws using [`ggdist`](https://mjskay.github.io/ggdist/index.html), an R package building on [`ggplot2`](https://ggplot2.tidyverse.org/) and ["designed for both frequentist and Bayesian uncertainty visualization"](https://mjskay.github.io/ggdist/index.html#:~:text=designed%20for%20both%20frequentist%20and%20Bayesian%20uncertainty%20visualization).

The below plot, @fig-showcase-probabilities uses `ggdist`'s [`stat_dotsinterval()`](https://mjskay.github.io/ggdist/articles/dotsinterval.html), which by default visualizes

* [an `rvar`'s median and central 66% and 95% intervals](https://mjskay.github.io/ggdist/reference/stat_dotsinterval.html#:~:text=point_interval%20%3D%20%22median_qi%22%2C%0A%20%20.width%20%3D%20c(0.66%2C%200.95)%2C) using a black dot and lines of varying thicknesses as when using `ggdist`'s [`stat_pointinterval()`](https://mjskay.github.io/ggdist/reference/stat_pointinterval.html#examples) and
* an `rvar`'s draws using grey dots as when using `ggdist`'s [`stat_dots()`](https://mjskay.github.io/ggdist/reference/stat_dots.html#examples):

```{r}
#| label: fig-showcase-probabilities
#| fig-cap: Probabilities of death for the two patient groups.
r0 = rvar(p0)
r1 = rvar(p1)
ggplot(data.frame(
    rv_name=c("control", "treatment"), rv=c(r0, r1)
)) +
    aes(xdist=rv, y=rv_name) +
    labs(x="probabilities of death", y="patient group") +
    stat_dotsinterval()
```

`rvar`s make it easy to compute functions of random variables, such as

* differences, e.g. $p_0 - p_1$: `r0 - r1` computes an `rvar` which prints as `r r0 - r1`, indicating the **sample mean** and the **sample standard deviation** of the difference of the probabilities of death,
* products, e.g. $p_0 \, p_1$: `r0 * r1` computes an `rvar` which prints as `r r0 * r1` which in this case
has no great interpretation, or
* the odds ratios needed in task 3.b).

Below, in @fig-showcase-odds-ratios, we compute the odds ratios using the `rvar`s and visualize its median, central intervals and draws, as above in @fig-showcase-probabilities:
```{r}
#| label: fig-showcase-odds-ratios
#| fig-cap: Odds ratios of the two patient groups.
rodds_ratio = (r1/(1-r1))/(r0/(1-r0))
ggplot(data.frame(
    rv=c(rodds_ratio)
)) +
    aes(xdist=rv) +
    labs(x="odds ratio", y="relative amount of draws") +
    stat_dotsinterval()
```

You can use @fig-showcase-odds-ratios to visually check whether the answers you computed for 3.b) make sense.



*This block showcases advanced tools, which you will be allowed and expected to use after this assignment.*
**For now, you should solve the assignment without the tools showcased herein.**



:::
::::



## (c)


Write your answers and code here!


# Inference for the difference between normal means (3 points)

Loading the library and the data.
``` {r}
data("windshieldy2")
# The new data are now stored in the variable `windshieldy2`.
# The below displays the first few rows of the new data:
head(windshieldy2)
```



## (a)


Write your answers here!



## (b)


Write your answers and code here!
```{r}
# Useful functions: mean(), length(), sqrt(), sum(),
# rtnew() (from aaltobda), quantile() and hist().
```



## (c)


Write your answers here!



:::: {.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()
```

:::
::::
