diff --git a/..rds b/..rds new file mode 100644 index 0000000..1a2d260 Binary files /dev/null and b/..rds differ diff --git a/classes-ASHgithubb.Rproj b/classes-ASHgithubb.Rproj new file mode 100644 index 0000000..8e3c2eb --- /dev/null +++ b/classes-ASHgithubb.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX diff --git a/exercises/class-1-exercises.Rmd b/exercises/class-1-exercises.Rmd index 087db7c..cc0215d 100644 --- a/exercises/class-1-exercises.Rmd +++ b/exercises/class-1-exercises.Rmd @@ -8,10 +8,6 @@ date: "13-02-24" knitr::opts_chunk$set(echo = TRUE) ``` -```{r} -library(rethinking) -``` - # Exercises for Class 1 ## Supporting code diff --git a/exercises/class-1-exercises_Anne.Rmd b/exercises/class-1-exercises_Anne.Rmd new file mode 100644 index 0000000..89f0d7c --- /dev/null +++ b/exercises/class-1-exercises_Anne.Rmd @@ -0,0 +1,626 @@ +--- +title: "Exercises for class 1" +output: html_document +date: "13-02-24" +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) + +#Answers: https://sr2-solutions.wjakethompson.com/bayesian-inference + +#git push origin main +#git remote add upstream https://github.com/methods-4-f24/methods-4-classes +#git fetch upstream +#git merge upstream/main +#git add . +#git commit -m "comment" +#git push origin main + +``` + +```{r} +install.packages(c("coda","mvtnorm","devtools", "usethis")) +devtools::install_github("rmcelreath/rethinking",ref="Experimental") +``` + +```{r} +pacman::p_load(dagitty) +library(rethinking) +``` + +# Exercises for Class 1 + +## Supporting code + +Below is code from Chapter 2 of the book that you can use to do the exercises. Not all the code you require is here though. Make sure to not rush it and take time to play around with each function. :)) + +```{r} +# Grid Approximation + +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +# define prior +prior = rep(1, 20) +#prior <- ifelse( p_grid < 0.5 , 0 , 1 ) +#prior <- exp( -5*abs( p_grid - 0.5 ) ) + +# compute likelihood at each value in grid +likelihood = dbinom(6, size = 9, prob = p_grid) + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + +``` + +```{r} +## R code 2.6 +library(rethinking) +globe.qa <- quap( + alist( + W ~ dbinom( W+L ,p) , # binomial likelihood + p ~ dunif(0,1) # uniform prior + ) , + data=list(W=6,L=3) ) + +# display summary of quadratic approximation +precis( globe.qa ) + +## R code 2.7 +# analytical calculation +W <- 6 +L <- 3 +curve( dbeta( x , W+1 , L+1 ) , from=0 , to=1 ) +# quadratic approximation +curve( dnorm( x , 0.67 , 0.16 ) , lty=2 , add=TRUE ) +``` +The function precis presents a brief summary of the quadratic approximation. In this case, +it shows the posterior mean value of p = 0.67, which it calls the “Mean.” The curvature is +labeled “StdDev” This stands for standard deviation. This value is the standard deviation of +the posterior distribution, while the mean value is its peak. Finally, the last two values in the +precis output show the 89% percentile interval, which you’ll learn more about in the next +chapter. You can read this kind of approximation like: Assuming the posterior is Gaussian, it +is maximized at 0.67, and its standard deviation is 0.16. + +The blue (full line) curve is the analytical posterior and the black (dotted) curve is the quadratic approximation. The black curve +does alright on its left side, but looks pretty bad on its right side. It even assigns positive + +```{r} +## R code 2.8 +n_samples <- 1000 +p <- rep( NA , n_samples ) +p[1] <- 0.5 +W <- 6 +L <- 3 +for ( i in 2:n_samples ) { + p_new <- rnorm( 1 , p[i-1] , 0.1 ) + if ( p_new < 0 ) p_new <- abs( p_new ) + if ( p_new > 1 ) p_new <- 2 - p_new + q0 <- dbinom( W , W+L , p[i-1] ) + q1 <- dbinom( W , W+L , p_new ) + p[i] <- ifelse( runif(1) < q1/q0 , p_new , p[i-1] ) +} + +## R code 2.9 +dens( p , xlim=c(0,1) ) +curve( dbeta( x , W+1 , L+1 ) , lty=2 , add=TRUE ) + +``` + + + +## Exercises + +### Easy + +2E1. Which of the expressions below correspond to the statement: the probability of rain on Monday? +(1) Pr(rain) +(2) Pr(rain|Monday) +(3) Pr(Monday|rain) +(4) Pr(rain, Monday) / Pr(Monday) + +```{r} +#2 og #4 - kan også skrives sådan: Pr(Monday|rain) Pr(rain) / Pr(Monday) +``` + +2E2. Which of the following statements corresponds to the expression: Pr(Monday|rain)? +(1) The probability of rain on Monday. +(2) The probability of rain, given that it is Monday. +(3) The probability that it is Monday, given that it is raining. +(4) The probability that it is Monday and that it is raining. + +```{r} +#3 +``` + +2E3. Which of the following expressions below correspond to the statement: the probability that it is Monday, given that it is raining? +(1) Pr(Monday|rain) +(2) Pr(rain|Monday) +(3) Pr(rain|Monday) Pr(Monday) +(4) Pr(rain|Monday) Pr(Monday) / Pr(rain) +(5) Pr(Monday|rain) Pr(rain) / Pr(Monday) + +```{r} +#1 og #4 +``` + +2E4. The Bayesian statistician Bruno de Finetti (1906–1985) began his 1973 book on probability theory with the dedication: “PROBABILITY DOES NOT EXIST.” The capitals appeared in the original, so I imagine de Finetti wanted us to shout this statement. What he meant is that probability is a device for describing uncertainty from the perspective of an observer with limited knowledge; it has no objective reality. Discuss the globe tossing example from the chapter, in light of this statement. What does it mean to say “the probability of water is 0.7”? + +```{r} + +#With the current knowledge, the probability is 0.7, however, this does not apply in the real world, as we for sure do not have all possible evidence. + +``` + +### Medium + +2M1. Recall the globe tossing model from the chapter. Compute and plot the grid approximate posterior distribution for each of the following sets of observations. In each case, assume a uniform prior for p. (1) W, W, W (2) W, W, W, L (3) L, W, W, L, W, W, W + +```{r} +# Grid Approximation nr. 1 +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +# define prior +prior = rep(1, 20) + +# compute likelihood at each value in grid +likelihood = dbinom(3, size = 3, prob = p_grid) + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + + +# Grid Approximation nr. 2 +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +# define prior +prior = rep(1, 20) + +# compute likelihood at each value in grid +likelihood = dbinom(3, size = 4, prob = p_grid) #3 succes out of 4 + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + + +# Grid Approximation nr. 3 +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +# define prior +prior = rep(1, 20) + +# compute likelihood at each value in grid +likelihood = dbinom(5, size = 7, prob = p_grid) #5 succes of 7 + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") +``` +And here, the tidyverse way: +```{r} +library(tidyverse) + + +dist <- tibble(p_grid = seq(from = 0, to = 1, length.out = 20), + prior = rep(1, times = 20)) %>% + mutate(likelihood_1 = dbinom(3, size = 3, prob = p_grid), + likelihood_2 = dbinom(3, size = 4, prob = p_grid), + likelihood_3 = dbinom(5, size = 7, prob = p_grid), + across(starts_with("likelihood"), ~ .x * prior), + across(starts_with("likelihood"), ~ .x / sum(.x))) %>% + pivot_longer(cols = starts_with("likelihood"), names_to = "pattern", + values_to = "posterior") %>% + separate(pattern, c(NA, "pattern"), sep = "_", convert = TRUE) %>% + mutate(obs = case_when(pattern == 1L ~ "W, W, W", + pattern == 2L ~ "W, W, W, L", + pattern == 3L ~ "L, W, W, L, W, W, W")) + +ggplot(dist, aes(x = p_grid, y = posterior)) + + facet_wrap(vars(fct_inorder(obs)), nrow = 1) + + geom_line() + + geom_point() + + labs(x = "Proportion Water (p)", y = "Posterior Density") +``` + +2M2. Now assume a prior for p that is equal to zero when p < 0.5 and is a positive constant when p ≥ 0.5. Again compute and plot the grid approximate posterior distribution for each of the sets of observations in the problem just above. + +```{r} +# In the book, McElreath says that we can do much better than choosing a uniform prior for estimating the coverage of water. Is this what we are doing here? How does that impact our estimation? + +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +prob = p_grid + +# compute likelihood at each value in grid +likelihood = dbinom(3, size = 3, prob = p_grid) + +# define prior +prior = ifelse(prob<0.5,0,1) + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +prob = p_grid + +# compute likelihood at each value in grid +likelihood = dbinom(3, size = 4, prob = p_grid) + +# define prior +prior = ifelse(prob<0.5,0,1) + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + + +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +prob = p_grid + +# compute likelihood at each value in grid +likelihood = dbinom(5, size = 7, prob = p_grid) + +# define prior +prior = ifelse(prob<0.5,0,1) + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + + +#The tidyverse one +dist <- tibble(p_grid = seq(from = 0, to = 1, length.out = 20)) %>% + mutate(prior = case_when(p_grid < 0.5 ~ 0L, + TRUE ~ 1L), + likelihood_1 = dbinom(3, size = 3, prob = p_grid), + likelihood_2 = dbinom(3, size = 4, prob = p_grid), + likelihood_3 = dbinom(5, size = 7, prob = p_grid), + across(starts_with("likelihood"), ~ .x * prior), + across(starts_with("likelihood"), ~ .x / sum(.x))) %>% + pivot_longer(cols = starts_with("likelihood"), names_to = "pattern", + values_to = "posterior") %>% + separate(pattern, c(NA, "pattern"), sep = "_", convert = TRUE) %>% + mutate(obs = case_when(pattern == 1L ~ "W, W, W", + pattern == 2L ~ "W, W, W, L", + pattern == 3L ~ "L, W, W, L, W, W, W")) + +ggplot(dist, aes(x = p_grid, y = posterior)) + + facet_wrap(vars(fct_inorder(obs)), nrow = 1) + + geom_line() + + geom_point() + + labs(x = "Proportion Water (p)", y = "Posterior Density") + +``` + +2M3. Suppose there are two globes, one for Earth and one for Mars. The Earth globe is 70% covered in water. The Mars globe is 100% land. Further suppose that one of these globes—you don’t know which—was tossed in the air and produced a “land” observatiion. Assume that each globe was equally likely to be tossed. Show that the posterior probability that the globe was the Earth, conditional on seeing “land” (Pr(Earth|land)), is 0.23. + +```{r} + +# Use the Bayes formula! + +# probability of land, given Earth +p_le <- 0.3 + +# probability of land, given Mars +p_lm <- 1.0 + +# probability of Earth +p_e <- 0.5 + +# probability of land +p_l <- (p_e * p_le) + ((1 - p_e) * p_lm) + +# probability of Earth, given land (using Bayes' Theorem) +p_el <- (p_le * p_e) / p_l +p_el + + + +``` + +For the exercises below, I highly suggest you to grab a piece of paper and try to solve the problems in a 'visual' way. + +2M4. Suppose you have a deck with only three cards. Each card has only two sides, and each side is either black or white. One card has two black sides. The second card has one black and one white side. The third card has two white sides. Now suppose all three cards are placed in a bag and shuffled. Someone reaches into the bag and pulls out a card and places it flat on a table. A black side is shown facing up, but you don’t know the color of the side facing down. Show that the probability that the other side is also black is 2/3. Use the counting method (Section 2 of the chapter) to approach this problem. This means counting up the ways that each card could produce the observed data (a black side faceing up on the table). + +```{r} + +# Write your probability here. + +inner_tree <- 3/6 +middle_tree <- 2/3 + +# Use the Bayes formula! + +p_black_black <- 2/3 + +p_black_white <- 1/3 + +p_black <- 0.5 + +p_white <- (p_black * p_black_black) + ((1 - p_black) * p_black_white) + +# probability of Earth, given land (using Bayes' Theorem) +p_black_black1 <- (p_black_black * p_black) / p_white +p_black_black1 + +#Shown w. right terms +card_bb_likelihood <- 2 +card_bw_likelihood <- 1 +card_ww_likelihood <- 0 + +likelihood <- c(card_bb_likelihood, card_bw_likelihood, card_ww_likelihood) +prior <- c(1, 1, 1) +posterior <- likelihood * prior +posterior <- posterior / sum(posterior) + +posterior[1] + +``` + +2M5. Now suppose there are four cards: B/B, B/W, W/W, and another B/B. Again suppose a card is drawn from the bag and a black side appears face up. Again calculate the probability that the other side is black. + +```{r} + +p_black_black <- 4/5 + +p_black_white <- 1/3 + +p_black <- 5/8 + +p_white <- (p_black * p_black_black) + ((1 - p_black) * p_black_white) + +# probability of Earth, given land (using Bayes' Theorem) +p_black_black1 <- (p_black_black * p_black) / p_white +p_black_black1 + +#shown w. right terms +bb_likelihood <- 4 +bw_likelihood <- 1 +ww_likelihood <- 0 + +bb_prior <- 1 +bw_prior <- 1 +ww_prior <- 1 + +posterior_bb <- (bb_likelihood*bb_prior)/(bb_likelihood*bb_prior+bw_likelihood*bw_prior+ww_likelihood*ww_prior) +posterior_bb +``` + +2M6. Imagine that black ink is heavy, and so cards with black sides are heavier than cards with white sides. As a result, it’s less likely that a card with black sides is pulled from the bag. So again assume there are three cards: B/B, B/W, and W/W. After experimenting a number of times, you conclude that for every way to pull the B/B card from the bag, there are 2 ways to pull the B/W card and 3 ways to pull the W/W card. Again suppose that a card is pulled and a black side appears face up. Show that the probability the other side is black is now 0.5. Use the counting method, as before. + +```{r} +bb_likelihood <- 2/3 +bw_likelihood <- 1/3 +ww_likelihood <- 0 + +bb_prior <- 1/6 +bw_prior <- 2/6 +ww_prior <- 3/6 + +bb_posterior <- (bb_likelihood*bb_prior)/(bb_likelihood*bb_prior+bw_likelihood*bw_prior+ww_likelihood*ww_prior) +bb_posterior + +``` +2M7. Assume again the original card problem, with a single card showing a black side face up. Before looking at the other side, we draw another card from the bag and lay it face up on the table. The face that is shown on the new card is white. Show that the probability that the first card, the one showing a black side, has black on its other side is now 0.75. Use the counting method, if you can. Hint: Treat this like the sequence of globe tosses, countng all the ways to see each observation, for each possible first card. + +```{r} +# probability of black-black, given black +bb_b_likelihood <- 6/9 + +# probability of black-white, given black +bw_b_likelihood <- 3/9 + +# probability of black +b_prior <- 3/5 #due to one whit card is removed + +# probability of white (marginalization) +p_w <- (bb_b_likelihood * b_prior) + ((1 - b_prior) * bw_b_likelihood) + +# probability of Earth, given land (using Bayes' Theorem) +bb_b_posterior <- (bb_b_likelihood * b_prior) / p_w + +bb_b_posterior +``` +This task is best solved by drawing up the 'garden of forking paths' from the beginning of the chapter. Essentially, we count how many ways **each starting card** could have produced observation *W* on the second draw. + +Let's start with BB. + +```{r, warning=FALSE} + +library(dagitty) + +# Create a dagitty object +garden <- dagitty("dag { + BB -> B1 + BB -> B2 + B1 -> BW + B1 -> WW + B2 -> BW + B2 -> WW +}") + +# Plot it +plot(garden) + +``` + +As we can see above, if our first card was BB, there are 6 ways the second draw could have been *W* - 3 for one side of BB (B1) and 3 for the other side of BB (B2). + +Now let's assume the first card drawn was BW and count the ways for the second draw to be *W*. + +```{r, warning=FALSE} + +# Create our garden +garden <- dagitty("dag { + BW -> B + B -> WW +}") + +# Plot it +plot(garden) + +``` + +If the first card was BW, then the only way to get *W* on the second draw is by pulling the WW card out - either of it's side. + +Now we table our findings. + +| Cards | Ways to produce observation *W* on second draw | +|-------|------------------------------------------------| +| BB | 6 | +| BW | 2 | +| WW | 0 | + +Total ways to produce observation *W* on second draw: 8 + +Therefore, to get the probability that the first card drawn was BB is **6/8** or **0.75**. + + + + 2H1. Suppose there are two species of panda bear. Both are equally common in the wild and live in the sample places. They look exactly alike and eat the same food, and there is yet no genetic assay capable of telling them apart. They differ however in their family sizes. Species A gives birth to twins 10% of the time, otherwise birthing a single infant. Species births twins 20% of the time, ottherwise birthing singleton infants. Assume these numbers are known with certainty, from many years of field research. + Now suppose you are managing a captive panda breeding program. You have a newe female panda of unknown species, and she has just given birth to twins. What is the probability that her next birth will also be twins? + +# After first birth, likelihood of species A and B is equal to the rate the +# species give birth to twins +a_likelihood <- 0.1 +b_likelihood <- 0.2 + +# Next calculate the posterior probability that the panda belongs to each +# species, assume species are equally likely +likelihood <- c(a_likelihood, b_likelihood) +prior <- c(1, 1) +posterior <- likelihood * prior +posterior <- posterior / sum(posterior) +posterior +#> [1] 0.333 0.667 + +# The probability the next birth is twins is the probability the panda belongs +# to each species times the likelihood each species gives birth to twins +(posterior[1] * a_likelihood) + (posterior[2] * b_likelihood) +#> [1] 0.167 + + 2H2. Recall all the facts from the problem above. Now compute the probability that the panda we have is from species A, asssuming we have observed only the first birth at that it was twins. + +# probability of species A +p_a <- 0.5 + +# probability of twins, given species A +p_ta <- 0.1 + +# probability of twins, given species B +p_tb <- 0.2 + +# probability of twins +p_t <- (p_a * p_ta) + ((1 - p_a) * p_tb) + +# probability of species A, given twins (using Bayes' Theorem) +# (note this is equivalent to `posterior[1]` above) +p_at <- (p_ta * p_a) / p_t +p_at +#> [1] 0.333 + + 2H3. Continuing on from the previous problem, suppose the same panda mother has a second birth and that it is not twins, but a singleton infant. Compute the posterior probability that this panda is species A. + +# likelihood for each species is Pr(twins) * Pr(singleton) +a_likelihood <- 0.1 * (1 - 0.1) +b_likelihood <- 0.2 * (1 - 0.2) + +# compute posterior probabilities +likelihood <- c(a_likelihood, b_likelihood) +prior <- c(1, 1) +posterior <- likelihood * prior +posterior <- posterior / sum(posterior) + +posterior[1] +#> [1] 0.36 + + 2H4. A common boast of Bayesian statisticians is that Bayesian inference makes it easy to use all of the data, even if the data are of different types. + So suppose now that a veterinarian comes along who has a new genetic test that she claims can identify the species of our mother panda. But the test, like all tests, is imperfect. This is the information you have about the test: + + The probability it correctly identifies a species A panda is 0.8. + The probability it correctly identifies a species B panda is 0.65. + + The vet administers the test to your panda and tells you that the test is positive for species A. First ignore your previous information from the births and compute the posterior probability that your panda is species A. Then redo your calculation, now using the birth data as well. + +# use Bayes' Theorem to determine the probability of species A, given a positive +# test +p_ap <- (0.8 * 0.5) / ((0.5 * 0.8) + (0.5 * 0.35)) +p_ap +#> [1] 0.696 + + +# Now include test data with observed births +# likelihood for each species is Pr(twins) * Pr(singleton) +a_likelihood <- 0.1 * (1 - 0.1) +b_likelihood <- 0.2 * (1 - 0.2) + +# compute posterior probabilities, using test result as prior +likelihood <- c(a_likelihood, b_likelihood) +prior <- c(p_ap, (1 - p_ap)) +posterior <- likelihood * prior +posterior <- posterior / sum(posterior) + +posterior[1] +#> [1] 0.563 \ No newline at end of file diff --git a/exercises/class-1-solutions.Rmd b/exercises/class-1-solutions.Rmd new file mode 100644 index 0000000..f7fcc7b --- /dev/null +++ b/exercises/class-1-solutions.Rmd @@ -0,0 +1,226 @@ +--- +title: "" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +# Solutions to exercises for Class 1 + +## Supporting code + +Below is code from Chapter 2 of the book that you can use to do the exercises. Not all the code you require is here though. Make sure to not rush it and take time to play around with each function. :)) + +```{r} +# Grid Approximation + +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +# define prior +prior = rep(1, 20) + +# compute likelihood at each value in grid +likelihood = dbinom(3, size = 4, prob = p_grid) + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + +``` + +## Exercises + +### Medium + +#### 2M1. + +```{r} + +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +# define prior +prior = rep(1, 20) + +# compute likelihood at each value in grid +likelihood = dbinom(5, size = 7, prob = p_grid) +# change to dbinom(3, size = 3 +# change to dbinom(3, size = 4 +# change to dbinom(5, size = 7 + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + +``` + +#### 2M2. + +```{r} + +# define grid +p_grid = seq(from = 0, to = 1, length.out=20) + +# step prior, which codes our assumption that all probabilities before 0.5 are impossible +prior = ifelse(p_grid < 0.5, 0 , 1) + +# compute likelihood at each value in grid +likelihood = dbinom(5, size = 7, prob = p_grid) +# change to dbinom(3, size = 3 +# change to dbinom(3, size = 4 +# change to dbinom(5, size = 7 + +# compute product of likelihood and prior +unstd.posterior = likelihood * prior + +# standardize the posterior so it sums to 1 +posterior = unstd.posterior / sum(unstd.posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +mtext("20 points") + +``` + +#### 2M3. + +```{r} + +# Pr(Earth|Land) = Pr(Land|Earth) * Pr(Earth) / Pr(Land) + +answer = 0.3 * 0.5 / ((1 + 0.3) / 2) + +answer + +``` + +#### 2M4. + +Using the table similar to one on page 25. + +*B* - black, *W* - white + +| Cards | Ways to produce observation *B* | +|-------|---------------------------------| +| BB | 2 | +| BW | 1 | +| WW | 0 | + +Total ways to produce observation *B*: 3 + +Since BB card can produce *B* in two different ways (via each side) out of total possible ways, the probability that the called pulled out was indeed BB is **2/3** or **0.66**. + +#### 2M5. + +Here we just add another BB card. + +*B* - black, *W* - white + +| Cards | Ways to produce observation *B* | +|-------|---------------------------------| +| BB | 2 | +| BB | 2 | +| BW | 1 | +| WW | 0 | + +Total ways to produce observation *B*: 5 + +Now that we have two BB cards, we simply sum the number ways the can produce observation *B* and divide it by the total number of possible observations of *B*. So **4/5** or **0.8**. + +#### 2M6. + +Here we add another column to our table - our prior. + +*B* - black, *W* - white + +| Cards | Ways to produce observation *B* | Prior | Updated count | +|-------|---------------------------------|-------|---------------| +| BB | 2 | 1 | 2 | +| BW | 1 | 2 | 2 | +| WW | 0 | 3 | 0 | + +Total ways to produce observation *B*: 4 + +So both BW and BB cards are **equally likely** to produce observation *B*. + +#### 2M7. + +This task is best solved by drawing up the 'garden of forking paths' from the beginning of the chapter. Essentially, we count how many ways **each starting card** could have produced observation *W* on the second draw. + +Let's start with BB. + +```{r, warning=FALSE} + +library(dagitty) + +# Create a dagitty object +garden <- dagitty("dag { + BB -> B1 + BB -> B2 + B1 -> BW + B1 -> WW + B2 -> BW + B2 -> WW +}") + +# Plot it +plot(garden) + +``` + +As we can see above, if our first card was BB, there are 6 ways the second draw could have been *W* - 3 for one side of BB (B1) and 3 for the other side of BB (B2). + +Now let's assume the first card drawn was BW and count the ways for the second draw to be *W*. + +```{r, warning=FALSE} + +# Create our garden +garden <- dagitty("dag { + BW -> B + B -> WW +}") + +# Plot it +plot(garden) + +``` + +If the first card was BW, then the only way to get *W* on the second draw is by pulling the WW card out - either of it's side. + +Now we table our findings. + +| Cards | Ways to produce observation *W* on second draw | +|-------|------------------------------------------------| +| BB | 6 | +| BW | 2 | +| WW | 0 | + +Total ways to produce observation *W* on second draw: 8 + +Therefore, to get the probability that the first card drawn was BB is **6/8** or **0.75**. + + + + + + + + diff --git a/exercises/class-1-solutions.html b/exercises/class-1-solutions.html new file mode 100644 index 0000000..213cae6 --- /dev/null +++ b/exercises/class-1-solutions.html @@ -0,0 +1,672 @@ + + + + + + + + + + + + + +class-1-solutions.knit + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + +
+

