Multiparameter models

class: center, middle, inverse, title-slide

# Multiparameter models
### Math 315: Bayesian Statistics

---

class: middle, clear

## 1. Introduction to multiparameter models

---

.left-column[
## Example
]

.right-column[

.large[
- Interest is in analyzing the average height of an adult for the Dobe area !Kung San, a foraging population
]

]

---
class: middle, inverse

# Noninformative analysis

---

# Sampling from the joint posterior

.large[
**Target:** `$p(\mu |y_1, \ldots, y_m)$`

**1.** Sample `$\sigma^2$` from `$\sigma^2 | y_1, \ldots, y_m \sim \text{InverseGamma}\left( \frac{n-1}{2},\ \frac{(n-1)s^2}{2} \right)$`

```r
n <- nrow(adults)
s2 <- var(adults$height)
sigma2s <- 1 / rgamma(1e4, (n-1)/2, (n-1) * s2 / 2) 
```

**2.** Sample `$\mu$` from `$\mu | \sigma^2, y_1, \ldots, y_m \sim \mathcal{N}\left( \bar{y},\ \frac{\sigma^2}{n} \right)$`

```r
ybar <- mean(adults$height)
mus <- rnorm(1e4, ybar, sqrt(sigma2s / n))
```
]

---
class: clear

```r
post_samples <- data.frame(mu = mus, sigma2 = sigma2s)
ggplot(data = post_samples) +
  geom_point(mapping = aes(x = mu, y = sigma2), alpha = 0.1) +
  theme_minimal() +
  labs(x = expression(mu), y = expression(sigma^2))
```

---
class: clear

```r
ggplot(data = post_samples) +
  geom_point(mapping = aes(x = mu, y = sigma2), alpha = 0.1) +
  geom_density_2d(mapping = aes(x = mu, y = sigma2), color = "orange", size = 0.8) +
  theme_minimal() +
  labs(x = expression(mu), y = expression(sigma^2))
```

---
class: clear

```r
ggplot(data = post_samples) +
  geom_density_2d(mapping = aes(x = mu, y = sigma2), color = "orange") +
  theme_minimal() +
  labs(x = expression(mu), y = expression(sigma^2))
```

---
class: clear

---

## Marginal posteriors

.pull-left[

`$\boldsymbol \mu$`

![](10-intro-multiparameter_files/figure-html/unnamed-chunk-8-1.svg)

```r
quantile(mus, probs = c(0.05, 0.9))
```

```
##       5%      90% 
## 153.9200 155.1236
```

]

.pull-right[

`$\boldsymbol \sigma$`

![](10-intro-multiparameter_files/figure-html/unnamed-chunk-10-1.svg)

```r
quantile(sigma2s, probs = c(0.05, 0.9))
```

```
##       5%      90% 
## 53.09276 66.33424
```

```r
quantile(sqrt(sigma2s), probs = c(0.05, 0.9))
```

```
##       5%      90% 
## 7.286478 8.144584
```
]

---
class: middle, inverse

# Informative analysis

---

# McElreath's analysis

.pull-left[

### Model

`$y_i \sim \mathcal{N}(\mu, \sigma)$`

`$\mu \sim \mathcal{N}(178,\ 20)$`

`$\sigma \sim \text{Unif}(0, 50)$`

]

.pull-right[

![](10-intro-multiparameter_files/figure-html/unnamed-chunk-13-1.svg)

.code70[

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    14.8   157.5   178.4   178.6   199.6   352.0
```
]

]

---

# Joint posterior via grid approximation

```r
# Create grid over the coordinate plane
param_grid <- expand.grid(
  mu = seq(from = 118, to = 238, length.out = 1000), 
  sigma = seq(from = 0, to = 50, length.out = 1000)
)

# Calculate joint log prior for each point on grid
logprior <- dnorm(param_grid[, "mu"], 178, 20, log = TRUE) +
      dunif(param_grid[, "sigma"], 0, 50, log = TRUE)

# Calculate log likelihood for each point on grid
ll <- numeric(length = nrow(param_grid))
for(i in 1:nrow(param_grid)) {
  ll[i] <- sum(dnorm(adults$height, mean = param_grid[i, "mu"], 
                     sd = param_grid[i, "sigma"], log = TRUE))
}

# Calculate log posterior, then exponentiate
logposterior <- logprior + ll
unstd_posterior <- exp(logposterior - max(logposterior)) # numeric stability
posterior <- unstd_posterior / sum(unstd_posterior)
```

---

# Joint posterior via grid approximation

```r
posterior_draws <- dplyr::sample_n(
  param_grid, size = 1e4, replace = TRUE, weight = posterior
)
```

---

# Marginal posteriors via grid approximation

.pull-left[

### `$\mu$`

![](10-intro-multiparameter_files/figure-html/unnamed-chunk-18-1.svg)

```r
quantile(posterior_draws$mu, probs = c(0.05, 0.95))
```

```
##       5%      95% 
## 153.9159 155.2372
```

]

.pull-right[

### `$\sigma$`

![](10-intro-multiparameter_files/figure-html/unnamed-chunk-20-1.svg)

```r
quantile(posterior_draws$sigma, probs = c(0.05, 0.95))
```

```
##       5%      95% 
## 7.307307 8.258258
```

<br>

]

---

# Vectorization in R

The grid approximation on a 1000 `$\times$` 1000 grid took about 30 seconds

Vectorizing the log-likelihood calculation improves speed to 14 seconds

```r
# vectorized log likelihood function
llnorm <- function(x, mu, sigma) { 
  sum(dnorm(x, mean = mu, sd = sigma, log = TRUE))
}

# use all of the x vector for each mu-sigma combo
llnorm <- Vectorize(llnorm, vectorize.args = c("mu", "sigma"))

# Calculate log likelihood for each point on grid
ll <- llnorm(x = adults$height, mu = param_grid$mu, sigma = param_grid$sigma)

# Everything else is the same...
```