Grid approximation example

Load packages

library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)

Calculations

Grid approximation is easy in this example, or indeed in any case with a single parameter in a bounded range. Simply calculate posteriors along a grid from 0 to 1, for a sequence of 6 W in 9 tosses.

• One key to this is that dbinom() is vectorized, so we can calculate likelihood for a grid of prob values all at once.
• The prior is U[0, 1], constant over \(\Pr(W)\), so prior=1 suffices. For other priors, calculate values over the same grid.

This borrows the tidyverse code ideas from Kurz, \S 2.4.3, https://bookdown.org/content/4857/small-worlds-and-large-worlds.html#grid-approximation.

set.seed(1224)
n_points <- 10
(
  d <-
  tibble(p_grid = seq(from=0, to=1, length.out=n_points),      # define grid
         prior  = 1) |>                                        # define prior
  mutate(likelihood = dbinom(6, size = 9, prob = p_grid)) |>   # compute likelihood at each value
  mutate(unstd_posterior = likelihood * prior) |>              # product of likelihood and prior
  mutate(posterior = unstd_posterior / sum(unstd_posterior))   # standardize posterior, so it sums to 1
)

## # A tibble: 10 x 5
##    p_grid prior likelihood unstd_posterior posterior
##     <dbl> <dbl>      <dbl>           <dbl>     <dbl>
##  1  0         1   0               0         0       
##  2  0.111     1   0.000111        0.000111  0.000123
##  3  0.222     1   0.00476         0.00476   0.00528 
##  4  0.333     1   0.0341          0.0341    0.0379  
##  5  0.444     1   0.111           0.111     0.123   
##  6  0.556     1   0.217           0.217     0.241   
##  7  0.667     1   0.273           0.273     0.303   
##  8  0.778     1   0.204           0.204     0.227   
##  9  0.889     1   0.0568          0.0568    0.0631  
## 10  1         1   0               0         0

Find the maximum posterior (~ MLE)

From the table, you can see that the posterior is a maximum at \(p = 0.667\) (unsurprising, with 6 / 9 waters). In the plot below, I add a dashed line at the p_grid value corresponding to the maximum posterior. Finding this in code:

(maxrow <- which.max(d$posterior))

## [1] 7

d[maxrow, ]

## # A tibble: 1 x 5
##   p_grid prior likelihood unstd_posterior posterior
##    <dbl> <dbl>      <dbl>           <dbl>     <dbl>
## 1  0.667     1      0.273           0.273     0.303

ggplot it

p1 <-
  d |> 
  ggplot(aes(x = p_grid, y = posterior)) +
  geom_point(size=3) +
  geom_line() +
  geom_vline(aes(xintercept = p_grid[which.max(posterior)]), 
             linetype="dashed") +
  labs(subtitle = "20 points",
       x = "probability of water",
       y = "posterior probability") +
  theme_bw(base_size = 16) +
  theme(panel.grid = element_blank())
p1

Generalize this to various number of points

Make a dataset, stacking the observations for each number of points. I’ll use n_points = c(5, 10, 15, 20) here.

The rest uses some fancy tidyr and purrr tools, perhaps worth explaining step by step. In the first step, create a tibble with one row for each number of points and a list of that number of p_grid values.

tibble(n_points = c(5, 10, 15, 20)) |> 
  mutate(p_grid = map(n_points, ~seq(from = 0, to = 1, length.out = .)))

## # A tibble: 4 x 2
##   n_points p_grid    
##      <dbl> <list>    
## 1        5 <dbl [5]> 
## 2       10 <dbl [10]>
## 3       15 <dbl [15]>
## 4       20 <dbl [20]>

Then unnest the lists and add the prior. This gives a tibble with sum(c(5, 10, 15, 20)) = 50 rows.

tibble(n_points = c(5, 10, 15, 20)) |> 
  mutate(p_grid = map(n_points, ~seq(from = 0, to = 1, length.out = .))) |> 
  unnest(p_grid) |> 
  expand(nesting(n_points, p_grid),
         prior = 1) 