Solutions to exercises for Class 1

+
+

Supporting code

+

Below is code from Chapter 2 of the book that you can use to do the +exercises. Not all the code you require is here though. Make sure to not +rush it and take time to play around with each function. :))

+
# Grid Approximation
+
+# define grid
+p_grid = seq(from = 0, to = 1, length.out=20)
+
+# define prior
+prior = rep(1, 20)
+
+# compute likelihood at each value in grid
+likelihood = dbinom(6, size = 9, prob = p_grid)
+
+# compute product of likelihood and prior
+unstd.posterior = likelihood * prior
+
+# standardize the posterior so it sums to 1
+posterior = unstd.posterior / sum(unstd.posterior)
+
+# plot the distribution
+plot(p_grid, posterior, type = "b",
+     xlab = "probability of water", ylab = "posterior_probability")
+mtext("20 points")
+

+
+
+

Exercises

+
+

Medium

+
+

2M1.

+
# define grid
+p_grid = seq(from = 0, to = 1, length.out=20)
+
+# define prior
+prior = rep(1, 20)
+
+# compute likelihood at each value in grid
+likelihood = dbinom(5, size = 7, prob = p_grid)
+# change to dbinom(3, size = 3
+# change to dbinom(3, size = 4
+# change to dbinom(5, size = 7
+
+# compute product of likelihood and prior
+unstd.posterior = likelihood * prior
+
+# standardize the posterior so it sums to 1
+posterior = unstd.posterior / sum(unstd.posterior)
+
+# plot the distribution
+plot(p_grid, posterior, type = "b",
+     xlab = "probability of water", ylab = "posterior_probability")
+mtext("20 points")
+

