.

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)
+
+```

exercises/class-5-exercises.Rmd

@@ -0,0 +1,88 @@
+---
+title: "class-5-exercises"
+output: html_document
+date: "2024-04-01"
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+```{r}
+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?
+
+#### 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 *daggity* package to illustrate it. 
+
+```{r}
+
+# ilustrate your dag
+# What are the conditional independencies of your DAG?
+
+```
+
+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
+
+```
+
+#### 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.