Statistical Rethinking (Code)

Code from Statistical Rethinking modified by R Pruim is shown below. Differences to the oringal include:

• a preference for putting data into containers (data frames, mostly), rather than working with lose vectors.
• use of lattice and ggplot2 rather than base graphics
• use of tidyverse for data transformation
• better (in my opinion) naming conventions

R code 2.1

ways <- c(0, 3, 8, 9, 0)
ways / sum(ways)

## [1] 0.00 0.15 0.40 0.45 0.00

R code 2.2

ways <- c(0, 3, 8, 9, 0)
ways / sum(ways)

## [1] 0.00 0.15 0.40 0.45 0.00

dbinom(6, size = 9, prob = 0.5)

## [1] 0.1640625

R code 2.3

# define grid
Water9 <-
  data_frame(p =  seq(from = 0, to = 1, length.out = 20)) %>% 
  mutate(
    unstd.prior = 1,                        # recycles automatically
    prior = unstd.prior / sum(unstd.prior), # standardize the prior
    likelihood =  dbinom(6, size = 9, prob = p),
    unstd.posterior = likelihood * prior,
    posterior = unstd.posterior / sum(unstd.posterior),
    posterior.dens = posterior / 0.001
  )  

R code 2.4

xyplot(posterior ~ p, data = Water9,
     type = "b",
     xlab = "probability of water",
     ylab = "posterior probability",
     main = "20 points")

xyplot(prior + likelihood + posterior ~ p, data = Water9, type = "b", 
     auto.key = list(points = TRUE, lines = TRUE),
     xlab = "probability of water",
     ylab = "posterior probability",
     main = "20-point grid")

ggplot(Water9, aes(x = p)) +
  geom_line(aes(y = prior, color = "prior")) +
  geom_line(aes(y = likelihood, color = "likelihood")) +
  geom_line(aes(y = posterior, color = "posterior")) +
  geom_point(aes(y = prior, color = "prior")) +
  geom_point(aes(y = likelihood, color = "likelihood")) +
  geom_point(aes(y = posterior, color = "posterior")) +
  labs(title = "20-point grid") 

R code 2.5

Water9a <-
  data_frame(p =  seq(from = 0, to = 1, length.out = 20)) %>% 
  mutate(
    unstd.prior = ifelse(p < 0.5, 0, 1),
    prior = unstd.prior / sum(unstd.prior), # standardize the prior
    likelihood =  dbinom(6, size = 9, prob = p),
    unstd.posterior = likelihood * prior,
    posterior = unstd.posterior / sum(unstd.posterior),
    posterior.dens = posterior / 0.001
  )  

Water9b <-
  data_frame(p =  seq(from = 0, to = 1, length.out = 20)) %>% 
  mutate(
    unstd.prior =  exp(-5 * abs(p - 0.5)),
    prior = unstd.prior / sum(unstd.prior), # standardize the prior
    likelihood =  dbinom(6, size = 9, prob = p),
    unstd.posterior = likelihood * prior,
    posterior = unstd.posterior / sum(unstd.posterior),
    posterior.dens = posterior / 0.001
  )  

R code 2.6

library(rethinking)
globe.qa <- 
  map(alist(w ~ dbinom(9, p),  # binomial likelihood
            p ~ dunif(0, 1)),  # uniform prior
      data = list(w = 6))

# display summary of quadratic approximation
precis(globe.qa)

##   Mean StdDev 5.5% 94.5%
## p 0.67   0.16 0.42  0.92

R code 2.7

# analytical calculation
w <- 6
n <- 9
plotDist("beta", params = list(shape1 = w + 1, shape2 = n - w + 1))

# quadratic approximation
plotDist("norm", mean = 0.67, sd = 0.16, add = TRUE, col = 2, lty = 2)