+
+
+

2M2.

+
# define grid
+p_grid = seq(from = 0, to = 1, length.out=20)
+
+# step prior, which codes our assumption that all probabilities before 0.5 are impossible
+prior = ifelse(p_grid < 0.5, 0 , 1)
+
+# compute likelihood at each value in grid
+likelihood = dbinom(5, size = 7, prob = p_grid)
+# change to dbinom(3, size = 3
+# change to dbinom(3, size = 4
+# change to dbinom(5, size = 7
+
+# compute product of likelihood and prior
+unstd.posterior = likelihood * prior
+
+# standardize the posterior so it sums to 1
+posterior = unstd.posterior / sum(unstd.posterior)
+
+# plot the distribution
+plot(p_grid, posterior, type = "b",
+     xlab = "probability of water", ylab = "posterior_probability")
+mtext("20 points")
+

+
+
+

2M3.

+
# Pr(Earth|Land) = Pr(Land|Earth) * Pr(Earth) / Pr(Land)
+
+answer = 0.3 * 0.5 / ((1 + 0.3) / 2)
+
+answer
+
## [1] 0.2307692
+
+
+

2M4.

+

Using the table similar to one on page 25.

+

B - black, W - white

+ + + + + + + + + + + + + + + + + + + + + +
CardsWays to produce observation B
BB2
BW1
WW0
+