## # A tibble: 50 x 3
##    n_points p_grid prior
##       <dbl>  <dbl> <dbl>
##  1        5  0         1
##  2        5  0.25      1
##  3        5  0.5       1
##  4        5  0.75      1
##  5        5  1         1
##  6       10  0         1
##  7       10  0.111     1
##  8       10  0.222     1
##  9       10  0.333     1
## 10       10  0.444     1
## # ... with 40 more rows

The remaining step is to calculate the likelihood and posterior for each row.

dat <- tibble(n_points = c(5, 10, 15, 20)) |> 
  mutate(p_grid = map(n_points, ~seq(from = 0, to = 1, length.out = .))) |> 
  unnest(p_grid) |> 
  expand(nesting(n_points, p_grid),
         prior = 1) |>
  mutate(likelihood = dbinom(6, size = 9, prob = p_grid)) |> 
  mutate(posterior = likelihood * prior / sum(likelihood * prior))
dat

## # A tibble: 50 x 5
##    n_points p_grid prior likelihood posterior
##       <dbl>  <dbl> <dbl>      <dbl>     <dbl>
##  1        5  0         1   0        0        
##  2        5  0.25      1   0.00865  0.00188  
##  3        5  0.5       1   0.164    0.0356   
##  4        5  0.75      1   0.234    0.0507   
##  5        5  1         1   0        0        
##  6       10  0         1   0        0        
##  7       10  0.111     1   0.000111 0.0000241
##  8       10  0.222     1   0.00476  0.00103  
##  9       10  0.333     1   0.0341   0.00741  
## 10       10  0.444     1   0.111    0.0241   
## # ... with 40 more rows

Make the plot

Plot posterior ~ p_grid, faceting by n_points

Define a tiny function to be used to label the facets

appender <- function(string, suffix = "points") paste(string, suffix)

ggplot(dat, aes(x = p_grid, y = posterior)) +
  geom_line() +
  geom_point(size=2) +
  geom_vline(aes(xintercept = p_grid[which.max(posterior)]), linetype="dashed") +
  labs(x = "probability of water",
       y = "posterior probability") +
  theme(panel.grid = element_blank()) +
  facet_wrap(~n_points, labeller = as_labeller(appender)) +
  theme_bw(base_size = 16)

Extensions

Different priors

McElrath suggests trying out other, non-uniform priors. Here’s one way to do that. Kurz has a fancier version that uses these in a grid approximation. Here, I just plot the priors.

n_points <- 10
p_grid <- seq(from=0, to=1, length.out=n_points)
prior1 <- ifelse(p_grid < 0.5, 0, 1)
prior2 <- exp(-5 * abs( p_grid - 0.5))

Make this stuff into a stacked tibble

pdat <- bind_rows(
          bind_cols(prior = "step", p_grid = p_grid, prob = prior1),
          bind_cols(prior = "exp",  p_grid = p_grid, prob = prior2),
)
pdat

## # A tibble: 20 x 3
##    prior p_grid   prob
##    <chr>  <dbl>  <dbl>
##  1 step   0     0     
##  2 step   0.111 0     
##  3 step   0.222 0     
##  4 step   0.333 0     
##  5 step   0.444 0     
##  6 step   0.556 1     
##  7 step   0.667 1     
##  8 step   0.778 1     
##  9 step   0.889 1     
## 10 step   1     1     
## 11 exp    0     0.0821
## 12 exp    0.111 0.143 
## 13 exp    0.222 0.249 
## 14 exp    0.333 0.435 
## 15 exp    0.444 0.757 
## 16 exp    0.556 0.757 
## 17 exp    0.667 0.435 
## 18 exp    0.778 0.249 
## 19 exp    0.889 0.143 
## 20 exp    1     0.0821

Plot them, overlaid using color=prior

ggplot(pdat, aes(p_grid, prob, color=prior)) +
  geom_point(size=4) + 
  geom_line(size = 1.5) +
  theme_bw(base_size = 16) +
  theme(legend.position = c(0.15, .80))

• The ifelse prior might correspond to a strong belief that the world must have at least 50% water.
• The exp prior: ???