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..b631070 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/exercises/class-8-exercises.Rmd b/exercises/class-8-exercises.Rmd new file mode 100644 index 0000000..a5b0929 --- /dev/null +++ b/exercises/class-8-exercises.Rmd @@ -0,0 +1,55 @@ +--- +title: "Class 8 Exercises" +output: html_document +date: "2024-04-24" +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +```{r} +library(rethinking) +``` + +## Exercises + +### Easy + +Do all the easy exercises from Chapter 9: **9E1** through **9E6**. + +### Medium + +#### 9M1 + +```{r} + +# Your solution here + +``` + +#### 9M2 + +```{r} + +# Your solution herre + +``` + +### Hard + +#### 9H1 + +```{r} + +# Your solution here + +``` + +#### 9H2 + +```{r} + +# Your solution herre + +``` diff --git a/exercises/class-8-solutions.Rmd b/exercises/class-8-solutions.Rmd new file mode 100644 index 0000000..f19b89f --- /dev/null +++ b/exercises/class-8-solutions.Rmd @@ -0,0 +1,237 @@ +--- +title: "Class 8 Solutions" +output: html_document +date: "2024-04-24" +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +```{r} +library(rethinking) +``` + +## Exercises + +### Easy + +#### 9E1 + +Only (3) is required. + +#### 9E2 + +Gibbs sampling requires that we use special priors that are *conjugate* with the likelihood. This means that holding all the other parameters constant, it is possible to derive analytical solutions for the posterior distribution of each parameter. These conditional distributions are used to make smart proposals for jumps in the Markov chain. Gibbs sampling is limited both by the necessity to use conjugate priors, as well as its tendency to get stuck in small regions of the posterior when the posterior distribution has either highly correlated parameters or high dimension. + +#### 9E3 + +Hamiltonian Monte Carlo cannot handle discrete parameters. This is because it requires a smooth surface to glide its imaginary particle over while sampling from the posterior distribution. + +#### 9E4 + +The effect number of samples n_eff is an estimate of the number of completely independent samples that would hold equivalent information about the posterior distribution. It is always smaller than the actual number of samples, because samples from a Markov chain tend to sequentially correlated or *autocorrelated*. As autocorrelation rises, `n_eff` gets smaller. At the limit of perfect autocorrelation, for example, all samples would have the same value and `n_eff` would be equal to 1, no matter the actual number of samples drawn. + +#### 9E5 + +`Rhat` should approach 1. How close should it get? People disagree, but it is common to judge that any value less than 1.1 indicates convergence. But like all heuristic indicators, `Rhat` can be fooled. + + +#### 9E6 + +A healthy Markov chain should be both *stationary* and *well-mixing*. The first is necessary for inference. The second is desirable because it means the chain is more efficient. A chain that is both of these things should resemble horizontal noise. + +A chain that is malfunctioning, as the problem asks, would not be stationary. This means it is not converging to the target distribution, the posterior distribution. Examples were provided in the chapter. A virtue of Hamiltonian Monte Carlo is that it makes such chains very obvious: they tend to be rather flat wandering trends. Sometimes they are perfectly flat. + +The best test of convergence is always to compare multiple chains. So the best sketch of a malfunctioning trace plot would be one that shows multiple chains wandering into different regions of the parameter space. Figure 8.7 in the chapter, left side, provides an example. + +### Medium + +#### 9M1 + +First, we load and preprocess the data. + +```{r} +data(rugged) +d <- rugged +d$log_gdp <- log(d$rgdppc_2000) +dd <- d[ complete.cases(d$rgdppc_2000) , ] +dd$log_gdp_std <- dd$log_gdp / mean(dd$log_gdp) +dd$rugged_std <- dd$rugged / max(dd$rugged) +dd$cid <- ifelse( dd$cont_africa==1 , 1 , 2 ) +dat_slim <- list( + log_gdp_std = dd$log_gdp_std, + rugged_std = dd$rugged_std, + cid = as.integer( dd$cid ) ) +``` + +Then we take Model 9.1 and place a uniform prior on $\sigma$. + +```{r} +# new model with uniform prior on sigma +m9.1_unif <- ulam( + alist( + log_gdp_std ~ dnorm( mu , sigma ) , + mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) , + a[cid] ~ dnorm( 1 , 0.1 ) , + b[cid] ~ dnorm( 0 , 0.3 ) , + sigma ~ dunif( 0 , 1 ) + ) , data=dat_slim , chains=4 , cores=4 ) + + +``` + +We then also make a version with an exponential prior on $\sigma$. + +```{r} +m9.1_exp <- ulam( + alist( + log_gdp_std ~ dnorm( mu , sigma ) , + mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) , + a[cid] ~ dnorm( 1 , 0.1 ) , + b[cid] ~ dnorm( 0 , 0.3 ) , + sigma ~ dexp( 1 ) + ) , data=dat_slim , chains=4 , cores=4 ) +``` + +Before looking at the posterior distributions for each model, it helps to visualize each prior. + +```{r} +curve( dexp(x,1) , from=0 , to=7 , + xlab="sigma" , ylab="Density" , ylim=c(0,1) ) +curve( dunif(x,0,1) , add=TRUE , col="red" ) +mtext( "priors" ) +``` + +Now let’s compare the posterior distributions of $\sigma$ for both models. + +```{r} +post <- extract.samples( m9.1_exp ) +dens( post$sigma , xlab="sigma" ) +post <- extract.samples( m9.1_unif ) +dens( post$sigma , add=TRUE , col="red" ) +mtext( "posterior" ) +``` + +The posterior distributions are almost identical. Why? Because there is a lot of data to inform $\sigma$. Can you find a prior that won’t wash out? + + +#### 9M2 + +Model code using using same data as in the previous problem: + +```{r} +m9M2 <- ulam( + alist( + log_gdp_std ~ dnorm( mu , sigma ) , + mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) , + a[cid] ~ dnorm( 1 , 0.1 ) , + b[cid] ~ dexp( 0.3 ) , + sigma ~ dexp( 1 ) + ) , data=dat_slim , chains=4 , cores=4 ) +precis( m9M2 , 2 ) +``` + +The big difference here is `b[2]`. In the original model, the mass of this parameter is almost entirely below zero. Now it cannot go below zero, because the prior is not defined below zero. So instead the mass presses up against zero tightly. + +### Hard + +#### 9H1 + +The code provided in the PDF version of the book doesn't work out of the box. Here is modified code that works: + +```{r} +mp <- ulam( + alist( + a ~ dnorm(0,1), + b ~ dcauchy(0,1) + ), + data=list(y=1), + iter=1e4, warmup=100) +``` + +What this code does is sample from the priors.There is no likelihood, and that’s okay. The posterior distribution is then just a merger of the priors. What is tricky about this problem though is that the Cauchy prior for the parameter `b` will not produce the kind of trace plot you might expect from a good Markov chain. This is because Cauchy is a very long tailed distribution, so it’ll occasionally make distance leaps out into the tail. We’ll look at the precis output, then the trace plot, so you can see what we mean. + + + +```{r} +precis(mp) +``` + +```{r} +traceplot( mp, n_col=2 , lwd=2 ) +``` + +The trace plot might look a little weird to you, because the trace for `b` has some big spikes in it. That’s how a Cauchy behaves, though. It has thick tails, so needs to occasionally sample way out. The trace plot for `a` is typical Gaussian in shape. Since the posterior distribution does often tend towards Gaussian for many parameters, it’s possible to get too used to expecting every trace to look like the one on the left. But you have to think about the influence of priors in this case. The trace on the right is just fine. + +#### 9H2 + +First,load and prepare the data. + +```{r} +data(WaffleDivorce) +d <- WaffleDivorce +d$D <- standardize( d$Divorce ) +d$M <- standardize( d$Marriage ) +d$A <- standardize( d$MedianAgeMarriage ) +d_trim <- list(D=d$D,M=d$M,A=d$A) +``` + +Now to fit the models over again, this time using `ulam`. Note that we need to add `log_lik=TRUE` to get the terms needed to compute PSIS or WAIC. + +```{r} +m5.1_stan <- ulam( + alist( + D ~ dnorm( mu , sigma ) , + mu <- a + bA * A , + a ~ dnorm( 0 , 0.2 ) , + bA ~ dnorm( 0 , 0.5 ) , + sigma ~ dexp( 1 ) + ) , data=d_trim , chains=4 , cores=4 , log_lik=TRUE ) +``` + + +```{r} +m5.2_stan <- ulam( +alist( + D ~ dnorm( mu , sigma ) , + mu <- a + bM * M , + a ~ dnorm( 0 , 0.2 ) , + bM ~ dnorm( 0 , 0.5 ) , + sigma ~ dexp( 1 ) + ) , data=d_trim , chains=4 , cores=4 , log_lik=TRUE ) +``` + + +```{r} +m5.3_stan <- ulam( + alist( + D ~ dnorm( mu , sigma ) , + mu <- a + bM*M + bA*A , + a ~ dnorm( 0 , 0.2 ) , + bM ~ dnorm( 0 , 0.5 ) , + bA ~ dnorm( 0 , 0.5 ) , + sigma ~ dexp( 1 ) + ) , data=d_trim , chains=4 , cores=4 , log_lik=TRUE ) +``` + +Now to compare the models: + +```{r} +compare( m5.1_stan , m5.2_stan , m5.3_stan , func=PSIS ) +``` + +```{r} +compare( m5.1_stan , m5.2_stan , m5.3_stan , func=WAIC ) +``` + +The model with only age-at-marriage comes out on top, although the model with both predictors does nearly as well. In fact, the PSIS/WAIC of both models is nearly identical. I’d call this is a tie, because even though one model does a bit better than the other, the difference between them is of no consequence. How can we explain this? Well, look at the marginal posterior for `m5.3_stan`: + +```{r} +precis(m5.3_stan) +``` + +While this model includes marriage rate as a predictor, it estimates very little expected influence for it, as well as substantial uncertainty about the direction of any influence it might have. So models `m5.3_stan` and `m5.1_stan` make practically the same predictions. After accounting for the larger penalty for `m5.3_stan` — 4.7 instead of 3.7 — the two models rank almost the same. This makes sense, because you already learned back in Chapter 5 that marriage rate probably gets its correlation with divorce rate through a correlation with age at marriage. So even though including marriage rate in a model doesn’t really aid in prediction, there is enough evidence here that the parameter bR can be estimated well enough, and including marriage rate doesn’t hurt prediction either. + +Or at least that’s what PSIS/WAIC expects. Only the future will tell which model is actually better for forecasting. PSIS/WAIC is not an oracle. It’s a golem. +