Exercise 3: Write your own PPL!

Goals

  • Introduce Monte Carlo integration in continuous spaces.
  • Implement importance sampling.
  • Build a universal probabilistic programming language in <30 lines of code.

Q.1: functions on the unit interval

For this question, use Simple Monte Carlo. The main twist compared to week one is that you will use a continuous random variable.

  1. Write a function mc_estimate that takes a function \(f:[0,1]\to\mathbb{R}\) and outputs a Monte Carlo estimate of \(\int_0^1 f(x)\mathrm{d}x\) using \(n=10000\) independent samples from \({\mathrm{Unif}}(0,1)\).

  2. Consider the function \(f:[0,1] \to [0,\infty)\) given by \[ f(x) = e^{-x^2}. \] It is possible to show that \[ \int_0^1 f(x) \, \mathrm{d}x = \frac{\sqrt{\pi}}{2} \text{erf}(1) \approx 0.7468241, \tag{1}\] where \(\text{erf}(\cdot)\) is the error function. Test your implementation of mc_estimate by checking that it produces an answer close to the value in Equation 1.

  3. Approximate the following integral using mc_estimate: \[ \int_0^1 \sin(\cos(\sin(x))) \mathrm{d}x. \]

Q.2: implementing SNIS for simPPLe

In this question, you will write the function posterior that we used in the PPL introduction.

  1. First, install the package distr, which allows us to work with distributions as objects—a necessary ingredient of every PPL. Load or install it using
if (!require(distr)){
  install.packages("distr")
  require(distr)
}
  1. Read this short tutorial on distr. Nothing to submit for this item.
  2. Read the “scaffold code”, and use distr and two of the functions below to create a fair coin, flip it, and to compute the probability of that flip:
ex03_scaffold.R
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) 
}
  1. Complete the implementation of the function posterior:
posterior = function(ppl_function, number_of_iterations) {
  numerator = 0.0
  denominator = 0.0
  for (i in 1:number_of_iterations) {
    weight <<- 1.0
    # update numerator and denominator
  }
  return(numerator/denominator)
}
  1. Test your program by checking that you can approximate the posterior probability of the fair coin obtained in exercise 1, Q.2.