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. You can download the qmd-files (full, simple) 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 to check your functions at the end of the notebook:
Loading required package: markmyassignmentassignment_path = paste("https://github.com/avehtari/BDA_course_Aalto/",
"blob/master/assignments/tests/assignment3.yml", sep="")
set_assignment(assignment_path)Assignment set:
assignment3: Bayesian Data Analysis: Assignment 3
The assignment contain the following (6) tasks:
- mu_point_est
- mu_interval
- mu_pred_interval
- mu_pred_point_est
- posterior_odds_ratio_point_est
- posterior_odds_ratio_intervalThe following installs and loads the aaltobda package:
if(!require(aaltobda)){
install.packages("aaltobda", repos = c("https://avehtari.github.io/BDA_course_Aalto/", getOption("repos")))
library(aaltobda)
}Loading required package: aaltobdaThe following installs and loads the latex2exp package, which allows us to use LaTeX in plots:
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, which allows us to use its rvar Random Variable Datatype:
Loading required package: posteriorThis is posterior version 1.4.0
Attaching package: 'posterior'The following object is masked from 'package:aaltobda':
mcse_quantileThe following objects are masked from 'package:stats':
mad, sd, varThe following objects are masked from 'package:base':
%in%, matchThe following installs and loads the ggdist package for advanced plotting functions:
Loading required package: ggplot2ggplot2::theme_set(theme_minimal(base_size = 14))
if(!require(ggdist)){
install.packages("ggdist")
library(ggdist)
}Loading required package: ggdistThis 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.
Loading the library and the data.
data("windshieldy1")
# The data are now stored in the variable `windshieldy1`.
# The below displays the data:
windshieldy1[1] 13.357 14.928 14.896 15.297 14.820 12.067 14.824 13.865 17.447The below data is only for the tests, you need to change to the full data windshieldy1 when reporting your results.
Write your answers here!
Write your answers and code here!
Keep the below name and format for the functions to work with markmyassignment:
# 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\).
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)$)')
)
Write your answers and code here!
Keep the below name and format for the functions to work with markmyassignment:
# 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}\).
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)$)')
)
Write your answers here!
Write your answers and code here!
The below data is only for the tests:
Keep the below name and format for the functions to work with markmyassignment:
# 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)
}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” 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, such that rvar(p0) prints as 0.05 ± 0.021 and rvar(p1) prints as 0.1 ± 0.029.
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.” The code below sets up an R data.frame() with the draws in p0 and p1 wrapped in an rvar, and uses that data frame to visualize the draws using ggdist, an R package building on ggplot2 and “designed for both frequentist and Bayesian uncertainty visualization”.
The below plot, Figure 3 uses ggdist’s stat_dotsinterval(), which by default visualizes
rvar’s median and central 66% and 95% intervals using a black dot and lines of varying thicknesses as when using ggdist’s stat_pointinterval() andrvar’s draws using grey dots as when using ggdist’s stat_dots():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()
rvars make it easy to compute functions of random variables, such as
r0 - r1 computes an rvar which prints as -0.05 ± 0.037, indicating the sample mean and the sample standard deviation of the difference of the probabilities of death,r0 * r1 computes an rvar which prints as 0.005 ± 0.0026 which in this case has no great interpretation, orBelow, in Figure 4, we compute the odds ratios using the rvars and visualize its median, central intervals and draws, as above in Figure 3:
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 Figure 4 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.
Write your answers and code here!
Loading the library and the data.
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)[1] 15.980 14.206 16.011 17.250 15.993 15.722Write your answers here!
Write your answers and code here!
Write your answers here!
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:
✔ | F W S OK | Context
⠏ | 0 | task-1-subtask-1-tests
⠏ | 0 | mu_point_est()
✖ | 1 3 | mu_point_est()
────────────────────────────────────────────────────────────────────────────────
Failure ('test-task-1-subtask-1-tests.R:17:3'): mu_point_est()
mu_point_est(data = 1:10) not equivalent to 5.5.
1/1 mismatches
[1] 14.5 - 5.5 == 9
Error: Incorrect result for vector 1:10
────────────────────────────────────────────────────────────────────────────────
⠏ | 0 | task-2-subtask-1-tests
⠏ | 0 | mu_interval()
✖ | 2 3 | mu_interval()
────────────────────────────────────────────────────────────────────────────────
Failure ('test-task-2-subtask-1-tests.R:15:3'): mu_interval()
mu_interval(data = 1:10, prob = 0.95) not equivalent to c(3.3, 7.7).
2/2 mismatches (average diff: 9)
[1] 13.3 - 3.3 == 10
[2] 15.7 - 7.7 == 8
Error: Incorrect result for 1:10 using a 95% interval.
Failure ('test-task-2-subtask-1-tests.R:16:3'): mu_interval()
mu_interval(data = 1:10, prob = 0.8) not equivalent to c(4.2, 6.8).
2/2 mismatches (average diff: 9)
[1] 13.3 - 4.2 == 9.1
[2] 15.7 - 6.8 == 8.9
Error: Incorrect result for 1:10 using a 80% interval.
────────────────────────────────────────────────────────────────────────────────
⠏ | 0 | task-3-subtask-1-tests
⠏ | 0 | mu_pred_interval()
✖ | 2 3 | mu_pred_interval()
────────────────────────────────────────────────────────────────────────────────
Failure ('test-task-3-subtask-1-tests.R:16:3'): mu_pred_interval()
mu_pred_interval(1:10, 0.95) not equivalent to c(-1.7, 12.7).
2/2 mismatches (average diff: 9)
[1] 11.8 - -1.7 == 13.5
[2] 17.2 - 12.7 == 4.5
Error: Incorrect result for 1:10, (95%)
Failure ('test-task-3-subtask-1-tests.R:18:3'): mu_pred_interval()
mu_pred_interval(1:10, 0.8) not equivalent to c(1.1, 9.9).
2/2 mismatches (average diff: 9)
[1] 11.8 - 1.1 == 10.7
[2] 17.2 - 9.9 == 7.3
Error: Incorrect result for 1:10, (80%)
────────────────────────────────────────────────────────────────────────────────
⠏ | 0 | task-4-subtask-1-tests
⠏ | 0 | mu_pred_point_est()
✖ | 1 3 | mu_pred_point_est()
────────────────────────────────────────────────────────────────────────────────
Failure ('test-task-4-subtask-1-tests.R:16:3'): mu_pred_point_est()
mu_pred_point_est(data = 1:10) not equivalent to 5.5.
1/1 mismatches
[1] 14.5 - 5.5 == 9
Error: Incorrect result for vector 1:10
────────────────────────────────────────────────────────────────────────────────
⠏ | 0 | task-5-subtask-1-tests
⠏ | 0 | posterior_odds_ratio_point_est()
✖ | 3 3 | posterior_odds_ratio_point_est()
────────────────────────────────────────────────────────────────────────────────
Failure ('test-task-5-subtask-1-tests.R:18:3'): posterior_odds_ratio_point_est()
posterior_odds_ratio_point_est(p1, p0) not equivalent to 0.5307909.
1/1 mismatches
[1] 2.65 - 0.531 == 2.12
Error: Incorrect result for p1 = rbeta(100000, 5, 95), p0 = rbeta(100000, 10, 90).
Failure ('test-task-5-subtask-1-tests.R:24:3'): posterior_odds_ratio_point_est()
posterior_odds_ratio_point_est(p0, p1) not equivalent to 1.876519.
1/1 mismatches
[1] 2.65 - 1.88 == 0.774
Error: Incorrect result for p0 = rbeta(100000, 5, 5), p1 = rbeta(100000, 5, 6).
Failure ('test-task-5-subtask-1-tests.R:25:3'): posterior_odds_ratio_point_est()
posterior_odds_ratio_point_est(p1, p0) not equivalent to 1.24902.
1/1 mismatches
[1] 2.65 - 1.25 == 1.4
Error: Incorrect result for p0 = rbeta(100000, 5, 6), p1 = rbeta(100000, 5, 5).
────────────────────────────────────────────────────────────────────────────────
⠏ | 0 | task-6-subtask-1-tests
⠏ | 0 | posterior_odds_ratio_interval()
✖ | 6 2 | posterior_odds_ratio_interval()
────────────────────────────────────────────────────────────────────────────────
Failure ('test-task-6-subtask-1-tests.R:17:3'): posterior_odds_ratio_interval()
posterior_odds_ratio_interval(p0, p1, 0.9) not equivalent to c(0.875367, 6.05911).
2/2 mismatches (average diff: 0.719)
[1] 0.68 - 0.875 == -0.196
[2] 7.30 - 6.059 == 1.242
Error: Incorrect result for p0 = rbeta(100000, 5, 95), p1 = rbeta(100000, 10, 90) and 90%.
Failure ('test-task-6-subtask-1-tests.R:19:3'): posterior_odds_ratio_interval()
posterior_odds_ratio_interval(p1, p0, 0.9) not equivalent to c(0.1650407, 1.142378).
2/2 mismatches (average diff: 3.34)
[1] 0.68 - 0.165 == 0.515
[2] 7.30 - 1.142 == 6.159
Error: Incorrect result for p1 = rbeta(100000, 5, 95), p0 = rbeta(100000, 10, 90) and 90%.
Failure ('test-task-6-subtask-1-tests.R:21:3'): posterior_odds_ratio_interval()
posterior_odds_ratio_interval(p1, p0, 0.8) not equivalent to c(0.2086461, 0.9392956).
2/2 mismatches (average diff: 3.42)
[1] 0.68 - 0.209 == 0.471
[2] 7.30 - 0.939 == 6.362
Error: Incorrect result for p1 = rbeta(100000, 5, 95), p0 = rbeta(100000, 10, 90) and 80%.
Failure ('test-task-6-subtask-1-tests.R:27:3'): posterior_odds_ratio_interval()
posterior_odds_ratio_interval(p0, p1, 0.9) not equivalent to c(0.2714472, 5.5970131).
2/2 mismatches (average diff: 1.06)
[1] 0.68 - 0.271 == 0.408
[2] 7.30 - 5.597 == 1.705
Error: Incorrect result for p0 = rbeta(100000, 5, 5), p1 = rbeta(100000, 5, 6) and 90%.
Failure ('test-task-6-subtask-1-tests.R:29:3'): posterior_odds_ratio_interval()
posterior_odds_ratio_interval(p1, p0, 0.9) not equivalent to c(0.1786667, 3.6839583).
2/2 mismatches (average diff: 2.06)
[1] 0.68 - 0.179 == 0.501
[2] 7.30 - 3.684 == 3.618
Error: Incorrect result for p0 = rbeta(100000, 5, 6), p1 = rbeta(100000, 5, 5) and 90%.
Failure ('test-task-6-subtask-1-tests.R:31:3'): posterior_odds_ratio_interval()
posterior_odds_ratio_interval(p1, p0, 0.8) not equivalent to c(0.252972, 2.633517).
2/2 mismatches (average diff: 2.55)
[1] 0.68 - 0.253 == 0.427
[2] 7.30 - 2.634 == 4.668
Error: Incorrect result for p0 = rbeta(100000, 5, 6), p1 = rbeta(100000, 5, 5) and 80%.
────────────────────────────────────────────────────────────────────────────────
Maximum number of failures exceeded; quitting at end of file.
ℹ Increase this number with (e.g.) `testthat::set_max_fails(Inf)`
══ Results ═════════════════════════════════════════════════════════════════════
Duration: 0.5 s
[ FAIL 15 | WARN 0 | SKIP 0 | PASS 17 ]
══ Terminated early ════════════════════════════════════════════════════════════