What is a PPL?

Outline

Topics

  • What is a Probabilistic Programming Language (PPL)?
  • How to use a PPL.

Rationale

PPLs are the main way you will compute posterior distributions in this course. This week you will write your own! This will give you insight on the strengths and limitations of these approaches.

PPL in a nutshell

  • It is often easier to do forward sampling than to compute a posterior:
    • Forward sampling: just go down a single path in the decision tree.
    • Computing a posterior: need to sum over all paths compatible with observed data.
  • PPLs allow you to:
    • code your model as if you were going to do forward sampling,
    • and the PPL will magically figure out how to approximate the posterior!

PPL: an example using simPPLe

We created possibly the simplest possible PPL for this course. Let’s call it simPPLe. The exercise this week will be to understand simPPLe by filling in a couple of key lines of code.

First I will show you what simPPLe can do.

Let us start with something we are familiar with: our bag of coins problem…

Mathematical description

Recall our “bag of coins” model can be written as:

\[ \begin{align*} X &\sim {\mathrm{Unif}}\{0, 1, 2\} \\ Y_i | X &\sim {\mathrm{Bern}}(X/2) \end{align*} \tag{1}\]

PPL description of the same model

Here is how to code up the “bag of coins” model in simPPLe for the purpose of computing \(\mathbb{P}(X = 1 | Y = (0, 0, 0, 0))\):

ex03_ppl.R
coin_flips = rep(0, 4) # "dataset" of four identical coin flips = (0, 0, 0, 0) 

# simPPLe's description of our "bag of coins" example
my_first_probabilistic_program = function() {
  
  # Similar to forward sampling, but use 'observe' when the variable is observed
  coin_index = simulate(DiscreteDistribution(supp = 0:2))
  prob_heads = coin_index/2
  for (i in seq_along(coin_flips)) { 
    observe(coin_flips[i], Bern(prob_heads)) 
  }
  
  # return the test function g(x, y)
  return(ifelse(coin_index == 1, 1, 0))
}

After solving this week’s exercise, you will be able to approximate this probability as follows:

source("../exercises/ex03_scaffold.R")
source("../exercises/ex03_ppl.R")
source("../../solutions/sol03_posterior.R") 

posterior(my_first_probabilistic_program, 10000)
[1] 0.05665221

Compare this to the solution of question 2.2 in exercise 1.

Extension: predicting the next draw

Recall the painful calculation we did last week to get the predictive.

Here is how to do it in simPPLe:

predict_next_flip <- function() {
  coin_index = simulate(DiscreteDistribution(supp = 0:2))
  prob_heads = coin_index/2
  for (i in seq_along(coin_flips)) { 
    observe(coin_flips[i], Bern(prob_heads)) 
  }
  next_flip = simulate(Bern(prob_heads))
  return(ifelse(next_flip == 0, 1, 0))
}

posterior(predict_next_flip, 10000)
[1] 0.9695877