simPPLe

Outline

Topics

Rationale

We will use simPPLe this week to do probabilistic modelling.

Background

If you missed Tuesday’s lecture last week, make sure to read the introduction to PPLs.

Setting up simPPLe on your computer

  • Let us combine the “scaffold” files and the answer of last week’s Q2 into one file. Copy paste the code below into a file called simple.R:
simple.R
###### solution

posterior = function(ppl_function, number_of_iterations) {
  numerator = 0.0
  denominator = 0.0
  for (i in 1:number_of_iterations) {
    weight <<- 1.0
    g_i = ppl_function()
    numerator = numerator + weight * g_i
    denominator = denominator + weight
  }
  return(numerator/denominator)
}

###### contents of the scaffold

suppressPackageStartupMessages(library(distr))

## Utilities to make the distr library a bit nicer to use

p <- function(distribution, realization) {
  d(distribution)(realization) # return the PMF or density 
}

Bern = function(probability_to_get_one) {
  DiscreteDistribution(supp = 0:1, prob = c(1-probability_to_get_one, probability_to_get_one))
}

## Key functions called by simPPLe programs

# Use simulate(distribution) for unobserved random variables
simulate <- function(distribution) {
  r(distribution)(1) # sample once from the given distribution
}

# Use observe(realization, distribution) for observed random variables
observe = function(realization, distribution) {
  # `<<-` lets us modify variables that live in the global scope from inside a function
  weight <<- weight * p(distribution, realization) 
}
  • This way, you can now load simPPLe by loading simple.R into your session (note: in the code below replace ../../solutions/ by the path to the file you just created)
source("../../solutions/simple.R")

# define your model as a function, e.g. my_function

# call: posterior(my_function, 100) to approximate posterior E[my_function | observations]

Testing

I will start this Tuesday’s lecture with the following simPPLe test (a continuous version of the rocket model): (note: in the code below replace ../../solutions/ by the path to the file you just created)

source("../../solutions/simple.R")

set.seed(1)

data = rep(0, 4) 

# simPPLe's description of our "bag of coins" example
beta_binomial = function() {
  
  # Similar to forward sampling, but use 'observe' when the variable is observed
  p = simulate(Beta(1, 1))
  for (i in seq_along(data)) { 
    observe(data[i], Bern(p)) 
  }
  
  # return the test function, here the parameter p
  return(p)
}

posterior(beta_binomial, 100)

What is the output?

  1. \(0.0010394\)
  2. \(0.0194394\)
  3. \(0.0826293\)
  4. \(0.1728016\)
  5. None of the above

\(0.1728016\).