```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(mosaic)
require(rethinking)

knitr::opts_template$set(
  original = list(eval = params$original, include = params$original),
  updated = list(eval = params$updated, include = params$updated))
```


```{r, opts.label = "original", results = "asis", echo = FALSE}
cat("The original code from *Statistical Rethinking* is shown below.")
cat("\n\n")
```

```{r, opts.label = "updated", results = "asis", echo = FALSE}
cat(paste(
"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",
sep = "\n"
))
cat("\n\n")
```

## R code 2.1

```{r, chunk2.1}
ways <- c(0, 3, 8, 9, 0)
ways / sum(ways)
```

## R code 2.2

```{r, chunk2.1}
```{r, chunk2.2}
dbinom(6, size = 9, prob = 0.5)
```

## R code 2.3

```{r, chunk2.3, opts.label = "original"}
# 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)

```


```{r, chunk2.3a, opts.label = "updated"}
# 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

```{r, chunk2.4, opts.label = "original"}
plot(p_grid,
     posterior,
     type = "b",
     xlab = "probability of water",
     ylab = "posterior probability")
mtext("20 points")
```




```{r, chunk2.4a, opts.label = "updated"}
xyplot(posterior ~ p, data = Water9,
     type = "b",
     xlab = "probability of water",
     ylab = "posterior probability",
     main = "20 points")
```

```{r, chunk2.4b, opts.label = "updated"}
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

  
```{r, chunk2.5, opts.label = "original"}
prior <- ifelse(p_grid < 0.5, 0, 1)
prior <- exp(-5 * abs(p_grid - 0.5))
```


```{r, chunk2.5a, opts.label = "updated"}
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

```{r, chunk2.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)
```

## R code 2.7

```{r, chunk2.7, opts.label = "original"}
# analytical calculation
w <- 6
n <- 9
curve(dbeta(x, w + 1, n - w + 1), from = 0, to = 1)
# quadratic approximation
curve(dnorm(x, 0.67, 0.16), lty = 2, add = TRUE)
```

```{r, chunk2.7a, opts.label = "updated"}
# 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)
```