probabilistic-programming-from-scratch.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "matplotlib.rcParams['figure.figsize'] = (10, 6)\n",
    "import random\n",
    "random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probabilistic programming from scratch\n",
    "\n",
    "> This article contains highlights from a series of [three interactive video tutorials](https://www.safaribooksonline.com/search/?query=%22Probabilistic%20Programming%20from%20Scratch%22&extended_publisher_data=true&highlight=true&is_academic_institution_account=false&source=user&include_assessments=false&include_courses=true&include_orioles=true&include_playlists=true&publishers=O%27Reilly%20Media%2C%20Inc.&field=title&sort=relevance&utm_source=oreilly&utm_medium=newsite&utm_campaign=probabilistic-programming-from-scratch-top-cta-orioles-link) on probabilistic programming from scratch published on O'Reilly Safari (login required).\n",
    ">\n",
    "> If you're interested in the business case for probabilistic programming the [Fast Forward Labs report](http://blog.fastforwardlabs.com/2017/01/18/new-research-on-probabilistic-programming.html) discusses it in detail, and compares modern industrial strength systems like Stan and PyMC3. Please [get in touch](http://www.fastforwardlabs.com/#contact) if you're interested in working with us.\n",
    ">\n",
    "> This article is [available as a Jupyter Notebook](https://github.com/fastforwardlabs/probabilistic-programming-from-scratch/blob/master/probabilistic-programming-from-scratch.ipynb).\n",
    "\n",
    "Real-world data is almost always incomplete or inaccurate in some way. This means that the uncertain conclusions we draw from it are only meaningful if we can answer the question: how uncertain?\n",
    "\n",
    "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.\n",
    "\n",
    "There are many probabilistic programming systems. Perhaps the most advanced is [Stan](http://mc-stan.org/), and the most accessible to non-statistician programmers is [PyMC3](https://pymc-devs.github.io/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](http://blog.fastforwardlabs.com/2017/01/18/new-research-on-probabilistic-programming.html).\n",
    "\n",
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "This \"from scratch\" approach follows in the footsteps of Joel Grus's book [Data Science from Scratch](http://shop.oreilly.com/product/0636920033400.do), and Jake\n",
    "VanderPlas's PyCon talk [Statistics for Hackers](https://www.youtube.com/watch?v=Iq9DzN6mvYA). 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](https://arxiv.org/abs/1603.00788) 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.\n",
    "\n",
    "## An A/B test\n",
    "\n",
    "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?\n",
    "\n",
    "But what if I tell you it was a very small test?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_visitors_a = 100  # number of visitors shown layout A\n",
    "n_conv_a = 4        # number of vistors shown layout A who converted\n",
    "\n",
    "n_visitors_b = 40\n",
    "n_conv_b = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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?\n",
    "\n",
    "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.\n",
    "\n",
    "## A simple algorithm for Bayesian inference\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def posterior_sampler(data, prior_sampler, simulate):\n",
    "    '''Yield samples from the posterior by Approximate Bayesian Computation.'''\n",
    "    for p in prior_sampler:\n",
    "        if simulate(p) == data:\n",
    "            yield p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function turns the _prior_ distribution into the _posterior_. What does that mean?\n",
    "\n",
    "I talk about these distributions in more detail in [the video tutorials](https://www.safaribooksonline.com/search/?query=%22Probabilistic%20Programming%20from%20Scratch%22&extended_publisher_data=true&highlight=true&is_academic_institution_account=false&source=user&include_assessments=false&include_courses=true&include_orioles=true&include_playlists=true&publishers=O%27Reilly%20Media%2C%20Inc.&field=title&sort=relevance&utm_source=oreilly&utm_medium=newsite&utm_campaign=probabilistic-programming-from-scratch-top-cta-orioles-link), 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.\n",
    "\n",
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "After we run a lot of simulations using lots 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.\n",
    "\n",
    "## Using the sampler\n",
    "\n",
    "And that's the whole algorithm. Abstracted away in that little `posterior_sampler` generator function, it's an extremely lightweight probabilistic programming system.\n",
    "\n",
    "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`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def simulate_conversion(p, n_visitors):\n",
    "    '''Return number of vistors who convert, given conversion fraction p.'''\n",
    "    outcomes = (random.random() < p for i in range(n_visitors))\n",
    "    return sum(outcomes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what happens when we run this function a few times to simulate 100 visitors converting with probability 0.1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 11, 12)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simulate_conversion(0.1, 100), simulate_conversion(0.1, 100), simulate_conversion(0.1, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. \n",
    "\n",
    "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\n",
    "outcome for each visitor in memory simultaneously. (If Python generators are new to you, I recommend you [watch](https://www.youtube.com/watch?v=EnSu9hHGq5o) or [read](https://nedbatchelder.com/text/iter.html) 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 video tutorials.)\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def uniform_prior_sampler():\n",
    "    '''Yield random numbers in interval (0, 1).'''\n",
    "    while True:\n",
    "        yield random.random()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's write a quick `take` function that will allow us to peek at a few samples from this generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "def take(n, iterable):\n",
    "    \"Return first n items of the iterable as a list.\"\n",
    "    return list(itertools.islice(iterable, n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this to draw 3 samples from our \"uniform\" prior sample like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.8943997782145745, 0.9196888444316101, 0.6267420468068673]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "take(3, uniform_prior_sampler())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to run `posterior_sampler` to create an object that will\n",
    "yield up samples from the posterior distribution for layout A's conversion\n",
    "fraction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "posterior_a_sampler = posterior_sampler(\n",
    "    data=n_conv_a,\n",
    "    prior_sampler=uniform_prior_sampler(),\n",
    "    simulate=lambda p: simulate_conversion(p, n_visitors_a)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(The `lambda` function here bakes in the number of visitors for layout A into the simulate argument. This is \"partial eveluation\". We could also have used `functools.partial` for this.)\n",
    "\n",
    "To get a few samples from the posterior, we can use `take` again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.03810094949124154, 0.067278098392767, 0.05786540496383208]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "take(3, posterior_a_sampler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAFpCAYAAACmgZ0NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFh1JREFUeJzt3X2MZmd5H+DfXS8QlVjFxIMDtpeB1kUyUTDN1DSFtBDA\nGJti2tLEVpuaxNEmtEiNGqndlAoi0j82qpI0qaO4W7CAKjVWPgxW12C2JK2DxNfaNdh8xcZZhBfH\nXjDlI1Aik7t/7HE7jGc8s/POPPN1XdLRnPOc55xzjx7P6ufnnPe81d0BAGBz/aWtLgAAYC8QugAA\nBhC6AAAGELoAAAYQugAABhC6AAAGELoAAAYQugAABhC6AAAGELoAAAbYt9UFLOfss8/u+fn5rS4D\nAGBVt99++5e6e261ftsydM3Pz+fYsWNbXQYAwKqq6vNr6ef2IgDAAEIXAMAAQhcAwABCFwDAAEIX\nAMAAQhcAwABCFwDAAEIXAMAAQhcAwABCFwDAAEIXAMAAQhcAwABCFwDAAPtW61BV5yd5Z5JzknSS\nw93961X11CQ3JplPcjzJj3X3V5Y5/uok/3ba/Hfd/Y6NKZ3lzB88clr9jx+6fJMqAQAWWzV0JXkk\nyc939x1VdWaS26vqaJLXJflAdx+qqoNJDib514sPnILZm5Ms5FRgu72qbl4unLE1TiekCWgAsH6r\n3l7s7ge6+45p/etJPp3k3CRXJHl01uodSV6zzOGvSHK0ux+egtbRJJduROEAADvJaT3TVVXzSZ6f\n5CNJzunuB6Zdf5pTtx+XOjfJFxZt3z+1AQDsKWsOXVX1vUl+L8nPdffXFu/r7s6p24frVlUHqupY\nVR07efLkLKcCANh21hS6quoJORW4fru7f39qfrCqnj7tf3qSh5Y59ESS8xdtnze1PUZ3H+7uhe5e\nmJubW2v9AAA7wqqhq6oqyduSfLq7f3XRrpuTXD2tX53kPcscfmuSS6rqrKo6K8klUxsAwJ6ylpmu\nFyb5iSQ/WlV3TstlSQ4leXlV3ZPkZdN2qmqhqt6aJN39cJJfSvKxaXnL1AYAsKes+sqI7v5gklph\n90uX6X8syU8v2r4+yfXrLRAAYDfwRnoAgAGELgCAAYQuAIABhC4AgAGELgCAAYQuAIABhC4AgAGE\nLgCAAYQuAIABhC4AgAGELgCAAYQuAIABhC4AgAGELgCAAYQuAIABhC4AgAGELgCAAYQuAIAB9m11\nAaxu/uCRrS4BAJiR0MWanU74O37o8k2sBAB2HrcXAQAGELoAAAYQugAABhC6AAAGELoAAAYQugAA\nBhC6AAAGELoAAAYQugAABlj1jfRVdX2SVyV5qLt/YGq7Mclzpi5PSfK/u/uiZY49nuTrSb6T5JHu\nXtigugEAdpS1fA3Q25Ncm+SdjzZ0948/ul5Vv5Lkq49z/Eu6+0vrLRAAYDdYNXR1921VNb/cvqqq\nJD+W5Ec3tiwAgN1l1me6fiTJg919zwr7O8n7q+r2qjow47UAAHastdxefDxXJbnhcfa/qLtPVNXT\nkhytqs90923LdZxC2YEk2b9//4xlAQBsL+ue6aqqfUn+QZIbV+rT3Semnw8luSnJxY/T93B3L3T3\nwtzc3HrLAgDYlma5vfiyJJ/p7vuX21lVT66qMx9dT3JJkrtnuB4AwI61auiqqhuSfCjJc6rq/qq6\nZtp1ZZbcWqyqZ1TVLdPmOUk+WFUfT/LRJEe6+30bVzoAwM6xlk8vXrVC++uWaftiksum9fuSPG/G\n+gAAdoVZH6SHZc0fPHJa/Y8funyTKgGA7cHXAAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwg\ndAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAMsG+r\nC9ir5g8e2eoSAICBzHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAM\nIHQBAAwgdAEADLBq6Kqq66vqoaq6e1HbL1bViaq6c1ouW+HYS6vqs1V1b1Ud3MjCAQB2krXMdL09\nyaXLtP9ad180Lbcs3VlVZyT5zSSvTHJhkquq6sJZigUA2KlWDV3dfVuSh9dx7ouT3Nvd93X3nyd5\nV5Ir1nEeAIAdb5Znut5QVZ+Ybj+etcz+c5N8YdH2/VMbAMCes97Q9VtJ/mqSi5I8kORXZi2kqg5U\n1bGqOnby5MlZTwcAsK2sK3R194Pd/Z3u/osk/zmnbiUudSLJ+Yu2z5vaVjrn4e5e6O6Fubm59ZQF\nALBtrSt0VdXTF23+/SR3L9PtY0kuqKpnVdUTk1yZ5Ob1XA8AYKfbt1qHqrohyYuTnF1V9yd5c5IX\nV9VFSTrJ8SQ/M/V9RpK3dvdl3f1IVb0hya1JzkhyfXd/clN+CwCAbW7V0NXdVy3T/LYV+n4xyWWL\ntm9J8pjXSQAA7DXeSA8AMIDQBQAwgNAFADDAqs90wQjzB4+sue/xQ5dvYiUAsDnMdAEADCB0AQAM\nIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0\nAQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEADCB0AQAMIHQBAAwgdAEA\nDLBq6Kqq66vqoaq6e1Hbv6+qz1TVJ6rqpqp6ygrHHq+qu6rqzqo6tpGFAwDsJGuZ6Xp7kkuXtB1N\n8gPd/YNJ/jjJLzzO8S/p7ou6e2F9JQIA7Hyrhq7uvi3Jw0va3t/dj0ybH05y3ibUBgCwa2zEM10/\nleS9K+zrJO+vqtur6sAGXAsAYEfaN8vBVfXGJI8k+e0Vuryou09U1dOSHK2qz0wzZ8ud60CSA0my\nf//+WcoCANh21j3TVVWvS/KqJP+4u3u5Pt19Yvr5UJKbkly80vm6+3B3L3T3wtzc3HrLAgDYltYV\nuqrq0iT/Ksmru/ubK/R5clWd+eh6kkuS3L1cXwCA3W4tr4y4IcmHkjynqu6vqmuSXJvkzJy6ZXhn\nVV039X1GVd0yHXpOkg9W1ceTfDTJke5+36b8FgAA29yqz3R191XLNL9thb5fTHLZtH5fkufNVB0A\nwC4x04P0sBXmDx45rf7HD12+SZUAwNr5GiAAgAGELgCAAYQuAIABhC4AgAGELgCAAYQuAIABhC4A\ngAGELgCAAYQuAIABhC4AgAGELgCAAYQuAIABhC4AgAGELgCAAYQuAIABhC4AgAH2bXUBu8X8wSNb\nXQIAsI2Z6QIAGEDoAgAYQOgCABhA6AIAGEDoAgAYQOgCABhA6AIAGMB7utj1TucdascPXb6JlQCw\nl5npAgAYQOgCABhA6AIAGGBNoauqrq+qh6rq7kVtT62qo1V1z/TzrBWOvXrqc09VXb1RhQMA7CRr\nnel6e5JLl7QdTPKB7r4gyQem7e9SVU9N8uYkL0hycZI3rxTOAAB2szWFru6+LcnDS5qvSPKOaf0d\nSV6zzKGvSHK0ux/u7q8kOZrHhjcAgF1vlme6zunuB6b1P01yzjJ9zk3yhUXb909tAAB7yoY8SN/d\nnaRnOUdVHaiqY1V17OTJkxtRFgDAtjFL6Hqwqp6eJNPPh5bpcyLJ+Yu2z5vaHqO7D3f3QncvzM3N\nzVAWAMD2M0voujnJo59GvDrJe5bpc2uSS6rqrOkB+kumNgCAPWWtr4y4IcmHkjynqu6vqmuSHEry\n8qq6J8nLpu1U1UJVvTVJuvvhJL+U5GPT8papDQBgT1nTdy9291Ur7HrpMn2PJfnpRdvXJ7l+XdUB\nAOwS3kgPADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUA\nMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA\n0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwwL71HlhVz0ly46KmZyd5U3f/\nh0V9XpzkPUn+ZGr6/e5+y3qvCZtt/uCRNfc9fujyTawEgN1m3aGruz+b5KIkqaozkpxIctMyXf+o\nu1+13usAAOwGG3V78aVJPtfdn9+g8wEA7CobFbquTHLDCvt+uKo+XlXvrarnbtD1AAB2lJlDV1U9\nMcmrk/zOMrvvSPLM7n5ekv+Y5N2Pc54DVXWsqo6dPHly1rIAALaVjZjpemWSO7r7waU7uvtr3f2N\naf2WJE+oqrOXO0l3H+7uhe5emJub24CyAAC2j40IXVdlhVuLVfX9VVXT+sXT9b68AdcEANhR1v3p\nxSSpqicneXmSn1nU9rNJ0t3XJXltktdX1SNJvpXkyu7uWa4JALATzRS6uvvPknzfkrbrFq1fm+Ta\nWa4BALAbeCM9AMAAM8107Xan83ZyAIDHY6YLAGAAoQsAYAChCwBgAKELAGAAoQsAYAChCwBgAKEL\nAGAAoQsAYAChCwBgAKELAGAAoQsAYAChCwBgAKELAGAAoQsAYAChCwBgAKELAGAAoQsAYAChCwBg\nAKELAGAAoQsAYAChCwBgAKELAGCAfVtdAOxU8wePnFb/44cu36RKANgJzHQBAAwgdAEADCB0AQAM\nIHQBAAwwc+iqquNVdVdV3VlVx5bZX1X1G1V1b1V9oqr+xqzXBADYaTbq04sv6e4vrbDvlUkumJYX\nJPmt6SfsKafzaUefdATYfUbcXrwiyTv7lA8neUpVPX3AdQEAto2NCF2d5P1VdXtVHVhm/7lJvrBo\n+/6pDQBgz9iI24sv6u4TVfW0JEer6jPdfdvpnmQKbAeSZP/+/RtQFgDA9jHzTFd3n5h+PpTkpiQX\nL+lyIsn5i7bPm9qWnudwdy9098Lc3NysZQEAbCszha6qenJVnfnoepJLkty9pNvNSf7p9CnGv5Xk\nq939wCzXBQDYaWa9vXhOkpuq6tFz/dfufl9V/WySdPd1SW5JclmSe5N8M8lPznhNAIAdZ6bQ1d33\nJXneMu3XLVrvJP98lusAAOx03kgPADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUA\nMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA0AUAMIDQBQAwgNAFADCA\n0AUAMMC+rS4AeKz5g0dOq//xQ5dvUiUAbBQzXQAAAwhdAAADCF0AAAPsqWe6Tvc5GQCAjWKmCwBg\nAKELAGAAoQsAYAChCwBggHWHrqo6v6r+sKo+VVWfrKp/sUyfF1fVV6vqzml502zlAgDsTLN8evGR\nJD/f3XdU1ZlJbq+qo939qSX9/qi7XzXDdQAAdrx1z3R19wPdfce0/vUkn05y7kYVBgCwm2zIM11V\nNZ/k+Uk+sszuH66qj1fVe6vquRtxPQCAnWbml6NW1fcm+b0kP9fdX1uy+44kz+zub1TVZUneneSC\nFc5zIMmBJNm/f/+sZQEAbCvV3es/uOoJSf5bklu7+1fX0P94koXu/tLj9VtYWOhjx46tu66VeCM9\nJMcPXb7VJQDsKlV1e3cvrNZvlk8vVpK3Jfn0SoGrqr5/6pequni63pfXe00AgJ1qltuLL0zyE0nu\nqqo7p7Z/k2R/knT3dUlem+T1VfVIkm8lubJnmVoDANih1h26uvuDSWqVPtcmuXa91wAA2C28kR4A\nYAChCwBgAKELAGAAoQsAYAChCwBggJnfSA/sLKfzkmAvUgXYOGa6AAAGELoAAAYQugAABhC6AAAG\nELoAAAYQugAABhC6AAAGELoAAAbwclRgRafzItXEy1QBHo+ZLgCAAYQuAIABhC4AgAGELgCAAYQu\nAIABfHoR2DCn82lHn3QE9hozXQAAAwhdAAAD7Pjbi6f78kYAgK1gpgsAYIAdP9MF7Eweugf2GjNd\nAAADCF0AAAMIXQAAAwhdAAADzPQgfVVdmuTXk5yR5K3dfWjJ/icleWeSH0ry5SQ/3t3HZ7kmsPds\np1fDeKgfWK91z3RV1RlJfjPJK5NcmOSqqrpwSbdrknylu/9akl9L8svrvR4AwE42y0zXxUnu7e77\nkqSq3pXkiiSfWtTniiS/OK3/bpJrq6q6u2e4LsCW2axXXWzmbJ7ZOdgeZgld5yb5wqLt+5O8YKU+\n3f1IVX01yfcl+dLjnfiuE1/dVrcTANZju/w7tl0C3enWISyy22ybl6NW1YEkB6bNb3/+l19191bW\nw8zOzirhmh3BOO4OmzaOtYkPjWzmuXcof4/b1zPX0mmW0HUiyfmLts+b2pbrc39V7UvyV3LqgfrH\n6O7DSQ4nSVUd6+6FGWpjixnD3cE47g7GcXcwjjvfLK+M+FiSC6rqWVX1xCRXJrl5SZ+bk1w9rb82\nyR94ngsA2IvWPdM1PaP1hiS35tQrI67v7k9W1VuSHOvum5O8Lcl/qap7kzycU8EMAGDPmemZru6+\nJcktS9retGj9/yT5R+s49eFZ6mJbMIa7g3HcHYzj7mAcd7hytw8AYPP5GiAAgAGGhq6qurSqPltV\n91bVwWX2P6mqbpz2f6Sq5hft+4Wp/bNV9YqRdfPd1juOVTVfVd+qqjun5brRtfP/rWEc/05V3VFV\nj1TVa5fsu7qq7pmWq5ceyzgzjuN3Fv09Lv0gFAOtYRz/ZVV9qqo+UVUfqKpnLtrn73Gn6O4hS049\nbP+5JM9O8sQkH09y4ZI+/yzJddP6lUlunNYvnPo/KcmzpvOcMap2y4aN43ySu7f6d7CseRznk/xg\nTn1/6msXtT81yX3Tz7Om9bO2+nfai8ss4zjt+8ZW/w6WNY/jS5L85Wn99Yv+XfX3uIOWkTNd/+9r\ng7r7z5M8+rVBi12R5B3T+u8meWlV1dT+ru7+dnf/SZJ7p/Mx3izjyPax6jh29/Hu/kSSv1hy7CuS\nHO3uh7v7K0mOJrl0RNE8xizjyPaxlnH8w+7+5rT54Zx6N2bi73FHGRm6lvvaoHNX6tPdjyR59GuD\n1nIsY8wyjknyrKr6X1X1P6vqRza7WFY0y9+Uv8ftY9ax+J6qOlZVH66q12xsaZyG0x3Ha5K8d53H\nsoW2zdcAsSc8kGR/d3+5qn4oybur6rnd/bWtLgz2qGd294mqenaSP6iqu7r7c1tdFCurqn+SZCHJ\n393qWjh9I2e6Tudrg7Lka4PWcixjrHscp9vDX06S7r49p55h+OubXjHLmeVvyt/j9jHTWHT3ienn\nfUn+R5Lnb2RxrNmaxrGqXpbkjUle3d3fPp1j2R5Ghq5Zvjbo5iRXTp+Ke1aSC5J8dFDdfLd1j2NV\nzVXVGUky/Z/1BTn10CfjrWUcV3Jrkkuq6qyqOivJJVMb4617HKfxe9K0fnaSFyb51KZVyuNZdRyr\n6vlJ/lNOBa6HFu3y97iTjHxqP8llSf44p2Y43ji1vSWn/iNKku9J8js59aD8R5M8e9Gxb5yO+2yS\nV271JxD28rLecUzyD5N8MsmdSe5I8ve2+nfZy8saxvFv5tTzIX+WUzPOn1x07E9N43tvkp/c6t9l\nLy/rHcckfzvJXTn1Sbm7klyz1b/LXl7WMI7/PcmD07+fdya5edGx/h53yOKN9AAAA3gjPQDAAEIX\nAMAAQhcAwABCFwDAAEIXAMAAQhcAwABCFwDAAEIXAMAA/xdXBusnB2By9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110ae61d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a_samples = take(10000, posterior_a_sampler)\n",
    "\n",
    "abbins = [i/200.0 for i in range(50)]          # 50 bins between 0 and 0.25\n",
    "plt.hist(a_samples, bins=abbins, normed=True)  # normed=True gives a probability density function\n",
    "plt.xlim(0, max(abbins));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the posterior distribution for the conversion fraction of layout A. It's a full picture of our knowledge after the experiment.\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0202"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(a > 0.1 for a in a_samples)/len(a_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we're using the fact that the sum of an iterable of booleans is equal the number of true elements.\n",
    "\n",
    "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%).\n",
    "\n",
    "## Completing the A/B test\n",
    "\n",
    "To complete the A/B test, however, we also need the posterior for layout B.\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normal_prior_sampler(mu=0.06, sigma=0.02):\n",
    "    '''Yield stream of samples from N(mu, sigma) in interval (0, 1).'''\n",
    "    while True:\n",
    "        x = random.normalvariate(mu, sigma)\n",
    "        if 0 <= x <= 1:\n",
    "            yield x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAF1CAYAAADSlV/tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHPZJREFUeJzt3X2QZXV95/H3JzPo+MBGHDpEGWAgohW0ArKzGDdEJSgO\nxIiurBkqIkasyYOuS2I2BdH4lGQxa8U1CdkgKyxouSORRGBXEFGIwgajAwEZFIUQ1B4JDDMKopAw\n43f/6DN6abqnH+7tX9/b/X5V3erz8Dvn9+1zuNMffufcc1NVSJIkaWH92GIXIEmStBwYuiRJkhow\ndEmSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ekgUhyRZJTF3D/leQZ3fQ5SX5/QPs9MMmDSVZ083+b\n5A2D2He3vwU9LpJGx8rFLkDScEpyF7AfsAv4HnAF8KaqenCq9lV1fKvaqurXZ9Ou+x3eUFWf3sO+\nvgE8eRB1JXkn8Iyqek3P/psdF0nDzZEuSXvyS1X1ZOBIYB3wtskNMmHe/5YkWbT/+VvMviUtP4Yu\nSTOqqq1MjHQ9B354Ce6Pkvw/4PvAIb2X5ZL8WJK3Jfl6knuTfCjJj3fr1naXCk9L8g3g6qn6TPJf\nktyd5FtJXj9p3QVJ/rCb3jfJ/03ynSQ7klzb9f9h4EDg/3SXD393qr57lvUGsJ9K8oUkDyS5NMlT\nu75elGR8Ui13JXlxkvXA7wG/3PV3c8+xmstxOTXJN5Lcl+StPf0clWRzV9M9Sd43j1MpaREZuiTN\nKMkBwAnAP/QsPgXYCOwNfH3SJq/rXscAhzBx+e7sSW1eCPw08NIp+lsP/A7wEuBQ4MV7KO8twDgw\nxsTl0N8DqqpOAb5BN1pXVf9tNn13Xgu8HngasBP4sz30DxMdfhL4r8BFXX+HT9Hsdcx8XI4GngUc\nC7w9yU93y/8U+NOq+jfATwF/NVNNkoaLoUvSnlyS5DvAdcBnmQgVu11QVbdW1c6qemTSdr8CvK+q\n7uzuATsT2DBpNOmdVfW9qnpoin5fDfyvqtpSVd8D3rmHGh9hIhwdVFWPVNW1NfOXyu6pb4AP9/T9\n+8Crd99o36fZHJd3VdVDVXUzcDOwO7w9Ajwjyb5V9WBVfX4A9UhqyNAlaU9eUVVPqaqDquo3J4WU\nb+5hu6fz6NGvrzPxwZ395rB97/rJI2m93gvcAXwqyZ1JzthD29n0PXn914G9gH1nsd+ZzOa4/HPP\n9Pf50U3+pwHPBG5L8sUkLxtAPZIaMnRJmq89jSZ9CzioZ/5AJi7T3TPL7e8GDpi0/dRFVH23qt5S\nVYcALwd+O8mxM/Qx00jY5L4fAe5j4lOcT9y9ohv9GpvDfmdzXKZUVbdX1cnATwB/DFyc5EkzbSdp\neBi6JC2ETcBvJTk4yZP50b1OO2e5/V8Br0tyWJInAu+YrmGSlyV5RpIA9zPxiIsfdKvvYeLeqbl6\nTU/f7wYurqpdwNeAVUl+McleTHya8/E9290DrN3DpznnfVySvCbJWFX9APhOt/gHe9pG0nAxdEla\nCOcDHwY+B/wT8DDwn2a7cVVdAbyfiU823sE0n3DsHAp8GngQuB74H1V1TbfuLOBt3Scbf2cO9X8Y\nuICJS32rgDd3dd0P/CbwQWArEyNfvZ9m/Fj3c3uSG6fYbz/HZT1wa5IHmbipfsMe7kmTNIQy8/2m\nkiRJ6pcjXZIkSQ0YuiRJkhowdEmSJDVg6JIkSWrA0CVJktTAypmbtLfvvvvW2rVrF7sMSZKkGd1w\nww33VdXYTO2GMnStXbuWzZs3L3YZkiRJM0qyp68q+yEvL0qSJDVg6JIkSWrA0CVJktTAUN7TJUmS\nlo9HHnmE8fFxHn744cUuZY9WrVrFmjVr2Guvvea1/YyhK8kBwIeA/YACzq2qP03yVOAiYC1wF/Dq\nqvr2FNufCrytm/3DqrpwXpVKkqQlaXx8nL333pu1a9eSZLHLmVJVsX37dsbHxzn44IPntY/ZXF7c\nCbylqg4DfhZ4Y5LDgDOAz1TVocBnuvlH6YLZO4DnAUcB70iyz7wqlSRJS9LDDz/M6tWrhzZwASRh\n9erVfY3GzRi6quruqrqxm/4u8BVgf+BEYPeo1YXAK6bY/KXAVVW1oxsFuwpYP+9qJUnSkjTMgWu3\nfmuc0430SdYCzwX+Htivqu7uVv0zE5cfJ9sf+GbP/Hi3TJIkaahccsklJOG2225bkP3P+kb6JE8G\n/ho4vaoe6E17VVVJqp9CkmwENgIceOCB/exKkiSNsLVnfGKg+7vrPb84q3abNm3i6KOPZtOmTbzr\nXe8aaA0wy5GuJHsxEbg+UlV/0y2+J8nTuvVPA+6dYtOtwAE982u6ZY9RVedW1bqqWjc2NuOT9CVJ\nkgbmwQcf5LrrruO8887jox/96IL0MWPoysSQ1nnAV6rqfT2rLgNO7aZPBS6dYvMrgeOS7NPdQH9c\nt0ySJGloXHrppaxfv55nPvOZrF69mhtuuGHgfcxmpOvngFOAX0hyU/c6AXgP8JIktwMv7uZJsi7J\nBwGqagfwB8AXu9e7u2WSJElDY9OmTWzYsAGADRs2sGnTpoH3MeM9XVV1HTDd7frHTtF+M/CGnvnz\ngfPnW6AkSdJC2rFjB1dffTW33HILSdi1axdJeO973zvQT1X6RHqNnmvOmlv7Y85cmDokSUvCxRdf\nzCmnnMIHPvCBHy574QtfyLXXXssLXvCCgfXjdy9KkqRlbdOmTbzyla981LJXvepVA7/E6EiXJEka\nKrN9xMOgXHPNNY9Z9uY3v3ng/Ri6NBzmeslQkqQR4+VFSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDo\nkiRJasDQJUmSlr0VK1ZwxBFHcPjhh3PkkUfyd3/3dwPvw0dGSJKk4TLoxwjN4ptJnvCEJ3DTTTcB\ncOWVV3LmmWfy2c9+dqBlONIlSZLU44EHHmCfffYZ+H4d6ZIkScveQw89xBFHHMHDDz/M3XffzdVX\nXz3wPgxdkiRp2eu9vHj99dfz2te+li1btpBkYH14eVGSJKnH85//fO677z62bds20P060qWF4Xcp\nSpJG1G233cauXbtYvXr1QPdr6JIkScve7nu6AKqKCy+8kBUrVgy0D0OXJEkaLrN4xMOg7dq1a8H7\nMHRp6ZvLpc5FeKNLkpYHb6SXJElqwNAlSZLUgKFLkiQtuqpa7BJm1G+Nhi5JkrSoVq1axfbt24c6\neFUV27dvZ9WqVfPehzfSS5KkRbVmzRrGx8cH/jDSQVu1ahVr1qyZ9/aGLkmStKj22msvDj744MUu\nY8F5eVGSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1MOMjI5KcD7wMuLeqntMtuwh4\nVtfkKcB3quqIKba9C/gusAvYWVXrBlS3JEnSSJnNc7ouAM4GPrR7QVX98u7pJH8C3L+H7Y+pqvvm\nW6AkSdJSMGPoqqrPJVk71bokAV4N/MJgy5IkSVpa+r2n6+eBe6rq9mnWF/CpJDck2binHSXZmGRz\nks3D/jUAkiRJc9Vv6DoZ2LSH9UdX1ZHA8cAbk7xguoZVdW5VrauqdWNjY32WJUmSNFzmHbqSrAT+\nA3DRdG2qamv3817g48BR8+1PkiRplPUz0vVi4LaqGp9qZZInJdl79zRwHLClj/4kSZJG1oyhK8km\n4HrgWUnGk5zWrdrApEuLSZ6e5PJudj/guiQ3A18APlFVnxxc6ZIkSaNjNp9ePHma5a+bYtm3gBO6\n6TuBw/usT5IkaUnwifSSJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA4Yu\nSZKkBgxdkiRJDcz4RHppWbnmrNm3PebMhatDkrTkGLo0e3MJJJIk6VG8vChJktSAoUuSJKkBQ5ck\nSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKk\nBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgMzhq4k5ye5N8mWnmXv\nTLI1yU3d64Rptl2f5KtJ7khyxiALlyRJGiWzGem6AFg/xfL/XlVHdK/LJ69MsgL4C+B44DDg5CSH\n9VOsJEnSqJoxdFXV54Ad89j3UcAdVXVnVf0r8FHgxHnsR5IkaeT1c0/Xm5J8qbv8uM8U6/cHvtkz\nP94tkyRJWnbmG7r+Evgp4AjgbuBP+i0kycYkm5Ns3rZtW7+7kyRJGirzCl1VdU9V7aqqHwD/k4lL\niZNtBQ7omV/TLZtun+dW1bqqWjc2NjafsiRJkobWvEJXkqf1zL4S2DJFsy8ChyY5OMnjgA3AZfPp\nT5IkadStnKlBkk3Ai4B9k4wD7wBelOQIoIC7gF/r2j4d+GBVnVBVO5O8CbgSWAGcX1W3LshvIUmS\nNORmDF1VdfIUi8+bpu23gBN65i8HHvM4CUmSpOXGJ9JLkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOX\nJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmS\npAYMXZIkSQ0YuiRJkhpYudgFaJFdc9ZiVzC65nrsjjlzYeqQJI0ER7okSZIaMHRJkiQ1YOiSJElq\nwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1ICh\nS5IkqQFDlyRJUgMzhq4k5ye5N8mWnmXvTXJbki8l+XiSp0yz7V1JbklyU5LNgyxckiRplMxmpOsC\nYP2kZVcBz6mqnwG+Bpy5h+2Pqaojqmrd/EqUJEkafTOGrqr6HLBj0rJPVdXObvbzwJoFqE2SJGnJ\nGMQ9Xa8HrphmXQGfSnJDko172kmSjUk2J9m8bdu2AZQlSZI0PPoKXUneCuwEPjJNk6Or6kjgeOCN\nSV4w3b6q6tyqWldV68bGxvopS5IkaejMO3QleR3wMuBXqqqmalNVW7uf9wIfB46ab3+SJEmjbF6h\nK8l64HeBl1fV96dp86Qke++eBo4DtkzVVpIkaambzSMjNgHXA89KMp7kNOBsYG/gqu5xEOd0bZ+e\n5PJu0/2A65LcDHwB+ERVfXJBfgtJkqQht3KmBlV18hSLz5um7beAE7rpO4HD+6pOkiRpifCJ9JIk\nSQ0YuiRJkhowdEmSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNGLokSZIa\nMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDo\nkiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAYMXZIkSQ2sXOwCNGDXnLXYFUiSpCkYuqRW5hKI\njzlz4eqQJC0KLy9KkiQ1YOiSJElqwNAlSZLUwKxCV5Lzk9ybZEvPsqcmuSrJ7d3PfabZ9tSuze1J\nTh1U4ZIkSaNktiNdFwDrJy07A/hMVR0KfKabf5QkTwXeATwPOAp4x3ThTJIkaSmbVeiqqs8BOyYt\nPhG4sJu+EHjFFJu+FLiqqnZU1beBq3hseJMkSVry+rmna7+qurub/mdgvyna7A98s2d+vFsmSZK0\nrAzkRvqqKqD62UeSjUk2J9m8bdu2QZQlSZI0NPoJXfckeRpA9/PeKdpsBQ7omV/TLXuMqjq3qtZV\n1bqxsbE+ypIkSRo+/YSuy4Ddn0Y8Fbh0ijZXAscl2ae7gf64bpkkSdKyMttHRmwCrgeelWQ8yWnA\ne4CXJLkdeHE3T5J1ST4IUFU7gD8Avti93t0tkyRJWlZm9d2LVXXyNKuOnaLtZuANPfPnA+fPqzpJ\nkqQlwifSS5IkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAYMXZIkSQ0YuiRJkhowdEmS\nJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1YOiSJElq\nwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1ICh\nS5IkqQFDlyRJUgMrF7sAzcI1Zy12BZIkqU+OdEmSJDUw75GuJM8CLupZdAjw9qp6f0+bFwGXAv/U\nLfqbqnr3fPuUlo25jm4ec+bC1CFJGph5h66q+ipwBECSFcBW4ONTNL22ql42334kSZKWgkFdXjwW\n+Meq+vqA9idJkrSkDCp0bQA2TbPu+UluTnJFkmdPt4MkG5NsTrJ527ZtAypLkiRpOPQdupI8Dng5\n8LEpVt8IHFRVhwN/Dlwy3X6q6tyqWldV68bGxvotS5IkaagMYqTreODGqrpn8oqqeqCqHuymLwf2\nSrLvAPqUJEkaKYMIXSczzaXFJD+ZJN30UV1/2wfQpyRJ0kjp6+GoSZ4EvAT4tZ5lvw5QVecAJwG/\nkWQn8BCwoaqqnz4lSZJGUV+hq6q+B6yetOycnumzgbP76UOSJGkp8In0kiRJDRi6JEmSGjB0SZIk\nNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAYMXZIkSQ0YuiRJkhowdEmSJDVg6JIkSWrA\n0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFL\nkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNbBysQtYtq45a7ErkCRJDRm6pKVgLiH+\nmDMXrg5J0rS8vChJktRA36EryV1JbklyU5LNU6xPkj9LckeSLyU5st8+JUmSRs2gLi8eU1X3TbPu\neODQ7vU84C+7n5IkSctGi8uLJwIfqgmfB56S5GkN+pUkSRoagwhdBXwqyQ1JNk6xfn/gmz3z490y\nSZKkZWMQlxePrqqtSX4CuCrJbVX1ubnupAtsGwEOPPDAAZQlSZI0PPoe6aqqrd3Pe4GPA0dNarIV\nOKBnfk23bPJ+zq2qdVW1bmxsrN+yJEmShkpfoSvJk5LsvXsaOA7YMqnZZcBru08x/ixwf1Xd3U+/\nkiRJo6bfy4v7AR9Psntf/7uqPpnk1wGq6hzgcuAE4A7g+8Cv9tmnJEnSyOkrdFXVncDhUyw/p2e6\ngDf2048kSdKo84n0kiRJDRi6JEmSGjB0SZIkNWDokiRJamBQ370oSZKWgbVnfGJO7e96zy8uUCWj\nx9AlNfL+z3xt1m1PP/aZC1iJevkHpD9zOX4eOy13hi4NhVEMJHOpWUuHIWNp8DxqMRi6FslSDxnD\nUrOGl3/0+uPxG06OnGpPRj50LeQ/PHN988zF6SN/5JcmR68ebSHfA3PhHzItpoV8HxieH20hj/Uw\nHL9l9ad/WP6AzNVCBgFHpJafUX0fLHWeFy0m//trY1mFLi1Pjl5Jo8lRoKXBQPcjQxm6btl6vydJ\nkgbEy7NLgyG0P8Nw/IYydGn0ObrUn7kePy8Ta7nzf9Q1CgxdGjkGusea2zH5mQWrY1j4B7g/Hj9p\nYRi6ljkDjKRWDHNa7gxdkrRADBmSehm6JElaYgz8/Vmo4/djC7JXSZIkPYojXdIyc/rKi2fd9v07\nT1rASiRpeXGkS5IkqQFDlyRJUgNeXhyQuVyykSRJy48jXZIkSQ0YuiRJkhowdEmSJDVg6JIkSWrA\n0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNzDt0JTkgyTVJvpzk1iT/eYo2L0pyf5Kbutfb\n+ytXkiRpNPXzNUA7gbdU1Y1J9gZuSHJVVX15Urtrq+plffQjSZI08uY90lVVd1fVjd30d4GvAPsP\nqjBJkqSlZCD3dCVZCzwX+PspVj8/yc1Jrkjy7D3sY2OSzUk27/r+/YMoS5IkaWj0HbqSPBn4a+D0\nqnpg0uobgYOq6nDgz4FLpttPVZ1bVeuqat2KJ/54v2VJkiQNlb5CV5K9mAhcH6mqv5m8vqoeqKoH\nu+nLgb2S7NtPn5IkSaOon08vBjgP+EpVvW+aNj/ZtSPJUV1/2+fbpyRJ0qjq59OLPwecAtyS5KZu\n2e8BBwJU1TnAScBvJNkJPARsqKrqo09JkqSRNO/QVVXXAZmhzdnA2fPtQ9LiOn3lxXNq//6dJy1Q\nJZI0+voZ6Vry5voHR5IkaTp+DZAkSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAl\nSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJamDlYhfQ0ukrL17s\nEqQlbS7vsffvPGkBK5Gk4eNIlyRJUgOGLkmSpAYMXZIkSQ0YuiRJkhowdEmSJDVg6JIkSWrA0CVJ\nktSAoUuSJKmBkX84qg88lUaTD1KVtNw40iVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGhvJG+p/g294g\nL+mH5vrvgTfeSxpGfY10JVmf5KtJ7khyxhTrH5/kom793ydZ209/kiRJo2reI11JVgB/AbwEGAe+\nmOSyqvpyT7PTgG9X1TOSbAD+GPjlfgqWpJn4OApJw6ifka6jgDuq6s6q+lfgo8CJk9qcCFzYTV8M\nHJskffQpSZI0kvq5p2t/4Js98+PA86ZrU1U7k9wPrAbu66NfSRqYhbx/1FE0Sb2G5kb6JBuBjd3s\nv/zWH12wZTHrUd/2xXC9FHge+3LBYhewm+dxafA8Dq+DZtOon9C1FTigZ35Nt2yqNuNJVgI/Dmyf\namdVdS5wLkCSzVW1ro/atMg8h0uD53Fp8DwuDZ7H0dfPPV1fBA5NcnCSxwEbgMsmtbkMOLWbPgm4\nuqqqjz4lSZJG0rxHurp7tN4EXAmsAM6vqluTvBvYXFWXAecBH05yB7CDiWAmSZK07PR1T1dVXQ5c\nPmnZ23umHwb+4zx2fW4/dWkoeA6XBs/j0uB5XBo8jyMuXu2TJElaeH73oiRJUgNNQ1c/XxuU5Mxu\n+VeTvLRl3Xq0+Z7HJGuTPJTkpu51Tuva9SOzOI8vSHJjkp1JTpq07tQkt3evUydvq3b6PI+7et6P\nkz8IpYZmcR5/O8mXk3wpyWeSHNSzzvfjqKiqJi8mbrb/R+AQ4HHAzcBhk9r8JnBON70BuKibPqxr\n/3jg4G4/K1rV7mtg53EtsGWxfwdfsz6Pa4GfAT4EnNSz/KnAnd3PfbrpfRb7d1qOr37OY7fuwcX+\nHXzN+jweAzyxm/6Nnn9XfT+O0KvlSFc/Xxt0IvDRqvqXqvon4I5uf2rPr39aGmY8j1V1V1V9CfjB\npG1fClxVVTuq6tvAVcD6FkXrMfo5jxoeszmP11TV97vZzzPxbEzw/ThSWoauqb42aP/p2lTVTmD3\n1wbNZlu10c95BDg4yT8k+WySn1/oYjWtft5Tvh+HR7/nYlWSzUk+n+QVgy1NczDX83gacMU8t9Ui\nGpqvAdKycDdwYFVtT/JvgUuSPLuqHljswqRl6qCq2prkEODqJLdU1T8udlGaXpLXAOuAFy52LZq7\nliNdc/naICZ9bdBstlUb8z6P3eXh7QBVdQMT9zA8c8Er1lT6eU/5fhwefZ2Lqtra/bwT+FvguYMs\nTrM2q/OY5MXAW4GXV9W/zGVbDYeWoaufrw26DNjQfSruYOBQ4AuN6tajzfs8JhlLsgKg+z/rQ5m4\n6VPtzeY8TudK4Lgk+yTZBziuW6b25n0eu/P3+G56X+DngC8vWKXakxnPY5LnAh9gInDd27PK9+Mo\naXnXPnAC8DUmRjje2i17NxP/EQGsAj7GxI3yXwAO6dn2rd12XwWOX+xPICzn13zPI/Aq4FbgJuBG\n4JcW+3dZzq9ZnMd/x8T9Id9jYsT51p5tX9+d3zuAX13s32U5v+Z7HoF/D9zCxCflbgFOW+zfZTm/\nZnEePw3c0/37eRNwWc+2vh9H5OUT6SVJkhrwifSSJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLU\ngKFLkiSpAUOXJElSA4YuSZKkBv4/wfBn5gpEzlwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110c7cfd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(take(100000, uniform_prior_sampler()), bins=abbins, label='A', normed=True)\n",
    "plt.hist(take(100000, normal_prior_sampler()), bins=abbins, label='B', alpha=0.5, normed=True)\n",
    "plt.title('Prior distributions')\n",
    "plt.xlim(0, max(abbins))\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The guesses for layout B are concentrated around 6%, which captures our prior knowledge.\n",
    "\n",
    "Now we have the prior sampler for layout B we can make its posterior sampler and take some samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "posterior_b_sampler = posterior_sampler(\n",
    "    data=n_conv_b,\n",
    "    prior_sampler=normal_prior_sampler(),\n",
    "    simulate=lambda p: simulate_conversion(p, n_visitors_b)\n",
    ")\n",
    "b_samples = take(10000, posterior_b_sampler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the two sets of posterior samples directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAF1CAYAAADMXG9eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF7NJREFUeJzt3X+w5WV9H/D3R5a4RBiDy4ai/LgY3WSiE9GuWqNBmVWD\nWouO1uIkio0t+dHE7qTTyW7MjNoZA45NAkk6VQyOGJNFS6MwxRRxWX8QNXGxVEARFKEurvxYtIgB\nK+vTP+53yXWzu/fc+5yz5/54vWbO3HO/3+f7fT7nPhx483yf8z3VWgsAAIvzqGkXAACwnAlTAAAd\nhCkAgA7CFABAB2EKAKCDMAUA0EGYAiaiqh6oqidO6NwvqKpdc36/qapeMKZz/1JVfWzO762qnjSO\ncw/nm9jfBZgOYQqWuaq6vaoeHP4jfVdVva+qju4438wQINb01NVaO7q1dlvPORbQ11Naa584VJtR\nX1dr7S9aay8eR11V9Ymq+jf7nf+w/V2Aw0OYgpXh5a21o5M8I8nGJL83rUJ6Q1jv8cu1b2D5EqZg\nBWmt3Znkr5M8NUmq6vFVdUVV3VdVX62qf7uvbVU9q6p2VtX9w4zWHw67PjX8/M4w2/Wcof2vVNWX\nq+rbVXVVVZ0y51ytqv5dVd2a5NY52540PH9sVb2/qu6pqjuq6veq6lHDvjdU1d9U1R9V1Z4kb93/\ndVXVUcOM27er6ktJnrnf/tur6oULfV0H6nvYdu1+Jby0qm6rqnur6p1zan9rVX1gTh2PzH5V1duT\n/EKSPx36+9NF/F2urar/PLzur1fVS+b09Yahpu8O+37pQP9MAJPn/8JgBamqk5K8NMlfDZsuTXJj\nkscn+ZkkV1fV11pr1yS5MMmFrbU/Hy4LPnU45vQkX0/yE621h4fznpXkd5O8PLNhaUuSbUl+fk73\nr0jy7CQPHqC0P0ny2CRPTLIuyceS7E5y8bD/2UOtxyc58gDHvyXJTw2Px2Q2MB7MQl7XTx+g7391\ngHO+MrMzfkcn+XiSryT5s0PUkNbam6vquUk+0Fo7WNtR/i6XJDkuyblJLq6qJyT58SR/nOSZrbWv\nVNUJSR53qHqAyTEzBSvDR6rqO0muTfLJJL8/BKvnJvmd1tpDrbXrMxsAXj8c84MkT6qq41prD7TW\nPneI8/9akvNaa18egsjvJzlt7uzUsP++1tqPhKmqOiLJ2Um2tta+21q7PckfJHndnGbfbK39SWvt\n4f2PH7wmyduH838js0HiYBbyukbpO0neMfT9f5JckOS185xzXiP+Xe5orb2ntbY3s6HqhMyGviT5\nYZKnVtVRrbXdrbWbemsCFkeYgpXhFa21n2itndJa+40hFDw+yX2tte/OaXdHkicMz9+YZEOSm6vq\n81X1zw9x/lOSXFhV3xlC231Jas65kuQbBzn2uMzO+NxxkDoOdew+j9+vzR0Ha5iFva5R+t6/zR1D\nPb1G+bt8a9+T1trfD0+Pbq19L7MzaL+WZHdVXVlVPzOGmoBFEKZg5fpmksdV1TFztp2c5M4kaa3d\n2lp7bZKfTPKOJJdV1WOStAOc6xtJfnUIbPseR7XWPjOnzYGOS5J7MztbNHcW65E65jl2n91JTtrv\n+ANa4Osape8coO9vDs+/l9lLbvv8kwWce5S/y0G11q5qrb0os7NVNyd5zyjHAeMnTMEKNVwO+0yS\n86pqbVX9XGZnbT6QJFX1y1W1vrX2wyTfGQ77YZJ7hp9z74X0riRbq+opw7GPrap/OWIde5N8KMnb\nq+qY4dLgb++rY0QfGvo/tqpOTPJbB2u4wNc1qv849H1Skn+f5IPD9uuTnF5VJ1fVY5Ns3e+4uw7W\nX8/fpaqOr6qzhpD4/SQPDK8NmAJhCla21yaZyexMyoeTvKW19vFh35lJbqqqBzK7aPvs1tqDw+Wk\ntyf5m+Gy3j9rrX04s7M8l1bV/Zld1P6SjO63MjuLc1tm13X9ZZL3LuD4t2X2EtjXM7tI+88P0Xbk\n17WA/i9Pcl1mw9OVGRaIt9auzmyw+uKw/3/sd9yFSV49fBrvQOu8Fvt3eVRmg9c3M3vJ9flJfn0B\nrwcYo2ptlBluAAAOxMwUAEAHYQoAoIMwBQDQQZgCAOggTAEAdDis38133HHHtZmZmcPZJQDAolx3\n3XX3ttbWz9fusIapmZmZ7Ny583B2CQCwKFV1qK+ueoTLfAAAHYQpAIAOwhQAQIfDumYKAFg9fvCD\nH2TXrl156KGHpl3KIa1duzYnnnhijjzyyEUdL0wBABOxa9euHHPMMZmZmUlVTbucA2qtZc+ePdm1\na1dOPfXURZ3DZT4AYCIeeuihrFu3bskGqSSpqqxbt65r9kyYAgAmZikHqX16axSmAIAV7SMf+Uiq\nKjfffPNEzm/NFABwWMxsuXKs57v9/JeN1G7btm153vOel23btuVtb3vbWGtIzEwBACvYAw88kGuv\nvTYXX3xxLr300on0IUwBACvW5ZdfnjPPPDMbNmzIunXrct111429D2EKAFixtm3blrPPPjtJcvbZ\nZ2fbtm1j78OaKQBgRbrvvvtyzTXX5IYbbkhVZe/evamqvPOd7xzrpwyFKSZrx3kLa3/G1snUAcCq\nc9lll+V1r3td3v3udz+y7fnPf34+/elP5/TTTx9bPy7zAQAr0rZt2/LKV77yR7a96lWvGvulPjNT\nLNxCZ5sAIKPfymBcduzY8Y+2velNbxp7P8IUS8tCgppLggAsAS7zAQB0EKYAADoIUwAAHYQpAIAO\nwhQAQAdhCgBYsY444oicdtppedrTnpZnPOMZ+cxnPjP2PtwaAQA4PMZ9n8IRbpFz1FFH5frrr0+S\nXHXVVdm6dWs++clPjrUMM1MAwKpw//3359hjjx37ec1MAQAr1oMPPpjTTjstDz30UHbv3p1rrrlm\n7H0IUwDAijX3Mt9nP/vZvP71r8+NN96YqhpbHy7zAQCrwnOe85zce++9ueeee8Z6XmEKAFgVbr75\n5uzduzfr1q0b63ld5gMAVqx9a6aSpLWWSy65JEccccRY+xCmAIDDY4RbGYzb3r17J96Hy3wAAB2E\nKQCADsIUAEAHYQoAmJjW2rRLmFdvjcIUADARa9euzZ49e5Z0oGqtZc+ePVm7du2izzHvp/mq6qQk\n709yfJKW5KLW2oVV9bgkH0wyk+T2JK9prX170ZUwr5ktV47c9vbzXzbBSgBgfieeeGJ27do19ptk\njtvatWtz4oknLvr4UW6N8HCS/9Ba+0JVHZPkuqq6OskbkmxvrZ1fVVuSbEnyO4uuhLESvACYtiOP\nPDKnnnrqtMuYuHkv87XWdrfWvjA8/26SLyd5QpKzklwyNLskySsmVSQAwFK1oDVTVTWT5OlJ/jbJ\n8a213cOub2X2MuCBjjm3qnZW1c6lPs0HALBQI4epqjo6yX9Psrm1dv/cfW12ZdkBV5e11i5qrW1s\nrW1cv359V7EAAEvNSGGqqo7MbJD6i9baXw2b76qqE4b9JyS5ezIlAgAsXfOGqaqqJBcn+XJr7Q/n\n7LoiyTnD83OSXD7+8gAAlrZRPs333CSvS3JDVV0/bPvdJOcn+VBVvTHJHUleM5kSAQCWrnnDVGvt\n2iR1kN2bxlsOU7HjvGlXAADLljugAwB0EKYAADqMsmYKlqaFXJ48Y+vk6gBgVTMzBQDQQZgCAOgg\nTAEAdBCmAAA6CFMAAB2EKQCADsIUAEAHYQoAoIMwBQDQQZgCAOggTAEAdBCmAAA6CFMAAB3WTLuA\n1Wxmy5XTLiFJcsH2WxbUfvOmDROqBACWHzNTAAAdhCkAgA7CFABAB2EKAKCDMAUA0EGYAgDoIEwB\nAHQQpgAAOghTAAAdhCkAgA7CFABAB2EKAKCDMAUA0EGYAgDoIEwBAHQQpgAAOghTAAAdhCkAgA5r\npl0Ay88F228Zue3mTRsmWAkATJ+ZKQCADsIUAEAHYQoAoIMwBQDQQZgCAOggTAEAdBCmAAA6CFMA\nAB2EKQCADsIUAEAHYQoAoIMwBQDQQZgCAOggTAEAdBCmAAA6rJl2ASvNzJYrp10CAHAYCVOsDjvO\nW1j7M7ZOpg4AVhyX+QAAOghTAAAdhCkAgA7CFABAB2EKAKCDMAUA0MGtEVaozWsum3YJALAqmJkC\nAOggTAEAdJg3TFXVe6vq7qq6cc62t1bVnVV1/fB46WTLBABYmkaZmXpfkjMPsP2PWmunDY+Pjrcs\nAIDlYd4w1Vr7VJL7DkMtAADLTs+aqd+sqi8OlwGPPVijqjq3qnZW1c577rmnozsAgKVnsWHqvyb5\nqSSnJdmd5A8O1rC1dlFrbWNrbeP69esX2R0AwNK0qDDVWrurtba3tfbDJO9J8qzxlgUAsDwsKkxV\n1Qlzfn1lkhsP1hYAYCWb9w7oVbUtyQuSHFdVu5K8JckLquq0JC3J7Ul+dYI1AgAsWfOGqdbaaw+w\n+eIJ1AIAsOy4AzoAQAdhCgCgw7yX+aDHBdtvWVD7zZs2TKgSAJgMM1MAAB2EKQCADsIUAEAHYQoA\noIMwBQDQQZgCAOggTAEAdBCmAAA6CFMAAB2EKQCADsIUAEAHYQoAoIMwBQDQQZgCAOggTAEAdBCm\nAAA6CFMAAB2EKQCADsIUAEAHYQoAoIMwBQDQQZgCAOggTAEAdBCmAAA6CFMAAB2EKQCADsIUAEAH\nYQoAoIMwBQDQQZgCAOggTAEAdFgz7QJgSdpx3uhtz9g6uToAWPKEKZaUC7bfMnLbzZs2TLASABiN\ny3wAAB2EKQCADsIUAEAHYQoAoIMF6MvI5jWXTbsEAGA/ZqYAADoIUwAAHYQpAIAOwhQAQAdhCgCg\ngzAFANBBmAIA6CBMAQB0EKYAADq4A/oIZrZcOe0SAIAlyswUAEAHYQoAoIMwBQDQQZgCAOggTAEA\ndBCmAAA6CFMAAB3cZ4pl64Ltt4zcdvOmDROsBIDVzMwUAEAHYQoAoIMwBQDQYd4wVVXvraq7q+rG\nOdseV1VXV9Wtw89jJ1smAMDSNMrM1PuSnLnfti1JtrfWnpxk+/A7AMCqM2+Yaq19Ksl9+20+K8kl\nw/NLkrxizHUBACwLi10zdXxrbffw/FtJjj9Yw6o6t6p2VtXOe+65Z5HdAQAsTd0L0FtrLUk7xP6L\nWmsbW2sb169f39sdAMCSstgwdVdVnZAkw8+7x1cSAMDysdgwdUWSc4bn5yS5fDzlAAAsL6PcGmFb\nks8m+emq2lVVb0xyfpIXVdWtSV44/A4AsOrM+918rbXXHmTXpjHXAgCw7LgDOgBAB2EKAKCDMAUA\n0EGYAgDoIEwBAHQQpgAAOghTAAAdhCkAgA7CFABAB2EKAKCDMAUA0EGYAgDoMO8XHQPz2HHewtqf\nsXUydQAwFWamAAA6CFMAAB2EKQCADsIUAEAHYQoAoIMwBQDQQZgCAOggTAEAdBCmAAA6CFMAAB2E\nKQCADr6bb4o2r7ls2iUAAJ3MTAEAdDAzxapwwfZbFtR+86YNE6oEgJXGzBQAQAdhCgCggzAFANBB\nmAIA6CBMAQB0EKYAADoIUwAAHYQpAIAOwhQAQIdVeQf0mS1XTrsEAGCFMDMFANBBmAIA6CBMAQB0\nEKYAADoIUwAAHYQpAIAOwhQAQAdhCgCggzAFANBBmAIA6CBMAQB0EKYAADoIUwAAHYQpAIAOwhQA\nQIc10y4AVp0d543e9oytk6sDgLEwMwUA0EGYAgDoIEwBAHQQpgAAOliADgdwwfZbRm67edOGCVYC\nwFJnZgoAoIMwBQDQQZgCAOggTAEAdOhagF5Vtyf5bpK9SR5urW0cR1EAAMvFOD7Nd0Zr7d4xnAeW\npYV88i/x6T+AlcZlPgCADr1hqiX5WFVdV1XnjqMgAIDlpPcy3/Naa3dW1U8mubqqbm6tfWpugyFk\nnZskJ598cmd3AABLS9fMVGvtzuHn3Uk+nORZB2hzUWttY2tt4/r163u6AwBYchYdpqrqMVV1zL7n\nSV6c5MZxFQYAsBz0XOY7PsmHq2rfef6ytfY/x1IVAMAysegw1Vq7LcnTxljLirB5zWXTLgEAOIzc\nGgEAoIMwBQDQQZgCAOggTAEAdBCmAAA6CFMAAB2EKQCADsIUAEAHYQoAoEPP18kAk7bjvNHbnrF1\ncnUAcFBmpgAAOghTAAAdhCkAgA7CFABAB2EKAKCDMAUA0EGYAgDoIEwBAHRw0044zC7YfsvIbTdv\n2jDBSgAYBzNTAAAdhCkAgA7CFABAhxWzZmpmy5XTLgEAWIXMTAEAdBCmAAA6CFMAAB2EKQCADsIU\nAECHFfNpvknavOayaZcAACxRwhSsFDvOW1j7M7ZOpg6AVUaYgiXM9/gBLH3WTAEAdBCmAAA6CFMA\nAB2EKQCADsIUAEAHYQoAoIMwBQDQQZgCAOjgpp2wQizkBp9JsvmMCRUCsMqYmQIA6CBMAQB0EKYA\nADoIUwAAHYQpAIAOwhQAQAdhCgCgg/tMwSo1s+XKkdvefv7LJlgJwPJmZgoAoIMwBQDQYVVe5tu8\n5rJplwAArBBmpgAAOghTAAAdVuVlPmBhfPIP4OCEKVilFrJ28IKHXz3BSgCWN5f5AAA6LNmZqYVc\nVgAAmJYlG6aApWMhlwRntizs3NZYAcvdiglT7h0FAEyDNVMAAB2EKQCADsIUAECHrjVTVXVmkguT\nHJHkz1pr54+lKmDZWvj6xdEXoE/yU74WwgOLtegwVVVHJPkvSV6UZFeSz1fVFa21L42rOGDlcxsU\nYLnrmZl6VpKvttZuS5KqujTJWUmEKWBkS+VO7JP8yhwzarCy9YSpJyT5xpzfdyV5dl85/8CtDoD9\nTfLfCwsJast1Ns13LMJkTPw+U1V1bpJzh1+/X1U3TrpPJu64JPdOuwi6Gccf8b5pF7Ao9Y7JjGO9\nY9xn5BC8F5euU0Zp1BOm7kxy0pzfTxy2/YjW2kVJLkqSqtrZWtvY0SdLgHFcGYzjymAclz9juPz1\n3Brh80meXFWnVtWPJTk7yRXjKQsAYHlY9MxUa+3hqvrNJFdl9tYI722t3TS2ygAAloGuNVOttY8m\n+egCDrmopz+WDOO4MhjHlcE4Ln/GcJmr1tq0awAAWLZ8nQwAQIexhamqOrOqvlJVX62qLQfY/+iq\n+uCw/2+rambOvq3D9q9U1S+OqyYWZrFjWFUzVfVgVV0/PN51uGvnH4wwjqdX1Req6uGqevV++86p\nqluHxzmHr2r21zmOe+e8H30waIpGGMffrqovVdUXq2p7VZ0yZ5/343LRWut+ZHYB+teSPDHJjyX5\n30l+dr82v5HkXcPzs5N8cHj+s0P7Ryc5dTjPEeOoy+OwjeFMkhun/Ro8Rh7HmSQ/l+T9SV49Z/vj\nktw2/Dx2eH7stF/Tanz0jOOw74FpvwaPkcfxjCQ/Pjz/9Tn/XvV+XEaPcc1MPfLVMq21/5dk31fL\nzHVWkkuG55cl2VRVNWy/tLX2/dba15N8dTgfh1fPGLJ0zDuOrbXbW2tfTPLD/Y79xSRXt9bua619\nO8nVSc48HEXzj/SMI0vHKOO4o7X298Ovn8vsPRsT78dlZVxh6kBfLfOEg7VprT2c5P8mWTfisUxe\nzxgmyalV9b+q6pNV9QuTLpaD6nk/eS8uHb1jsbaqdlbV56rqFeMtjQVY6Di+MclfL/JYpmjiXyfD\nqrA7ycmttT1V9U+TfKSqntJau3/ahcEqdUpr7c6qemKSa6rqhtba16ZdFAdXVb+cZGOS50+7FhZu\nXDNTo3y1zCNtqmpNkscm2TPisUzeosdwuES7J0laa9dldo3AholXzIH0vJ+8F5eOrrFord05/Lwt\nySeSPH2cxTGykcaxql6Y5M1J/kVr7fsLOZalYVxhapSvlrkiyb5PI7w6yTWttTZsP3v4pNipSZ6c\n5O/GVBejW/QYVtX6qjoiSYb/E35yZhdLcvj1fM3TVUleXFXHVtWxSV48bOPwW/Q4DuP36OH5cUme\nm+RLE6uUQ5l3HKvq6UnendkgdfecXd6Py8m4VrIneWmSWzI7K/HmYdt/yuw/IEmyNsl/y+wC879L\n8sQ5x755OO4rSV4y7VX5q/Wx2DFM8qokNyW5PskXkrx82q9lNT9GGMdnZnb9xfcyOzt805xjf2UY\n368m+dfTfi2r+bHYcUzy80luyOwnx25I8sZpv5bV/BhhHD+e5K7h35/XJ7lizrHej8vk4Q7oAAAd\n3AEdAKCDMAUA0EGYAgDoIEwBAHQQpgAAOghTAAAdhCkAgA7CFABAh/8PbZgxRTlbU90AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110c6af28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(a_samples, bins=abbins, label='A', normed=True)\n",
    "plt.hist(b_samples, bins=abbins, label='B', alpha=0.5, normed=True)\n",
    "plt.title('Posterior distributions')\n",
    "plt.xlim(0, max(abbins));\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "And we're finally in a position to answer the original question: how confident can we be that layout B is better?\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.659"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(b > a for a, b in zip(a_samples, b_samples))/len(a_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's it! Given the data and our prior beliefs, we are 66% sure that layout B is better than A.\n",
    "\n",
    "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.\n",
    "\n",
    "Approximate Bayesian Computation is slow and in some ways crude. There are often [better options](http://blog.fastforwardlabs.com/2017/01/30/the-algorithms-behind-probabilistic-programming.html). But it's not cheating. In the tutorials, I show how to speed it up for more realistic problems.\n",
    "\n",
    "## Next steps\n",
    "\n",
    "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. \n",
    "\n",
    "Log in to O'Reilly Safari and check out the three tutorials for a deeper video and text discussion, interactive code examples, exercises, and ideas for where to go next:\n",
    "\n",
    " - [Part 1: A/B Testing with Approximate Bayesian Computation](https://www.safaribooksonline.com/oriole/probabilistic-programming-from-scratch-1-a-b-testing-with-approximate-bayesian-computation)\n",
    " - [Part 2: Bayes' Theorem and Online Learning](https://www.safaribooksonline.com/oriole/probabilistic-programming-from-scratch-2-bayes-theorem-and-online-learning)\n",
    " - [Part 3: Improving Algorithmic Performance and PyMC3](https://www.safaribooksonline.com/oriole/probabilistic-programming-from-scratch-3-improving-algorithmic-performance-and-online-learning)\n",
    "\n",
    "The [Fast Forward Labs report](http://blog.fastforwardlabs.com/2017/01/18/new-research-on-probabilistic-programming.html), meanwhile, discusses the business case for probabilistic programming in detail, and compares modern industrial strength systems like Stan and PyMC3. Please [get in touch](http://www.fastforwardlabs.com/#contact) if you're interested in working with us."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}