Total ways to produce observation B: 3

+

Since BB card can produce B in two different ways (via each +side) out of total possible ways, the probability that the called pulled +out was indeed BB is 2/3 or 0.66.

+
+
+

2M5.

+

Here we just add another BB card.

+

B - black, W - white

+ + + + + + + + + + + + + + + + + + + + + + + + + +
CardsWays to produce observation B
BB2
BB2
BW1
WW0
+

Total ways to produce observation B: 5

+

Now that we have two BB cards, we simply sum the number ways the can +produce observation B and divide it by the total number of +possible observations of B. So 4/5 or +0.8.

+
+
+

2M6.

+

Here we add another column to our table - our prior.

+

B - black, W - white

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
CardsWays to produce observation BPriorUpdated count
BB212
BW122
WW030
+

Total ways to produce observation B: 4

+

So both BW and BB cards are equally likely to +produce observation B.

+
+
+

2M7.

+

This task is best solved by drawing up the ‘garden of forking paths’ +from the beginning of the chapter. Essentially, we count how many ways +each starting card could have produced observation +W on the second draw.

+

Let’s start with BB.

+
library(dagitty)
+
+# Create a dagitty object
+garden <- dagitty("dag {
+  BB -> B1
+  BB -> B2
+  B1 -> BW
+  B1 -> WW
+  B2 -> BW
+  B2 -> WW
+}")
+
+# Plot it
+plot(garden)
+
## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.
+

