Probabilistic Programming

The first thing that makes a language a probabilistic programming language (PPL) is a set of primitives for drawing random numbers. At this point, a PPL looks like any old imperative language with a rand call. Here's an incredibly boring webppl program:

var b = flip(0.5); b ? "yes" : "no"

This boring program just uses the outcome of a fair coin toss to return one string or another. It works exactly like an ordinary program with access to a flip function for producing random Booleans. Functions like flip are sometimes called elementary random primitives, and they're the source of all randomness in these programs.

If you run this little program in the webppl editor, you'll see that it does exactly what you expect: sometimes it outputs “yes,” and sometimes it prints “no.”

Things get slightly more interesting when we realize that webppl can represent entire distributions, not just individual values. The webppl language has an Enumerate operation, which prints out all the probabilities in a distribution defined by a function:

var roll = function () { var die1 = randomInteger(6) + 1; var die2 = randomInteger(6) + 1; return die1 + die2; }

var dist = Enumerate(roll); print(dist); viz.auto(dist);

You get a printout with all the possible die values between 2 and 12 and their associated probabilities. That viz.auto call also instructs webppl's in-browser interface to display a pretty graph.

This may not look all that surprising, since you could imagine writing Enumerate in your favorite language by just running the roll function over and over. But in fact, Enumerate is doing something a bit more powerful. It's not sampling executions to get an approximation of the distribution; it's actually enumerating every possible execution of the function to get an exact distribution. This begins to reveal the point of a probabilistic programming language: the tools that analyze PPL programs are the important part, not actually executing the programs directly.

This is enough to code up the math for our paper-recommender model. We can write functions to encode the relevance piece and the class registration piece, and we can test it out by randomly generating “student profiles.”

// Class attendance model. var attendance = function(i_pl, i_stats, busy) { var attendance = function (interest, busy) { if (interest) { return busy ? flip(0.3) : flip(0.8); } else { return flip(0.1); } } var a_4110 = attendance(i_pl, busy); var a_4780 = attendance(i_stats, busy); var a_4242 = attendance(i_pl && i_stats, busy);

return {cs4110: a_4110, cs4780: a_4780, cs4242: a_4242}; }

// Relevance of our three papers. var relevance = function(i_pl, i_stats) { var rel1 = i_pl && i_stats; var rel2 = i_pl; var rel3 = i_stats;

return {paper1: rel1, paper2: rel2, paper3: rel3}; }

// A combined model. var model = function() { // Some even random priors for our "student profile." var i_pl = flip(0.5); var i_stats = flip(0.5); var busy = flip(0.5);

return [relevance(i_pl, i_stats), attendance(i_pl, i_stats, busy)]; }

var dist = Enumerate(model); viz.auto(dist);

Running this will show the distribution over all the observed data. This isn't terribly useful, but it is interesting. We can see, for example, that if we know nothing else about a student, our model says they're quite likely to take none of the classes and to be interested in none of the papers.

The next important piece of a PPL is a conditioning construct. Conditioning lets you determine how much weight to give to a given execution in a program. Crucially, you can choose the weight of an execution partway through the run—after doing some of the computation. You can even mark an execution as completely irrelevant, effectively filtering the executions to a subset.

This conditioning operation is particularly useful for encoding observations. If you know something about the outcome of a process, you can use conditioning to communicate that the observation must be true (or is likely to be true). Let's return to our dice example for a contrived example. Say that we saw the first die, and it's a 4. We can condition on this information to give us the distribution on total values given that one of the dice is a 4:

var roll = function () { var die1 = randomInteger(6) + 1; var die2 = randomInteger(6) + 1;

// Only keep executions where at least one die is a 4. if (!(die1 === 4 || die2 === 4)) { factor(-Infinity); }

return die1 + die2; }

Unsurprisingly, the distribution is flat except for the outcome 8, which is less likely because only one 4-containing roll can produce that total. Outcomes like 2 and 12 have probability zero because they cannot occur when one of the dice comes up 4.

What if, on the other hand, we don't get to see either of the dice, but someone told us that the total on the dice was 10. What does this tell us about the values of the dice themselves? We can encode this observation by conditioning on the outcome of the roll.

var roll_condition = function () { var die1 = randomInteger(6) + 1; var die2 = randomInteger(6) + 1;

// Discard any executions that don't sum to 10. var out = die1 + die2; if (out !== 10) { factor(-Infinity); }

// Return the values on the dice. return die1 + "+" + die2; }

The results probably don't surprise you, but we can use the same principle with our recommender.

Let's use the same philosophy now to actually produce recommendations. It's simple: we just need to condition on the class registration of the person we're interested in. Here's an example that describes me: I attend my own class, CS 4110, and the fictional PPL class, CS 4242, but not the ML class, 4780.

// A model query that describes my class attendance. var rec = function() { var i_pl = flip(0.5); var i_stats = flip(0.5); var busy = flip(0.5);

// Require my class attendance. var att = attendance(i_pl, i_stats, busy); require(att.cs4242 && att.cs4110 && !att.cs4780);

return relevance(i_pl, i_stats); }

(In this example, we define require to wrap factor(-Infinity) and completely eliminate program executions that don't satisfy a condition.) Calling Enumerate on this rec function finally gives us something useful: a distribution over paper relevance! It's a little easier to understand if we just look at one paper at a time:

return relevance(i_pl, i_stats).paper1;

Suddenly, this is pretty nifty! By telling the enumerator which executions are relevant to us, it can tell us what it knows about data under those conditions. To me, and I think to most programmers, this already feels like a much more natural tool for expressing probabilistic models. Rather than carefully writing down which random variables depend on which others, you use the flow of program execution to build up those dependencies.

The point of this section: Writing generative models feels very comfortable and straightforward, and a programming language is a great way to write down a generative algorithm. We've successfully made our job easy by shifting the burden of doing inference on these simple models to the compiler and tools.

That Enumerate operation might not look like much, but it's actually an implementation of the central problem that PPLs were meant to solve: statistical inference. Especially in the presence of conditioning, inference on a general probabilistic program is a hard problem. Think about how Enumerate must work: it needs to know the outcome of a program for every single valuation of the program's random primitive draws. It's not hard to see how this grows exponentially. And it's even less clear if you introduce floating-point numbers: if you draw a number between 0.0 and 1.0, that's a lot of possible valuations—and it's not representable as a histogram.

Enumerate just won't do for nontrivial problems. It's for this reason that efficient inference algorithms for PPLs are a huge focus in the community. Here are a couple of other inference algorithms that don't involve exploring every possible execution of a program.

The second most obvious inference algorithm uses sampling. The idea is to run the program a large number of times, drawing different random values for each random primitive on each execution. Apply the program's conditioning to weight each sample and total them all up. The interaction with weighting makes this strategy rejection sampling, so called because you reject some executions when you reach a conditioning statement.

The webppl language has this strategy built in. It's called ParticleFilter, and we can run it on our examples from above:

var sampled = ParticleFilter(rec('paper1'), 1000);

2.5.2. MCMC

Rejection sampling works for small examples, but it runs into trouble in the presence of conditioning. It can waste a lot of work taking samples that don't matter (i.e., they're destined to be rejected when they hit a factor call). Smarter sampling strategies exist, the most prominent of which are Markov chain Monte Carlo methods.

The webppl standard library provides MCMC algorithms. Making MCMC correct and efficient is a popular problem in PPL research, and a full discussion is currently out of scope of this lecture.

To be as starry-eyed as possible, the promise of PP is nothing less than the democratization of machine learning. In case you're not convinced yet, here are a sampling of popular applications of probabilistic programming.