# Probabilistic programming from scratch

Working with uncertainty in real-world data.

Working with uncertainty in real-world data.

Real-world data is almost always incomplete or innaccurate in some way. This means that the uncertain conclusions we draw from it are only meaningful if we can answer the question: how uncertain?

One way to do this is using Bayesian inference. But, while Bayesian inference is conceptually simple, it can be analytically and computationally difficult in practice. Probabilistic programming is a paradigm that abstracts away some of this complexity.

There are many probabilistic programming systems. Perhaps the most advanced is Stan, and the most accessible to non-statistician programmers is PyMC3. At Fast Forward Labs, we recently shared with our clients a detailed report on the technology and uses of probabilistic programming in startups and enterprises.

But in this article, rather than use either of these advanced comprehensive systems, we’re going to build our own extremely simple system from scratch.

We’ll write clear, functional Python 3. We’ll use generators to build up a pipeline that will allow us to answer concrete questions. We won’t use any libraries (except for random number generation and plotting). And I’ll go easy on the mathematics. The code will be slow compared to Stan and PyMC3, but hopefully you’ll understand every line.

This “from scratch” approach follows in the footsteps of Joel Grus’s book Data Science from Scratch, and Jake VanderPlas’s PyCon talk Statistics for Hackers. I recommend both. In his talk, Jake said, “if you can write a for loop, you can do statistical analysis.” That isn’t always true (good luck implementing ADVI with a for loop). But when it is true, we can focus on the big, fundamental ideas without getting lost in algebraic or computational details.

Let’s take a specific data analysis problem: a simple A/B test for a website. Suppose our site has two layouts. During our test, 4% of visitors to layout A convert (i.e., buy something, sign up for the mailing list, whatever), and 5% to layout B convert. Clearly, layout B is better, so we should use that layout, right?

But what if I tell you it was a very small test?

`n_visitors_a`

`=`

`100`

`# number of visitors shown layout A`

`n_conv_a`

`=`

`4`

`# number of vistors shown layout A who converted`

`n_visitors_b`

`=`

`40`

`n_conv_b`

`=`

`2`

Are you still sure B is better? And what if it’s going to cost us $1 million to change the layout if we get this decsion wrong. Are you sure *enough*? If not, how much more data would you need?

To answer these questions, we need to quantify exactly how confident we are that layout B is better, given the slice of data we do have.

We can do that using Bayesian inference. Bayesian inference is a method for updating your knowledge about the world with the information you learn during an experiment. It derives from a simple equation called Bayes’s Rule. In its most advanced and efficient forms, it can be used to solve huge problems. But we’re going use a specific, simple inference algorithm called Approximate Bayesian Computation (ABC), which is barely a couple of lines of Python:

`def`

`posterior_sampler`

`(`

`data`

`,`

`prior_sampler`

`,`

`simulate`

`):`

`'''Yield samples from the posterior by Approximate Bayesian Computation.'''`

`for`

`p`

`in`

`prior_sampler`

`:`

`if`

`simulate`

`(`

`p`

`)`

`==`

`data`

`:`

`yield`

`p`

This function turns the *prior* distribution into the *posterior*. What does that mean?

I talk about these distributions in more detail in the Orioles, but for this article, the rough idea is sufficient: samples from the prior distribution are our best guesses of the values of the unknown parameter of our system. In the case of an A/B test, this is the conversion fraction of a layout. These guesses are made before we do the experiment.

Samples from the posterior distribution, meanwhile, are guesses of the same parameters made after the experiment, in the light of the data we gathered. Once you have the posterior, you can answer concrete questions about the implications of the data, such as how likely it is that layout B is better, given our data.

In the function above, the `prior_sampler`

argument is a generator that yields samples from the prior. We run the function `simulate`

on a single sample from the prior. This simulates the experiment, assuming the prior sample is correct. We then compare the simulated outcome to the real outcome (`data`

). If these agree, the sample from the prior can be used as a sample from the posterior. If not, we try again with another sample from the prior, until they do agree.

After we run a lot of simulations using a lot of samples from the prior, we’ll have a good idea of which values were most likely to produce the observations we in fact saw, and are therefore most likely to be correct.

And that’s the whole algorithm. Abstracted away in that little `posterior_sampler`

generator function, it’s an extremely lightweight probabilistic programming system.

Let’s use it to finish our A/B test, starting with layout A. We need to prepare three arguments: `data`

, `prior_sampler`

and `simulate`

. We already have the `data`

for our A/B test. Let’s now write a function that simulates the conversion of `n_visitors`

visitors to a website with known probability `p`

:

`import`

`random`

`def`

`simulate_conversion`

`(`

`p`

`,`

`n_visitors`

`):`

`'''Return number of vistors who convert, given conversion fraction p.'''`

`outcomes`

`=`

`(`

`random`

`.`

`random`