+

As we can see above, if our first card was BB, there are 6 ways the +second draw could have been W - 3 for one side of BB (B1) and 3 +for the other side of BB (B2).

+

Now let’s assume the first card drawn was BW and count the ways for +the second draw to be W.

+
# Create our garden
+garden <- dagitty("dag {
+  BW -> B
+  B -> WW
+}")
+
+# Plot it
+plot(garden)
+
## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.
+

+

If the first card was BW, then the only way to get W on the +second draw is by pulling the WW card out - either of it’s side.

+

Now we table our findings.

+ + + + + + + + + + + + + + + + + + + + + +
CardsWays to produce observation W on second draw
BB6
BW2
WW0
+

Total ways to produce observation W on second draw: 8

+

Therefore, to get the probability that the first card drawn was BB is +6/8 or 0.75.

+
+
+
+
+ + + + +
+ + + + + + + + + + + + + + + diff --git a/exercises/class-2-exercises.Rmd b/exercises/class-2-exercises.Rmd new file mode 100644 index 0000000..f672786 --- /dev/null +++ b/exercises/class-2-exercises.Rmd @@ -0,0 +1,273 @@ +--- +title: "Chapter 3" +output: html_document +date: "2024-01-08" +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +set.seed(100) + +pacman::p_load(dagitty) +``` + +```{r} +library(rethinking) +``` + + +# Exercises for Class 2 + +What's good my young Bayesians? Welcome to the exercises notebook for your 2nd class of Methods 4. As you might have noticed, a lot of Chapter 3 deals with summarizing the posterior. Summarizing the posterior is helpful, but as both Chris and Richard have re-iterated multiple times, you should **always report the full posterior distribution**. + +Therefore, I am providing the solutions to the *Easy* exercises aimed at summarizing the posterior, so you don't waste time scanning the book for the right lines of code. Play around with them, but don't wait too long to get to the **juice** of this notebook (and Chapter 3) - **simulation of the model's implied observations**. + +## Exercises + +### Easy. + +These problems use the samples form the posterior distribution for the globe tossing example. This code will give you a specific set of samples, so that you can check your answers correctly. + +```{r} + +p_grid = seq(from = 0, to = 1, length.out=1000) +prior = rep(1, 1000) +likelihood = dbinom(6, size = 9, prob = p_grid) #The mean of the amount of succes is 6/9 = 66% +posterior = likelihood * prior +posterior = posterior / sum(posterior) +samples = sample(p_grid, prob = posterior, size = 1e4, replace = TRUE) + +# Let's also visualize our samples, just to know what we are working with +dens(samples) + +``` + +#### 3E1. + +How much posterior probability lies below p = 0.2? + +```{r} + +sum(samples < 0.2) / 10000 + +``` + +#### 3E2. + +How much posterior probability lies above p = 0.8? + +```{r} + +sum(samples > 0.8) / 10000 + + +``` + +#### 3E3. + +How much posterior probability lies above p = 0.2 and below p = 0.8? + +```{r} + +sum(samples > 0.2 & samples < 0.8) / 10000 + +``` + +#### 3E4. + +20% of the posterior probability lies below which value of p? + +```{r} + +quantile(samples, 0.2) + +``` + +#### 3E5. + +20% of the posterior probability lies above which value of p? + +```{r} + +quantile(samples, 0.8) + +``` + +#### 3E6. + +Which values of p contain the narrowest interval equal to 66% of the posterior probability? + +```{r} + +HPDI(samples, prob = 0.66) + +``` + +#### 3E7. + +Which values of p contain the 66% of the posterior probability, assuming equal posterior probability both below and above the interval? + +```{r} + +PI(samples, prob = 0.66) + +``` + +### Sampling to simulate prediction + +Section **3.3** of this Chapter is, imo, the most important one. Simulating predictions is something you'll frequently doing going forward, so it's worth taking time to unpack the mechanics of it. + +```{r} +# First, play around with this function. What does it do? What does the output mean? + +rbinom(1, size = 9, prob = 0.7) #It is 70% likely that the function outputs "succes" for each of the 9 trials. This is done once. + +# Now let's scale it up. +dummy_w = rbinom(1e4, size = 100, prob = 0.2) #A trial results in succes 20% of the time. 100 trials are done. The count of succes trials is noted and put in the histogram. This is done 10000 times. + +simplehist(dummy_w) + +# What happens if we increase the size? The histogram becomes ore detailed. +# What about probability? The peak of the histogram is the probability/percentage of the size +``` + +Now let's do a posterior predictive check. Take a minute to contemplate what are we doing here. How does the resulting distribution compare to our sampled posterior? + +```{r} +w = rbinom(1e4, size = 9, prob = samples) +simplehist(w) + +# Compare to sampled posterior here - choose one of the summarization methods above! Is it more or less conservative? + +#? +``` + +The resulting distribution is for predictions, but it incorporates all of the uncertainty embodied in the posterior distribution for the parameter p. As a result, it is honest. While the model does a good job of predicting the data - the most likely observation is indeed the observed data - predictions are still quite spread out. + +### Medium + +### 3M1 & 3M2. + +Suppose the globe tossing data had turned out to be 8 water in 15 tosses. Construct the posterior distribution, using grid approximation. Use the same flat prior as before. + +Draw 10,000 samples from the grid approximation from above. Then use the samples to calculate the 90% HDPI for p. + +```{r} + +p_grid = seq(from = 0, to = 1, length.out=1000) +prior = rep(1, 1000) +likelihood = dbinom(8, size = 15, prob = p_grid) +posterior = likelihood * prior +posterior = posterior / sum(posterior) +samples = sample(p_grid, prob = posterior, size = 1e4, replace = TRUE) + +dens(samples) + +``` + + +```{r} + +#Posterier distribution +w1 = rbinom(1e4, size = 15, prob = samples) +simplehist(w1) + +#Which values of p contain the narrowest interval equal to 90% HDPI of the posterior probability? + +HPDI(samples, prob = 0.90) + +``` +Between p=0.35 and p=0.73 in the density plot, 90% of trials lie. + +#### 3M3. + +Construct a posterior predictive check for this model and data. This means simulate the distribution of samples, averaging over the posterior uncertainty in p. What is the probability of observing 8 water in 15 tosses? + +```{r} +#the samples are simulated in the above task +#Averaging is done by taking the sum and deviding by the amount +#dbinom is used to get the probability for this specific result instead of generating a lot of outcomes + +w2 <- sum(dbinom(8 , size=15 , prob=samples))/10000 +w2 + +``` + +#### 3M4. + +Using the posterior distribution constructed from the new (8/15) data, now calculate the probability of observing 6 water in 9 tosses. + +```{r} +#the samples are simulated in the above task +#Averaging is done by taking the sum and deviding by the amount +#dbinom is used to get the probability for this specific result instead of generating a lot of outcomes + +w3 <- sum(dbinom(6 , size=9 , prob=samples))/10000 +w3 + +``` + +#### 3M5. + +Start over at 3M1, but now use a prior that is zero below p=0.5 and a constant above p=0.5. This corresponds to prior information that a majority of the Earth’s surface is water. Repeat each problem above and compare the inferences. What difference does the better prior make? If it helps, compare inferences (using both priors) to the true value p = 0.7. + +```{r} + +# Your code here. + +``` + + +### Hard + +**Introduction**. The practice problems here use the data below. These data indicate the gender (male = 1, female = 0) of officially reported first and second born children in 100 two-children families. + +```{r} +birth1 <- c(1,0,0,0,1,1,0,1,0,1,0,0,1,1,0,1,1,0,0,0,1,0,0,0,1,0, +0,0,0,1,1,1,0,1,0,1,1,1,0,1,0,1,1,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0, +1,1,0,1,0,0,1,0,0,0,1,0,0,1,1,1,1,0,1,0,1,1,1,1,1,0,0,1,0,1,1,0, +1,0,1,1,1,0,1,1,1,1) +birth2 <- c(0,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0, +1,1,1,0,1,1,1,0,1,0,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1, +1,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0,0,0,1,0,0,0,1,1,0,0,1,0,0,1,1, +0,0,0,1,1,1,0,0,0,0) +``` + +So, for example, the first family in the data reported a boy (1) and then a girl (0). The second family reported a girl (0) and then a boy (1). The third family reported two girls. You can load these wo vectors into R's memory by typing: + +```{r} + +#library(rethinking) +data(homeworkch3) + +``` + +#### 3H1. + +Using grid approximation, compute the posterior distribution for the probability of a birth being a boy. Assume a uniform prior probability. Which parameter value maximizes the posterior probability? + +```{r} + +# Your code here. + +``` + +#### 3H2. + +Using the sample function, draw 10,000 random parameter values from the posterior distribution you calculated above. Use these samples to estimate the 50%, 89% and 97% highest posterior density intervals. + +```{r} + +# Your code here. + +``` + +#### 3H3. + +Use rbinom to simulate 10,000 replicates of 200 births. You should end up with 10,000 numbers, each one a count of boys out of 200 births. Compare the distribution of predicted numbers of boys to the actual count in the data (111 boys out of 200 births). There are many good ways to visualize the simulations, but the dens command (part of rethinking package) is probably the easiest way in this case. Does it look like the model fits the data well? That is, does the distribution of the predictions include the actual observation as a central, likely outcome? + +```{r} + +# Your code here. + +``` diff --git a/exercises/class-2-solutions.Rmd b/exercises/class-2-solutions.Rmd new file mode 100644 index 0000000..9aa76a8 --- /dev/null +++ b/exercises/class-2-solutions.Rmd @@ -0,0 +1,261 @@ +--- +title: "Chapter 3" +output: html_document +date: "2024-01-08" +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +set.seed(100) +``` + +# Exercises for Class 2 + +What's good my young Bayesians? Welcome to the exercises notebook for your 2nd class of Methods 4. As you might have noticed, a lot of Chapter 3 deals with summarizing the posterior. Summarizing the posterior is helpful, but as both Richard and Chris have re-iterated multiple times, you should **always report the full posterior distribution**. + +Therefore, I am providing the solutions to the *Easy* exercises aimed at summarizing the posterior, so you don't waste time scanning the book for the right lines of code. Play around with them, but don't wait too long to get to the **juice** of this notebook (and Chapter 3) - **simulation of the model's implied observations**. + +## Exercises + +### Easy. + +These problems use the samples form the posterior distribution for the globe tossing example. This code will give you a specific set of samples, so that you can check your answers correctly. + +```{r} + +p_grid = seq(from = 0, to = 1, length.out=1000) +prior = rep(1, 1000) +likelihood = dbinom(6, size = 9, prob = p_grid) +posterior = likelihood * prior +posterior = posterior / sum(posterior) +samples = sample(p_grid, prob = posterior, size = 1e4, replace = TRUE) + +# Let's also visualize our samples, just to know what we are working with +dens(samples) + +``` + +#### 3E1. + +How much posterior probability lies below p = 0.2? + +```{r} + +sum(samples < 0.2) / 10000 + +``` + +#### 3E2. + +How much posterior probability lies above p = 0.8? + +```{r} + +sum(samples > 0.8) / 10000 + + +``` + +#### 3E3. + +How much posterior probability lies above p = 0.2 and p = 0.8? + +```{r} + +sum(samples > 0.2 & samples < 0.8) / 10000 + +``` + +#### 3E4. + +20% of the posterior probability lies below which value of p? + +```{r} + +quantile(samples, 0.2) + +``` + +#### 3E5. + +20% of the posterior probability lies above which value of p? + +```{r} + +quantile(samples, 0.8) + +``` + +#### 3E6. + +Which values of p contain the narrowest interval equal to 66% of the posterior probability? + +```{r} + +HPDI(samples, prob = 0.66) + +``` + +#### 3E7. + +Which values of p contain the 66% of the posterior probability, assuming equal posterior probability both below and above the interval? + +```{r} + +PI(samples, prob = 0.66) + +``` + +### Sampling to simulate prediction + +Section **3.3** of this Chapter is, imo, the most important one. Simulating predictions is something you'll freqeuently doing going forward, so it's worth taking time to unpack the mechanics of it. + +```{r} +# First, play around with this function. What does it do? What does the output mean? +rbinom(1, size = 9, prob = 0.7) + +# Now let's scale it up. +dummy_w = rbinom(1e4, size = 9, prob = 0.1) +simplehist(dummy_w) + +# What happens if we increase the size? +# What about probability? +``` + +Now let's do a posterior predictive check. Take a minute to contemplate what are we doing here. How does the resulting distribution compare to our sampled posterior? + +```{r} +w = rbinom(1e4, size = 9, prob = samples) +hist(w) + +``` + +The resulting distribution is for predictions, but it incorporates all of the uncertainty embodied in the posterior distribution for the parameter p. As a result, it is honest. While the model does a good job of predicting the data - the most likely observation is indeed the observed data - predictions are still quite spread out. + +### Medium + +### 3M1 & 3M2. + +Suppose the globe tossing data had turned out to be 8 water in 15 tosses. Construct the posterior distribution, using grid approximation. Use the same flat prior as before. + +Draw 10,000 samples from the grid approximation from above. Then use the samples to calculate the 90% HDPI for p. + +```{r} + +p_grid = seq(from = 0, to = 1, length.out=1000) +prior = rep(1, 1000) +likelihood = dbinom(8, size = 15, prob = p_grid) +posterior = likelihood * prior +posterior = posterior / sum(posterior) +samples = sample(p_grid, prob = posterior, size = 1e4, replace = TRUE) + +dens(samples) + +# plot the distribution +# plot(p_grid, posterior, type = "b", +# xlab = "probability of water", ylab = "posterior_probability") + +``` + +```{r} +HPDI(samples, prob = 0.9) +``` + +#### 3M3. + +Construct a posterior predictive check for this model and data. This means simulate the distribution of samples, averaging over the posterior uncertainty in p. What is the probability of observing 8 water in 15 tosses? + +```{r} +w = rbinom(1e4, size = 15, prob = samples) +simplehist(w) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +``` + +#### 3M4. + +Using the posterior distribution constructed from the new (8/15) data, now calculate the probability of observing 6 water in 9 tosses. + +```{r} + +p_grid = seq(from = 0, to = 1, length.out=1000) +prior = posterior +likelihood = dbinom(6, size = 9, prob = p_grid) +posterior = likelihood * prior +posterior = posterior / sum(posterior) + +# plot the distribution +plot(p_grid, posterior, type = "b", + xlab = "probability of water", ylab = "posterior_probability") +``` + +#### 3M5. + +Start over at 3M1, but now use a prior that is zero below p=0.5 and a constant above p=0.5. This corresponds to prior information that a majority of the Earth’s surface is water. Repeat each problem above and compare the inferences. What difference does the better prior make? If it helps, compare inferences (using both priors) to the true value p = 0.7. + +### Hard + +**Introduction**. The practice problems here use the data below. These data indicate the gender (male = 1, female = 0) of officially reported first and second born children in 100 two-children families. + +```{r} +birth1 <- c(1,0,0,0,1,1,0,1,0,1,0,0,1,1,0,1,1,0,0,0,1,0,0,0,1,0, +0,0,0,1,1,1,0,1,0,1,1,1,0,1,0,1,1,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0, +1,1,0,1,0,0,1,0,0,0,1,0,0,1,1,1,1,0,1,0,1,1,1,1,1,0,0,1,0,1,1,0, +1,0,1,1,1,0,1,1,1,1) +birth2 <- c(0,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0, +1,1,1,0,1,1,1,0,1,0,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1, +1,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0,0,0,1,0,0,0,1,1,0,0,1,0,0,1,1, +0,0,0,1,1,1,0,0,0,0) +``` + +So, for example, the first family in the data reported a boy (1) and then a girl (0). The second family reported a girl (0) and then a boy (1). The third family reported two girls. You can load these wo vectors into R's memory by typing: + +```{r} + +#library(rethinking) +data(homeworkch3) + +``` + +3H1. Using grid approximation, compute the posterior distribution for the probability of a birth being a boy. Assume a uniform prior probability. Which parameter value maximizes the posterior probability? + +```{r} + +# Okay, so since we are talking about any birth, we should find two things - how many boys were born and how many observation there were in total. + +observations = length(birth1) + length(birth2) +boys_born = sum(birth1) + sum(birth2) + +p_grid = seq(from = 0, to = 1, length.out=1000) +prior = rep(1, 1000) +likelihood = dbinom(boys_born, size = observations, prob = p_grid) +posterior = likelihood * prior +posterior = posterior / sum(posterior) + +# Find maximum a posterior (MAP) estimate. +p_grid[which.max(posterior)] +``` + +3H2. Using the sample function, draw 10,000 random parameter values from the posterior distribution you calculated above. Use these samples to estimate the 50%, 89% and 97% highest posterior density intervals. + +```{r} + +set.seed(100) +samples = sample(p_grid, prob = posterior, size = 1e4, replace = TRUE) + +HPDI(samples, prob = 0.5) +HPDI(samples, prob = 0.89) +HPDI(samples, prob = 0.97) + +``` + +3H3. Use rbinom to simulate 10,000 replicates of 200 births. You should end up with 10,000 numbers, each one a count of boys out of 200 births. Compare the distribution of predicted numbers of boys to the actual count in the data (111 boys out of 200 births). There are many good ways to visualize the simulations, but the dens command (part of rethinking package) is probably the easiest way in this case. Does it look like the model fits the data well? That is, does the distribution of the predictions include the actual observation as a central, likely outcome? + +```{r} + +w = rbinom(n = 1e4, size = 200, prob = samples) +dens(w) + +``` diff --git a/exercises/class-5-exercises.Rmd b/exercises/class-5-exercises.Rmd new file mode 100644 index 0000000..41eecc4 --- /dev/null +++ b/exercises/class-5-exercises.Rmd @@ -0,0 +1,163 @@ +--- +title: "class-5-exercises" +output: html_document +date: "2024-04-01" +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +```{r} +pacman::p_load(dagitty, tidyverse, brms, here, tidybayes, geomtextpath, gganimate, splines, ggplot) + +library(rethinking) +library(dagitty) +``` + +Welcome back to your favorite course of all time. This week we start getting into causal inference - arguably the most important part of all of your Methods courses. Why? Because we finally start practicing *doing science*, not just statistics. + +The exercises for this week are versions of exercises from Chapter 5, modified to ease you into our next portfolio which we will start on next week. + +## Exercises + +### Easy. + +Do this by discussing the exercises in pairs. No need to code anything. :)) + +### Medium + +Throughout this course, there has been a hard emphasis on clarifying and formalizing your assumptions about the world in your statistical models. Up until now we have been doing it in a form of priors for our parameters. Today we begin formalizing our *causal assumptions* about the *generative model* of the data, in the form of DAGs. + +#### 5M0. + +Conceptual question. Let's say you have bought a farm in Vestjylland. The previous owner was a **very** detail oriented and meticulous farmer and has logged all of the relevant data, such as soil moisture and air temperature for the period of 3 years. + +You are far from a good farmer. You only know how to code. In the morning, the televised weather forecast forecasts substantial heat wave for the next 3 days. You don't know what to do, should you water your crops extra, given this new information? You don't trust your common sense, so you throw in all of the data (soil moisture and air temperature) into your machine learning model. The model predicts that soil moisture for the forecasted values of air temperature will be normal. So you decide not to give the soil extra water. Was that a good decision? Why? + +Bad decision, cause the model does not include amount of watering. The heat wave leads to extra watering, which leads to more moisture. Or maybe these are too extreme values for the model to predict realistically. + +#### 5M1. + +Invent your own example of a spurious correlation. An outcome variable should be correlated with both predictor variables. But when both predictors are entered in the same model, the correlation between the outcome and one of the predictors should mostly vanish (or at least be greatly reduced). + +See if you can come up with a cogsci-inspired phenomenon. Anxiety? Bliss? Make a DAG and use the *dagitty* package to illustrate it. + + +# In this (simulated) example, I'll predict ice cream sales from the temperature and the number of shark attacks. + +```{r} + +# ilustrate your dag +# What are the conditional independencies of your DAG? + +# temperature -> shark attacks (more people swimming) -> ice creams + +library(dagitty) +dag_new <- dagitty( "dag{ + temp -> ice; + temp -> attacks; + attacks -> ice }" ) +coordinates(dag_new) <- list( x=c(A=0,D=1,M=2) , y=c(A=0,D=1,M=0) ) +drawdag( dag_new ) + +``` + +Now see if you can play god and generate the data. This will require you to think even deeper about your variables - what is the scale of each variable and how they interact. Formulate your assumptions in natural language and perhaphs ask ChatGPT to help you with the data simulation process to avoid spending a lot of time looking for the right code. Just make sure that the output matches your desired generative structure. :)) + + +```{r} + +# generate your data here +set.seed(123) +n <- 100 +temp <- rnorm(n) +shark <- rnorm(n, temp) +ice_cream <- rnorm(n, temp) + +spur_exp <- tibble(ice_cream, temp, shark) %>% + mutate(across(everything(), standardize)) + +#d = ice, m=sharks, a=temp +m_temp <- quap( + alist( + ice_cream ~ dnorm( mu , sigma ) , + mu <- a + b_temp * temp , + a ~ dnorm( 0 , 0.2 ) , + b_temp ~ dnorm( 0 , 0.5 ) , + sigma ~ dexp( 1 ) + ) , data = spur_exp ) + +m_shark <- quap( + alist( + ice_cream ~ dnorm( mu , sigma ) , + mu <- a + b_shark * shark , + a ~ dnorm( 0 , 0.2 ) , + b_shark ~ dnorm( 0 , 0.5 ) , + sigma ~ dexp( 1 ) + ) , data = spur_exp ) + +m_both <- quap( + alist( + ice_cream ~ dnorm( mu , sigma ) , + mu <- a + b_shark*shark + b_temp*temp , + a ~ dnorm( 0 , 0.2 ) , + b_shark ~ dnorm( 0 , 0.5 ) , + b_temp ~ dnorm( 0 , 0.5 ) , + sigma ~ dexp( 1 ) + ) , data = spur_exp ) + + +precis( m_temp ) +precis( m_shark ) +precis( m_both ) + +plot( coeftab(m_temp,m_shark,m_both), par=c("b_temp","b_shark") ) + +``` + + +```{r} +ggplot(spur_exp, aes(temp, ice_cream))+ + geom_point()+ + geom_smooth(formula = y ~ x, method = "lm") + +ggplot(spur_exp, aes(temp, shark))+ + geom_point()+ + geom_smooth(formula = y ~ x, method = "lm") +``` + + + +#### 5M2. + +Invent your own example of a masked relationship. An outcome variable should be correlated with both predictor variables, but in opposite directions. And the two predictor variables should be correlated with one another. + +No need to illustrate or code anything here. Just think it up. Try another cog-sci related phenomenon. + +#### 5M4. + +In the divorce data, States with high numbers of Mormons (members of The Church of Jesus Christ of Latter-day Saints, LDS) have much lower divorce rates than the regression models expected. Find a list of LDS population by State and use those numbers as a predictor variable, predicting divorce rate using marriage rate, median age at marriage, and percent LDS population (possibly stan- dardized). You may want to consider transformations of the raw percent LDS variable. + +You don't need to find any data tho, we got you covered. + +```{r} + +data(WaffleDivorce) + +d <- WaffleDivorce + +d$pct_LDS <- c(0.75, 4.53, 6.18, 1, 2.01, 2.82, 0.43, 0.55, 0.38, + 0.75, 0.82, 5.18, 26.35, 0.44, 0.66, 0.87, 1.25, 0.77, 0.64, 0.81, + 0.72, 0.39, 0.44, 0.58, 0.72, 1.14, 4.78, 1.29, 0.61, 0.37, 3.34, + 0.41, 0.82, 1.48, 0.52, 1.2, 3.85, 0.4, 0.37, 0.83, 1.27, 0.75, + 1.21, 67.97, 0.74, 1.13, 3.99, 0.92, 0.44, 11.5 ) + +d$L <- standardize( d$pct_LDS ) +d$A <- standardize( d$MedianAgeMarriage ) +d$M <- standardize( d$Marriage ) +d$D <- standardize( d$Divorce ) + +``` + +Feel free to proceed other exercises from Chapter 5 if you have completed the exercises above. \ No newline at end of file