# simPPLe

## Outline

### Topics

- Last week’s Q2 exercise solution.

### 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
= function(ppl_function, number_of_iterations) {
posterior = 0.0
numerator = 0.0
denominator for (i in 1:number_of_iterations) {
<<- 1.0
weight = ppl_function()
g_i = numerator + weight * g_i
numerator = denominator + weight
denominator
}return(numerator/denominator)
}
###### contents of the scaffold
suppressPackageStartupMessages(library(distr))
## Utilities to make the distr library a bit nicer to use
<- function(distribution, realization) {
p d(distribution)(realization) # return the PMF or density
}
= function(probability_to_get_one) {
Bern 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
<- function(distribution) {
simulate r(distribution)(1) # sample once from the given distribution
}
# Use observe(realization, distribution) for observed random variables
= function(realization, distribution) {
observe # `<<-` lets us modify variables that live in the global scope from inside a function
<<- weight * p(distribution, realization)
weight }
```

- 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)
= rep(0, 4)
data
# simPPLe's description of our "bag of coins" example
= function() {
beta_binomial
# Similar to forward sampling, but use 'observe' when the variable is observed
= simulate(Beta(1, 1))
p 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?**

Click for choices

- \(0.0010394\)
- \(0.0194394\)
- \(0.0826293\)
- \(0.1728016\)
- None of the above

Click for answer

\(0.1728016\).