diff --git a/lab/lab_week_05.Rmd b/lab/lab_week_05.Rmd new file mode 100644 index 0000000..99be6b2 --- /dev/null +++ b/lab/lab_week_05.Rmd @@ -0,0 +1,71 @@ +--- +title: "EEEB UN3005/GR5005 \nLab - Week 05 - 24 and 26 February 2020" +author: "USE YOUR NAME HERE" +output: pdf_document +fontsize: 12pt +--- + +```{r setup, include = FALSE} + +knitr::opts_chunk$set(echo = TRUE) + +library(rethinking) +``` + + +# Statistical Distributions and Summary Statistics + + +## Exercise 1: Grid Approximation, Our Old Friend + +Imagine that the globe tossing example from the *Statistical Rethinking* text and class resulted in 8 water observations out of 15 globe tosses. + +With this set of data, use grid approximation (with 101 grid points) to construct the posterior for *p* (the probability of water). Assume a flat prior. + +Plot the posterior distribution. + +```{r} + +``` + + +## Exercise 2: Sampling From a Grid-Approximate Posterior + +Now generate 10,000 samples from the posterior distribution of *p*. Call these samples `post.samples`. Visualize `post.samples` using the `dens()` function. + +For your own understanding, re-run your sampling and plotting code multiple times to observe the effects of sampling variation. + +```{r} + +``` + + +## Exercise 3: Summarizing Samples + +Return the mean, median, and mode (using `chainmode()`) of `post.samples`. Then calculate the 80%, 90%, and 99% highest posterior density intervals of `post.samples`. + +```{r} + +``` + + +## Exercise 4: Model Predictions + +Using `post.samples`, generate 10,000 simulated model predictions (you can call these `preds`) for a binomial trial of size 15. Visualize the model predictions using the `simplehist()` function. + +Based on these posterior predictions, what is the probability of observing 8 waters in 15 globe tosses? + +```{r} + +``` + + +## Exercise 5: More Model Predictions + +Using the *same* posterior samples (i.e., `post.samples`), generate 10,000 posterior predictions for a binomial trial of size 9. + +Using these new predictions, calculate the posterior probability of observing 8 waters in 9 tosses. + +```{r} + +``` diff --git a/lecture/lecture_week_04.R b/lecture/lecture_week_04.R index afb411c..e37282b 100644 --- a/lecture/lecture_week_04.R +++ b/lecture/lecture_week_04.R @@ -1,24 +1,27 @@ -# Demonstrate grid approximation of the posterior for a globe tossing problem -# with 6 water observations out of 9 globe tosses +# Demonstrate grid approximation of the posterior for a +# globe tossing problem with 6 water observations out of +# 9 globe tosses -#============================================================================== +#========================================================== # 1) Define the grid to be used to compute the posterior # In this case, let's say we choose 20 grid points -# Since the parameter of interest (proportion of water on the globe) is -# bounded by 0 and 1, our grid should have those limits as well +# Since the parameter of interest (proportion of water on +# the globe) is bounded by 0 and 1, our grid should have +# those limits as well p_grid <- seq(from = 0, to = 1, length.out = 20) p_grid -#============================================================================== +#========================================================== -# 2) Compute/define the value of the prior at each parameter value on the grid +# 2) Compute/define the value of the prior at each +# parameter value on the grid # In this case, simply choose a flat prior prior <- rep(1, 20) @@ -26,42 +29,46 @@ prior <- rep(1, 20) # You could visualize this prior plot(p_grid, prior, pch = 19, type = "b") -#============================================================================== +#========================================================== -# 3) Compute the likelihood at each parameter value on the grid +# 3) Compute the likelihood at each parameter value on +# the grid -# This requires the use of a likelihood function applied to the observed data, -# evaluated at each potential parameter value on the grid +# This requires the use of a likelihood function applied +# to the observed data, evaluated at each potential +# parameter value on the grid likelihood <- dbinom(6, size = 9, prob = p_grid) # You could also visualize the likelihood plot(p_grid, likelihood, pch = 19, type = "b") -#============================================================================== +#========================================================== -# 4) Compute the unstandardized posterior at each parameter value on the grid +# 4) Compute the unstandardized posterior at each +# parameter value on the grid -# The unstandardized posterior is simply the product of the likelihood and -# prior +# The unstandardized posterior is simply the product of +# the likelihood and prior unstd.posterior <- likelihood * prior # Again, could visualize plot(p_grid, unstd.posterior, pch = 19, type = "b") -# Note that this unstandardized posterior is not a proper probability -# distribution since it does not add to 1 +# Note that this unstandardized posterior is not a proper +# probability distribution since it does not add to 1 sum(unstd.posterior) -#============================================================================== +#========================================================== -# 4) Standardize the posterior +# 5) Standardize the posterior posterior <- unstd.posterior/sum(unstd.posterior) -# This standardized posterior is a now proper probability distribution +# This standardized posterior is a now proper probability +# distribution sum(posterior) # Visualize the posterior diff --git a/lecture/lecture_week_05.R b/lecture/lecture_week_05.R new file mode 100644 index 0000000..5d909b6 --- /dev/null +++ b/lecture/lecture_week_05.R @@ -0,0 +1,292 @@ + + +library(rethinking) + +#========================================================== + + +# Probability mass function for the binomial distribution + +# Parameters: +# +# 1) number of trials = 10 +# (Note: This is a parameter that we normally will NOT +# need to estimate in our statistical models. Rather, it's +# something we observe. But it is a parameter in the sense +# that it's a numeric characteristic used to define +# the binomial distribution.) +# +# 2) probability of success = 0.3 + +# If we have 10 binomial trials, then only 11 different +# outcomes are possible. We could have anywhere from 0 +# to 10 successes + +# Let's get the probability masses associated with each +# of these 11 outcomes, given our arbitrarily chosen +# parameter values +prob.masses <- dbinom(0:10, size = 10, prob = 0.3) + +prob.masses + +# We can plot these probability masses +plot( + x = 0:10, y = prob.masses, + xlab = "Number of successful binomial trials", + ylab = "Probability", + pch = 19 +) + +# Note that the probability mass of all outcomes +# sums to 1 +sum(prob.masses) +# In other words, if we conduct 10 binomial trials, one +# of these 11 outcomes MUST be observed. But some outcomes +# are more likely than others + + +# Random generation function for the binomial distribution + +# 10 simulated draws (observations) from the binomial +# distribution +# Note that each of these observations is composed of +# 10 binomial trials. The "size" and "prob" parameters +# have not changed from the examples above +draws.10 <- rbinom(10, size = 10, prob = 0.3) + +simplehist( + draws.10, xlim = c(0, 10), + xlab = "Number of successes observed" +) + +# 100 simulated draws (observations) +draws.100 <- rbinom(100, size = 10, prob = 0.3) + +simplehist( + draws.100, xlim = c(0, 10), + xlab = "Number of successes observed" +) + +# 10,000 simulated draws (observations) +draws.10000 <- rbinom(10000, size = 10, prob = 0.3) + +simplehist( + draws.10000, xlim = c(0, 10), + xlab = "Number of successes observed" +) + +#========================================================== + + +# Probability density function for the normal distribution + +# Parameters: +# +# 1) mean = 2.5 +# +# 2) standard deviation = 1.5 + +# Plot the probability density function for the +# normal distribution +curve( + dnorm(x, mean = 2.5, sd = 1.5), + from = -10, to = 10, + ylab = "Probability density" +) + +# 10 simulated draws (observations) from the +# normal distribution +rnorm(10, mean = 2.5, sd = 1.5) + + +# Probability density function for the Cauchy distribution + +# Parameters: +# +# 1) location = 2.5 +# +# 2) scale = 1.5 + +# Plot the probability density function for the +# Cauchy distribution +curve( + dcauchy(x, location = 2.5, scale = 1.5), + from = -10, to = 10, + ylab = "Probability density" +) + +# 10 simulated draws (observations) from the +# Cauchy distribution +rcauchy(10, location = 2.5, scale = 1.5) + +#========================================================== + + +# Sampling from a posterior distribution + +# First, let's generate a posterior distribution using +# grid approximation (with 101 points). The data observed +# are 6 successes out of 9 binomial trials +p_grid <- seq(from = 0, to = 1, length.out = 101) +prior <- rep(1, 101) +likelihood <- dbinom(6, size = 9, prob = p_grid) +unstd.posterior <- likelihood * prior +posterior <- unstd.posterior / sum(unstd.posterior) + +# Plot the posterior +plot( + p_grid, posterior, + pch = 19, type = "b", + xlab = "Proportion of water on globe", + ylab = "Posterior probability" +) + +# Generate 10,000 samples of p (the proportion of water +# parameter) from the grid-approximate posterior +samples <- sample( + p_grid, size = 10000, + prob = posterior, replace = TRUE +) +# Reading the arguments within "samples()" would go +# something like this: I want 10,000 samples (with +# replacement) from the values in "p_grid", taking values +# according to the probability assigned to them by the +# "posterior" + +# Plot the posterior samples +plot(samples, pch = 19, col = alpha("black", 0.2)) +dens(samples, xlim = c(0, 1), + xlab = "Proportion of water on globe") + + +# Grid-approximate posterior for 5 successes in 5 +# binomial trials (leads to an extremely skewed +# posterior distribution for p) +p_grid2 <- seq(from = 0, to = 1, length.out = 101) +prior2 <- rep(1, 101) +likelihood2 <- dbinom(5, size = 5, prob = p_grid2) +unstd.posterior2 <- likelihood2 * prior2 +posterior2 <- unstd.posterior2 / sum(unstd.posterior2) + +# Plot the posterior +plot( + p_grid2, posterior2, + pch = 19, type = "b", + xlab = "Proportion of water on globe", + ylab = "Posterior probability" +) + +# Generate 10,000 samples of p (the proportion of water +# parameter) from the grid-approximate posterior +samples2 <- sample( + p_grid2, size = 10000, + prob = posterior2, replace = TRUE +) + +# Plot the posterior samples +plot(samples2, pch = 19, col = alpha("black", 0.2)) +dens(samples2, xlim = c(0, 1), + xlab = "Proportion of water on globe") + +#========================================================== + + +# Calculate percentile and highest posterior density +# intervals from our posterior samples + +# 50% PI +PI(samples2, prob = 0.5) + +dens(samples2, xlim = c(0, 1), xaxs = "i", + xlab = "Proportion of water on globe", + main = "50% PI shaded") +shade( + density(samples2, adjust = 0.5), + PI(samples2, prob = 0.5), + col = "darkgrey" +) + +# 50% HPDI +HPDI(samples2, prob = 0.5) + +dens(samples2, xlim = c(0, 1), xaxs = "i", + xlab = "Proportion of water on globe", + main = "50% HPDI shaded") +shade( + density(samples2, adjust = 0.5), + HPDI(samples2, prob = 0.5), + col = "darkgrey" +) + + +# 70% PI +PI(samples2, prob = 0.7) + +dens(samples2, xlim = c(0, 1), xaxs = "i", + xlab = "Proportion of water on globe", + main = "70% PI shaded") +shade( + density(samples2, adjust = 0.5), + PI(samples2, prob = 0.7), + col = "darkgrey" +) + +# 70% HPDI +HPDI(samples2, prob = 0.7) + +dens(samples2, xlim = c(0, 1), xaxs = "i", + xlab = "Proportion of water on globe", + main = "70% HPDI shaded") +shade( + density(samples2, adjust = 0.5), + HPDI(samples2, prob = 0.7), + col = "darkgrey" +) + + +# With a less skewed posterior, the PI and HPDI are +# often going to be very, very similar + +# 70% PI +PI(samples, prob = 0.7) + +dens(samples, xlim = c(0, 1), + xlab = "Proportion of water on globe", + main = "70% PI shaded") +shade( + density(samples, adjust = 0.5), + PI(samples, prob = 0.7), + col = "darkgrey" +) + +# 70% HPDI +HPDI(samples, prob = 0.7) + +dens(samples, xlim = c(0, 1), + xlab = "Proportion of water on globe", + main = "70% HPDI shaded") +shade( + density(samples, adjust = 0.5), + HPDI(samples, prob = 0.7), + col = "darkgrey" +) + +#========================================================== + + +# Generate a posterior predictive distribution for +# binomial trials of size 10, given the posterior in +# samples2 +preds <- rbinom(10000, size = 10, prob = samples2) + +plot(preds, pch = 19, col = alpha("black", 0.2)) +simplehist(preds, xlim = c(0, 10)) + + +# Note that even if we fix the probability of success value +# there would still be variation in our outcomes because +# the data-generating process is subject to outcome +# uncertainty +test.outcomes <- rbinom(10000, size = 10, prob = 0.8) + +simplehist(test.outcomes, xlim = c(0, 10)) diff --git a/lecture/lecture_week_05.pptx b/lecture/lecture_week_05.pptx new file mode 100644 index 0000000..3be3d23 Binary files /dev/null and b/lecture/lecture_week_05.pptx differ