`()`

`<`

`p`

`for`

`i`

`in`

`range`

`(`

`n_visitors`

`))`

`return`

`sum`

`(`

`outcomes`

`)`

Here’s what happens when we run this function a few times to simulate 100 visitors converting with probability 0.1:

`simulate_conversion`

`(`

`0.1`

`,`

`100`

`),`

`simulate_conversion`

`(`

`0.1`

`,`

`100`

`),`

`simulate_conversion`

`(`

`0.1`

`,`

`100`

`)`

`(`

`5`

`,`

`11`

`,`

`12`

`)`

Effectively, this function runs a fake A/B test in which we already know the conversion fraction. Here’s how it works. The `random`

function returns a random floating point number between 0 and 1. That number will be smaller than `p`

with probability equal to `p`

. We use a generator expression to do this `n_visitors`

times. `outcomes`

is then an iterable of booleans, `True`

if that trial was a success (i.e., the ith visitor converted) and `False`

if not. Finally, we take advantage of the fact that the `sum`

of an iterable of booleans is equal to the number of `True`

elements it contains.

We use a generator comprehension in this function because `n_visitors`

could potentially be very large, and we want to do these calculations lazily and avoid holding the

outcome for each visitor in memory simultaneously. (If Python generators are new to you, I recommend you watch or read Ned Batchelder’s fantastic PyCon 2013 talk, “How To Loop Like a Native” for an introduction. If you’re not using generators, you’re not using the full power of Python. I talk about them much more in the Orioles.

Finally we need `prior_sampler`

. This should be a generator that yields a large (potentially infinite) number of guesses for the conversion fraction of a layout. Suppose we’ve never used layout A. Logically, we know that it must be somewhere between 0 and 100%, but other than that, we have no idea. This generator captures that idea:

`def`

`uniform_prior_sampler`

`():`

`'''Yield random numbers in interval (0, 1).'''`

`while`

`True`

`:`

`yield`

`random`

`.`

`random`

`()`

Let’s write a quick `take`

function that will allow us to peek at a few samples from this generator.

`import`

`itertools`

`def`

`take`

`(`

`n`

`,`

`iterable`

`):`

`"Return first n items of the iterable as a list."`

`return`

`list`

`(`

`itertools`

`.`

`islice`

`(`

`iterable`

`,`

`n`

`))`

We can use this to draw three samples from our “uniform” prior sample like this:

`take`

`(`

`3`

`,`

`uniform_prior_sampler`

`())`

`[`

`0.8943997782145745`

`,`

`0.9196888444316101`

`,`

`0.6267420468068673`

`]`

Now we’re ready to run `posterior_sampler`

to create an object that will

yield up samples from the posterior distribution for layout A’s conversion

fraction:

`posterior_a_sampler`

`=`

`posterior_sampler`

`(`

`data`

`=`

`n_conv_a`

`,`

`prior_sampler`

`=`

`uniform_prior_sampler`

`(),`

`simulate`

`=`

`lambda`

`p`

`:`

`simulate_conversion`

`(`

`p`

`,`

`n_visitors_a`

`)`

`)`

(The `lambda`

function here bakes in the number of visitors for layout A into the simulate argument. This is “partial evaluation.” We could also have used `functools.partial`

for this.)

To get a few samples from the posterior, we can use `take`

again:

`take`

`(`

`3`

`,`

`posterior_a_sampler`

`)`

`[`

`0.03810094949124154`

`,`

`0.067278098392767`

`,`

`0.05786540496383208`

`]`

These are three guesses for the unknown conversion

fraction of layout A. They are all around 4% percent. That’s consistent with the

fact that we had four conversions from 100 visitors. But note they are not

*exactly* 4%. We can build up a picture of their distribution by getting

thousands of samples from `posterior_sampler`

and plotting a histogram.

`a_samples`

`=`

`take`

`(`

`10000`

`,`

`posterior_a_sampler`

`)`

`abbins`

`=`

`[`

`i`

`/`

`200.0`

`for`

`i`

`in`

`range`

`(`

`50`

`)]`

`# 50 bins between 0 and 0.25`

`plt`

`.`

`hist`

`(`

`a_samples`

`,`

`bins`

`=`

`abbins`

`,`

`normed`

`=`

`True`

`)`

`# normed=True gives a probability density function`

`plt`

`.`

`xlim`

`(`

`0`

`,`

`max`

`(`

`abbins`

`));`

This is the posterior distribution for the conversion fraction of layout A. It’s a full picture of our knowledge after the experiment.

The posterior distribution is the object we use to answer concrete questions about our knowledge after an experiment. We can, for example, use this one to say how likely it is that layout A’s conversion fraction is greater than 10%. We can do this by measuring the fraction of `a_samples`

for which that is true:

`sum`

`(`

`a`

`>`

`0.1`

`for`

`a`

`in`

`a_samples`

`)`

`/`

`len`

`(`

`a_samples`

`)`

` ``0.0202`

Again, we’re using the fact that the sum of an iterable of booleans is equal the number of true elements.

This tells us there’s a 2% chance layout A’s conversion fraction is greater than 10% (or equivalently, we can be 98% certain it’s less than 10%).

To complete the A/B test, however, we also need the posterior for layout B.

Let’s suppose that layout B is *not* brand new. We do have a rough idea of its conversion fraction before we conduct the experiment. Perhaps we think it’s 6%, plus or minus a couple of per cent. In this case, `uniform_prior_sampler`

is not appropriate, but something like this would work well:

`def`

`normal_prior_sampler`

`(`

`mu`

`=`

`0.06`

`,`

`sigma`

`=`

`0.02`

`):`

`'''Yield stream of samples from N(mu, sigma) in interval (0, 1).'''`

`while`

`True`

`:`

`x`

`=`

`random`

`.`

`normalvariate`

`(`

`mu`

`,`

`sigma`

`)`

`if`

`0`

`<=`

`x`

`<=`

`1`

`:`

`yield`

`x`

Mathematically, this function samples from a normal distribution with a known mean and standard deviation, and support in the interval (0, 1). Conceptually, though, it’s easier to plot a histogram, and compare this prior to the uniform prior we used for layout A.

`plt`

`.`

`hist`

`(`

`take`

`(`

`100000`

`,`

`uniform_prior_sampler`

`()),`

`bins`

`=`

`abbins`

`,`

`label`

`=`

`'A'`

`,`

`normed`

`=`

`True`

`)`

`plt`

`.`

`hist`

`(`

`take`

`(`

`100000`

`,`

`normal_prior_sampler`

`()),`

`bins`

`=`

`abbins`

`,`

`label`

`=`

`'B'`

`,`

`alpha`

`=`

`0.5`

`,`

`normed`

`=`

`True`

`)`

`plt`

`.`

`title`

`(`

`'Prior distributions'`

`)`

`plt`

`.`

`xlim`

`(`

`0`

`,`

`max`

`(`

`abbins`

`))`

`plt`

`.`

`legend`

`();`

The guesses for layout B are concentrated around 6%, which captures our prior knowledge.

Now that we have the prior sampler for layout B, we can make its posterior sampler and take some samples.

`posterior_b_sampler`

`=`

`posterior_sampler`

`(`

`data`

`=`

`n_conv_b`

`,`

`prior_sampler`

`=`

`normal_prior_sampler`

`(),`

`simulate`

`=`

`lambda`

`p`

`:`

`simulate_conversion`

`(`

`p`

`,`

`n_visitors_b`

`)`

`)`

`b_samples`

`=`

`take`

`(`

`10000`

`,`

`posterior_b_sampler`

`)`

Let’s visualize the two sets of posterior samples directly.

`plt`

`.`

`hist`

`(`

`a_samples`

`,`

`bins`

`=`

`abbins`

`,`

`label`

`=`

`'A'`

`,`

`normed`

`=`

`True`

`)`

`plt`

`.`

`hist`

`(`

`b_samples`

`,`

`bins`

`=`

`abbins`

`,`

`label`

`=`

`'B'`

`,`

`alpha`

`=`

`0.5`

`,`

`normed`

`=`

`True`

`)`

`plt`

`.`

`title`

`(`

`'Posterior distributions'`

`)`

`plt`

`.`

`xlim`

`(`

`0`

`,`

`max`

`(`

`abbins`

`));`

`plt`

`.`

`legend`

`();`

Note the histogram for layout B is to the right of that for layout A. This is telling us that layout B is probably better than layout A.

And we’re finally in a position to answer the original question: how confident can we be that layout B is better?

In this simple case, we can do this by making pairwise comparisons between the two lists of estimates and measuring the fraction of pairs for which the estimated conversion fraction of B is bigger than that for A.

`sum`

`(`

`b`

`>`

`a`

`for`

`a`

`,`

`b`

`in`

`zip`

`(`

`a_samples`

`,`

`b_samples`

`))`

`/`

`len`

`(`

`a_samples`

`)`

` ``0.659`

And that’s it! Given the data and our prior beliefs, we are 66% sure that layout B is better than A.

If a mistake is going to cost us $1 million, this probably isn’t confident enough. We need to run a bigger A/B test before making this business decision. If we collect more data, however, we’ll quickly run into one of the limitations of ABC: it’s slow. Change the data and rerun this code to see for yourself.

Approximate Bayesian Computation is indeed slow and in some ways crude. There are often better options. But it’s not cheating. In the Orioles, I show how to speed it up for more realistic problems.

This article hopefully helped you to form a mental model of the relationship between the prior and the posterior distributions. This should leave you in good shape to look deeper into probabilistic programming and inference.