source("../exercises/ex03_scaffold.R")
source("../exercises/ex03_ppl.R")
source("../../solutions/sol03_posterior.R")
posterior(my_first_probabilistic_program, 100)
[1] 0.03713188
PPLs is 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 limitation of these approaches.
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 coin problem…
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}\]
Here is how to code up the “bag of coin” model in simPPLe for the purpose of computing \(\mathbb{P}(X = 1 | Y = (0, 0, 0, 0))\):
ex03_ppl.R
= rep(0, 4) # "dataset" of four identical coin flips = (0, 0, 0, 0)
coin_flips
# simPPLe's description of our "bag of coins" example
= function() {
my_first_probabilistic_program
# Similar to forward sampling, but use 'observe' when the variable is observed
= simulate(DiscreteDistribution(supp = 0:2))
coin_index for (i in seq_along(coin_flips)) {
= coin_index/2
prob_heads observe(coin_flips[i], Bern(1 - 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 compute this probability as follows:
source("../exercises/ex03_scaffold.R")
source("../exercises/ex03_ppl.R")
source("../../solutions/sol03_posterior.R")
posterior(my_first_probabilistic_program, 100)
[1] 0.03713188
Recall the painful calculation we did last week to get the predictive.
Here is how to do it in simPPLe:
<- function() {
predict_next_flip = simulate(DiscreteDistribution(supp = 0:2))
coin_index = coin_index/2
prob_heads for (i in seq_along(coin_flips)) {
observe(coin_flips[i], Bern(1 - prob_heads))
}= simulate(Bern(1 - prob_heads))
next_flip = 1 - next_flip
same_as_observed return(same_as_observed)
}
posterior(predict_next_flip, 1000)
[1] 0.9690606