Homework 2 solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer all questions and submit them either as an IPython notebook, LaTeX document, or Markdown document. Each question is worth 25 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pymc as pm\n",
    "import pylab as plt\n",
    "import pandas as pd\n",
    "\n",
    "# Set seed\n",
    "np.random.seed(10011)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1\n",
    "\n",
    "The goal of this problem is to investigate the role of the proposal distribution in a Metropolis-Hastings algorithm designed to simulate from the posterior distribution of the mixture parameter $\\delta$. \n",
    "\n",
    "1. Simulate 200 realizations from the mixture distribution:\n",
    "    $$y_i \\sim \\delta N(7, 0.5^2) + (1-\\delta) N(10, 0.5^2)$$\n",
    "    with $\\delta = 0.7$. Plot a histogram of these data. \n",
    "2. Implement an ABC procedure to simulate from the posterior distribution of $\\delta$, using your data from part (1). \n",
    "3. Implement a random walk M-H algorithm with proposal $\\delta^{\\prime} = \\delta^{(i)} + \\epsilon$ with $\\epsilon \\sim Unif(−1,1)$. \n",
    "4. Reparameterize the problem letting $U = \\log\\left[\\frac{\\delta}{1 - \\delta}\\right]$ and $u^{\\prime} = u^{(i)} + \\epsilon$. Implement a random walk chain in U-space. \n",
    "5. Compare the estimates and convergence behavior of the three algorithms.\n",
    "\n",
    "In part (1), you are asked to simulate data from a distribution with $\\delta$ known. For parts (2)–(4), assume $\\delta$ is unknown with prior $\\delta \\sim Unif( 0,1)$. For parts (2)–(4), provide an appropriate plot and a table summarizing the output of the algorithm. To facilitate comparisons, use the same number of iterations, random seed, starting values, and burn-in period for all implementations of the algorithm. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Part 1**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Sample data\n",
    "n = 200\n",
    "delta = 0.7\n",
    "\n",
    "y = np.where(np.random.random(n) < delta, np.random.normal(7, 0.5, n), np.random.normal(10, 0.5, n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 14.,  59.,  46.,  17.,   0.,   5.,  12.,  31.,  13.,   3.]),\n",
       " array([  5.82439521,   6.38620264,   6.94801007,   7.5098175 ,\n",
       "          8.07162492,   8.63343235,   9.19523978,   9.75704721,\n",
       "         10.31885463,  10.88066206,  11.44246949]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEBCAYAAABlki5mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADWtJREFUeJzt3X+s3Xddx/Hnax3t3OxWV5aWCH9s9Jb/DI4/xI250F6N\nNGA0xMw/DCgmbBadkNjUMON6EyMhRQtOJpGYoH/ghBAxJpXetjasGSAJi5uJQn+4uIWljW3TesG7\n3o2+/eN+l1zYvfece3e/5/SzPh/JTe75nu/t552lfe57P/d7zk1VIUlq13XjHkCS9OoYcklqnCGX\npMYZcklqnCGXpMYZcklq3MCQJ3ljkmNJjif50+7YZPf4eJId/Y8pSVpKBt1HnuQx4M+r6mvd4+uA\n48Bkd8oh4N7yhnRJGotlr8iTrAPe/HLEOxPAiaqarapZ4DSwrccZJUnLuH7A87cBNyT5MnAz8Ahw\nBriY5EB3ziVgM3CytyklSUsaFPLzzIf6vcA64Angt4BNwG4gwKPAuR5nlCQtY9mQV9WLSZ4DtlbV\nd5NcBk4B2xecNlFVp5b6M44ePereuSStws6dOzPMeYOuyAH2Ap9Ncgvwhar6vyRTwOHu+X2D/oA7\n77xzmFkkSZ0nn3xy6HMHhryqngV2/cixaWB6xZNJktbcMFfkWgNnZi5zdmZuJGtt2bierRs3jGQt\nSeNnyEfk7Mwcew4u+aOENbV/1zZDLl1DfIm+JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXO\nkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS\n4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDVuYMiTfC7J15McS/K+7thkkuPdx47+x5Qk\nLeX6Ic4p4L6qehYgyXXAFDDZPX8oybGqqp5mlCQtY9itlSz4fAI4UVWzVTULnAa2rflkkqShDHNF\nPgN8PskF4CPArcDFJAe65y8Bm4GT/YwoSVrOwJBX1YMASd4K7Af2ApuA3cxfqT8KnOtxRknSMlZy\n18oLwIvAKWD7guMTVXVqTaeSJA1t4BV5kseANwDfA3ZX1ZUkU8Dh7pR9/Y0nSRpkmK2VX1vk2DQw\n3ctEkqQV8QVBktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4\nQy5JjRvmF0uoMevXhaeenxnJWls2rmfrxg0jWUvS4gz5a9CF2ZeYOvLMSNbav2ubIZfGzK0VSWqc\nIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZek\nxhlySWqcIZekxg0V8iQbkvx3kg91jyeTHO8+dvQ7oiRpOcP+hqAHgG8BlSTAFDDZPXcoybGqqj4G\nlCQtb+AVeZIbgZ8H/hEIMAGcqKrZqpoFTgPbep1SkrSkYa7IHwT+AtjSPd4MXExyoHt8qTt2cu3H\nkyQNsuwVeZJbgHdU1VeYvxoHOA9sAj4KPNR9fq7PISVJSxt0RX43cEOSvwNu784/DmxfcM5EVZ3q\naT5J0gDLhryqDgIHAZK8H7ipqp5OMgUc7k7b1+uEkqRlDXvXClX1Nws+nwame5lIkrQiviBIkhpn\nyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWp\ncYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZc\nkhpnyCWpcYZckhp3/aATkvwxcBdwBfhgVf1Xkkng4e6Uh6vqX3qcUdIaODNzmbMzcyNZa8vG9Wzd\nuGEka2mIkFfVHwIkuRvYm+QBYAqY7E45lORYVVV/Y0p6tc7OzLHn4KmRrLV/1zZDPkIr2Vp5O/Cf\nwARwoqpmq2oWOA1s62M4SdJgA6/IAZI8DrweuAfYDlxMcqB7+hKwGTjZy4SSpGUNdUVeVT8H/Abw\nt8B5YBPwUeCh7vNzPc0nSRpgJVsrZ5i/gj/F/FX5yyaqajQbb5KkVxjmrpW/Z35bZQ74naq6kmQK\nONydsq+/8SRJgwxz18p9ixybBqZ7mUiStCK+IEiSGmfIJalxhlySGmfIJalxhlySGmfIJalxhlyS\nGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfI\nJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGjcw5Ek+k+RYkq8muaM7\nNpnkePexo/8xJUlLuX7QCVX1AEAX7D1JdgNTwGR3yqEkx6qq+htTkrSUlWytzABzwARwoqpmq2oW\nOA1s62M4SdJgA6/IF/gA8ClgM3AxyYHu+KXu2Mk1nk2SNIShQp7kPcB3qurbSbYDm4DdQIBHgXP9\njShJWs4wP+x8G3BvVX2yO3Qa2L7glImqOtXHcJKkwYa5Iv8i8FySY8DTVfV7SaaAw93z+/oaTpI0\n2DB3rdyxyLFpYLqXiSRJK+ILgiSpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhq3kjfN\nkl5h/brw1PMzI1lry8b1bN24YSRrSS0x5HpVLsy+xNSRZ0ay1v5d2wy5tAi3ViSpcYZckhpnyCWp\ncYZckhpnyCWpcYZckhpnyCWpcYZckhp3Tb8g6Pz353jhpRrJWldGs4yka9A1HfJvfXeGTzz+bO/r\nbPnx9Xz4HW/qfR1J1ya3ViSpcYZckhp3TW+tSOqH74o5WoZc0przXTFHy60VSWqcIZekxrm1Io3R\nmZnLnJ2ZG8lacz+4MpJ1NHqGXBqjszNz7Dl4aiRrPTx5+0jW0ei5tSJJjRsY8iT3JPlmkv0Ljk0m\nOd597Oh3REnScobZWtkAfAy4CyDJdcAUMNk9fyjJsary3UQkaQwGXpFX1RHgwoJDE8CJqpqtqlng\nNLCtp/kkSQOs5oedtwIXkxzoHl8CNgMn12wqSdLQVhPy88AmYDcQ4FHg3FoOJUka3rB3rWTB56eB\n7QseT1TVaO6fkiS9wsAr8iR7gXcBW5PcXFX3J5kCDnen7OtxPknSAANDXlUfBz7+I8emgem+hpIk\nDc8XBElS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXO\nkEtS4wy5JDXOkEtS41bzq96k17wzM5c5OzPX+zpzP7jS+xp67TPk0iLOzsyx52D/v8Hw4cnbe19D\nr31urUhS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXO2w8lNW39uvDU8zMjWWvLxvVs3bhhJGuthCGX\n1LQLsy8xdeSZkay1f9e2qzLkbq1IUuMMuSQ1zpBLUuMMuSQ1btUhTzKZ5Hj3sWMth5IkDW9Vd60k\nuQ6YAia7Q4eSHKuqWrPJJElDWe0V+QRwoqpmq2oWOA1sW7uxJEnDWu195LcCF5Mc6B5fAjYDJ9dk\nKknS0FYb8vPAJmA3EOBR4NxaDTUqb7ntRu7/mZ/sfZ2b1l9H0vsykq5RWc22dpJ1wOPM75EHOFxV\ndy927tGjR903l6RV2Llz51CXgKsKOUCSXwD+qHs4VVWHV/UHSZJelVWHXJJ0dfAFQZLUOEMuSY0z\n5JLUuN5CnuRzSb6e5FiS9/e1Tp+SvLGb/3iSPxv3PMNKcnM398sfl8Y900oleV+Sf03yRJJ3jnue\nlUpyf/f3/3CSiXHPM0iSe5J8M8n+BceaeRuOJeZ/xbGr1RLzf6b79/vVJHcs9/V9/mKJAu6rqmd7\nXKNvnwAeqqqvjXuQlaiq/wXeCZDkp4DfHe9Eq/L7wE8DNwGHgJ8d7zjDS3Ij8JtV9fYkrwf+EvjV\nMY81yAbgY8Bd0OTbcPzQ/Mscu1q9YtaqegCg+5/oHuC3l/rivrdWmn0ZTHev/Jtbi/giHgQeGfcQ\nq/AfwL3Au4FvjHmWlQrwuiQbgIvA1iSvG/NMy6qqI8CFBYeaehuOReZf9NjVasCsM8Dccl/f5xX5\nDPD5JBeAj1TVqR7X6sNtwA1JvgzcDDxSVf8w5plWJMlm4E1V9fS4Z1mFaeDDwHrg02OeZUWq6vtJ\n/gT4Z+b/HfwE86+E/p+xDrYyvg3H1eMDwKeWO6G3kFfVgwBJ3grsB36lr7V6cp75v7zvBdYBTyT5\nSnd10ooPAn817iFWqtsPfHdV/VL3+PEkR1r6b19VXwK+BJDkyapqKeLwGnkbjtYleQ/wnar69nLn\njeKulReAF0ewzpqqqheB54CtVTUHXB7zSCuS5HrmtyWa+i6is47uIiNJgB9j/mcuzUmyC/i3cc8x\npIVboaeB7QseTzTwXfViW7ktbe/+0KxJ3gbcW1WfHPSFvV2RJ3kMeAPz31p+qK91erYX+GySW4Av\ntHRFCPwy8E9VdWXcg6xUVZ1M8o0kB5m/2Ph0Vb0w7rlWIslfA28Bvgf8+pjHGSjJXuBdzO/n31xV\n9yeZAl5+6419YxtuCEvM/wfALy48Nt4pl7bY/MAXgeeSHAP+/eVdjkW//ur9IbQkaRi+IEiSGmfI\nJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalx/w+6SA561yKfVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1067e0908>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Part 2**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7.9231652131648911, 1.5150267222378637)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(y), np.std(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_samples = 100\n",
    "epsilons = [0.05, 0.01]\n",
    "\n",
    "# List for holding accepted samples\n",
    "samples = []\n",
    "    \n",
    "while len(samples) < n_samples:\n",
    "        \n",
    "    # Simulate from prior\n",
    "    d = np.random.random()\n",
    "    \n",
    "    # Simulate data\n",
    "    y_sim = np.where(np.random.random(n) < d, np.random.normal(7, 0.5, n), np.random.normal(10, 0.5, n))\n",
    "    \n",
    "    if (np.abs(np.mean(y_sim) - np.mean(y)) < epsilons[0]) & (np.abs(np.std(y_sim) - np.std(y)) < epsilons[1]):\n",
    "        \n",
    "        samples.append(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  1.,   4.,   8.,  13.,  23.,  16.,  19.,  10.,   4.,   2.]),\n",
       " array([ 0.57736545,  0.59805126,  0.61873708,  0.63942289,  0.66010871,\n",
       "         0.68079452,  0.70148034,  0.72216616,  0.74285197,  0.76353779,\n",
       "         0.7842236 ]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEBCAYAAABysL6vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC8hJREFUeJzt3V+onHV+x/HPV9Oku5BU1i1G+gd0o71bxL0oWBaphtKV\nLlRKsVcV9satgttCg/SmmqtFAnXLsiIstF6J24VakAr+I6x263bbCnpRWjUISsUUN2hPacxR8+vF\nGcshJDnnzHlmPPme1wsGPM/MPL9fHh/emfPLPDM1xggA/Vz2WU8AgMUQeICmBB6gKYEHaErgAZoS\neICmNgx8VT1SVcer6kdVde1s26NV9dJs+52LnyYAW7VnoweMMb6ZJFV1S5IjSf4oyUhyxxjjrcVO\nD4B5bWWJZiXJmXU/18RzAWBCWwn8N5I8MvvvlSSPVdWTVXVo+mkBsF21mY8qqKqvJ/nSGOM752y/\nIcn9Y4zbFzQ/AOa04Rp8VX0lyc1jjD89z90fJvnoQs99/vnnfdANwBxuvfXWbS+Dbxj4JD9M8nZV\nHU/y6hjjW1X1gyQHs7ZUc8/FnnzjjTdud44Au8rLL788yX428y6aa8+z7Y5JRgdgYVzoBNCUwAM0\nJfAATQk8QFMCD9CUwAM0JfAATQk8QFMCD9CUwAM0JfAATQk8QFMCD9CUwAM0tZnPg4cd4d2VMzm5\nsrq08a7avzcH9+9b2ngwNYHnknFyZTVHnnpjaeMdu+2QwHNJs0QD0JTAAzQl8ABNCTxAUwIP0JTA\nAzQl8ABNCTxAUwIP0JTAAzQl8ABNCTxAUwIP0JTAAzQl8ABNCTxAUwIP0JTAAzQl8ABNCTxAUwIP\n0NRFA19Vj1TV8ar6UVVdO9t2uKpenN1uWc40AdiqPRe7c4zxzSSZhfxIVd2d5GiSw7OHPF1Vx8cY\nY7HTBGCrNrtEs5JkNcl1SV4bY5weY5xOciLJoUVNDoD5XfQV/DrfSPKXSa5M8n5VPTTb/sFs2+sL\nmBsA27Bh4Kvq60n+Y4zx71V1fZIrktydpJI8nOS9xU4RgHls9I+sX0ly8xjjO7NNJ5Jcv+4h140x\n3ljU5ACY30av4H+Y5O2qOp7k1THGt6rqaJJnZ/c/sMjJATC/jd5Fc+15tj2T5JmFzQiASbjQCaAp\ngQdoSuABmhJ4gKYEHqCpzV7JCizQuytncnJldWnjXbV/bw7u37e08fhsCDzsACdXVnPkqeVdM3js\ntkMCvwtYogFoSuABmhJ4gKYEHqApgQdoSuABmhJ4gKYEHqApgQdoSuABmhJ4gKYEHqApgQdoSuAB\nmhJ4gKYEHqApgQdoSuABmhJ4gKYEHqApgQdoSuABmhJ4gKYEHqApgQdoas9nPQHYqfZeXnnlnZWl\njLX6ydmljMPuIvBwAadOf5yjz725lLHuP3zNUsZhd7FEA9CUwAM0tWHgq+qrVfXTqjq2btujVfVS\nVR2vqjsXO0UA5rGZNfh9Sb6d5KZ120aSO8YYby1kVgBs24av4McYzyU5dZ67avrpADCVedfgV5I8\nVlVPVtWhKScEwDTmepvkGOPeJKmqG5IcS3L7lJMCYPs2+wr+QssxHyb5aKK5ADChDV/BV9V9Sb6W\n5GBVHRhj3FVVP0hyMGtLNfcseI4AzGHDwI8xHkzy4Dnb7ljYjACYhAudAJoSeICmBB6gKYEHaErg\nAZoSeICmBB6gKYEHaErgAZoSeICmBB6gKYEHaErgAZoSeICmBB6gqbm+sg+S5N2VMzm5srq08VY/\nObu0saADgWduJ1dWc+SpN5Y23v2Hr1naWNCBJRqApgQeoCmBB2hK4AGaEniApgQeoCmBB2hK4AGa\nEniApgQeoCmBB2hK4AGaEniApgQeoCmBB2hK4AGaEniApgQeoCmBB2hK4AGaumjgq+qrVfXTqjq2\nbtvhqnpxdrtl8VMEYB57Nrh/X5JvJ7kpSarqsiRHkxye3f90VR0fY4zFTRGAeVz0FfwY47kkp9Zt\nui7Ja2OM02OM00lOJDm0wPkBMKeNXsGf6wtJ3q+qh2Y/f5DkyiSvTzorALZtq4H/WZIrktydpJI8\nnOS9qScFwPZtJvC17r9PJLl+3c/XjTHemHZKwKLtvbzyyjsrSxvvqv17c3D/vqWNx5qLBr6q7kvy\ntSQHq+rAGOOuqjqa5NnZQx5Y8PyABTh1+uMcfe7NpY137LZDAv8ZuGjgxxgPJnnwnG3PJHlmkZMC\nYPtc6ATQlMADNCXwAE0JPEBTAg/QlMADNCXwAE0JPEBTAg/QlMADNCXwAE0JPEBTAg/QlMADNLXV\nb3RiB3t35UxOrqwubbzVT84ubSxg6wS+kZMrqzny1PK+YOv+w9csbSxg6yzRADQl8ABNCTxAUwIP\n0JTAAzQl8ABNCTxAUwIP0JTAAzQl8ABNCTxAUwIP0JTAAzQl8ABNCTxAUwIP0JTAAzQl8ABNCTxA\nUwIP0JTAAzQ1d+Cr6tGqeqmqjlfVnVNOCoDt27ON544kd4wx3ppqMgBMZ7tLNDXJLACY3HYCv5Lk\nsap6sqoOTTUhAKYx9xLNGOPeJKmqG5IcS3L7VJMCYPumeBfNh0k+mmA/AExo7lfwVfV4kquztlRz\nz2QzAmAS21mi+YMpJwLAtFzoBNCUwAM0JfAATQk8QFMCD9DUdj6Lhk14d+VMTq6sLmWs1U/OLmUc\n2Kq9l1deeWdlKWNdtX9vDu7ft5SxdjqBX7CTK6s58tQbSxnr/sPXLGUc2KpTpz/O0efeXMpYx247\nJPAzlmgAmhJ4gKYEHqApgQdoSuABmhJ4gKYEHqApgQdoSuABmhJ4gKYEHqApgQdoSuABmhJ4gKYE\nHqApgQdoyhd+AK0s89ujkp39DVICD7SyzG+PSnb2N0hZogFoSuABmhJ4gKYEHqApgQdoSuABmtp1\nb5P84PTHGRlLGWvPZbWUcQDOZ9cF/q/+5Z3863/+91LGuuvXfykH9u26QwzsELuuPu+f/ij/9T8f\nLWWsDz8+mwM78/oHYBewBg/QlMADNDV34KvqcFW9OLvdMuWkANi+udbgq+qyJEeTHJ5terqqjo8x\nlvP2FAA2NO8r+OuSvDbGOD3GOJ3kRJJD000LgO2a9100X0jyflU9NPv5gyRXJnl9klkBsG3zBv5n\nSa5IcneSSvJwkvemmtQi/favfTFfvnr/Usa6/oufz/unP17KWADnqnmWzavq8iQvZG0NvpI8O8b4\njXMf9/zzz1uTB5jDrbfeuu1L4ecKfJJU1W8l+fPZj0fHGM9udzIATGfuwAOws7nQCaApgQdoSuAB\nmpor8Fv5mIKqerSqXqqq41V150bbLzVbPBa/PPvzvlhVfzHPPnayiY7FrjovqurA7M/66e2Dre5j\np5voWOyq82L22D+sqn+qqh9X1W/Os4+MMbZ0y9pfCj9O8rnZ7YXM/rH2Ao//6yS/utntl9JtjmPx\neJKbtrOPnXqb4ljs1vNi3fO+nOT7u/m8ON+x2K3nRZJXk1ye5ECSl+bZxzyv4Of5mIILvZ/zUv/K\no00fi9m1A18aY/zjvPvY4aY4Fv//kAXNcVnm/X96b5LvbnMfO80Ux+JTu+28+LckNyf5nSQ/mWcf\n81zJutWPKVhJ8lhVnUryJ2OMNzbYfinZyrH4xSQ/X1V/l7W/kb87xnhii/vYyaY4FsnuOy+SJFV1\nZZJfGWO8Ou8+dqgpjkWyO8+LZ5L8cZK9Sb431z7m+DXj+qz9uvS5JJ9P8miSQ5t43g1Jntjs9kvh\ntpVjkeTnkvxD1n7l2pvkn2fPm+t47rTbFMdiN54X657zZ0l+bzv72Im3KY7Fbjwvklyb5G/X/fzC\nPL2YZ4nmxGyQT103Nve36YdJzvddeRfafinY9LEYY3yU5O0kB8cYq0nObHUfO9wUx2K9XXFeJElV\n7cnar+FPrNu8686L5ILHYr3dcl5cntkKS1VV1oI+triPrS/RjDE+qaqjST79aIIHPr2vqn4/yf+O\nMf5+3bbHk1ydtV+x7tlo+6Vkq8ciyX1Jvl9Vv5Dkb8baGloutI9LyYTHYjeeF7+b5MkxxtnN7ONS\nMsWxmD12V50XY4zXq+onVfVU1v5h9XtjjA9nj930eeGjCgCacqETQFMCD9CUwAM0JfAATQk8QFMC\nD9CUwAM0JfAATf0fkdBiTxP/lQAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c01ea20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Part 3**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 14.,  59.,  46.,  17.,   0.,   5.,  12.,  31.,  13.,   3.]),\n",
       " array([  5.82439521,   6.38620264,   6.94801007,   7.5098175 ,\n",
       "          8.07162492,   8.63343235,   9.19523978,   9.75704721,\n",
       "         10.31885463,  10.88066206,  11.44246949]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEBCAYAAABlki5mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADWtJREFUeJzt3X+s3Xddx/Hnax3t3OxWV5aWCH9s9Jb/DI4/xI250F6N\nNGA0xMw/DCgmbBadkNjUMON6EyMhRQtOJpGYoH/ghBAxJpXetjasGSAJi5uJQn+4uIWljW3TesG7\n3o2+/eN+l1zYvfece3e/5/SzPh/JTe75nu/t552lfe57P/d7zk1VIUlq13XjHkCS9OoYcklqnCGX\npMYZcklqnCGXpMYZcklq3MCQJ3ljkmNJjif50+7YZPf4eJId/Y8pSVpKBt1HnuQx4M+r6mvd4+uA\n48Bkd8oh4N7yhnRJGotlr8iTrAPe/HLEOxPAiaqarapZ4DSwrccZJUnLuH7A87cBNyT5MnAz8Ahw\nBriY5EB3ziVgM3CytyklSUsaFPLzzIf6vcA64Angt4BNwG4gwKPAuR5nlCQtY9mQV9WLSZ4DtlbV\nd5NcBk4B2xecNlFVp5b6M44ePereuSStws6dOzPMeYOuyAH2Ap9Ncgvwhar6vyRTwOHu+X2D/oA7\n77xzmFkkSZ0nn3xy6HMHhryqngV2/cixaWB6xZNJktbcMFfkWgNnZi5zdmZuJGtt2bierRs3jGQt\nSeNnyEfk7Mwcew4u+aOENbV/1zZDLl1DfIm+JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXO\nkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS\n4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDVuYMiTfC7J15McS/K+7thkkuPdx47+x5Qk\nLeX6Ic4p4L6qehYgyXXAFDDZPX8oybGqqp5mlCQtY9itlSz4fAI4UVWzVTULnAa2rflkkqShDHNF\nPgN8PskF4CPArcDFJAe65y8Bm4GT/YwoSVrOwJBX1YMASd4K7Af2ApuA3cxfqT8KnOtxRknSMlZy\n18oLwIvAKWD7guMTVXVqTaeSJA1t4BV5kseANwDfA3ZX1ZUkU8Dh7pR9/Y0nSRpkmK2VX1vk2DQw\n3ctEkqQV8QVBktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4Qy5JjTPkktQ4\nQy5JjRvmF0uoMevXhaeenxnJWls2rmfrxg0jWUvS4gz5a9CF2ZeYOvLMSNbav2ubIZfGzK0VSWqc\nIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZekxhlySWqcIZek\nxhlySWqcIZekxg0V8iQbkvx3kg91jyeTHO8+dvQ7oiRpOcP+hqAHgG8BlSTAFDDZPXcoybGqqj4G\nlCQtb+AVeZIbgZ8H/hEIMAGcqKrZqpoFTgPbep1SkrSkYa7IHwT+AtjSPd4MXExyoHt8qTt2cu3H\nkyQNsuwVeZJbgHdU1VeYvxoHOA9sAj4KPNR9fq7PISVJSxt0RX43cEOSvwNu784/DmxfcM5EVZ3q\naT5J0gDLhryqDgIHAZK8H7ipqp5OMgUc7k7b1+uEkqRlDXvXClX1Nws+nwame5lIkrQiviBIkhpn\nyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWp\ncYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZc\nkhpnyCWpcYZckhp3/aATkvwxcBdwBfhgVf1Xkkng4e6Uh6vqX3qcUdIaODNzmbMzcyNZa8vG9Wzd\nuGEka2mIkFfVHwIkuRvYm+QBYAqY7E45lORYVVV/Y0p6tc7OzLHn4KmRrLV/1zZDPkIr2Vp5O/Cf\nwARwoqpmq2oWOA1s62M4SdJgA6/IAZI8DrweuAfYDlxMcqB7+hKwGTjZy4SSpGUNdUVeVT8H/Abw\nt8B5YBPwUeCh7vNzPc0nSRpgJVsrZ5i/gj/F/FX5yyaqajQbb5KkVxjmrpW/Z35bZQ74naq6kmQK\nONydsq+/8SRJgwxz18p9ixybBqZ7mUiStCK+IEiSGmfIJalxhlySGmfIJalxhlySGmfIJalxhlyS\nGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfI\nJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalxhlySGjcw5Ek+k+RYkq8muaM7\nNpnkePexo/8xJUlLuX7QCVX1AEAX7D1JdgNTwGR3yqEkx6qq+htTkrSUlWytzABzwARwoqpmq2oW\nOA1s62M4SdJgA6/IF/gA8ClgM3AxyYHu+KXu2Mk1nk2SNIShQp7kPcB3qurbSbYDm4DdQIBHgXP9\njShJWs4wP+x8G3BvVX2yO3Qa2L7glImqOtXHcJKkwYa5Iv8i8FySY8DTVfV7SaaAw93z+/oaTpI0\n2DB3rdyxyLFpYLqXiSRJK+ILgiSpcYZckhpnyCWpcYZckhpnyCWpcYZckhpnyCWpcYZckhq3kjfN\nkl5h/brw1PMzI1lry8b1bN24YSRrSS0x5HpVLsy+xNSRZ0ay1v5d2wy5tAi3ViSpcYZckhpnyCWp\ncYZckhpnyCWpcYZckhpnyCWpcYZckhp3Tb8g6Pz353jhpRrJWldGs4yka9A1HfJvfXeGTzz+bO/r\nbPnx9Xz4HW/qfR1J1ya3ViSpcYZckhp3TW+tSOqH74o5WoZc0przXTFHy60VSWqcIZekxrm1Io3R\nmZnLnJ2ZG8lacz+4MpJ1NHqGXBqjszNz7Dl4aiRrPTx5+0jW0ei5tSJJjRsY8iT3JPlmkv0Ljk0m\nOd597Oh3REnScobZWtkAfAy4CyDJdcAUMNk9fyjJsary3UQkaQwGXpFX1RHgwoJDE8CJqpqtqlng\nNLCtp/kkSQOs5oedtwIXkxzoHl8CNgMn12wqSdLQVhPy88AmYDcQ4FHg3FoOJUka3rB3rWTB56eB\n7QseT1TVaO6fkiS9wsAr8iR7gXcBW5PcXFX3J5kCDnen7OtxPknSAANDXlUfBz7+I8emgem+hpIk\nDc8XBElS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXO\nkEtS4wy5JDXOkEtS41bzq96k17wzM5c5OzPX+zpzP7jS+xp67TPk0iLOzsyx52D/v8Hw4cnbe19D\nr31urUhS4wy5JDXOkEtS4wy5JDXOkEtS4wy5JDXO2w8lNW39uvDU8zMjWWvLxvVs3bhhJGuthCGX\n1LQLsy8xdeSZkay1f9e2qzLkbq1IUuMMuSQ1zpBLUuMMuSQ1btUhTzKZ5Hj3sWMth5IkDW9Vd60k\nuQ6YAia7Q4eSHKuqWrPJJElDWe0V+QRwoqpmq2oWOA1sW7uxJEnDWu195LcCF5Mc6B5fAjYDJ9dk\nKknS0FYb8vPAJmA3EOBR4NxaDTUqb7ntRu7/mZ/sfZ2b1l9H0vsykq5RWc22dpJ1wOPM75EHOFxV\ndy927tGjR903l6RV2Llz51CXgKsKOUCSXwD+qHs4VVWHV/UHSZJelVWHXJJ0dfAFQZLUOEMuSY0z\n5JLUuN5CnuRzSb6e5FiS9/e1Tp+SvLGb/3iSPxv3PMNKcnM398sfl8Y900oleV+Sf03yRJJ3jnue\nlUpyf/f3/3CSiXHPM0iSe5J8M8n+BceaeRuOJeZ/xbGr1RLzf6b79/vVJHcs9/V9/mKJAu6rqmd7\nXKNvnwAeqqqvjXuQlaiq/wXeCZDkp4DfHe9Eq/L7wE8DNwGHgJ8d7zjDS3Ij8JtV9fYkrwf+EvjV\nMY81yAbgY8Bd0OTbcPzQ/Mscu1q9YtaqegCg+5/oHuC3l/rivrdWmn0ZTHev/Jtbi/giHgQeGfcQ\nq/AfwL3Au4FvjHmWlQrwuiQbgIvA1iSvG/NMy6qqI8CFBYeaehuOReZf9NjVasCsM8Dccl/f5xX5\nDPD5JBeAj1TVqR7X6sNtwA1JvgzcDDxSVf8w5plWJMlm4E1V9fS4Z1mFaeDDwHrg02OeZUWq6vtJ\n/gT4Z+b/HfwE86+E/p+xDrYyvg3H1eMDwKeWO6G3kFfVgwBJ3grsB36lr7V6cp75v7zvBdYBTyT5\nSnd10ooPAn817iFWqtsPfHdV/VL3+PEkR1r6b19VXwK+BJDkyapqKeLwGnkbjtYleQ/wnar69nLn\njeKulReAF0ewzpqqqheB54CtVTUHXB7zSCuS5HrmtyWa+i6is47uIiNJgB9j/mcuzUmyC/i3cc8x\npIVboaeB7QseTzTwXfViW7ktbe/+0KxJ3gbcW1WfHPSFvV2RJ3kMeAPz31p+qK91erYX+GySW4Av\ntHRFCPwy8E9VdWXcg6xUVZ1M8o0kB5m/2Ph0Vb0w7rlWIslfA28Bvgf8+pjHGSjJXuBdzO/n31xV\n9yeZAl5+6419YxtuCEvM/wfALy48Nt4pl7bY/MAXgeeSHAP+/eVdjkW//ur9IbQkaRi+IEiSGmfI\nJalxhlySGmfIJalxhlySGmfIJalxhlySGmfIJalx/w+6SA561yKfVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c074160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several ways of solving this; here is the easiest (though not the most general)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "dnorm = norm.pdf\n",
    "\n",
    "def calc_posterior(d):\n",
    "    return np.log(d*dnorm(y, 7, 0.5) + (1-d)*dnorm(y, 10, 0.5)).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0\n",
      "Iteration 1000\n",
      "Iteration 2000\n",
      "Iteration 3000\n",
      "Iteration 4000\n",
      "Iteration 5000\n",
      "Iteration 6000\n",
      "Iteration 7000\n",
      "Iteration 8000\n",
      "Iteration 9000\n"
     ]
    }
   ],
   "source": [
    "n_iterations = 10000\n",
    "\n",
    "# Initialize trace for parameters\n",
    "trace = np.empty(n_iterations+1)\n",
    "\n",
    "# Set initial values\n",
    "trace[0] = 0.5\n",
    "# Calculate joint posterior for initial values\n",
    "current_log_prob = calc_posterior(trace[0])\n",
    "\n",
    "# Initialize acceptance counts\n",
    "accepted = 0\n",
    "\n",
    "for i in range(n_iterations):\n",
    "\n",
    "    if not i%1000: print('Iteration %i' % i)\n",
    "\n",
    "    # Grab current parameter values\n",
    "    d_current = trace[i]\n",
    "\n",
    "    # Propose new value for d\n",
    "    d_prop = d_current + np.random.uniform(-1, 1)\n",
    "    while (d_prop<0) or (d_prop>1):\n",
    "        d_prop = d_current + np.random.uniform(-1, 1)\n",
    "    \n",
    "    # Calculate log posterior with proposed value\n",
    "    proposed_log_prob = calc_posterior(d_prop)\n",
    "\n",
    "    # Log-acceptance rate\n",
    "    alpha = proposed_log_prob - current_log_prob\n",
    "\n",
    "    # Test proposed value\n",
    "    if np.log(np.random.random()) < alpha:\n",
    "        # Accept\n",
    "        trace[i+1] = d_prop\n",
    "        current_log_prob = proposed_log_prob\n",
    "        accepted += 1\n",
    "    else:\n",
    "        # Reject\n",
    "        trace[i+1] = trace[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  3.00000000e+00,   0.00000000e+00,   2.50000000e+01,\n",
       "          1.53000000e+02,   7.48000000e+02,   2.51400000e+03,\n",
       "          3.01300000e+03,   2.70600000e+03,   7.11000000e+02,\n",
       "          1.28000000e+02]),\n",
       " array([ 0.5       ,  0.52758049,  0.55516098,  0.58274147,  0.61032196,\n",
       "         0.63790245,  0.66548293,  0.69306342,  0.72064391,  0.7482244 ,\n",
       "         0.77580489]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEBCAYAAAB4wNK4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADqtJREFUeJzt3VGIXNd9x/HvT+vKTcGusFKiQFuIqtVjMQmEkOCYyksh\naQM1pShPNuTFqRycGmJModRaCAQjqAghIpCHmDwYu4GmUGrQStYSq6nTFALxQ5NIuxgcCBbIQmJL\n115b+vdhjspYHc3OzI5mruTvBy7M/O+Zc+7Za+3P587emVQVkiTtmvcBSJK6wUCQJAEGgiSpMRAk\nSYCBIElqDARJErBNICT5epIzSU4n2d9qzyV5Nclqkkf72i4lOdu2Q9vVJUndklHuQ0jyGeCRqnos\nyfeAZ6rqjb79u4CzwFIrnayqzw6qAw+WNz9IUufcNWK7TwG/6HueG/YvAueqahMgyXqSRXorkPfV\ngQPA+R0dtSRp6rYNhCSvAB8GHmilDeD5JJeAJ6tqDbgPuJzkeGtzBdhLLzgG1Q0ESeqYbQOhXfr5\nJPB94M+q6gmAJPcDx4CHgbeAPcAReiFwArhIb4UwqC5J6phRLxm9OaDt28C77fE6cLBv32JVrSVZ\nGFS/2SAvv/yy7y1I0gQeeuihGy/lj21oICR5kd7loi3gK632AvBRepeOHgeoqqtJloFT7aVHh9WH\n+fjHPz7uHCTpA+1nP/vZVPoZGghVdXhA7Ys3absCrIxalyR1izemSZIAA0GS1BgIkiTAQJAkNQaC\nJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANB\nktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUjM0EJJ8PcmZJKeT7G+1pSRn23aor+1YdUlSt9w1\nbGdV/R1Aks8ATyf5MrAMLLUmJ4EzSXaNWk+yWlU13WlIknZqaCD0+RTwC2AROFdVmwBJ1pMs0ltp\njFQHDgDnpzsNSdJObRsISV4BPgw8ABwELic53nZfAfYCGbNuIEhSx2wbCFX12SSfBL4PPAnsAY7Q\n+2V/ArhIbyUwTl36QHpz4x0ubGzNZeyP3LObfffcPZexdXsY9ZLRm63tGr1VwnWLVbWWZGGc+o6O\nWLqNXdjY4qmX5vNP4NjnDxgIGmpoICR5kd7loi3gK1V1LckycKo1OQpQVVfHqUuSume7vzI6PKC2\nAqzstC5J6hZvTJMkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmA\ngSBJagwESRJgIEiSGgNBkgSM/o1pkm5zuxfCz3+zMZex/frO24OBIH1AXNp8j+XTr89lbL++8/bg\nJSNJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJwDaBkOQ7SVaT/CjJ/lZ7Lsmrrf5oX9ulJGfb\ndmi7uiSpW4bemFZVXwZov8ifAv4aKOBwVb1xvV2SXcAysNRKJ4Ezg+pJVquqpjoLSdKOjXqn8gbw\nTt/z3LB/EThXVZsASdaTLNJbgbyvDhwAzu/oqCVJUzdqIHwJ+GZ7vAE8n+QS8GRVrQH3AZeTHG9t\nrgB76QXHoLqBoLl5c+MdLmxszWXsravX5jKuNIptAyHJF4BfVdUvAarqiVa/HzgGPAy8BewBjtAL\ngRPARXorhEF1aW4ubGzx1Etrcxn7maWPzWVcaRRDAyHJJ4AHq+prA3a/DbzbHq8DB/v2LVbVWpKF\nQfWdHLAk6dbYboXwA+DXSVaB16rqq0leBPbRu3T0OEBVXU2yDJxqrzs6rC5J6p7t/spo/4Da4Zu0\nXQFWRq1LkrrFG9MkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkS\nYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSp\nGRoISb6TZDXJj5Lsb7WlJGfbdqiv7Vh1SVK33DVsZ1V9GaD9In8qyRFgGVhqTU4CZ5LsGrWeZLWq\narrTkCTt1NBA6LMBbAGLwLmq2gRIsp5kkd5KY6Q6cAA4P91pSJJ2atRA+BLwTWAvcDnJ8Va/0moZ\ns24gSFLHbBsISb4A/KqqfpnkILAHOELvl/0J4CK9lcA4dUlSxwwNhCSfAB6sqq+10jpwsK/JYlWt\nJVkYpz6NA5ckTdd2K4QfAL9Osgq8VlVfTbIMnGr7jwJU1dVx6pKk7tnur4z2D6itACs7rUuSusUb\n0yRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQY\nCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkZGghJHkjy0yTH\n+mrPJXk1yWqSR/vqS0nOtu3QdnVJUrfctc3+u4FvAJ/uqxVwuKreuF5IsgtYBpZa6SRwZlA9yWpV\n1TQOXpI0PUNXCFV1Grg0YFdueL4InKuqzaraBNaTLA6qAwemcNySpCnbboUwyAbwfJJLwJNVtQbc\nB1xOcry1uQLspRccg+rnd3bYkqRpGzsQquoJgCT3A8eAh4G3gD3AEXohcAK4SG8FMqguSeqYUQLh\nxstD170NvNserwMH+/YtVtVakoVB9fEPU5J0qw0NhCRPA58D9iW5t6oeS/IisI/epaPHAarqapJl\n4FR76dFhdUlS9wwNhKp6Fnj2htrhm7RdAVZGrUuSusUb0yRJgIEgSWoMBEkSYCBIkhoDQZIEGAiS\npMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJ\nEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkZGghJHkjy0yTH+mpLSc627dCkdUlSt9y1zf67gW8AnwZI\nsgtYBpba/pPAmXHqSVarqqY3BUnSNAwNhKo6neTBvtIicK6qNgGSrCdZpLfSGKkOHADOT38qkqSd\n2G6FcKP7gMtJjrfnV4C9QMasGwiS1DHjBsJbwB7gCL1f9ieAi/RWAuPUJUkdM0ogpO/xOnCw7/li\nVa0lWRinPvnhSpJulaGBkORp4HPAviT3VtVjSZaBU63JUYCqujpOXZLUPdu9qfws8OwNtRVgZUDb\nseqSpG7xxjRJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiS\npGbc70OQpubNjXe4sLE183G3rl6b+ZjS7cBA0Nxc2NjiqZdm//UYzyx9bOZjSrcDA0HSLbd7Ifz8\nNxszH/cj9+xm3z13z3zc25WBIOmWu7T5HsunX5/5uMc+f8BAGINvKkuSAANBktQYCJIkwECQJDUG\ngiQJMBAkSY2BIEkCDARJUjNRICR5LsmrSVaTPNJqS0nOtu1QX9uBdUlSt0x6p3IBh6vqDYAku4Bl\nYKntPwmcGVRPslpVtYNjliTdAjv56Ir0PV4EzlXVJkCS9SSL9FYg76sDB4DzOxhXknQLTBoIG8Dz\nSS4BTwL3AZeTHG/7rwB76YXGoLqBIEkdM1EgVNUTAEnuB44BTwN7gCP0QuAEcJHeCmFQXZLUMTv9\ntNO3gXeBNeBgX32xqtaSLAyq73BMSdItMFEgJHkB+Cjw38CRqrqWZBk41ZocBaiqq4PqkqTumfSS\n0RcH1FaAlVHrkqRu8cY0SRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIa\nA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkpqJvkJTd443N97hwsbWXMbeunptLuNK\nGsxA+IC7sLHFUy+tzWXsZ5Y+NpdxJQ1mIEi6Y+1eCD//zcZcxv7IPbvZd8/dcxl7UgaCpDvWpc33\nWD79+lzGPvb5A7ddIPimsiQJmGEgJFlKcrZth2Y1riRpNDO5ZJRkF7AMLLXSySSrVVWzGF+StL1Z\nrRAWgXNVtVlVm8A6cGBGY0uSRjCrN5XvAy4nOd6eXwH2AudnNH6neS+ApC6YVSC8BewBjgABTgAX\nBzW8dm0OV5EC710rdi/M5z127wWQ1AWZxWX8JAvAK/TeQwhwqqo+c2O7l19+2fcUJGkCDz30UHba\nx0wCASDJnwJ/354uV9WpmQwsSRrJzAJBktRt3pgmSQIMBElSYyBIkoAZBcI4H1uR5LkkryZZTfLo\nJH3M2pTmN7DeBWPO7/fbHM4m+YdJ+pi1Kc2vk+dv1Lklubcd+/Xtyrh9zMOU5tfJcwdj/7f5SJL/\nSPLjJH8ySR9U1S3d6IXOj4EPte0V2pvZN2n/PeAPd9LHLLdpzG9Yfd7bBPN7Afj0HXz+/t/8unr+\nJv25A38MfPdOO3eD5tfVczfJ/IDXgAXgXuDVSfqYxQphko+tuPHvabv80RfTmN929XkaeX7tfpM/\nqqp/n7SPOZjG/P6vyS06xklN+nN/AvjWDvuYhWnM77qunTsYf37/BTwI/Dnwk0n6mMWdyuN+bMUG\n8HySS8CTVbU2QR+zNI35DavP2zjz+z3gt5P8M73/S/lWVf1wzD5mbRrzg26ev7F/7kn2An9QVa9N\n2scMTWN+0M1zB+PPbwX4G2A38O2J+pjBsucgvSXZh4DfAZ4DDozwuvuBH+6kjxkt63Y8v1Hqt8P8\ngN8C/o3esnU38J/tdXfE+bvZ/Lp6/ib5uQN/C/zlTvq4nebX1XM37vyA/cA/9T1/ZZJ/e7O4ZLTe\nDuq6xRotgd8G3t1hH7MwjfmNUp+XkedXVe8Cvwb2VdUW8M64fczBNObXr0vnb6yfe5K76F1u+GFf\n+Y44d3DT+fXr0rmD8ea3QLvikyT0AqDG7OPWXzKqqqtJloHrH1Vx9Pq+JH8F/E9V/Wtf7QXgo/SW\ncY9v18e8TWN+w+rzNu78gKeB7yb5XeAfq3fdkjvl/HHz+XXu/E0wt78A/qWqro3Sx7xNY36tbefO\nHYw3v6o6n+QnSV6i90byt6vq7dZ25PPnR1dIkgBvTJMkNQaCJAkwECRJjYEgSQIMBElSYyBIkgAD\nQZLUGAiSJAD+F/OCscf2/1tRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c304978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0\n",
      "Iteration 1000\n",
      "Iteration 2000\n",
      "Iteration 3000\n",
      "Iteration 4000\n",
      "Iteration 5000\n",
      "Iteration 6000\n",
      "Iteration 7000\n",
      "Iteration 8000\n",
      "Iteration 9000\n"
     ]
    }
   ],
   "source": [
    "logit = lambda p: np.log(p/(1-p))\n",
    "invlogit = lambda x: 1./(1 + np.exp(-x))\n",
    "\n",
    "n_iterations = 10000\n",
    "\n",
    "# Initialize trace for parameters\n",
    "trace = np.empty(n_iterations+1)\n",
    "\n",
    "# Set initial values\n",
    "trace[0] = 0.5\n",
    "# Calculate joint posterior for initial values\n",
    "current_log_prob = calc_posterior(trace[0])\n",
    "\n",
    "# Initialize acceptance counts\n",
    "accepted = 0\n",
    "\n",
    "for i in range(n_iterations):\n",
    "\n",
    "    if not i%1000: print('Iteration %i' % i)\n",
    "\n",
    "    # Grab current parameter values\n",
    "    d_current = trace[i]\n",
    "\n",
    "    # Propose new value for d\n",
    "    d_prop = invlogit(logit(d_current) + np.random.uniform(-1, 1))\n",
    "    \n",
    "    # Calculate log posterior with proposed value\n",
    "    proposed_log_prob = calc_posterior(d_prop)\n",
    "\n",
    "    # Log-acceptance rate\n",
    "    alpha = proposed_log_prob - current_log_prob\n",
    "\n",
    "    # Test proposed value\n",
    "    if np.log(np.random.random()) < alpha:\n",
    "        # Accept\n",
    "        trace[i+1] = d_prop\n",
    "        current_log_prob = proposed_log_prob\n",
    "        accepted += 1\n",
    "    else:\n",
    "        # Reject\n",
    "        trace[i+1] = trace[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10c510f28>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEBCAYAAACT92m7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4HcW5/7+jakm2JFvuvcmVYiB0CGCMKTGhB8hNQg1J\nCJDkl0ISklycyg1JfJMAySV0EgyYGlNiG7DB2GCaKwZsyb03SZatLs3vj3P2aM/uzO5s33PO+3ke\nPfbZs2d2dnZ23pm3DeOcgyAIgiD05EVdAYIgCCJ+kHAgCIIgTJBwIAiCIEyQcCAIgiBMkHAgCIIg\nTJBwIAiCIEzYCgfG2DTG2OLk31Sbc7/GGFvGGFvCGDvLTRkEQRBE9DCrOAfGWB6AxQCmJQ/NA3AG\nl/yIMbYKwDEAygDM45yf7LQMgiAIInrsVg7VANZxzps5580AagGMtTh/LYAzAMwA8K7LMgiCIIiI\nKbD5vg+AesbYrOTnBgBVANZLzp8P4LsAigDc67IMgiAIImLshMN+AJUAbgbAANwHYJ/oRMbYaAAz\nOOdfTH5+izH2mpMyCIIgiHhgJxxqAYzTfa7mnNdIzs3XymOMMQAlALjDMvD666+TLYIgCMIFZ599\nNvOrLEvhwDnvZIzNBLAgeehO7TvG2BUAmjjnLyfPXc8Ye5cx9goStox7OectyXOFZcg49thjHd4G\nQRBEbvPRRx/5Wp7dygGc8/lI2BKMx+cIjv3WSRkEQRBEPKEgOIIgCMIECQeCIAjCBAkHgiAIwgQJ\nB4IgCMIECQeCIAjCBAkHgiAIwgQJB4IgMo72zi40tnZEXY2shoQDQRAZxz1Lt+Gyx1dHXY2shoQD\nQRAZx67GtqirkPWQcCAIH+Gco62jK+pqEIRnSDgQhI+8vakBMx5ZGXU1CMIzJBwIwkd2N7ZGXQWC\n8AUSDkRWwjnHUyt3h3/d0K9IEMFAwoHISlo6uvDg+zuiroYym+qa0dzeGXU1CCJFzgmHlo4u3LVw\nU9TVIIg0bnr2Uzz8wU7Xv29s7cD9y7b7WCMi18k54bDzYCveqK2LuhqxoKmNZqp+40Wt1OrBy2nl\nzkN4ZvUeD1ePjufX7EFDCwW0xY2cEw5EgoaWDlz82CpwTlpyIlr+9u52vLmBJmxxg4RDjtLWSb74\nRG5w5b9W46PtB6OuRsZBwiFH0RYMtG7IDnzbVT4iglzA1jV34OPdh4O7QJZCwiFHSQkHkg7+ElF7\n0mO0hvq5c0g45ChdPg0nz6/Zgzmrwo8nILKLoMfuLpIOjiHhkKP49a7cv2w7/vGefTzBv9fuxaa6\nZn8u6oDOrtwYFDJdrRQ0H21vjLoKGQcJhxwlbC+le5Zuw8uf7A/tetrtvVZzILRr6mnr6MK6fU2R\nXJswc5jcth1DwiFH6YrEIB3+LL65PRqvrOc+3oNbXvjM0W8YTf+JGEHCIUfRZtYvrNnjSR/LYjKi\nrdjRiN0xyPGvtWSLj0LppU/2YfmOcNUih9s6sf9wu6vfrtvXhA+3ketopkPCIUfRDNL3v7cDh1oz\nf8n9o1dq8L9vb/G93I4ujgXrnavDWlxEO8vE7F+WbMU/Qk6N8evXN+Lq2Wtc/faXr23AT/5T63ON\nvJEblid/IeGQxWzY34yf/qdG+F2mO2+s29eE/U3pM1vRCsirbaVmXxPuflNd6PDkMJQJMYZLNtXj\nnqVbhd/VNbtbNbiFIvXjBwmHLOb9bQfxwTaxOsKvd9GJUsnP1/+WFz7D7xdtSjumH5B/80b6dyq8\nt7UBVz+RPlvOz3OnNnP5M9e40e69uHYv/r12X+pzc3snGls7cN6Dy2F08mrt6CKjbo5BwiFCOOeR\npWnmuqG6oaUDd7+5OZJ6eME4gHXqJN77LnTeq3cdNq1GChyO8iwpLsOeB3sR9r99YyM+2XMYFz26\nCnVNHeji5vJ+/fpGXP74Km+VJDIKJeHAGJvGGFuc/JtqcV45Y2yh7q9B990jjLF3ksev8aPymc5L\nn+zDRY9G88LpB9aPdx/GgvXWLp+zV+zCWxsFydEEY2d7Z5dQTeBmMv39l9ZL4yOMs+UgYhryHU7J\nNaGbSVqSRRvqU2qkVP0N5yzbehCdGXRPblhYewB1TWZ1WmcXR3sIesKwnQ7ssBUOjLE8ADMBTE/+\n3ckkLiqc84Oc87M452cB+A6Ap/VfA7gy+f2j3que+Wysawm0fNFD+nDbQdy7dJvjwevhD3bi8Q93\nmb8QlPOFh1fiuTV7VU61ZfWuQ1ix45DwO+P9dXKOuqZ2fPO5T11cSYxbtZJfawdVgedGrWTVB+IQ\nUfz7RZvwzuYG+xN94ncLN+OFj8399o+Lt+DLsz8O9Nqcc9z+itg+GBUqK4dqAOs4582c82YAtQDG\nKvzuNgB/NRwLRBPb2cWx77B3N8awo2nDmI0YmfvJPry4dq/w5d/roQ3vfnNzWubLzXUtOCCYhblB\nbqxM706b6lqwsa4ZGw6IVxqNrR246dlPpNcRdc78PLs6pOMloSEz1KCLc5z/0ArTee2dXfh0T7CJ\n5LzIBuN9uOW1mjrMX+dP4KTq/YhOq9nXlJP7TagIhz4A6hljsxhjswA0AKiy+gFjrArAMM65XmfS\nCOAJxthcxpiKcFHm2dV7hJL9ofd3KKfq3VTXLHwRgyQv6BgBi+L1L4FWjdnLXeRISv52wfoDaaqp\nPYfbcJXBuOv33RqbL48xy0FgZ2MbNjlcrWkDnVPXVDeD60uf7kv7LCvjtZo63Pbvdc4vYEDU/bjh\nXzfwECwuVinnP9vrXnBGv16KDyrCYT+ASgA/BXBH8v/7LH8B3ATgfv0BzvltnPNTAfwcwN3Oqwrs\nP9yOTwQzpnqJVH9y5W6hekOEcZb793e3ZXxengcsch6JBp5WjysZfZkibYzfrWm8Rh4L7uV2OnP0\nWyuz12FA2vQHljvas8M421ddKU1/YLkvq3anzSVLwc05x60vygXn3sNtmP7AcouKZPY77ycqwqEW\nwDjd52rOuVQ5xhgrADADwPOSU1oAuNI3/OGtzfiODzMmEcY+8dyavVKh44aWji5c8li68Vl14dDY\n2oFrn05fGf177V4canVev9TM0M1L4GLqv7+pHWslL/LP59X66t9uW5TN9yI1XxizYBUaWjrQlly9\nFOWrPQi7yY21zUG5aiYPLzu10iMf7MCtL6qlFvG79euac0895BZb4cA570TCIL0AwHwAd2rfMcau\nYIx9wfCTiwHM5ZynvWmMsScZY28C+AOAH3qst+8EPQQ0tnaY/MTzFEfbXY1t2HEwfXZ2z9JtWLzJ\nvbFO/3C0QcKN2sfuN/+7eAu+O3ed8NxlWw86GoTM1zZf3Ysh9VnDKvP6OWvR1in23rEjCKGiPbOK\nHgUA3EVhqyKr//QHlntOKLhs60F8tjeapIR2/TUeU4F4UKByEud8PhKCwXh8juDYM5IyrnJcu+7f\nor2L++KRIStDuLm7y57S0cXT/OPbOyUBRIr346d6S7uk0NVUV59VOw+hrCgPY6pKTb+VIaold/C9\nk3IBtZWXl5bb1tDqOrVIkINMWVE+gPR+4ZcRWOsW+i5n7Cvr9jZhXN9SuMWJ/HYulLVrcGHer7CD\nE1WJo1DKiCC4RRvqMePhlYEacP3cU/mCh1akqSjuWboNNz3r3r2yU+FtunfpVnzzObknjhE7efOD\nl9fj7+8Gl8+nPulX7ynpn+EzB3C4zd/ZdGoG7bCamaC6tnqd9PW3Wyn4ca8dnRwvfbIPvxcEY3LO\nfY/OrpekB9Hu5WAOeicZyQjhsONgK4DETE6EoxQO0o7sr+DRD767JNlCVWcxVnJLcz/9aHsjNhyQ\ne+IcbuvEzsZW3czK/rodXVy8otLhttV+/fom5Xo4waiO15fvVdXDOY+Vk4JXe43Vz/Vt5WeGWRkP\nvL8Df1myFa8JgjEXb6o32es0VIWG8VbnrRMHfXIAq3Y24vJ/rk4dc+rhli1khHDQ0ISEU/YcajO9\nSJxzbK5rxr/XpuuZ65rb8d7WhC7fL72xLPJRNrByztPST8tWDgzAf83+GHsO2XuL/OO97bjpme6V\nhX7Grs0gjaqJj3cfxq9e32gqa+OB5pRx1A3vbW1I6cu9tLBo5uv32K1Ph/Hkyt3K7s5WA29zeyee\n9mFrVae3uvFAs3KAYJxWPlapw41CQx5non5DB1vimUOqdn8Tvjc3GIccERklHNzylSc/xttJ4602\noLy39SC+/uynuGfptuRZic7zxPJd+Nm8DQCA376xCQtrBSkjFFCbUYvPWrb1IL76VLd3Un3Sw+Ln\n89LTIGvdvUkhP9P2hla0dvJum4OoNoLqvLfVHCfyjec+xdOr9wAAWh3mVOCc42fzNqSupTL71Z+y\nu7EtJThVDNJ+jnEbdcF19c3tlkGDVhOLFTsOWboZS8s0TXCc/X5TXYs0QNBIms3Bph5BEtS1LCd+\nPigRdh5s9X1b3A+3N0pdeIMgI4SDH93D6PZpTves/a+7Z6zZfdhVLn9VZDpf41JZU+0sEwzUQEJf\nq4qVWmn5jkb87Z1t5i8EiNRNKi+y1s5us7ne+Own+LoW5aywcvjPZ/LnV9/cbltn/SCiN3D+8JUa\n/Jcg8FJlLLMznXVxjhU7GpWEvhOs7tX4PlhhXKnWt3TgGYuV0L1Lt2FbgzvVzN7D7a42lHK7/8Uz\nq/egyQf7xvdeWufKzhin1OUZIRzs8Gqn5pzjfxZtti1r+gPLTS/G9AeWi9NgeKiTyNBqdV6zgorH\n2OdEM6ddjW14XpBbRlYvL7h9Zq0dXSmVlKkIwYu1pb5F+vWX/rVGOXsr5902osc+3InNNnpoq3fc\n7tY/3dOEH71Sgz+85W3zImMdrFRuLxrUq1aOAk+t2pNmsH19/QHcn1wJ/c+iTabJzYtr92LRhnp9\nzQCopY8pKXA3RBndkp1wyAfh0OWDmWZnY2ug7sp2ZIZwsJvduRC2+p/oXxq7F3dXo9nu0eYyXaXs\nWk4HzpUK2RyNNXRbZ6/YXbWlo0s5jbnIoO/U+8mJd5M2g120Qa5q5IZ/3aDNHhttghz9fIJGQ3u6\nId/mfN1zeL2mTji4irq0XSZgq9+6RWllp/u/aEW1qLbO9xm+qLhrnlqLB9/rXgFp9fq/d9NX9/ub\n2tNUnn4Re+GgGYfdomJUvlevSrHpiaqDqlaMcXaU1qmk13L2Ouw51G57l6t3pWc2deOB4kRopQy5\nxhVL8oDeXrBd54X241dqcMMcVZdckc1BvY5AwrtJzQOJO3pZrAYP23a0MgwpXkOEiuCcm9z8x0nJ\nTgdvzavOKobE7f7VMtyO5cbNnwDgtws3odEm/sWt6Ghq70pLNigStMZV0cwFG/ANHzMRa8RaOHR2\nJYyX7T64oCzaUC91y3zpk+5UUSaVjuHSbZ1dePh9uTHRmKOp2TAI64vza0bkRkVTaJGCwU1aDsAY\n8CZ+ZqKX9Lo5a1MCYkt9C/YZZ2syby3DLbR2cvxpsTNVzAtr9yp5IK3b14z5ijNdwHpwsAtY8yug\nzYjKW6Ql//NrZnxj0kPO6R3p96920r/1tb5+zlrl3xm9FqNi2ZYGqTrROJZo+DE+ioi1cNBueW9S\nz+8muvHTPYkAnuU7GrFsi/0qxO4SrR1dmL1yN77w8AphwjEtFuOj7Y3oEDw0pYWDgVk2A54rtZrF\nby59fLX8S49olzW2zQdJ3b+z1UkC1QBG0S3XNakJwqdWdu9loVJFP8ZWYxHthlWr00tYjSHG+jbp\nByLB79ImORYPTW/zMRFwtLIoLkqr90Pv78T7OgePbq/FcPnJqzVYvqMR339pPQDrZyr7LigbdryF\nQ/KuX6tJ6HinV3dnCm9q60wNmlYDilPjkqpnRHsnN+U7SpCo88/nb8DqneZNalReqtSg19FlabTT\nvFk4uPJ75ttqRfE84wpC+9S7JD1zyyMf7nRdh8M2S/y5glmh1reKC8x3khbDIXjx9M9NpqoxHm3t\n6Eqp9uy6mOx7k1ux00EhQk+YRz7cibc31tufqMAWF0Fpojt/8H2xR5Mb7ygjqll8P9zeiNU7D2Ft\nMtu08RG9XlOHFy2cRG598TNfszvoibVwMFKhG1BqDzTj1aSb4tOr9qgVwID7Xbq4qUbG6h9uJzcr\nV1T0vlrfvOm5T3CHIbZBj9P9jQH/jJjCchQK1wbloqQXyt1vuvfI6dKVZ8VaQZp34+PUz3DF6bG7\n27qlo1sYXaZbZS2qrcOjEiE3d+3e1OzQLUb1wXNrFPt9ErdDiF/BoEs2pwsHld4rUpk0uFR7Auld\nVJZRQGll6LoGZvT3uHSzWYD+a7l5B8aDLR343tx1+GxvU2AR7LEWDtZ6W3cFfrAt4dkjG1NkbqSP\nfLBD8DuB2sj42exDmsJubN9xsA21+9O9EPYcasMf39osrKsTjC+qCh1d7oSk1fG3N9nXQ3bFxRvr\ncf+y7a5eVO0+tPrc+IzeCG5d4p5D3cLjcFtnypY1Z3W3r7/xufviHMbT6/yU6qQo/eeR4UY1/tne\nJszV2QS3N7RYTgjc2ErcxES0dnTZppYBgIc/2IELbGxa+jq/bZFpWf++b61vCTwgLtbCIeAgRiFr\ndov3K9byI6UZXW10saLP+u5kvIfbX1mPxRvrUZQvfyw1+5ukeWFkiF6YNzc4Fw5b6ltSidGU1UoS\n2Wg2/KefuGbXIfz2jY1pvxHh1p9dS0kiKjvtmMKNilKMLLXY+9ht33WsRTJ+jlg6uL2+Pqbkujmf\nuOq7skt/tP0g5qxOF7IqWqVvPPcpfvSK/Upw9ordQttj2vXsL2dJUE813sLBgJobqDf2yVzo1JXs\nKTYcaDbZHVr0PvyGXrh8xyG8XnMAfcsKpZfVe7Ko6kbPfbB75uIkmlrjAZ2v9bpkHn5hhLTg/8ar\nrdklFr5GFm2oMwROyXETfau9r0LhILoRi1dwx8FWvL2xHjX7uld5Vi7PtjYHUT1sEJW5pb4FV/xz\ndSpgzY2xU1YPN153qiuHJTarSVkcTGcXx83PW28iZHQi+fGrZrXtdoUcbofbOv1LyKfwHrd0dKWr\n2HQ/8WuvdiNK+zlEhbEv6W0LIne/97Y2YMOBZlx19ED315S6BJi/F886u4+K8ufoB3RRl7B7f/T9\n6C9Ltir/TkNk85i9wqzT1KNvd+33Rs8ZI1zQXgBSkegilm9vTPmPO5H9P3hZujGhlFQ72Ax8qu36\nS8HqQYatK6vka5UB/G/vbMMNxw8GkMjY29DSgcbWTpT3KPBVraSfZKj3PbXzZr5m3ZayYto6u7DZ\nyjsKwJ0L7J/T85LV6Kuf7kO1h30sNH7z+kZcemR/y9xcRr4yew0OShwvcnLlYGlzELxAj3+0Cw+9\n79zrRY/Jx16BPYfasDmZZMvJC9irON/xtYR7M3vsHWKvKzHaC75+f5NSaL/eeMu5hVoJwO2vdg/y\nfniMpMoSHLNaOegPWotyl3WwXTmIT7B6zuuTey48/7Fue9vUPSb+42bvjA8k6UV+u7B7kFVNTmlU\nK7lWr0WkHZv19lbMXuk9m+6bG+vx1oY6rNmVsBmotINMMARJLIXDhQ+vcLW5h6VPNdxJWC4bzXTl\nfeO5T/H1ZJItlXprL6lKXh+/O4XKktkKbaD5bG8T7viP3JNKwzhbVNU7pzV3AINBt0HaXHhXABd0\no4YRlsO5cMD+p86jxbgo+kdyBWuz1hMe3XCgWfjNRou9Q2SIphJOZs8a0sV9CEJjsc4d19Nz1P0+\n4HAP18RSOLR2ctQ3t+OVT/fZnqt/uWURhBpBby7e0NKBX7+xyfKcfYfbcF7SBrBiR0L/7sRQJ5pV\n6n+9etchLN8uz7Uk23hIFf1gX7s/fYcwbvEJcBbgJl4heXv7tZ/X7m/CO1u0PTvk5wEJl2lIz3TO\n9AeWWz4fPesFO7Btrm9JEwQiNOGm3cc7SeO49nn2il3e09IoPAujIdb4m/X7m4WZbe3okuinrGq0\nW2HPE6c46REXP7pS/lsX0iGoKHo9sbY5WG1TqQ0e2gCrgswH3RrjkCc3Oqi4tmkrixrJ1ovvbG7A\nV4+R20zsBtjfvLERh1o78dJ1U2zr4gbL3cN03y1PPhe9KqOuuaO7yWzuI18gHfyaGP51yTZh7IMI\nvyJnH9SlXNlmt3pL3rrRx59zrpTtU2ZO0bL3PvzBTkzsX2rwipM/EJHBU1QN49aab9Ske9UZ+45q\ngkUjb/oUTBcmTcaJq37pEFNiuXKwY8+httQq4YBkL1gRXvahFbnPmaJ/HYxeNbr4hddr6tLqZmlr\nERzT97PDbV2BZlzVz/5MHV54fvf/F2+sT608llj4cwNIi+/w+25abSJK/VRPdHZxx7mqrMYMlYC0\na5/Wcgqln6uPUTHPPOXlCnP9CE7/jWHVbBRuy7YexHadZ5nT2e8vbQzVUbCrsRU/ebUGbZ1dWLVT\nbUWooapWEk0I/QpMtCK2wuGwxcDzlSc/xl2LNoVUE/mjU9kjwYgxqA0Afv/m5tTS3w6rhHlA8JOR\nA15Vc8lGs/P9XilKPRLA+yB20xRdSN6ysr3NgYQKx69cVU+v2uMokMx4bw0unp3srkWRy8aVgKht\nnSQvNGIXMBl2HAcDsGbXYXy4vREL1h9w5TUHAE+sMBu5t9rYT78311u0vQqxFQ6ipHZ6DigmTJOx\nW7Avgxj56sDYFVW65t+SqjLjvtCqM4FCUYCcIE1HUDjN2KHaRsb3WjUS21ldEmWytGP2dZGfaY9I\n123VhCt3NEq/332oTbiPulGdo2GssXHFtN1CqKlgfEdVWiitbX2eydhdPwjhob0Pmh3k+TV7pM8j\nrS6wvv0bnlFNWx8csRUOQfPuFrUdwJzwpGAGIKPNYJ9QXWLfZ7ONp5vMtf5hfvnMKc/VXlDhoO2i\nRiLSl+mCOvt0HcD5aueHr9RYqr2Mqhsg3Z5hde20cZmlu227GTft8r2JBmP9Mb+7qt096INB/UJz\nudb+/du727FCZfOtiNxxnZCzwsEuWEbD6hkaH/CrFvsV25WrH9StrvnZXrMhW3S+aIYZDd7fAm3G\n73bmV2+YyekFsWiBIr6Ku6HMmYLKsgJCmts7pf3Oaj2pbnGQ4yYUxcmui07YVNeMy//pb6p5FRul\nltn18Y90zi6SG7tr4abU/4NYGftNzgoHt/C0/7t/wMa+IXJbVEU0ZnYbJf3F6fis+g5YnfbQ+zuV\n9huW8cG2Rmyqa8bupBuvflAT3o+P761QOPg4KlpFqpvuzUX6GeMOglZYrVREx/wMdKyX2FOmP7Dc\ndZmXPLbK9hwtCaPeTT5P0rhv6IIFX/p0n6PJZBTE2pU1mzHOgkUGWCdoek4/XzgRTlZSdufb/lhH\nfUsHKnq4764/fqVG2ZguFvouJYbFfb2zuQGf7DmM65PpLjS+5zG1t8Zyg3ojXd2v1k+cqF9VJktp\naiUfu2qsvEIVK6Pi7RclGS8cwtbd/VHn1ufl2sbfivcQUOegi5xEoeDT8+Hc3oPDCr1gSFfhmSvo\n54rfyuD9xIpd+GxvE04dWeHfBXW8XpOe1kI1RiVIglIrBTwncoRVVX5usT+LjKCDd2XEVq2kOrMx\nzo7CxMv7ZPytfucoFW8HGVG+JKLNh9zk8xHRxTm+ZZNxUxV931KtnmxjGDusVlOaU8KtL65zVba+\nLKeE0U38yuaqRnykQ57FS7jMuJtfjImtcFBlgQe/aRWCmm1ZDZrLHHpSBbXBuFPKBWof5ZrZSLV7\nfdzj93C7dcDh1nrrDWWcYKWg8iNY0Uk9jd5KXjGmYlGpSlCxCLFaOcSoLl6wFQ6MsWmMscXJv6kW\n55Uzxhbq/hqclpFpeDFIW78jzsqNehMXX7C5Bz9nXPrNY0RNLYwI9hHtVr0mQYyaH7xsbRsRPVH9\nPCYT02CokCWywdrmwBjLAzATwLTkoXmMsYVcMBpxzg8COCv5u6MA3Oq0DL948eO9uNcmHsAPPNkc\nXH4nIig9rh/EXW5Jg/IiuWqw6F85Y7p1P8IojSUYY3kApG1k5Sdx6me5snKoBrCOc97MOW8GUAtg\nrEK5twH4q8cyXOMmDbAbPNkcfOzNaR6KMeuZqncZ1bsddGy5MCV4wAZvEfsPt0M/VoviZTzXxVCZ\nRkFOqTgN4kERRsbUMLDzVuoDoJ4xNiv5uQFAFQDpepIxVgVgGOdccxJ2XAaQCGpxi5VByE+8BJr5\nGoUb2dBqxhw4FJ+6CQm4eqLivTgcKF1AwNWz11gXE0A7iFK9+OWgYCY+/azJZbbZuGEnHPYDqARw\nMxIai/sA2G2ycBOA+z2WgYc/cL+jW1gpJETbgKpSXCBftDl9f+I0GzO+GKp1i+oWZLEPfrWpqByV\nHfRUmbN6j/1JCnh1pU7ALT4lCMp3Ik7vwCEP2Z/jhJ1wqAUwTve5mnMuTT3IGCsAMAPA6W7L8IOw\nVg5GBvQsUt5U5OPd8v0EXlbY5EhPcO6Bzpj7yT6UFTnf+hSw35M6U2kQrBL87J/P+CQc/MAUIS0Y\nsYNS+UbRe2Qq3DgJKi9YCgfOeSdjbCaABclDd2rfMcauANDEOX9Z95OLAczlnHeplBEUmaC/tkrR\n7XR2pV+qR21yeNKwx25MvGwjQ7Rp0ToPqVLijClbh+CcIlFWYY/MW7c/EBuKHVnhJWiBbYQ053w+\ngPmC43MEx55xUkZQRPXQVHbpCgIXKXNCI072ECf4VWuvabEzCZU2C6I3/DFg12OnBGdXCZeMD4IT\n4aux10FhQWxMr0LaVeMmHTIU0daYbvCy+2AUHDWwp3+FiSKks2TgtCJbVstZKRyimqzGYeUQN+Jc\nNytujMFmK1HgyZhqtDl4q0rskdscsuPOs1I4+OsmGs11naBX3XjdIc9v/PGCCR8vKcL1ZNowseGA\nexdyFZtDprWHG2jlEGM4gPH9SkO/blS6xiyZqMSGPObfCx7XWWQQ2kfTvQp3ggvgwjEjrs/cKVkp\nHMC5b7EOToLxGluj0S9nSV+MDQzZP8MN4v5UVg65QLx3aVAnK4UDh38h7FEN+E6IyhCercQtBUmm\noLJHSRQup0EhczYgb6WY42eUdNz90rOkL8YGEg3uMHbDhbV1pnO82DQyhWx5H7NSODy9ao+vs79b\nXjBvMhPM6SzJAAAgAElEQVSnAcQq2ppwQZwebgYhigbPRVpivv2nKlkpHIDg8ytlyeSAENCRpak8\ngibTYjqCYrYhS0CmkrXCgdTGhAiVbhGVKzRBxInsFQ6kG3DN1DG9cePxg6OuRlbQnCUqBiL3yFrh\nEEB+r5zh9jNHmPYHzhZoRUlkOv0C2k3PSNYOobRycA9jLGuNstkSvUrEF6u9WvwgrC0JslY4hLXh\nT1QMqygO9gI0iBKEKyb1Dz87QxBkrXAg9QFBEFGQLavTrBUOhzIgstkLWdL/iCxmxoS+UVchEkg4\nxJwseT4EkbmEuHovyo+PqoAS7xGREp9XgSAIPUGvHMLaXZGEAyEkU7f3JOJDmE4hHSHocr572jCl\n8zoVVw6PXjnJS3UCJ2uFg5t+me0eTgQRJmG9TgN6FoWi5+9VXKDkpmrMPHveuCrheYN6ufM4DEtr\nlbXCwQ1nj+0TdRWUCbp/5EqciCyg6K7zx4Rck3AZ2btH1FXwjWMG9wrlOvl5wNdPyJ3MAVknHAb1\nKnL922+cOMTHmkSL143ic0WtJJuFRb2nQ8+i/EDLLykM49VPb8Pzx4tn0J6vEvN5TKa+S1knHB69\ncnLUVYgFlD5EDdnGLFE33/BK/2b2N5881HSsvjn49NrGQfvUkRWBXCcsdXCurKY1on4HYkXcZyDO\nyKqb8Y07zxmV9lmmq4565eCFE4eVp/5fVpQvXIWEkYIhrBbM5GflhqKQZn4ZLxwm9S/zray4dDFj\nPUa40A/n2PuizNGD0vXTsgV/JjsnTNC9E4n9sM132doZfLbYsJQp2qMqjOlD83vF0ZcS74WPcQai\n6UgH9HRvx3BXj/TPlx3R33EZMX1PYsNxQxJCQhaw5Ea4xqHNf3/BWHx5yoDUZ9kAHcZuZWHHgn3r\n5KH41fTRgV4jBo84NDJeOBQGGBk5pqoEAPCrc4PtcEaMqo4ozFlur3nlUc4FmQonDivHvBum+Fbe\nlKSHi3zl4Lxf5cdAOhTkMdMkRzRzPRTArm0n6NRZgHdDrGpzalcZ1acHThwejF0DSEwY7FRYo7LI\nCyzjhcONMtcyH9/TwjxxM133uUH+XSRA+vcMZxkKADecEIzHV2lRvmfdcoFgtBHZHCb0K8WQcuc+\n6GHpgq0QtVBY3jJnju6d9tlr7IHoeQlJXidsg7FoYnpOtdkdnryVAmTyALldQerV4eJ5yLpWUYH4\nm8uPDGaWbMLF+lz/ouSal4WMBy6fmGoJ7V+RWukvF41HeY8Cx+XfINg9L/SUP4brRZnnx3jtoPuh\nH3Y2q7HGWLxo8NxxMHs2yVISDoyxaYyxxcm/qTbnDmWMLUye+yfd8UcYY+8kv7vGSSVvOcXsihcE\nss7Vr0xscwjLS6Kk0LnPe9wM0k6W25ce0S+QOuib5LihCbWSl9ntDz4/PO2zpobUM9BlFKwVF0yQ\nxwt4GYBPNKiFnGLsc276bXp5avfi58z8DMPqBwDKixP3YVwpiOrXFoahPyR5bzs9YozlAZgJYFry\n0DzG2EIun5L8AcAdnPOlhuMcwJWc8y2uayuun/B4lY1Ff0h5MbYfbLU8Jw7JFU8bWZkayJwQlmwY\nWlGMbQ3W7QgAlSUFQJ1amT2Lzd3S7/sZnFQbtXS4f5lLdS6i8288BocFeny3A5fV73qXyPu2l0lB\nnkebyVGD0gMvldVCEpz+Oqg+/8MzRqB3SSGq+5ZgZ+P+7usJLlgWcPAiEJ4NUmXlUA1gHee8mXPe\nDKAWwFjRiYyxfABjBIIhdYq7aqbzOd1gKSuwqlT+At126jDcJIiGdroSCGMA/sW0UULjaGWJWK5P\njOkuVDGQs5HQx6IfWtHeKW8xq37nqU96eEgje/cwrbBlAYaqGLv9+H7ivu1n3xLNeUsK8zCuX6nA\n0G/GLumeaHUZV1SEQx8A9YyxWYyxWQAaAMjWtf0A9GCMvcAYe4Mxdonuu0YATzDG5jLGhMJFFbvU\nAr8+d7SlPWDGxL44eYTZqyEumpjTR1XanvPTs0bi1wIvqj9/cTyA7hfrnovGh5KxMiOI4gG7bPrq\nvu6EfNzUiV4w3opsYNXG4+DuXb3gzw0tx4R+pehVLB6jtH0n/t/pw4XfxwkV4bAfQCWAnwK4I/n/\nfRbnNgC4DMB5AH7KGCsBAM75bZzzUwH8HMDdXiqtF86ix3bCsAr0ldgJrIjLe2Wsh6hefcuKMKRc\nrsfXdM/j+pXi+KHudcmP+5SOxM0k8jGfUxrnBzxyij2F3OF28xovNoegvWqcTlKcrOTHVpVIbYN+\noq+RaFV40vAK/OWi8Sm1pYwYOLbZolLFWgDjdJ+rOec1ohM55+0AtgIYyDlvAyBSRrcAaHdaUT0W\nK25vKPRFoxHSb0RLZ9ntWr07fUoLMLpPQnjceqpzg742mA/QJTIUrWgONKk9SjfePyrpkVV5+IqJ\n6N+zKNVmQTgTiIp0q1kpdDt6+FgHv5kgUQvJUJkkadx3yQTXKjw9TsYWKwFu17u8PJPYbPbDOe9E\nwiC9AMB8AHdq3zHGrmCMfcHwk9sB/IMxtgTAnKSdAoyxJxljbyJhsP6hl0qfOVo3SAU4GRQ9ArsZ\ngZGxDnWMM89RD7izuvXigjz8/dKJADwMNAa0qGI9qt10VB9vulavUepBeA1ZMb5fKaYM7ula7z7J\nyqXS4sHHZfULpA+AJw+vQG/J4C1Lg60qv/0Ufk0OggO9emO5Jk77OXDO53POT0v+LdAdn8M5f9lw\n7hbO+QWc81M553/WHb+Kc34G53wG53yzl0p/XuBu5gcqfVF/ThBaCqZYD1NlQmBadR9866R0Q35Y\ns9KvHect4ND4rMJoukReI3eMqSrBz84eKfxubJV8Bi66Ly/99MKJfd3/WBGZc0UUgk61rW47dRju\ndDCRy0QyQPNljZ8dyKhqEHkuBD4WCqSD7B7DfnmK8vMwsnf6CkDURiKHgbPHOBfo+vsTpaZwonZy\noqJwW6aehIqDuRaeVmWfPKJC6rljHNw4ugV4pY1qT1TXChfqQKfI2kh1xRtmBLLWvEPKi1Ntk6bJ\ncIkogDJqMl44+EkQftVOZ21eBq0eLnX0mm3CDaLXUnRsUHmxr3sUjK0qwQtfO0r5/DACFvXX+MHn\nh4Mx9+6cdr+S3c2BJvk+Db+YNspxPU4RePW5wkXzG/uLZmz3Gj+hijD7qaAfucmpFRMzkCUZIRys\nPDDCzuWuerUfnzkieb75F49dOUm6PaUT+4C+7FtOGYpnv3pk6rNDvxBHZ+txMvY5Fr42z9ZTojvJ\nT/0ygvcoyPMcc+DG82hYZbptRf98jI9KZTLBGPDwFT54jVn0k8L8PLFrueT2jU4hQak2f3PuGFM9\nVJ+I/jeVuqBFP1KIxCkILmdwPMtn6YOzHm12IypzYK9i/OvqI0zHX7zmKJQV5bvqgJMHlLk2PKuM\nsTL1gtdAJyHJMit6FAh3MQsSL6+u8bednKOhxTyT/+Ikez2+7cpBUlFVNVCv4nzlGXifUvsyvcbS\neEndcYTHLXFlWO6bYNN0+tdi8oAyvHTt0dbnx3AtkfHCwY91w9XJ/PfGspw8rnPHGeICWdo/Smje\nD6ZITIVCTFHUFpX3Evdg7MT6T099OSHw/FzLybKjen2VfLE5JAsReTQyxjC4vBh7D7ebrnfisIpU\nBP8dU0cKy7aTufoZ6PRkJtD5Nx6j7EHj9x7VohQqxysP+NxRjiujUDsvqL2pFb8b1dvaE49zoMiw\nSvPkyhonb6U4cO44cypcv+iVzOVjHJS7HPTYonyGW3WrB02poDKw2232Lnvh9WU70bDccPzgtM6t\nFyyyyE5Z8aKOen0IxrXpgtTIYcM5cMboSpyVNLYbn/X1n+tuB7Ps5qkyJKUr12PqWLmxX58sTv9s\nGVO7gpeBKM192aJ/HjNYLXeYlh1ZXeio0aMgD/NvPAaA80nDq9dPwRUB7GFywYQq/OSsEb6X64SM\nEQ7f/7y4oewGX+PAO7RC3d/diU6bMYbp+tVDauUQjk3ESV0L8hh6GgaKVDkOdWtc8P8ZPro/iqoz\nuLwYl7jYHS+9XPF9OlUt3jF1FH505kgA6fs5MKi1qzFZnYbdmKz6uL+v08+P7F2i29sjvYAfnuEt\nuFOzselJu4Lkhip6FKCypFConjTeYmlhPv5zwxRP+1//7rwxpmOythYdN146X7C5khWpSYHumGh7\nU86jD17MGOEgI48xYXCWhnEDlocUjWuXHNEPxwzpZXIZVH1eTmbyfrg7hrFhvB3h9eX46WdF6J+J\n6fEkb0EWOc65zSxW8XGfMiLdzVJbhRr758T+ZZ5adepYdys5kU0GQNoqXIPDez8/VjdWlGoTR90L\nKFIW6K84Lpnzyq4W+mra2YH07sZxIuOFA2CdgVUF2YMe2bsEf71ovMsy1dVKbgO00jb0UTc5KMdR\npJ0jOOnk4RVp9gtbA2f08iuFnzYHlWuYYmhs6iF7fpoHjWskBcvq0c9DZHraLdu0lXFQriwpVHp3\n/vCFaod16i60VNHuoq+alqJDVUZ9++ShadHuqXdWQRpELS+yQjgEMWkulnj+uPEkAhIRldJzFcu0\n+qEnLxvBj2V9V3985vTR+I1uma6fIclsF7lCYr/h7s8dPmwC06ekQKhvd6O6vEKQtXhMnxLTikI2\n61VJPW1XK33q/dNHVpo8084a0zstbYpWnr5rytRyTtCXp+p9p7pvg8p5DNZuxcb0IrHJrRQH/Bj8\nbzllKL57mnyANvLlYwamfdYCxVQfCzP86ySPu+r9WqmVrPo4M3yvMrg4fQSnjfQeNRoExhm8amTq\n0IpiHOnAZTKPMctZqm2aaQ5To//gjG69vldb1vkTzHahaz83CC9fNwWAfS4r7epHWbSJE118VVkh\nLp6cvgPgOdVVQluIX1ufarXTq55FRYvuwpgpwIi277xMOBgvc+Ek8+6HKhO0IAk+Nj4mHD+0HIMk\nbpGlAm8hoyT/04xxuPixVabzxlaVoGZ/s/zC6qtIx7gVmsaqDC4vwto9hx2VYadC+t7pw/HqZ/st\nz1FBdBU/2pIx4OXrjkZBHsND7++wFfriQcO6DfTf6r2GuG7uJ5cN3PSgZKsx0Uxy3g1THG+ByhhL\nueWOrirB7kP2+yHL8iIB/qjuxvczJyD0c3uSX5w9CgUu0qPbvXt3TB2JzfUtyisbq/cpKm1sRqwc\nguaEYT6lCNChbeVoNZB9XrKpj7rNIcFlR/RDb4uXVPjb5I/vOn+MmpHPcIo+GjtTKczPU9rdy/jd\nRQpBbIAxslbmHSU+7nUAZIyJPdh8Hmks28vDtbRI7+KCPFN7e4qM18EYcNqoSpw0vPv992sH6N6l\nhZgyuJcnA3rqNiNyNskI4dBDMLPXZwf1OpPM86kV9I9QiyrVUgAHksE1Weg1nxvsMO1Gd5sdO8Sd\nz3hU6YpnXViNX013ZpTVfNj1uHkc2jNUTbGhv4Y+T5BIUNx26jAMM7pZmxwHfOhE+nQaaS9OfLwF\nhlYUC9U2Whu6zSFm5JsnCjyiXA4mv79gLEYo5A5r70qIH+0qXzqqP7509ADhuVr2aeOTOWVktzAb\n36/U9cZQdsRaOPz4zBF49qtHYlAyF7/TvRHs0PzxjU17oyC/vJsuc9Gkvhhq0WGkj1QgSY4Q5Pe3\n7hLyGhu/UVs4MJtSVcowY7fntbFukwf0xPDe/iXwCxJt1njTCYNNExzjIDR5QFnaxkp2Y9SlR3Tr\nqP2JpfFf7+lX3rNTRlTimMG9cKaLzL5WnCZYubudaE4Z3Av/uHyi7XlGjcKNJwzBNZJ09CJV09GD\neuKqo7vtoeP7lWKyxd4fXoi1cJg6tk8qehlI2A00w+4I3QChV5E4QdaoTjJRWvWlb58yDBXFmW/W\n+fMXx+F3DttWhKitvn96NFGg+oFLKVLY6/UclmfnkXLqyEqhx1GcEIkGN3U+Zkgv/M8FnradV0aU\nrcCBR64tfu5uCAS73ou1cDBy3fGD8bdLJgBIqEOM6gKRioSZ/uMOrw9B+HtJoarXslK9Wg0tTu9l\nYv8yW+8ML9x1/hhhlKgq5QG6zeoNyU4nwnqdsf6nooFf1Usm7TfOqmODu/a/+Ih+uMqgFrn8yP4o\nLcwTzny/fuIQ07EwGFNVkjZDHyVZffYtKzKlghcZrMNUwln1uwn9g1k1ADngrcRN//FWjrora/Dd\nJ8h9kIPwpZYJyGOHlKOipAD7DnvaWlwJP9wgVYvQPx/js5oxoS8aWrujg00JDQXXkD0Tu2cl+/bc\ncVVoaOnAG7V11r/XFVBckIfWjm6z7eQBPTF5QLpHTkWPArxwjXUW0rCZXt0nLeVKnsVkRHM75hz4\nv0snYKSPakw3b2zvEnmQ79lj+2D+Ou9egSKyXjgETVVpIWqtXFmTiMbx7omlusfMKSMqsHRzg3oF\nJdx+5gg0t5t9M9yOnd8IaEYYVm4qKRz460Xj0LMoH39ZstX29MkDyvDxboNbMOdpqzwGhut8SE54\n+ZH9U+kc3PCtZNCZSDgUSQbPBy6biFlvb8FBSdqLbMPr3udGUoF8Dl40u5iToMIeMkqt5AY3aiUn\nA+RPzxqZSlUtLCv56IIY5NyWyAGcOLzCVwPfYEkMyVEDe+KJqyc7KistcEjFWB6Qq9/nhvbCCcPK\nMb5fGYZU9FCqzE0CIcmRGMivdbAPtkoX7FNa6LuRVmPGxL6485xRANInNgN6FeGu88f6qrbMKVw0\nTlRp07JeOJQWqex2pd76xjz4pUX56C3K7eRgUHP77K2q7WyHNu+9T1YXxhJ6XPm1E0zQJTg8f0Jf\n3Huxu5xWfvLb88biv3WbyJ80PGHTcvqyciR0w8ao+7RzeMIrR/PecjIbDGLiUd6jIJW0T+tLfnsL\nBskFE6rwrZOG4GqJm6g8+DBYSgvz0aMgD9WCFd+8G6Z01yPqxErIcuHw2JWT8McZ43wtc1SfEvz6\n3NH2J3rActBPO09+Yn+LpajxVwN7qSdXs1viWmE1hN0xdVTKk0X28shwskS33UTHok01nbXjF5cb\nPwoM0kjM1v/8xfGpA2HOGFWudV/SGQRwPqEJc/L7u/PG4JrjBuGSI/oL1Xc3nTAYXz1WfRXnJ72K\n8/Hva48WGpLD3vLYjqwWDgN7FStlbFV5JMUFealtA3sWKZhqTMEEChdxiKzIl649Gl8SbEAi27f6\nyqMH4Gdnj0yUaVPPYZU9hEFl0joq3nd+HktLp+wEP9KVewmKs8PNJDCQ7VcVufSIfrYTBtXa/T0p\nUETnB3WLxw0ttzTiXn7UAOGe1UHz6JWT8DUbofTq9VNMx+y6WT+LlbkXslo42CHKvHrbqcOE+vOC\nPIYnkvs+T+hfarl6+OX00akIaStvKa/R8bKfFRWY00IA8syw+XkMnx8VjO7afqZu873CNfzIpuDm\nGagObiqeX8YzOIAug79AqSQq3Q/PMr3//TdPGuotAl7XlqMzSBUFJDLfBsWgXsWm7UKNuEkN8p1T\nh2HOV/xPZ5PT3kqnjKxIzWw0VHYxy2PMMh/TicPKzYNz8qPeq8dK9WNE+PrHaxXqGSuPrtDqoHCO\nY62S4QdCG4Gg0DJd/EZZUT6GKaRnULm+kee/dpRyCmq36O/4xGHlWLb1YKDXc0MY7aCKXuDbTVyK\nCvJshY4bsmLlcNGkfviKhbEPEL/0eYz5O7NJXiQt+panfYXLdBGimoFSafBIFWi6nPIAOsCBbSEK\nwtCkBH0Jt/cginPQby1qtT+GV4O0mwHRy7M6cbh8YnXl0QNMQXVhERfBYEWYZonYrhycbCA/tm8p\nxnrw944KLYGYaZGhOIN2YsBSsRNccWR/yxTMQeBHX/dri9RRvXugSmKX0eNUHWiM/lZTM3FoWs+v\nHTcIA/3akc0FfsxK9Xf8hQlVOKe6D+5csMF0nur+GlETlfE4TFNULIXD0IpiHDfUnXFSRgw8w6T0\nKSlEY2un49/53T2DSG2g30hFe6F+MW0UfvnaRl/Kv/WUoSjvUYDfvLEp7fj546sc7ydxn0HFKMNJ\nX3rw8olKW22aXnrerX+2WxU7LtshBXnMmROC3feMobiAYdKAMqzf1+StciEThaNAVIGgsRQOD10x\nKdTr+TYJsNl9TcaNJwy2XfloRV80uS8GlScGm5h5vpl44LKJaSuR648fhD2H2nHayMpE2nAgLbGi\nCLt71HbQ+mDbQexqbMPKnYcAJITDjoOtqc8aMoMfs/jOC6p2An3XGVJejMEVxdgfQjoRt1g9F+Eq\nV3Dsq8cOisyl1C2RTDL10fVxUisxxqYB+O/kx//mnL9hce5QAI8ny32fc/7/nJYRFFZtarerWdAU\n5jMll1sgkYdISzAYc9lgSq2tN+IzlpjRlkv2KO5G7S6///kRmL1il0kYGOlRkGdKrPanGdW+vXR+\neA49/KXE5Ki+WS1FRVh7CqddM85L8QCJ+r5jo1ZijOUBmAlgWvLQPMbYQi6POvoDgDs450s9lBE6\nJYX5eC5GO5s5tTkE0ZBOn47TsbVvWSE6OoPtArIB37if8xEO9oZOK9/Vr9TJj/vSUJHJA8ocpcGP\nM1GMWlH1ArtpWzWAdZzzZgBgjNUCGAtgvfFExlg+gDF6weC0jCjpGdC+C1dPGYAjB/YMdDkYtZwd\n3acEQ4y7mNlw38UThL79fhFUe08Z3BMnDS/HM6v3iE/w8SZ6l6r1yciTE9ow60J/sxRESVcEqzT9\n09X69ZmjK7FoQ32g17XrfX0A1DPGZiU/NwCognhg7wegB2PsBQDlAP7KOX/eYRkZjajbDOxVjIG9\nilG7X93wJnrV/eySKmU5GVz/fqmaIVePvTopnvz+gurU//20UYjke6/iAkeG4HCxMrDFW1h5IWp9\nh3b9Cf3LIhcO+wFUArgZiTHrPgD7LM5tAHAZgHwASxhj/3FYBuECp5vRq6gr4rJ8djPOzJjQN23P\n5iC456LxGNvXvxiZsMbTCyf1zTgPoTgRZVqTsLETDrUA9GvCas55jehEznk7Y2wrgIGc8+2MsVan\nZRAJRGqCtk7z3gtuidoAL8Ov1+6208RpQvxkXD//4mruPGeUp30Z7NC36xcn9ZOeR9hT0aNAmqMs\nKPQTh9h4K3HOOxljMwEsSB66U/uOMXYFgCbO+cu6n9wO4B+MsQoAT+vsDMIywiTuelk71u+Tbyjk\n1OaQxav+yHEj4LTU2ET8KSnMx7+ulu/fkk3YKn455/MBzBccnyM4tgXABapl5BJKmVyTHD24JzYc\naMafLqy2PxnxDvDzSv+yIvRVdPMlxAQxF7Caj4Qx9zhtZCXe3hSszj1KZO37/dNHAAinjTPTKugC\nQQLWUBnQqwgvXyfeV9e4qjljdG+cMVotS+o1xw3KiJwwbhnQqwhPWOy0Z0UU6rM4CupAXJ0tvhvh\n457LMioz1KHBC2VF+aHmR8uJFn7w8omWu5GFRaFEQqkEMR01sCdW7TIHeP2Xx9QKMuI4yDnlxuOH\nYHp1S6jXzBV7pew+4+td5Zy4PEpKnxEgblMdOyXIzlSYT4YCp4yuKsm4/QQyhagHzgsn9U3tmRIE\nPz97FMotMuGGSZ5gThlG++eEcIg7mW4sJ5xR2SPzbShRB16O6lOCUX2CE/ynj4qHk8CfZlRjREiT\nWyMkHAgiRP597dGpVO1+Mb5fqXBb2CCJeuWQzfxiWrdrsz61S9hTSBIOMaBvhH7ThN9YD5t+CwYg\nkSvqxhPk6dZPHFaOy46g+IZM4bSR9qsW8lbKMNwstbPJgEfEkz6lhfjGSUN9LTNXDO+5TFZsE0oQ\ncaFHYTyMmEFDsiH7IeGQIfipClIpKhNf/jjEe4zrW4p/XjU56mqEQCb2EMIJpFbKMS4/sj96l8TT\nW2Zi/zJc6cGw+oUJfXHcEH+3l3VDfw/7PWcKpFYKl1tOGYrSkFelJBx8RBbkFiduCmCfaL8oK8rH\nDRaGVTvy8xiGVETj9pdrOM0ETHjDmDDx9FGVyjsFuiX+o1kG4XQjdoIgCDf0LSvCdccPDvQaJBwI\ngnBMFPtWE+FCaiVCDCmVCQtuOH4w9h5qj7oaRICQcCAIwjFnjekTdRWIgCG1EkEQBGGChEOGcOTA\nnuhTQgs9giDCgUabDOHqKQNx9ZRg9m4gCIIwQisHgiAIwgQJB4Igsp7yHNxW1CvUYgRBZDVPffkI\nVJK9zjHUYoQQinIgsoXepfHMJRZ3SK1EEARBmCDhQBAEQZgg4UAIoewZBJHbkHAghHSRdCCInIaE\nAyGE8vUTRG5D3kqEkDPH9EY7SQiCyFlIOBBCThtZidNGVkZdDYIgIkJJrcQYm8YYW5z8m2pz7iOM\nsXcYYwsZY9fYHScIgiDih+3KgTGWB2AmgGnJQ/MYYws5l1osOYArOedbFI8TBEEQMUNl5VANYB3n\nvJlz3gygFsBYm98wh8cJgiCIGKEiHPoAqGeMzWKMzQLQAKDK4vxGAE8wxuYyxsYqHCcIgiBihopB\nej+ASgA3IzHzvw/APtnJnPPbAIAxNgXA3QAusTpOEARBxA+VlUMtgHG6z9Wc8xqF37UAEO1ALjtO\nEARBxATblQPnvJMxNhPAguShO7XvGGNXAGjinL+sO/YkgEFIqJG+bXecIAiCiB9KcQ6c8/kA5guO\nzxEcu0pShvA4QRAEET8ofQZBEARhgoQDQRAEYYKEA0EQBGGChANBEARhgoQDQRAEYYKEA0EQBGGC\nhANBEARhgoQDQRAEYYKEA0EQBGGChANBEARhgoQDQRAEYYKEA0EQBGGChANBEARhgoQDQRAEYYKE\nA0EQBGGChANBEARhgoQDQRAEYYKEA0EQBGGChANBEARhgoQDQRAEYYKEA0EQBGGChANBEARhgoQD\nQRAEYYKEA0EQBGGChANBEARhgoQDQRAEYYKEA0EQBGGChANBEARhgoQDQRAEYcJWODDGpjHGFif/\nptqc+whj7B3G2ELG2DVuyiAIgiCip8DqS8ZYHoCZAKYlD81jjC3knHPJTziAKznnWzyUQRAEQUSM\n3Zlq/RAAAANWSURBVMqhGsA6znkz57wZQC2AsTa/YT6UQRAEQUSI5coBQB8A9YyxWcnPDQCqAKyX\nnN8I4AnG2AEA3+Oc17gogyAIgogYO+GwH0AlgJuRWBHcB2Cf7GTO+W0AwBibAuBuAJc4LYMgCIKI\nHjvhUAtgnO5zdXI1YEcLgHa3ZXz00UcKlyAIgiCCgtnZhRlj0wH8IvlxJud8QfL4FQCaOOcv6859\nEsAgJNRL3+acb7YqgyAIgogntsKBIAiCyD0oCI4gCIIwQcKBIAiCMEHCgSAIgjARG+GQKyk2GGN/\nT6YXeZMxNjp5THjvTo9nIoyxYsbYZsbYt5Ofc7ItGGNDk/1iMWPsj8ljudoWX2OMLWOMLWGMnZU8\nljNtwRg7nTH2HmPsbt0xX+7fUbtwziP/Q0JILQFQkvx7C0ljebb+AZgK4G9IxH6k3busTbKxrQB8\nB8Bz6I6Dycm2APAkgFN0n5XvOQvbYhWAfADlAN7JtX6BRKqhSwDc7VdfcNMudnEOYZFKsQEAjDEt\nxUY2R1E3AmiD4N4ZY9VIPEil48jQtmKMlQI4B8AcAD2Ro23BGMsHMIZzvlR3OCfbIslaAGcAGAjg\nXeRYW3DOX2OMnaE75Pn+3bRLXIRDLqbYuB7An5G4T9G9M4fHM7GtbgNwD4AByc+52hb9APRgjL2A\nxGz5rwB2ITfbAgDmA/gugEIkMirkar/QkI2PTu/fUbvERTjkVIoNxtiFAD7jnH/KGBsH8b3nOTye\nUTDGKgCcxjm/izF2bfKwrB9kdVsgcd8NAC5DQp2yBMANyMG2YAk73AzO+ReTn98CcAtysC10+PVe\nOGqXuAgHt2k6Mg7G2HEAzuCc/yB5SHjvSVWD8vFgax0IpyIxW54NYBQSfXExcrAtOOftjLGtAAZy\nzrczxloB1CAH2wIJ4VgAAIwxhoRuPBfbQp/d2pcxwnG7RG180RlhpgN4O/l3TtT1CfA+NwB4E8BC\nAH+2unenxzP1D8A1AG7O5bYAMBzAK0isGr6T423x02Rb/AfAtbnWFgBuB7AIwKcA/s/P+3fSLpQ+\ngyAIgjARmzgHgiAIIj6QcCAIgiBMkHAgCIIgTJBwIAiCIEyQcCAIgiBMkHAgCIIgTJBwIAiCIEyQ\ncCAIgiBM/H/svsGfnq/EhQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c483fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2\n",
    "\n",
    "Carlin (1992) considers a Bayesian approach to meta-analysis, and includes the following examples of 22 trials of beta-blockers to prevent mortality after myocardial infarction. These data are given below.\n",
    "\n",
    "In one possible random effects model we assume the true baseline mean (on a log-odds scale) $m_i$ in a trial $i$\n",
    "is drawn from some population distribution. Let $r^C_i$ denote number of events in the control group in trial $i$, and $r^T_i$ denote events under active treatment in trial $i$. Our model is:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "r^C_i &\\sim \\text{Binomial}\\left(p^C_i, n^C_i\\right) \\\\\n",
    "r^T_i &\\sim \\text{Binomial}\\left(p^T_i, n^T_i\\right) \\\\\n",
    "\\text{logit}\\left(p^C_i\\right) &= \\mu_i \\\\\n",
    "\\text{logit}\\left(p^T_i\\right) &= \\mu_i + \\delta \\\\\n",
    "\\mu_i &\\sim \\text{Normal}(m, s).\n",
    "\\end{aligned}$$\n",
    "\n",
    "In this case, we want to make inferences about the population effect $m$, and the predictive distribution for the effect $\\delta_{\\text{new}}$ in a new trial. \n",
    "\n",
    "This particular model uses a random effect for the population mean, and a fixed effect for the treatment effect. There are 3 other models you could fit to represent all possible combinations of fixed or random effects for these two parameters.\n",
    "\n",
    "Build all 4 models to estimate the treatment effect in PyMC and \n",
    "\n",
    "1. use convergence diagnostics to check for convergence in each model \n",
    "2. use posterior predictive checks to compare the fit of the models\n",
    "3. use DIC to compare the models as approximations of the true generating model\n",
    "\n",
    "Which model would you select and why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r_t_obs = [3, 7, 5, 102, 28, 4, 98, 60, 25, 138, 64, 45, 9, 57, 25, 33, 28, 8, 6, 32, 27, 22]\n",
    "n_t_obs = [38, 114, 69, 1533, 355, 59, 945, 632, 278,1916, 873, 263, 291, 858, 154, 207, 251, 151, 174, 209, 391, 680]\n",
    "r_c_obs = [3, 14, 11, 127, 27, 6, 152, 48, 37, 188, 52, 47, 16, 45, 31, 38, 12, 6, 3, 40, 43, 39]\n",
    "n_c_obs = [39, 116, 93, 1520, 365, 52, 939, 471, 282, 1921, 583, 266, 293, 883, 147, 213, 122, 154, 134, 218, 364, 674]\n",
    "N = len(n_c_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fixed():\n",
    "    \n",
    "    delta = pm.Normal('delta', 0, 1e-5, value=0)\n",
    "    mu = pm.Normal('mu', 0, 1e-5, value=1)\n",
    "    \n",
    "    theta_t = pm.Lambda('theta_t', lambda mu=mu, d=delta: pm.invlogit(mu + d)*np.ones(N))\n",
    "    theta_c = pm.Lambda('theta_c', lambda mu=mu: pm.invlogit(mu)*np.ones(N))\n",
    "    \n",
    "    treat = pm.Binomial('treat', n_t_obs, theta_t, value=r_t_obs, observed=True)\n",
    "    control = pm.Binomial('control', n_c_obs, theta_c, value=r_c_obs, observed=True)\n",
    "    \n",
    "    treat_pred = pm.Binomial('treat_pred', n_t_obs, theta_t, value=r_t_obs)\n",
    "    control_pred = pm.Binomial('control_pred', n_c_obs, theta_c, value=r_c_obs)\n",
    "    \n",
    "    return locals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 10000 of 10000 complete in 4.7 sec"
     ]
    }
   ],
   "source": [
    "M_fixed = pm.MCMC(fixed())\n",
    "M_fixed.sample(10000, 5000)\n",
    "M_fixed.sample(10000, 5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEgCAYAAADlpDdIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE0hJREFUeJzt3HuUJGV5x/Hvs0LwcFF3l6iIIKBoMnjBGAyixpXxkuXg\nhaNoEjkK61FEEw0xaNRE8MpJgABqjEqCoAiKSQAv8QYLCoqCEhBnJSDuyoKgxF3lIhCEJ39UDTTj\nTE93T8/0M8v3c06fna6qrnpq3ur69fv2uxOZiSRJVS0ZdQGSJHVjUEmSSjOoJEmlGVSSpNIMKklS\naQaVJKk0g2oAEfGAiDgyItZGxA0RcVJEbDllmwMj4taIWN8+romI5e26JRFxSrv8WxGxXbv82Ih4\nc5+1PDQiTm7rWB8RayLiHcM72xmPe0pEHD7Duo9ExAunWb4uIvbuss8t2nP4aUTcHRHzdn1GxNkd\nx1kfER8ccD9HRMQnh13fsEXEWe153hkR46OuR9Nrr6fJ+8Z1EXFZRLy0j9eviIj1Q6rlwIg4fxj7\nmqvNRl3AIvVG4BnAE4HbgBOAY4GDO7ZJ4DOZuWqa1z8f2D4zd4iII4A3RcTJwDjw5D5r+QJwGfDo\nzLw1Ih4E7NjnPgaR7eO3V2S+rstrZt5h5h3ADhHxKGDt3MrrLjOf03GcR2Xm3fN5vFHLzBcBRMRa\nZmkHjVQC/5GZrwSIiKcCZ0fEDzLzitGWNjr2qAazAjg9M2/OzN8AhwKviIhtO7aJ9jGdO4EHRMQD\ngC2A3wDHAX+TmXf1WkRE7AM8Ejg4M28FyMybMvMHU7Y7qe0Bfqjt2V0XEU/sWB8R8bcRcVXb6/lI\nRDxwyj5eHhE/aj/pfQLYeur5RcT32vW3R8SrZyj7CRFxXkT8LCIuiIhdpzu1Wc57VUT8sD2X0yNi\nWbftu+2qyzHOi4jDIuLj7e/kxxGxYso2l9C0/X4dveZHTNlm3/ZT8fqI+EpE7Dhl/U5tr+6ZEXFp\n2zYndaz/WER8YMprTouI93Y8XxoRx0XERERcHxHfj4hn9P3LaH6vl7fncU2/vXsNxX3uG5l5Ec2H\nqd/raycR+0XExRFxY0ScGhHRsW7W6yUiPgccA+wR944K/dGczmwuMtNHnw/g74ALgOU0F9VTgV8C\nz+7Y5s/aZeuAbwDPmrKPfwauBr4EvBY4c4A6DgfO6GG7k4CfA69tn28zZf2hwOXAdjQfXj4OfKhj\n/e/T9BxXtM+fBvwCeOcMxzsXWDXN8nXt+S5tf2/HAN+bZrudgLuBJdOs2w/4KfDY9vm7gS8M2I7d\njnMesAb4g47jnD9DG3xihv0/BbgJ2LN9/mqa3m9MU8OXgYe3yx7Usf4PgRuBzdvnDwZuoekFTm6z\nBfC8yfMA3gVMzFDTWmDvaZbv2rbxWMeyrRb6vXV/fwBHAJ9sf14CHABcCyzt8fUraD4Iv7+9LpYD\nG4Dxfq8X4FXTXfOjeNijGsyRwDk0AXQlcAhwB/DQyQ0y8zRg18zciebiOzMidulY/4bMfDTwYuBN\nNGHRr+1pblqzSZoe4MfaY988Zf3BwHsz8/pshsDeTnNTnfQy4POZeV77+guB/2KWns8MdRyVmRuz\neSe8E3hSOwTXq4OB4zPzyvb5u4EV0X7PN0QJHJ2Zl7TPz2f6IdVuPefXAJ/KzG8DZOa/0dwk9pxm\n2wMy84Z2u5vuKSLzu8B1wL7topcBF2TmTzq2uSMzv5r3Dl/+O/C42U/xPn5Fcw0/LyK2b/d7a5/7\n0HC8uB2ivY2m3ffKzI19vP6GzHx7e138guYD1w6TK/u4Xvp9f88bg2oAmXlXZh6embtl5q40N/Vt\ngJ9N2e7G9t/VNKG272/tDA4DzgAOiYgLI+KzEbG0x1J+RUc4zqLbhb4DcEw0k0PWAt8Cbuu4+T+M\npjc0VO2NcCPwu328bAfg0I5arwJuZ36+l+t8o95J/++XHYD9J2tt611GM1w71YYu+zkBOLD9+VXt\n83uLbBwaEd+IiAtoeutLoo/JKJn5c5qRgR2B1RFxSTjpYlTOyMydgdOBB2fmNZMrIuJpHUNx6yPi\nmz3s7z7X7jCul4VWtrBF5oU0n0a/3WWb36H5hHSPiHgk8ArgO8Aemfk0YAL4qx6P+x2aMeSt+q74\nvq4GXpmZO3c8lmXm9e36a4Gdp7xmMwb7Un7zyR+imQW5DLhm5s2nrfUdU2rdNjO/M0Ats5nrpIOr\ngX+ZUutDM/Ozfe7nUzS9xr2AxwBnTVn/l8BLgf0z8xk0Yda3zLwyM/86Mx9HMwrw+YjYYpB9aWDJ\nvR+QXg/sGhGH3bMy88LM3KHj8fQBjjGU62UhGVQDiIitIuIh7c+7AP9I833N7R3b7BjNZAki4rk0\nMwTPnLKro2luCL9pt5scRnogvTmL5vuaj0Yz24+I2DYipvbcZuvC/xNwXETc0/2fPL/WZ4B9ImLP\n9tPYy4AXzLLP6Y4ZwD9ExMMiYjOa39uX20/zve7jWODwaGZDTdbaaw+0n1q7Le+0gWaCyOYRsdmU\nIcgPAwdHxPPv2eEAtWbmL2na+jTg5Gwm8HR6JHAD8PNoJvQc1S7fnOlNe16dQ9M0H6zuAHqe3KOh\n6JxIcTPw58ARndf7XPdL79fLRpqg3KZ9328/xxoGZlAN5gnApRFxHfBV4IOZOfX/4bwZWB8RPwHe\nAjxvcigQICKeCTwsMz8DfAW4nqZn8XzgA/QgM++kmdJ+F3BVRFxL813K1CnuM04lb/dzEs3F+ul2\nttda4L0d66+m+RR2RlvjCprvqLqWN8OyT9EE9s+AXWi+x5nJuoi4dEqt59JMPjm+HfpYSzP5oy8R\ncTbNEGe2x5nu/1FNPYfpzulU4Nc0vc7Lab5znKz1CppAf2vH7/Vz0/RSeum5nUBzg/nXadYdA2xL\n86Hli8DHaD78zHRjOa393d0z47L9oHN6RFwbzf/DeQOwcppQ1Py6z3s1m1l/7wFOjYht+thHt2W9\nXi9fBi4BfgxcQfO+G4loZ3dIklSSPSpJUmkGlSSpNINKklSaQSVJKm1R/FHac845xxkfkhgfHy/x\n1xK8J82Pmdp3Ucz627hxY/0iJc2rpUuXlggp8J40H7q1r0N/kqTSDCpJUmkGlSSpNINKklSaQSVJ\nKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJpRlU\nkqTSDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklSa\nQSVJKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJ\npRlUkqTSDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmmbjboASYtb\nZrJmzRoAxsbGiIgRV6RNjT0qSQPLTA466P2sXLmOlSvXsWrVkWTmqMvSJiYWw0W1cePG+kVK9zPL\nli2ddZsNGzYO7XhLly4t01XznjR83drXoT9Js+ollAZ93TDDTJsmg0rSrGYKk8mhv9WrdwdgfPwy\nTjzxbX5PpaHqa+gvIi7OzD26rP84sAL4i8z84pR1r8nMEwYp0m62VNdCTaZw6G/T1q19hzqZIjMP\nAk6aYfVrh3ksSaokM5mYmGBiYsIJJUM2a1BFxNsi4uKIOAnYul32JxHxrYj4ZkS8pId9fAJ4XESc\nGxF/37H8oIg4KyK+HxFvnMN5SBoBZ/01/D3Mr65DfxGxHfCfwNOBrYDLgZ2BS4G9gDuA1cBzM/OO\n9jWHA9+dZujvt4YNI2LzzLwzIrYALsrMJ01Xh91sqZ5eJ1gMa7JE5aG/iYkJVq5cxy23HADA1luf\nwpe+tBO77bbbSOpbjOYy628HmtC5G7g5Im4EtgW2B77QbvMQ4BHA2gFq++OI2Be4BdhygNdLmkeD\nzvbrZR/O9lOvZguqHwNPjoglwDJgu8y8MSJ+CLwoM2/q41i/ExFL2tCbdDzwBGBH4E/7KVzS/Jst\nTJz11xgbG2Pvvc9k9erm+fj4ZYyN7TPaojYhs876i4i3APsDPwD2yMzHR8TTgfcBCfw0M1/Rsf3h\nwEuAkzPzmI7lRwOPB9Zm5iHtso/SBNV/A0/JzD2nq8GhP6kuZ/01/FNSc9Otff3LFJIWhepBpblZ\nsOnpkiQNm0ElSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQ\nSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwqSVJp\nBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIk\nlWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwq\nSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrN\noJIklWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZVJKk\n0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkEl\nSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZ\nVJKk0gwqSVJpBpUkqTSDSpJUmkElSSrNoJIklWZQSZJKM6gkSaUZVJKk0gwqSVJpBpUkqTSDSpJU\nmkElSSrNoJIklWZQSZJKM6gkSaVtNuoCFlJmsmbNGgDGxsaIiBFXJEmazf2mR5WZHHTQ+1m5ch0r\nV65j1aojycxRlyVJmkUshpv1xo0b51zksmVLZ1y3YcPGue5e0jxbunRpmSGQYdyTdF/d2neTHPrr\nFkq9bG9wSVIdm2RQTRc0k0N/q1fvDsD4+GWceOLb/J5Kkoq73wz9gZMppMXMob9N2/1u6G8mEcFu\nu+026jIkSX1YkFl/EXFgRJwSERMR8fqI+GFEPCoiLu7Y5uJu+xiGzGRiYoKJiQln/EkaKu8v82eh\npqcnsBY4EdgGOA3YfYGO3RTg9HRJ88T7y/xayKG/G9p/bwEevsDHZvnyZcDR9zw/66wDWL68+dlZ\nfpLmYs2aNaxevTu33HIAAOec0yzzq4bhGPV3VEsAImJLYMth7rifKepOT5ekukYdVBdFxNE0vayh\n9pOnho3T0yXNl7GxMfbe+0xWr26ej49fxtjYPqMtahPi9HRJi0L16eneX+amW/ver4JK0uJVPag0\nN93a937zR2klSYuTQSVJKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklSaQSVJKs2g\nkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJpRlUkqTS\nDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklSaQSVJ\nKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJpRlU\nkqTSDCpJUmkGlSSpNINKklSaQSVJKs2gkiSVZlBJkkozqCRJpRlUkqTSDCpJUmkGlSSpNINKklTa\nZqMuoBeXXHLJqEuQNHo5Pj4eoy4CvCfNkxnbNzJzoYuRJKlnDv1JkkozqCRJpRlUkqTSDCpJUmkG\n1QAi4iMRcW5EfD0idulh++dExPntY++FqHEuIuKZEXFRRBzV4/YnRcSF7e/kVfNd31wMcG6Lre36\nqncxtV11/V5bozjuQrZ3v/fJbhbF9PRqMvN1AO2N4DDgkJm2jYglwLuA57SLvhIR52bt6ZZbAEcC\ne/W4fQIvz8xr5q+koen53BZb2w1Y72Jqu+r6fd+M4rgL1t793CdnY49qbm4G/m+WbXYFrszM2zLz\nNuBq4DHzXtkcZObZwIY+X1bi/7fMps9zW2xtN2i9i6LtqhvwfTOK4y50e9/nPhkRK/rtzdmj6iIi\nngu8ZcriN2fm99ufVwHHz7KbZcAvI+LY9vmvgOXAVUMrdEA9nF+vbgZOjYgNwKGZ+aOhFDgHQzq3\nxdZ276H/esu1nebVKNp7FXB8RDwJOA54CPDAiDgQeF8btF0ZVF1k5teAr023LiJeAPxPZl4xy25+\nQdMwr6f5JPNh4H+HWeegup1fn/t5I0BE7A4cBew3133O1ZDObVG1XUQ8lj7rrdh2mj8L3d7T3Cef\nHRHPAnbKzJN73Y9DfwOIiKcAz8rM43rY/GrgsR3Pd10kn1oHGR64Hbhz2IXMg17PbbG13VzqXSxt\nV92ohlH7Pe68t3ef98mu7FEN5rPA+og4F7h88lMKQETsD/w6M78IkJl3RcS7uPfT7xELXWy/IuKt\nwErg4RHxoMw8uGPdfc6vXfZpYDuaYYU3LHS9/ejn3BZb281W72Jvu+q6XVujOG6B9p72PpmZXwe+\n3s+O/Ft/kqTSHPqTJJVmUEmSSjOoJEmlGVSSpNIMKklSaQaVJKk0g0qSVJpBJUkq7f8BfVnWd7Le\nRA0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c5a3a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.Matplot.summary_plot([M_fixed.delta, M_fixed.mu])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def random_mean():\n",
    "    \n",
    "    delta = pm.Normal('delta', 0, 1e-5, value=0)\n",
    "    \n",
    "    m_mu = pm.Normal('m_mu', 0, 1e-5, value=0)\n",
    "    tau_mu = pm.Gamma('tau_m', 1, 0.01, value=1)\n",
    "    mu = pm.Normal('mu', m_mu, tau_mu, value=[1]*N)\n",
    "    \n",
    "    theta_t = pm.Lambda('theta_t', lambda mu=mu, d=delta: pm.invlogit(mu + d))\n",
    "    theta_c = pm.Lambda('theta_c', lambda mu=mu: pm.invlogit(mu))\n",
    "    \n",
    "    treat = pm.Binomial('treat', n_t_obs, theta_t, value=r_t_obs, observed=True)\n",
    "    control = pm.Binomial('control', n_c_obs, theta_c, value=r_c_obs, observed=True)\n",
    "    \n",
    "    treat_pred = pm.Binomial('treat_pred', n_t_obs, theta_t, value=r_t_obs)\n",
    "    control_pred = pm.Binomial('control_pred', n_c_obs, theta_c, value=r_c_obs)\n",
    "    \n",
    "    return locals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 10000 of 10000 complete in 6.2 sec"
     ]
    }
   ],
   "source": [
    "M_random_mean = pm.MCMC(random_mean())\n",
    "M_random_mean.sample(10000, 5000)\n",
    "M_random_mean.sample(10000, 5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEgCAYAAAAOk4xLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2YXHV99/H3ZwPGQATZxAooVkFqTQJoLYpIeAqsBkhS\ni9Ryl0ISrwpIRUWBlMSC3hRCiA/4XL3bBJoiKiKbGOQhCUhEFNAiZMGKCSsosZVsgoRCpNnv/cc5\nszmZzD7M7sycM7Of13XtlcyZc2a+u/Pwm3Pm9/0cRQRmZmZF1JZ3AWZmZv3xIGVmZoXlQcrMzArL\ng5SZmRWWBykzMyssD1JmZlZYHqRGSNIYSVdKelzSbyUtlbRH2TqzJT0n6cn05wlJE9Lr2iQtS5f/\nUNJ+6fLPSPpolbX8kaRr0zqelPSIpPm1+237vd9lki7t57qvSJpZYXm3pOMHuM2x6e/wlKReSXV7\nrkpalbmfJyV9fpi3c5mkf6t1fbUmqTP9PV+UNC3vemyH9DlUeq/4jaSfSXpPFdsfK+nJGtUyW9La\nWtzWSOyWdwEt4HzgKOBQ4Hnga8BngLMz6wTwjYiYW2H7dwKviogDJF0GfEjStcA04M1V1vJd4GfA\nQRHxnKS9gNdUeRvDEenPrldEnDPANv3fYMQ24ABJfww8PrLyBhYRJ2Tu548joree95e3iJgFIOlx\nBnkcrOEC+HZEnAkg6a3AKknrIuLn+ZaWD+9JjdyxwDcj4tmI+F/gI8DfSJqYWUfpTyUvAmMkjQHG\nAv8LfBb4WERsH2oRkk4CXg2cHRHPAUTE7yNiXdl6S9M9vy+ke3S/kXRo5npJmifpsXRv5yuSXlp2\nG++V9Mv00951wPjy30/ST9LrX5D0vn7KPkTSXZL+S9IPJB1c6Vcb5PeeK+nR9Hf5pqT2gdYf6KYG\nuI+7JF0oaUn6N9kg6diydX5K8ti/O7O3vH/ZOqekn4yflHSbpNeUXf/adG9uqqQH08dmaeb6r0r6\nXNk2X5d0eebyPpI+K6lL0kZJD0k6quo/RvJ3fTj9PZ6odq/ehm2n94qIuI/kw9OfVnUj0rsl3S/p\nd5Kul6TMdYM+RyQtBz4FHK4dR4DeNqLfbLgiwj8j+AEWAD8AJpA8ud4KbAGOy6xzerqsG7gbOKbs\nNr4IrAe+B7wfuHkYdVwKfGcI6y0F/ht4f3r5ZWXXfwR4GNiP5EPMEuALmevfSLLHeGx6+e3AJuAf\n+7m/O4G5FZZ3p7/vPunf7VPATyqs91qgF2ircN27gaeAP0kvfxL47jAfx4Hu5y7gEeDPMveztp/H\n4Lp+bv8twO+BI9LL7yPZ61WFGm4F9k2X7ZW5/s+B3wG7p5f3BraS7P2V1hkLdJR+D+ATQFc/NT0O\nHF9h+cHpYzwps2zPRr+2RuMPcBnwb+n/24AzgF8D+wxx+2NJPvhekT4XJgA9wLRqnyPAWZWe543+\n8Z7UyF0JrCYZfH4BnAtsA/6otEJEfB04OCJeS/IkvFnSgZnrz4uIg4C/AD5EMlBU61Ukb1iDCZI9\nv6+m9/1s2fVnA5dHxMZIDntdQvKGWvJXwIqIuCvd/l7gFgbZ4+mnjqsjYnMkr4h/BA5LD7sN1dnA\nNRHxi/TyJ4FjlX6vV0MBLI6In6aX11L5MOpAe8x/B/x7RPwIICL+heTN4ogK654REb9N1/t9XxER\nDwC/AU5JF/0V8IOI+FVmnW0RcXvsOGR5I/CGwX/FnTxD8hzukPSq9Hafq/I2bPj+Ij0U+zzJY31k\nRGyuYvvfRsQl6XNhE8kHrANKV1bxHKn2NV0XHqRGKCK2R8SlETE5Ig4meUN/GfBfZev9Lv13DcmA\ndsouNwYXAt8BzpV0r6RvSdpniKU8Q2ZgHMRAT/gDgE8pmQjyOPBD4PnMG/8rSfaCaip9E9wMvKKK\nzQ4APpKp9THgBerzPVz2Bfsi1b92DgBOK9Wa1ttOcoi2XM8At/M1YHb6/7PSyzuKTHxE0t2SfkCy\nl96mKiaeRMR/kxwReA2wRtJP5QkWjfSdiHgd8E1g74h4onSFpLdnDr89KemeIdzeTs/XWjxHGqmQ\nRTW5mSSfQn80wDovIfmU1EfSq4G/AX4MHB4Rbwe6gA8P8X5/THL8eM+qK97ZeuDMiHhd5qc9Ijam\n1/8aeF3ZNrsxvC/gdy/9R8lsx3bgif5Xr1jr/LJaJ0bEj4dRy2BGOsFgPfDlslr/KCK+VeXt/DvJ\n3uKRwOuBzrLrPwi8BzgtIo4iGciqFhG/iIgLIuINJHv/KySNHc5tWVWCHR+IPgAcLOnCvisj7o2I\nAzI/7xjGfdTkOdIoHqRGSNKekl6e/v9AYBHJ9zMvZNZ5jZKJEUg6kWQm4M1lN7WY5M3gf9P1SoeO\nXsrQdJJ8P/PPSmb1IWmipPI9tsF24T8NfFZS3+5/6fdLfQM4SdIR6SeyvwJmDHKble5TwFWSXilp\nN5K/263pp/ih3sZngEuVzIAq1TrUPc9qah1oeVYPyWSQ3SXtVnbY8UvA2ZLe2XeDw6g1IraQPNZf\nB66NZLJO1quB3wL/rWTyztXp8t2prOLvlT0cTfKhahsw5Ik8NmzZSRPPAv8HuCz7HB/p7TL058hm\nkkHyZelr/VUjrGFYPEiN3CHAg5J+A9wOfD4iyvtsPgo8KelXwEVAR+nwH4CkqcArI+IbwG3ARpI9\nincCn2MIIuJFkmnr24HHJP2a5LuT8mns/U4XT29nKcmT9oZ0VtfjwOWZ69eTfBL7TlrjsSTfSQ1Y\nXj/L/p1ksP4v4ECS72360y3pwbJa7ySZaHJNeujjcZKJHlWRtIrksGak91OpT6r8d6j0O10P/A/J\n3ubDJN8xlmr9OclgfnHm77q8wt7JUPbYvkbyRvP/Klz3KWAiyQeWlcBXST749PcG8/X0b9c3szL9\nkPNNSb9W0nNzHjC9woBotbfT6zOS2X3/F7he0suquI2Blg31OXIr8FNgA/BzktdawymdxWFmZlY4\n3pMyM7PC8iBlZmaF5UHKzMwKy4OUmZkVVtMGzK5evdozPsxGoWnTpuWehOD3n9rr73Ed8ew+SfdH\nxOEDXL+EZJry30fEyrLr/i4ivlZxw0Fs3rzZT5IhiggeeeQRACZNmkQma5L29nZ6egYKODArjn32\n2Sf3AQr8/lNrAz2udT/cFxFzSEJNK8ll3v1oEhHMmXMFHR0b6OjYwNy5V+K2A7Paiwi6urro6ury\na6yGhjVISfqHNAZ+KclpGpD0LiUn7btH0qlDuI3rgDdIulPSxzPL5yg5KdtDks4fTn22Q1dXF8uX\nL+b558/k+efPpLPzarq6uvIuy6yllD4MTp/ezfTp3f4wWENVfyeVRr3MBN4G7Ak8nEb4XAUcSRKf\nskbSdyM5cV1FEXFmeqjwuLKrlkXEkrQT/z6GmLhgla1fvx6YutOyo4+eSk9PNaHKZjaQRx55hDVr\n3sTWrWcAsHp1smzy5Mk5V9b8hjNx4gDggTTm/VlJvyOJ2HgVyZlhAV4O7M/wzqh6dJo3txXYY7CV\nbWAHHXQQydk2pqRLusikHHHRRRflUJWZ2dAM53DfBuDNktrScML90hy6R4FZEXFcRBwaEUMZoF5S\nIR7+GuACKueSWZUmT57MjBnbGDeum3Hjupk5cxubNu3Yi5o3b16O1Zm1hkmTJnH88Q8yfvwyxo9f\nxrRpP2PSpEl5l9UShjW7T9JFwGnAOpLTSkyR9A7gn0iCDJ+KiL/JrH8pcCpJavOnMssXk3zEfzwi\nzk2X/TNJaOt/AG+JiEonhfPsmioMNLvPrJkUeXafX2fDN9Dj2rQBsx6kzEafIg9SNny5TkE3MzMb\nLg9SZmZWWB6kaqgZm/kWLlyYdwlmZv2q2SAlabakZZK6JH1A0qOS/riKdV+TXvdAZr37a1VfvTVr\nssOiRYvyLsGsJTTjh9RmUMuA2SDpi/oP4GXADcCbgF8Ncd03k5yOvCkf3VKyQ0lnJ3R1rWXKlCkD\nbGVmraD0IXXNmjcBMG1aJ//6r//gGX41UOsU9N+m/24F9h3k9qtZt9Da2/ehPNUBkmSHLKc8mLUm\nJ07UTxEHhjYASXvQVIkTAcxn12SH5JOUBygzs+rVe5Cq5tBdad370ibfrVVun5uens1EBLNnb2PV\nqm4ATjxxG0uWbPbuvtkokCRO3MyaNcnlJHHipHyLahFu5q2hZuw4X7hwoaORrGkUuZm3GV//RZFb\n4oSkOyssjog4fqS3XcRByszqq8iDlA3fQI9rXQ/3VTgNh5mZ2ZAN96SHf1frQszMzMoNt5nXp303\nY2gNnG7yHN38+I9M1YNU2WnfF6TLKp7yPZsYMVh6RDWJFWZFMJRThvu04qObH/+RG+75pO6PiMMz\nl3ePiBdLp3yPiMPK1yvfpsJtngW8HthC8l3ZS4EHI6Kz0vr+4rI2RvvsvqQRu7mNph68Zps40dXV\nxfTp3X1NvuPHL+N733utm3zLNGLiRK1O+d4yKRTNYtGiRYUcpFph8GiUZv9bjaZB1qo33EHgJZLa\nIqI3vXwNydl0XwP8dWa9Jk2PsLw1wxvXrnltP9slr20o61jrcpPvyA33cF/ptO8bIuID/Z3yXdKX\ngedI9ozeExH9pq2mh/vGpxdLe1KPRcRNldb34b7aaG9vp6enJ+8ymtZQGjjd5Fk7zXa4D/z4D4VP\nH2/98iBlzaQZBykbXG7NvJXUM4XCzMxaS8MHKadQFMtFF12UdwlmZv3y4T4zaxo+3NeaBnpca3b6\n+KGQdFna8PuXmWVTJd0n6eqydS+X9LikkxtZ41BkO8h7e3vdTW5mVXEKxdA1+nBfAPMi4pbMsrHA\nlcCRO60YsUDSi40sbiiyU4ojgr33vpLf//50wKeMNrPB+VTz1RlwkJI0GzgBeDPwReCDwDsj4glJ\nD0TEn6frDZgmUX6z2QsRsUrSMVVX3mA7N0wu7vvfc8/9bd//OzvPYMKE5P/N0OdjZo3nU81XZ7A9\nqQAeJ+l/ehlwA8mA9QRNctbcoap11/5At+cBzMxsaIZyuG9URBUNdeDY9XDfDZnDfc2XJjDas/vM\nGs0pFNUZyYBTy8ijpnlXl8SSJZf0dZC/8Y3X8eijjwIwadJJTTVAQXGz+8xaVfl7SDO+bzTScAap\n0mG++9J4pK2M4NCfpIuB6cC+kvaKiLOzVw/3dutJ0k7Hj30s2cyqUf4eYv0bcJCKiGsHuO7czMXL\nhnh/W4CPSxpbyuSLiKuAq8pXlHQ5cBJw1xBv28zMWoybeUc5Z/dZM3Ezb2sqTDOvmZlZNRo+SJWn\nTkj6Snoq+u9LOjCzXmETJwYSEaxbt47Ozk7WrVtX+G5yZ/eZ1ZbTJGqr4Yf7JF0K3F+WOoGk44HT\nst91pes+EBEry2+niLvbpenpK1dOYvv2NsaMWcXJJ+/LkiWXePaOWQ0U/XCfT3I5PIU6VUeqUkHP\nAn9odCG1sqN5d0caxfbts1m+HKdQmI0STpOovSI15s4lOQ19oY0kmaJ8Ww9aZmYDK8QgJWkG8J8R\n8fO8axnMQAOLD/eZjW5Ok6i93AcpSW8BjomIj+Vdy0iVOsm7urpYv349Bx10NpMnT/YAZTZKOE2i\n9oowBf1bwOHpDL/PlV3XdI+uJKZMmcKsWbOYMmVK4Z+gCxcuzLsEs5ZSSpPwB9TayGOQKqVO/CVA\nRBwYEcdExHERcX5ppTRxYibwXA41jhqLFi3KuwQzs345cWKUc+KENZOiT0G34SniFPSWFRGZ49GT\nvLtvZjYCRUicuFzSGkmrmj1xojS7r6NjAx0dG5g790p3nJuNMk6cqK089qQCmFdKnIiIBQCS3gFc\nDJxdWi7pxRzqG5ZKzbydnUkjr/uhzEaHXRMnOp04MUJFSpw4Ani00YWM1FCaeyutU5SBy9l9ZrXj\nxInaK8R3UpLuBiYCU/OuZTADD0oBzAempJe7gMspH5OLMkABPiuvmRVaIQapiDha0luB64BCfwc1\nWOLE7NnbWLWqG4ATT9zGkiWbvatvNko4caL2CjFIpX5LseqpmiSWLp3v2X1mo5QTJ2ov90FB0jdI\nDvX9Afj78qsbX9HIlLrNzWx08ntAbeUxSJUSJ8ZGxE0R8d5KK6WJEycBdzWyODMzK46G90lFxDUR\ncURE3DTIegsi4s8i4q4GlTYqObvPzIos92bedNlYSb+SdF5mWeGbeVuhac/ZfWZWZLk386bOAX6S\nXpesVPBm3lLT3h13HApAR4eb9sxGM0ei1UfuzbyS9gBOJDllx/ic6qlaV1cXy5fvnC7R1bWWKVOm\nDLCVmbUiJ03UT+6z+4DzgS8Ar8y7kGocffSufcelZUVq1jWz+nPSRP3ketJDSXsDR0XErTTZdPO7\n714LXAJcn/7MJ3O00szMaiCvQar0bv4O4KWSvk7yvdQcSU3x0WPy5MnMmLGNceO6GTeum5kzt7Fp\n0+am24tydp/ZyCVJEw8yfvwyxo9fliZNTMq7rJbQ8JMeSroUeCAiVpYtPwvYMyK+lFl2GXB/+bpQ\njJOO+YtSs8Yq8kkP/X4wfEU76eFOzbylhRFxbXalZmjmdWe5mZX4/aA+fPp4M2saRd6TsuEb6HHN\ndeKEmZnZQHJPnJC0VNK9ku5Mv5cqrVf4xAmA3t5eOjs76ezspLe3N+9yzCxnrZBEUyR57EmVEidu\nylx+b0Qcl/1eKj2t/NIc6huy3t5eDjnkTObMgTlz4NBDz2q6gcrZfWa1U2rqnT69m+nTu5k790oP\nVCOU1+G+8uOPhTjOXK2JEyewceNKYDYwm6ee+i4TJ07IuarqOLvPrHayTb1bt57B6tWH9c34s+Ep\nQuLEs8D1knqAj0TEL/MuaCADnz5+13WarW/KzKxIch+kIuJ8AElvAq4G3p1vRQPLDjqlw30bN54G\nwP7738hDD11LW5vno5iNRj59fO3lPkhlvAAUNvW8kra2Nh5++DpWrFgBwIwZHqDMRjOfPr72ch+k\nJN0A7Edy2O+88qsbX1F12tramDVrVt5lmFlBuKm3tnJPnIiIv660UjMkTrQCZ/eZWZE5ccLMmoYT\nJ1qTEyfMzKwpFSFx4tVp2sRaSZ/OrOfECTNrKk6bqL08vpMqJU7ckl5eDMyPiB/utFLEAkmFnu23\n6xT0szwF3WyU8ink6yOv2X0CkDQGOKh8gGoGOxp2d5zq6qmnZjNx4s7ruZnXbHTwKeTrI+8p6K8g\nOTPvzcBewOcj4js511TRUJImhrNd3oPYwoULmTdvXq41mJn1J+9BahPwDHAqMAa4R9KtEfF8vmXt\nqtJg0gqJE4sWLfIgZVYDTpuoj1wHqYh4UdKTwL4R8RtJ2/Ksp1pOnDCzEqdN1Efee1IAFwNfk7Q3\n8M2yvajCP8JOnDCzEqdN1F4REieeIEmW2IkTJ8zMzIkTo1x7ezs9PT15l2E2JE6caE1OnLB+ObvP\nzIos18QJSXulaROln2cy6zVF4kS1ipZQ4Zl9ZrXjxInaK0LixHEAkg4FPti3UhMkTlTLCRVmrcuJ\nE/WR17tjpUftfODzjS6kkVasWMHGjSuB2cBsnnrqu33T182suWUTJ7ZuPYPVqw/rm45uw1eIj/CS\nJgAHRMRDeddST3PmzK64bLhpFmZmra4QgxTwfuCreRdRb08/vYn99jsZWAosZf/9T+HppzflHo1k\nZiOXJE48yPjxyxg/flmaODEp77KaXu7NvJJ2A04BpuZdS70VMaHC2X1mteHEifpoeJ+UpEuBByJi\nZXr5PcDrI2JhhXUvA+4vrZvlPoXacJ+UNRP3SbWmgR7XIiRO3FhpJSdOmJmZEydGOe9JWTPxnlRr\nKtqeVMuLiMxx6Uk+Lm1mNky5Jk6kl8+U9GNJ90g6LrNeUyZOlBr6Ojo20NGxgblzr3Tnudko4cSJ\n2itC4sTHgDcDewK3AW+H5k2c6OrqYvnyxX2XOzuhq2stU6ZMybGq/jm7z6w2nDhRH3kd7ss+ao8A\nxwD7Aj/Kp5zaSJpyd51Jf/TRO5YVrSfK08/NaiObOAGwenWyzOeXGpkifCd1O/Bh4CXAF3OupQYC\nmA+U9py6gMsBFW6AMjMrulw7SSUdCJwSETMj4l3AhZLG5VnTSPT0bGbTps3MmLGNceO6GTeum5kz\nt7Fp02YPUGYtzokT9ZH3ntSYUg1KDtyOI9kVaVqSWLp0vmf3mY0yTpyoj1wHqYh4TNKPJN1Cslf3\nxYh4IbNKUz7Cknwc2mwU8mu/9vI43FdKnPhLgIi4IiJOioh3RcTS0kpp4sRM4Lkcahw1Fi7cJY3K\nzKwwnDgxyjlxwpqJEyda00CPaxGaec+WdK+kOyQdnFmvKZt5s9zYZ2Y2Mrk280raA5gTEUdImgh8\nGTgNmreZtyTb2Ld9ey9HHPFvfOtbV+R+ag4zawzHo9VG3s28AnaXNJbku6p9Je0eEU07OJVMmNAO\n7EieuPPOM5k4ETZt6vGT1azFOX2idvKe3fecpCuA7wHPAvsALwd+l2ddIzHYqeCTwSvh3imz1uT0\nidrJu0+KiPg28G0AST+NiKYdoCAZeJKBqv/kidJ6ReDsPjMrssJ8QSLpJODBvOuohZ6ezfT0bOHp\npz/KMcfcxrhx3ey5558ya9ZFbNrUU5gBCpzdZ1YPTp+ondz3pCT9C/AGYCtwRvnVja+odtra2rjp\nps9kvjw92cekzUYBp0/UThFOH/++Siu1yunj3YFuNjr5tV8bbuY1s6bhZt7WVKhmXjMzs6Gq2yBV\nIVliqqT7JF1dtt4JktamP8dnljdV4kRvby+dnZ10dnbS29ubdzlD5uw+s9pwwkx91O1wn6RLgftL\np4mXdALwMuDIiLgwXdYGrAVOSDe7DTgm0qLS23ggIlaW336Rdrd7e3s55JAz2bjxNAD23/9GHnro\n2qZIl3B2nzWToh7u27V592du3q1Cnof7+u44IlYB5e+GBwO/iIjnI+J5YD3w+jrXVHMrVqxg48aV\nwGxgNk899V1WrFiRc1Vm1ijZ5t2tW89g9erD+mb22cjkPQW9Hdgi6TPp5WeACcBj+ZVUG3PmzC5U\nP5SZWTPKe5DaRBKD9AGSva4vAU/nWtEwzJgxg/32O3mXw32el2I2OiTNuzezZk1yOWnePSnfolpE\nowep8uOO64E/yVw+OCJ+2cB6aqKtrY2HH76u7xDfjBnN8X2UmdWGm3frp2GDlKSLgekkSed7RcTZ\nEbFd0ieAO9LVLivfrFH1jVRbWxuzZs3Ku4yqObvPrDbcvFsf9Zzd9yHgdGBRRNw0jO1LiRMXRMRd\n5dcXaXafmTVGUWf32cgM9Lg6ccLMmoYHqdbkxAkzM2tKRUic6G950yRONGvahJnVjhMn6qOeEycC\nmFdKnADGAlcCR5atV3F5RCyQVPjTyO+aNnFW06RNmFlt+HTx9ZN34kS/y5vFxIkTmjptwtl9ZiPn\nxIn68cf9EUhOE7+rOXNm096+T7/XF8miRYvyLsHMrF95J040tZ6ezU0dLmtmteHEifrJO3FisOWF\n57QJM3PiRP3kmjgx0PLSZo2qbySaNW3CzGrHiRP1Uc9BagvwcUljI+KmiLgKuKp8pf6WZxIn7qpj\njWZmVmB1G6Qi4hrgmhFsvwBYULuKrBJn95lZkTkWqQAiInMse5KPZZv1w7FIrSmXWKQqEie+IulO\nSd+XdGBmedMkToxEqQmwo2MDHR0bmDv3SnermzUhJ07UR+6JExFxDoCk44ELgXPT5U2RODFSEya0\nA4v7Lnd2QlfXWqZMmZJfUWZWFSdO1E+9Z/ftlDgh6ZgB1n0W+EOd6ymU/pp9jz56at//fQp6s+LL\nJk4ArF6dLPNsv5ErUjPvXEYw0aLZ7BigApgPlPacuoDLAXmAMrNRrxBdp5JmAP8ZET/Pu5ZG6enZ\nTE/PZjZt2syMGdsYN66bceO6mTlzG5s2bW7YAOXsPrORSxInHmT8+GWMH78sTZyYlHdZLaGeZ+a9\nFHggIlZmlh0LnBwRF2aWvQU4PSI+NpTbKGml2TV5zu5rb2+np6dp831tlCny7D7P0h2+gR7X3BMn\ngG8BT0q6E3g4Is7Pbtao+vLkTnWz5ufXcX0UIXHiwF03deKEmZm5mXfU8+E+ayZFPtxnw1f0Zt7L\nJa2RtKpozbxuzjMzy1c9Z/eVmnlvSi+Xmnl3XiliQUQcD1wKXJxdDiytY339igjWrVvHqafOY/r0\nbqZP727ZJAhn95kVkz8kJ4rUzHsE8Gid6xlUqXN8+fLFwI6m2s7OM5gwofWaa+fNm5d3CWZWxgkW\nOxSimVfS3cBEsqNCDnY02C4ewjqtN2CZWTE4wWKHQgxSEXG0pLcC1wEN/w6qcjxR/0kQ4AHKzKwR\ninL6eIDfktOgmR1wsrvZEW/k8MNv55OffB+TJ78LaUse5ZnZKJMkWNzMmjXJ5STB4qR8i8pJ7s28\nkr5BcqjvD8Dfl2/WqPr67lBiyZJLMp3jnx6Vx4HNLD+7vg+dNGrfh+oZi/Qh4HRgUWaGXzXbl5p5\nL4iIu8qvd59CbSxcuNCTJ6xpuE+qNQ30uLqZd5RzM681Ew9SrSmXZl4zM7ORyj1xIr1urKRfSTov\nsyz3xIlybq4zM2us3BMnUucAP0m3STbOMXGiktKsv46ODXR0bGjZBAozGx5/iK2P3BMnJO0BnEhy\nyo7xda5n2Lq6utIUikRnJ3R1rWXKlCkDbGVmo4ETIuqnCM285wNfAF6ZdyGV7Gj03TUM4+ijk2XN\n3Njr7D6zkXNCRP3kOnFC0t7AURFxK4U/wWEAlwDXpz/zyRydbFqefm5mRZZ34sQ7gJdK+jrwOmA3\nSXdFRFeD6+pXaS8pIpg9exurVnUDcOKJ21iyZLN3583MCRF1lGviRETcAtySXn8WsGfZAFWYEUAS\nS5fOz3SAT/IAZWaAEyLqyYkTZtY03Mzbmpw4YWYtwYNUa3LihPVr4cKFeZdgZtav3BMnJC2VdK+k\nO9PvpUrLC5c4US+lJsB169axbt26hjYDLlq0qCH3Y2Zu+B2Oek6cKCVO3JJeLiVOHFlhvfdGxBM7\nLYxYIOnFOtZXCKUmwNWrD+OFF24FpjFu3EvcDGjWYtzwOzy5J06UrzfaTJjQzo7T1f8tAFu3Qmfn\nGTzyyA/VIsfzAAANL0lEQVTcDGjWItzwOzxFSJx4FrheUg/wkYj4Zd4FNULlU9bvbOrUo/r+38yp\nFmZmw5X7IBUR5wNIehNwNfDufCtqjGyTcOXDfT/zoQCzFuKG3+HJO3Ei6wWg5b+DKpdtAow4p29Z\no5oBnd1n1hhu+B2eXBMn0uU3APuRHPY7r3yzRtWXJ0m5HZd2dp9Z4+T5Wm9W9RyktgAflzQ2Im6K\niKuAq8pXioi/rrRxJnHirjrWaGZmBebECTNrGk6caE1OnDAzs6ZUhMSJV6dpE2slfTqzfFQkTkQE\n69ato7Ozk3Xr1rkL3azFOXWiOvXckyolTpQS0EuJE+UWA/MjYmpEXNC3ccQCYGkd68tdafr5ccc9\nxpw5cNxx/8ycOVc09Inr7D6zxim95qdP72b69G7mzr3SA9Ug6n24b6fECaBnpyulMcBBEfHDOtdR\nSBMmtLN8+WK2b58LzGb79mUsX35m3xTVRnB2n1njZFMntm49g9WrD2vo670Z5d3M+wqSM/PeDOwF\nfD4ivpNzTQ3Rf+LEoUyduuOSkybMbDTLe5DaBDwDnAqMAe6RdGtEPJ9vWfXX07O5b9d/5cpJbN/e\nxpgxqzj55H1ZsuQSN/mZtSCnTlQv18SJiHhR0pPAvhHxG0nbGlxPrkod6F1dXaxfv56DDjqbyZMn\ne4Aya1FOnahe7okTwMXA1yTtDXyzbC+q5R89SUyZMoUpU6bkXYqZNYBTJ6pThMSJJ0iSJXbixInG\ncHafmRWZEyfMrGk4caI1DfS45j1xYlSIiMwx6Ek+Bm1mNkS5Jk5I2itNmyj9PJO5riUSJ3p7ezn1\n1Hl0dGxw855Zi3OaRO3Vc0+qlDhxS3q5lDhxZN8KEb8HjgOQdCjwwcx1CyQ19fmldvRCfbVvWWfn\nGUyY4P4ns1ZTailZs+ZNAEyb1ukTl9ZAvQ/37ZQ4IemYAdY9H/hcnetpiKGcGj67jgcss+aXTZMA\nWL06WeaZfCNTiBR0SROAAyLiobxrqb0ALgGuT3/mp8t2GMqgVi/O7jOzIivEIAW8n+wxsSbX07OZ\nnp7NbNrUw8yZF7Lnnm9k3Lhujj32dp5++gJ6erb0rVP6yYuz+8xqI0mTeJDx45cxfvyyNE1iUt5l\nNb1cEycAJO0GnAJM3XX15rZzd/nrmDTpHB+fNmtRTpOojyIkTvwFsCIieitt1qj66sXd5Wajh1/v\ntVeExIkbK23sxAkzM3PixCjX3t5OT0/P4CuaFYATJ1rTQI9rEU4ff6akH0u6R9JxmeWFaeZt5QY9\nZ/eZWZHl2syb+hjwZmBP4Dbg7VCcZt5Sg94ddxwKQEdHazXozZs3L+8SzFqOo9BqpwjNvI8AxwD7\nAj+qcz1V2dG/tLhvWWcnTowws345eaK2ihAwezvwYeAlwBfzLKSaptpK63rgMjMnT9RWroOUpAOB\nUyJiZnr5bkmrinX6+CBJiSidlLALuJxKs+Pb2/fxQGVmVkN5N/OOKdWgZF94HOWZQQ1UaYCJCGbP\n3saqVd0AnHjiNpYs2exddzOrKEmeuJk1a5LLSfLELud1tSHKtZk3Ih6T9CNJt5DMNPxiRLyQ3axR\n9fVHEkuXzm/ZL0EXLlzoyRNmNeTkidqqW5+UpA8BpwOLIuKmYWxfaua9ICLuKr/efQq14T4paybu\nk2pNAz2ubuYd5TxIWTPxINWacmnmNTMzG6kiJE6cLeleSXdIOjizvKGJE62cKmFm1qzquSdVSpwo\nfR9VSpzoI2kPYE5EvJ3k+6sr+jaOWAAsrWN9pfth3bp1nHrqPKZP72b69G7mzr3SA5WZDZs/9NZO\n3okTAnaXNJYkNX1fSbtHREPikEqd4cuXLyZ7OqvOzjNGTaqEs/vMasuJE7WVazNvRDwn6Qrge8Cz\nwD7Ay4Hf1eP+KidKLK6wrPL6rThoefq5WW05caK2co9FiohvA98GkPTTiGjgANVXBUNJlSjdRisO\nVmZmRZR34sSOK6STgAfrdcf9pUmUdssj3sjhh9/OJz/5PiZPfhfSlnqVYmYtzIkTtZX76eMl/Qvw\nBmArcEb5ZnWuqawz/NM+bmxmI+LEidpy4oSZNQ0387amXJp5I+KaiDhiOANUuv2CiPizSgOU1c7C\nhQvzLsHMrF+ORRrlHItkzcR7Uq0plz2pCokTX5F0p6Tvp+eRKq13gqS16c/xmeUNTZzIS29vL52d\nnXR2dtLb25t3OWY2TG7grY96TpwoJU7cAhAR5wCkA9GFwLmS2oBPACek29wm6c5ILJDUkKbevPT2\n9nLIIWeyceNpAOy//1k89NC1tLU5UtGsmbiBt34aljiR8Szwh/T/BwO/KJ2JV9J64PXAY3WuqxAm\nTpwArOy7/NRTs5k40X1YZs3GDbz1k0cz71zgmvT/7cAWSZ9JLz8DTKCFB6mBm4p3XccDlpmNZg0d\npCTNAP4zIn6eLtpEEoP0AZK9ri8BTzeypkbLDjq7Hu67seGH+5zdZzZybuCtn3r2SV0KPBARK9PL\nbwFOj4iPZdYZA9xN8p2UgDsi4h393UZWq8yu6e3tZcWKFQDMmDHD30eZDaDIs/siItPAO8nfR1Vh\noMe1kXtS3wKelHQn8HBEnB8R2yV9ArgjXeeysm1a/lFua2tj1qxZeZdhZiMkyd9B1UE9B6ktwMcl\njY2ImyLiwEorRcTtwO3lyzOJE3fVsUYzMyswN/OaWdMo8uE+G75cmnnNzMxGqgiJE1Ml3Sfp6rLt\nmzJxotkSJJzdZ1YbTpyoj3rP7ru/lDiRWX48cFpEnJtePgF4GXBkRFxY4TaaZnZfEaaUV8vZfdZM\ninq4b9fEiZ85caIKec7uGyxxgohYJemYOtfREE6QMBudnDhRP3knTrSEwVIknCBhZjY8eSdOtITS\nwNOMh/vMbOScOFE/jTx9/FuAY7KJE9mrG1VHPbW1tfHww9dlEiQ8QJmNBj5lfP3kmjgBIOliYDqw\nr6S9IuLszDZN9yg3W4KEs/vMasOJE/VRz9l9HwJOBxYN5xTymcSJCyqdQr6Is/vMrL6KOrvPRmag\nx9WJE2bWNDxItaaiBMw2FScam5nlrwiJE/0tzy1xore3l1NPnUdHxwamT+9m7twr3UFuZoNy6kTt\n5Z44MdDyPBInIoIJE9p3Wb527Q/8pahZzop8uM+pE8OXZ8DsoIkTQ1jeUJUGKICpU4+ivX2fIZ3+\nvZk4u8+sNrKpE1u3nsHq1Yf1fWVgw5dHE89c4MtVLG+YHQNQAJcA16c/89NlrWfRokV5l2Bm1q9C\nJE4UJYmip2dz3y776tWH0dvbzdve9gtuvPFTtLVtybM0Mys4p07UR+6JE4MkUTTczp3jr2PSpHN8\nTNnMBuXUifrIPXFigOWQU+KEO8fNbDj83lF79RyktgAflzQ2Im6KiAMrrdTf8kzixF31K9HMzIqs\naRMnVq9e3ZyFm9mITJs2LfdjaH7/qb3+HtemHaTMzKz1+TwSZmZWWB6kzMyssDxImZlZYXmQMjOz\nwmqpQUrSWEm/knRe3rWU6y/tPU+STpC0Nv05Pu96Sor4typX8Ofaq9O/31pJn867nixJZ0r6saR7\nJB2Xdz21JGmqpPskXV3E+5S0VNK96XPjrDrXVbPXcKudT+oc4CcUMGgvIs6BvrT3C4FzB96iviS1\nAZ8ATkgX3SbpzijAdM+i/a36UdjnGrAYmB8RP8y7kAo+BrwZ2BO4DXh7vuXU1FjgSuDIgt5nAO+N\niCfqW1JtX8MtsyclaQ/gRKCTnJIqhqgQae/AwcAvIuL5iHgeWA+8PueayhXlb7WTIj/XJI0BDiro\nAAXwCHAMcArwo5xrqamIWAX0FPw+G/183ek1LOnYavfimm5PStKJwEVliz8KTAe+ALyy4UVl9Fdf\nRDyU/n8ucE1jq6qoHdgi6TPp5WeACcBj+ZW0i6L8rcqdTwGea/14BfBSSTcDewGfj4jv5FxT1u3A\nh4GXAF/MuZbR5lngekk9wEci4pcNuM+5wDWSDgM+C7yc5Pk5G/indJAdUNMNUhFxB3BHdpmkvYGp\nEXFV+svnplJ9JUVJe09tInnCfIDk09WXgKdzrSijYH+rPulz7aiIWJj3c60fm0g+cJwKjAHukXRr\nurecq/S7iVMiYmZ6+W5Jq4pQ22hQykWV9CbgauDd9by/Cq/h4yQdA7w2Iq4d6u003SDVj3eQjM5f\nB14H7JZ+v1KYM44VLe2d5PDen2QuH9ygT1aDKuDfKqvQz7WIeFHSk8C+EfEbSdvyriljDOl7jpJ4\n8HEU8zu9kcjj8G+19/kC8GI9Cimp5Wu45WKR0uOde0bEl/KuJUvSBuBJoJdd095zIakD+Mf04ifS\nvcDcFfFvVUmBn2uvAb4C7A18MyIKc8hU0iXAUSTfh98QEUvzrah2JF1M8rXDvsD3I+LsvO5T0mnA\n/0TEysy6NwD7kRz2Oy8iflXHumr2Gm65QcrMzFpHy8zuMzOz1uNByszMCsuDlJmZFZYHKTMzKywP\nUmZmVlgepMzMrLA8SJmZWWF5kDIzs8L6/7mPJuMy0AXjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d847668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.Matplot.summary_plot([M_random_mean.delta, M_random_mean.m_mu, M_random_mean.tau_mu, M_random_mean.mu])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def random_treat():\n",
    "    \n",
    "    mu = pm.Normal('mu', 0, 1e-5, value=1)\n",
    "    \n",
    "    m_delta = pm.Normal('m_delta', 0, 1e-5, value=0.5)\n",
    "    tau_delta = pm.Gamma('tau_delta', 1, 0.01, value=1)\n",
    "    delta = pm.Normal('delta', m_delta, tau_delta, size=N)\n",
    "    \n",
    "    theta_t = pm.Lambda('theta_t', lambda mu=mu, d=delta: pm.invlogit(mu + d))\n",
    "    theta_c = pm.Lambda('theta_c', lambda mu=mu: pm.invlogit(mu)*np.ones(N))\n",
    "    \n",
    "    treat = pm.Binomial('treat', n_t_obs, theta_t, value=r_t_obs, observed=True)\n",
    "    control = pm.Binomial('control', n_c_obs, theta_c, value=r_c_obs, observed=True)\n",
    "    \n",
    "    treat_pred = pm.Binomial('treat_pred', n_t_obs, theta_t, value=r_t_obs)\n",
    "    control_pred = pm.Binomial('control_pred', n_c_obs, theta_c, value=r_c_obs)\n",
    "    \n",
    "    return locals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 10000 of 10000 complete in 6.8 sec"
     ]
    }
   ],
   "source": [
    "M_random_treat = pm.MCMC(random_treat())\n",
    "M_random_treat.sample(10000, 5000)\n",
    "M_random_treat.sample(10000, 5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEgCAYAAADhUed1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucXHV9//HXOwFjZBHZpBWoUCVSaxJF688WlXAnbZAE\nL6XWnykk6aMi2mJViBGw4O9HZQmI0qql2ppgI17QyCYGL+SCIoqClkt2sWIgSgVbyYZL+IVIs5/f\nH+fMZnYyOzuzcznnTN7Px2Mf2XPmzJzvZGbnM+d7vt/3UURgZmaWN5OyboCZmVk1LlBmZpZLLlBm\nZpZLLlBmZpZLLlBmZpZLLlBmZpZLLlAdJmmypMslPSjpV5JWSnpOxTaLJD0l6aH05xeSpqW3TZK0\nKl3/PUmHpus/Kul9DbbltyVdl7bjIUmDki5q3bMdc7+rJF0yxm3XSlpQZf1WSSfVeMwp6XN4WNKw\npLa9tyWtL9vPQ5L+cYKPc6mkf2t1+1pNUn/6PJ+RdHLW7el26fui9Pf/S0l3S/rTBu5/gqSHWtSW\nRZJubcVjTcR+We14H3YecCzwcmAn8Gngo8A5ZdsE8MWIWFLl/n8M/E5EHC7pUuDdkq4DTgZe2WBb\nvgbcDcyIiKckPRc4osHHmIhIf/a+IeIdNe4z9gNG7AIOl/S7wIPNNa+2iDilbD+/GxHD7dxf1iLi\nDABJDzLO62AtEcBXIuIsAEl/CKyXtDkifpJt0zrLR1CddwLwpYh4MiL+B3gP8DZJ08u2UfpTzTPA\nZEmTgSnA/wAfA86PiN31NkLSacALgHMi4imAiHgiIjZXbLcyPeL7eHok90tJLy+7XZKWSbo/Pcq5\nVtKzKx7jLZJ+ln4j/CzQU/n8JP0ovf1pSX85RrNfJukWSf8l6buSjqr21MZ53ksk3Zc+ly9J6q21\nfa2HqrGPWyRdIGlF+n/ygKQTKrb5Mclr/8ayo+TDKrY5Pf32/JCkb0o6ouL2F6ZHcXMk3ZW+NivL\nbv+UpH+ouM/nJV1WtnywpI9JGpD0iKR7JB3b8H9G8v96b/o8ftHo0byNMurvPyJ+SPJl6PcbehDp\njZLukPRrSddLUtlt477uktYAHwFerT29OX/U1DNrVET4p4M/wMXAd4FpJG/CPwQeA04s2+at6bqt\nwHeA4yse4xPAFuDrwNuBGyfQjkuAr9ax3Urgv4G3p8sHVtz+HuBe4FCSLzwrgI+X3f5SkiPFE9Ll\n1wDbgL8bY3+bgCVV1m9Nn+/B6f/bR4AfVdnuhcAwMKnKbW8EHgZ+L13+P8DXJvg61trPLcAg8Adl\n+7l1jNfgs2M8/quAJ4Bj0uW/JDnaVZU2fAM4JF333LLb/xfwa2D/dPkgYAfJUV9pmynA3NLzAD4E\nDIzRpgeBk6qsPyp9jWeWrTug039b3fIDXAr8W/r7JGAh8J/AwXXe/wSSL7IfTl/facAQcHKjrztw\ndrX3bqd+fATVeZcDG0gKz0+Bc4FdwG+XNoiIzwNHRcQLSd6sN0o6suz2d0XEDOANwLtJikSjfofk\nw2o8QXLE96l0309W3H4OcFlEPBJJV9eFJB+mJX8GrI2IW9L7fx+4iXGOdMZox5URsT2Sv5y/A45O\nu9rqdQ5wTUT8NF3+P8AJSs/jtVAAV0XEj9PlW6nedVrrSPmvgM9FxO0AEfGvJB8qx1TZdmFE/Crd\n7omRRkTcCfwSOD1d9WfAdyPi52Xb7IqIb8WebsovAy8Z/ymO8jjJe3iupN9JH/epBh/DRntD2qW6\nk+T1e21EbG/g/r+KiAvT13cbyRemw0s3NvC6N/p32lIuUB0WEbsj4pKImBURR5F8mB8I/FfFdr9O\n/91IUsxO3+vB4ALgq8C5kr4v6QZJB9fZlMcpK4rjqPWHcTjwESWDPh4EvgfsLPvQfz7J0U9LpR+A\n24HfauBuhwPvKWvr/cDTtOe8W/kf9jM0/rd2OHBmqa1pe3tJumUrDdV4nE8Di9Lfz06X9zQy8R5J\n35H0XZKj80lqYJBJRPw3SU/AEcBGST+WB1M066sR8SLgS8BBEfGL0g2SXlPW5faQpNvqeLxR78FW\nvO6dkKvG7KMWkHz7vL3GNs8i+SY1QtILgLcBPwBeHRGvAQaAv61zvz8g6Vs+oOEWj7YFOCsiXlT2\n0xsRj6S3/yfwoor77MfETrbvX/pFyajGXuAXY29eta0XVbR1ekT8YAJtGU+zgwm2AP9U0dbfjogb\nGnycz5EcJb4WeDHQX3H73wB/CpwZEceSFLGGRcRPI+K9EfESkqP+tZKmTOSxjGDPF5x3AkdJumDk\nxojvR8ThZT+vm8A+WvK6t5sLVIdJOkDS89LfjwSWk5yPebpsmyOUDIJA0qkkI/5urHioq0g+CP4n\n3a7UXfRs6tNPcj7mn5WM3kPSdEmVR2rjHeJfDXxM0kj3QOn5pb4InCbpmPRb258B88d5zGr7FHCF\npOdL2o/k/+0b6bf3eh/jo8AlSkZFldpa7xFnI22ttb7cEMnAj/0l7VfR1fhJ4BxJfzzygBNoa0Q8\nRvJafx64LpKBOeVeAPwK+G8lA3WuTNfvT3VVn1d5FzTJF6pdQN2DdmyU8gESTwL/G7i0/H3b7ONS\n/+u+naRAHpj+/f5Ok21oiAtU570MuEvSL4FvAf8YEZXzaN4HPCTp58BSYG6pyw9A0hzg+RHxReCb\nwCMkRxJ/DPwDdYiIZ0iGpu8G7pf0nyTnSiqHqo85JDx9nJUkb+4vpKO3HgQuK7t9C8m3ta+mbTyB\n5BxUzeaNse5zJIX6v4AjSc7TjGWrpLsq2rqJZFDJNWnXyIMkgzoaImk9SVdmpPupNg+q8jlUe07X\nA/+P5CjzXpJziqW2/oSkkL+/7P91TZWjknqO1D5N8oH0L1Vu+wgwneTLyjrgUyRfesb6IPp8+n83\nMoIy/YLzJUn/qWT+zbuAeVWKodVn1N9cJKP4/i9wvaQDG3iMWuvqfd2/AfwYeAD4CcnfT8coHalh\nZmaWKz6CMjOzXHKBMjOzXHKBMjOzXHKBMjOzXHJYLLBhwwaPFDErsJNPPjmzxAN/fjRvrNcvN6P4\nJN0REa8eZ5sVJMOU/zoi1lXc9lcR8emqdxzH9u3b8/GfUENEMDg4CMDMmTORRG9vL0NDtUIEzLrf\nwQcfnGkcTxE+P/Ks1utXqC6+iFhMEl5aTUfH53dSRLB48YeZO/cB5s59gCVLLicvXyzMbG8RwcDA\nAAMDA/5bbUJLC5SSi1utSiPc36nksgZjhnlK+oCSOPiVJJdgKK3/EyUX47tN0pvr2O9ngZdI2iTp\ng2XrFyu52No9ks5r8ullZmBggDVrrmLnzrPYufMs+vuvZGBgIOtmmVkVpS+U8+ZtZd68rf5C2YSW\ndvFJOpsk7+sxkvNbzwbuiojK/C/SWJfVwOuAA4B7I+KFaVjhj9P1u4CNwKmRXJAOJVdivbNKF99e\nXYSS9o+IZ9LZ9z+MiKOrtTvvh+i9vWMl3AgIhoYaCTk26y556+IbGBhg3ryt7NixEICenlV8/esv\nZNasWZm0L+9qvX7tGCTxq/TfHcAhNfZxOEmhGQaelFSK8plOEsvytXT5ecBhTOwqqcel2XI7gOeM\nt3F+BXARMDtdHiBJE6p61XQzs66Q5Si+B4BXpkdMvSQXvIPkAmv3AWeUX9umDs+SNClGX377GpLs\nuyOAP29BmzOxbdt2Fi3axfr1WwE49dRdrFixHend1L4Shpl12syZMznppBvZuDFZPvnku5k587Rs\nG1VQnShQVbvPIuJRJZcU/gGwmaRbkIgISctIgjEDeDgi3lZx98sl/X5EfKRs3c3ATZIejIhz03W3\npT//TnIV10KSxMqVF+01is/M8kcSK1ZcWPb3epr/XicoN8PMs5T3c1BmNra8nYOyxnT6HNQokjZV\nWR0RcVK7921mZsXV9gIVESe2ex9mZtZ9Gp4HJemOOrZZIelBSa+vcluti8w10o5L0/lNbypbN0fS\nDyVdWbHtZWO1p8j6+voATwo0s+7UliSJDiU+BLAsIlaXrZsCXF6lPRfXaE9hVBai5cuXe1KgmbVM\n3r7s1tXFJ+kDwJtIJuCMSnwA/o6kWFwdEV8Z53FGEh+AjRHxf9P1i0kud/0i4F8ioq7LlpPMVB0R\nEeslHV/nfQulVIg2bnwFACefnMx9njatF7hqZLv+/oVMm5b87gm8Zlavap8xn/nMBzIdgThugUoT\nHxYAf0Sa+JCunwT0UZb4IOlrpcSHaiLirDTxofK81KqIWFFKfADqLVD7hD1JEqMLUfly7fu5WJlZ\nbYODg2zc+IqRBIwNG5J1WSZg1HME5cSHDho71mgs1VMmhoYea2m7zMw6rZ4CVbTEh0LPiBvrSGfv\nw++76e+HoaHHiDi/bFLgnyC5OJlZY/KYgDFugSpS4oOk9wPzgEMkPTcizim/eaKPmwfVZqdfccWk\nkdscRGlmzchjAkZhkySUBNG9FVheMZKv2raXAacB742IWypv90xws+JykkSx1Xr9ClugWslvMLPi\ncoEqtq65oq6Zme07ClugKpMkJF2bXlH325KOLNuusCkS5ZPmhoeHczWBzswmJm+TYfMsy+tBNauU\nJHETQES8A0DSScAFwLnp+oslPZNZKyeofNReRHDQQZfzxBNvBfIxgc7MGpfHybB5VuQCBdVH5j0J\n/KbTDWmlahNzn3rqL0Z+L6VFLF36fpYtW9bh1pnZROVxMmyeFb1AVbOEZF5V7jU+KXe05cuXs3z5\nFTW3cYKEmRVVVxUoSfOB/4iIn2TdlnrUKh57d/F9oayL724+85kPMG2aC5BZkeRxMmyedU2BkvQq\n4PiIOD/rtrRC5aS5l770s9x3331APibQmVnj8jgZNs+6pkABNwAPpUnp90bEeWW3FfIdUJkQ4X5q\ns+Jz8kv9ilygHgM+KGlKRKyOiCOrbVSWInFLJxtnZmbNKWyBiohrqGMwRHqxwovb36LOW7p0adZN\nMDNrG0cd4agSsyJz1FGxdWXUUZUkicskbZS0vmhJEpUzyz3T3My6RTOfZ4Xt4mPvJImLASS9Dng/\ncE5pfZ6TJKrNLB8eHmbTpleOLHumuZkVUbPJGUUuUFB9dN4xJBdSLITBwUHWrKm8lPs9wMtHlgcH\nv+tRP2ZWOM0mZxS9QI0i6Tskl6Kfk3VbmvPyUUtz5hwLeFKume1bCnsOqpqIOA5YBHw246bUbebM\nmSxYcD49Pavo6VnFGWdcwPz57xu1vG3bUNXi1NfXl0GLzczqkyRn3DXyeZYkZ8ys+/6FHcUn6RLg\nzohYV7H+CODTEfHH421bkvUonIgom1mevHjly2P11/b29jI0NNSZRprllEfx5Vvl51vl51mt169r\nuvgkfZGke+83wF9X3tz5FtWv2sxyn3Mys27QTHJGkQtUZZLEW6pt5CQJM7NiKmwXXysV9RDdXXxm\n7uIruq6cqGtmZt2tsAWqMkkiXTdF0s8lvatsXe6TJOpVPiN7eHiYxYsXO23CzHLNSRJ7vAP4UXpb\nslHOkyTqVTkj+8AD+3jiibdwww1bnTZhZrnkJInSL9JzgFNJrgvVk1mL2qRyRnbyb5I44bQJM8uj\nZpMkCtvFV8V5wMezboSZmbVG0Y+gAJB0EHBsRPRJWpR1e9ohmZF9Ixs3JssHHvgFnnjiLUj3pLOz\nP5BtA83MKlR+biWfVafVff+iF6jSuabXAc+W9HngRcB+km6JiIHsmtZaklix4sKRGdkvfel13Hdf\nkok7c+ZpPv9kZrlT+bnV6GdV0QuUANKBEjcBSDobOKCiOHXFp3fljOy1a9eybNmyDFtkZlZbM0kS\nRT4HVUqSeFP5yoi4LiI+WVpOkyQWAE91uH1tt3z58qybYGbWNk6SoLgzwZ0kYeYkiaLbJ8Jiu914\nicBmZt2msF18lUkSklZK+r6kTel5qNJ2hU+SGB4e5s1vXsbcuQ8wb95Wliy53OkRZjnWTHqC7VHk\nI6jKJIkA3hIRvxi1UcGTJCKC6dOnAZ8aWdffv5Bp07Jrk5mNrdn0BNujyAUK9h6d1xXvgN7eg+vc\n8pKRbX05eLN8aDY9wfYobBdfFU8C10taK+nFWTemtQK4ELg+/bkoXXfpyBb1FzUzs2Io+hHUiIg4\nD0DSK4ArgTdm26KJKz8aKnUXbNhwNMPDW/mjP/opX/7yR5g06bEMW2hmY2k2PcH26JoCVeZpoLDn\nnCqNnon9ImbOfIf7ss1yrNn0BNujawqUpC8Ah5J09b2r8ubOt6h1mpmJbWad57/Z1ihygSolSUyJ\niNUR8efVNkqTJE4Dbulk48zMrDmFHSQREddExDERsXqc7S6OiD+IiFs61LSO6evry7oJZmZtU9gC\nVWWi7gvSSbq3Srq6bLvCT9StVJoEuHz5ck8CNLOuVdgCxZ6JuqUjqKuAiyJiTkS8d2SjiIuBlRm0\nry1Ko/rmzn0AwKkSZgXgZImJKXKBgnTwg6TJwIyI+F7G7Wm7gYEB1qy5ip07zwKgv/9KBga65rJX\nZl2n9KVy3rytjiprUNELVMlvkVyw8EZJGyUVdg7UeLZs2bLXuuOOOyiDlphZPcqTJXbsWMiGDUeP\nDEG32rqlQG0DHgfeDPwJcKGkqdk2qT1mzJjBnlQJSFIlXuYkCTPrOl1RoCLiGeAh4JCI+A2wK+Mm\ntc2sWbOYP38XU6duZb/9TmLBgl1s27bdWXxmOZUkS9xFT88qenpWpckSM7NuViEUeR5UpfcDn5Z0\nEPCliNhZdluhJ+qWk8TKlRelXQTH+dpQZjnnZImJK3KBqpyo+wuSCbmjdONEXc9SNysW/81OjC/5\nji/ZbFZkvuR7sdV6/briHJSZmXWfwhao8iQJSc9NUyRKP4+Xbdd1SRKQXAa+v7+f/v5+hoeHs26O\nmY3Bk3QnrsjnoCov+X4igKSXA38zslHBL/lezfDwMC972Vk88sj+wBs47LCzueee65g0qbDfN8y6\nki//3pyif6JVe5XPA/6x0w3ppLVr1/LII+uAfmARDz/8NdauXZt1s8ysgifpNqfIR1B7kTQNODwi\n7sm6Le20ePGimus8J8rMukHRj6AqvR34VNaNaLdHH93GoYeWTqmt5LDDTufRR7cxNOQJu2Z54km6\nzemaIyhJ+wGnA3Oybku7TZo0iXvv/SzTp09nxQqYP9/nn8zyyJN0m9M1BQp4A7A2IqoNaeu6d0Sp\nIJ1xxhkZt8TMavEk3YkrcoGqTJL4crWNujFJomTp0qVZN8HMrG2cJIFngpsVmZMkis1JEmZmVjiF\nLVDlSRLp8lmSfiDpNkknlm3XdUkSTpEwyzenR7RGkc9BVSZJnA+8EjgA+CbwGui+JIk9KRJnAjhF\nwixnnB7ROkUuUDB6dN4gcDxwCHB7Ns1pv+nTfwmsG1l++OFFrF270qP5zHKiPD0CYMOGZJ1H8jWu\n6AWq3LeAvwWeBXwi47a0RXJZ9/JLu18KXOoUCTPrSl1RoCQdCZweEQvS5e9IWl9xVd3CGxraXtHF\n9yEOO+xOd/GZ5UiSHnEjGzcmy0l6xF7XUrU6dEWBAiaTPhclHb1TSc5RdZ1SisTatWtZvBgXJ7Oc\ncXpE63RFgYqI+yXdLukmkpGJn4iIp8s26ap3x6RJk0bOObk4meWP0yNao8gFqjJJ4sPVNurmJAkz\ns27mJAmKOxO8t7eXoaGhrJthliknSRSbkyS6lLP4zKybFfYIStKlwJuASyNitaRzgEXADuCdEXF/\nut1lwNuAv46IddUeqwjfgCKCwcHBkVnpkpg5c6ZPvto+z0dQ+Vb67AKqfmbVev2KfA5qJElC0nOA\nxRFxjKTpwD8BZ0J3JEmUZqZv2HA0Tz/9DeBkpk59lmeom1muNZuqUeQCBXtG5wnYX9IUksETh0ja\nPyIKXZhKBgcHWbPmqnTpLwDYsQP6+xcyOPhdjxYys1xqNlWj6AUKgIh4StKHga8DT5LELTwP+HWm\nDWtSkhwBcOyY28yZc6zTI8ysK3VFgQKIiK8AXwGQ9OOIKHRxgj2xRWN38d3NZz7zAbpsmpeZdYlm\nUzW6pkCVSDoNuCvrdrRS+cz0iHeMrFu79l6ffzKz3Go2VaNrCpSkfwVeQjKKb2HlzZ1vUWtVm5k+\nZ84cli1bllGLzMzG10yqRpELVGWSxF9W28hJEmZmxVTYeVCtVNR5DE6SMPM8qKLr1nlQ+4TxJrmZ\nmXWr3EcdSbpU0j2S3pQuz5H0Q0lXVmx3iqRb05+TytZfJulBSa/vdNubVRq9N3fuA8yd+wBLllyO\nj3jN8i8iGBgYYGBgwH+zTSjCEdRIYkS6PAW4HHhtaQNJk4APAaekq74paVMkCpskMTAwUDZBF/r7\nYWDgVmbPng04i88sj5pNT7A9cn8ElRp5ZSNiPVB54uUo4KcRsTO9iu4W4MUdbF9bbNmyZa91xx03\nZ2QCr0fwmeVPeXrCjh0L2bDh6JFuemtMEY6g6tELPCbpo+ny48A04P7smtS8xYvPBi4EZqdrBoDL\n6IJR82Zm4yrKEdR4tpFEG10IXJT+/mimLWqBbdu2M3/+LqZO3crUqVtZsGAX27Ztd7SRWY4l6Ql3\n0dOzip6eVWl6wsysm1VIRT2CqjyE2AL8XtnyURHxsw62py0ksXLlRR7FZ1YgzaYn2B6FK1CS3g/M\nI0ksf25EnBMRuyV9CLg53ezSyrt1so2t1MwsbDPLhv9uW6MIXXylxIg3AUTEFRFxQkT8fkScU9oo\nIr4VEcemP6VCVUqSWAA81fGWt1lfX1/WTTAzaxsnSVDcmeBOkjBzkkTR1Xr9cn8E1cBE3bHWF3ai\nbokn/ZnZvij3BYo9E3VXp8ulibqVqq6PiIuBlW1rXZs5TcKsuPzlsjlFGSQxaqKupOMrNxhrfdFN\nm9YLjJ0mYWb55ESJ5hWlQO2T9lzyfbTjjpvT4ZaYWaPKEyUANmxI1nl0X/1coHJsaGg7vb3PI5l7\nXC1N4pKsmmZm1nZFOAdVzVjHyF137FwrTWJo6N1ZN8/MxuBEieYV7giq2kTdWutLd8ugqS3hNAmz\nYnKiRPNyPw9K0ruBtwLLy0byNXL/0iXf3xsRt1TbxvMYzIrL86CKrdbrl/sC1Ql+g5kVlwtUsRV6\noq6Zme2bcl+gGkiSuFbSJknflnRk2frCJ0mMpTyLr3JCoCcImlkeNPNZlPsuPkmXAHeULvku6RTg\nQOC1EXFBle1PAs6MiHMrHuPOiFhXbR9FPUQvZfFVTgg86aS7AMomCN7tCYLWtdzFl197T1be+7Oo\n1utXlFF84yZJlHkS+E37m5Qfg4ODrFmzJ21izZqFo27v71/ItGn4Qodm1lHNTlYuSoFqxBLgmqwb\n0W6llInk32Mbuo8LlZkVQVcVKEnzgf+IiJ9k3ZZ2S1Imkn/dxWdmeZRMVr6RjRuT5WSy8ml137+o\nBWqvT1pJrwKOj4jzM2hPpqpNCAQ8QdDMMtXsZOXCFagaiRE3AA9J2gTcGxHnld+t0+3shKVLl478\nXu0S0w6lNLOsVftsqvu+BRjF5yQJMxuTR/EVm5MkxuE3mFlxuUAVm5MkzMyscHJfoBpIkrhM0kZJ\n64uUJBERbN68mf7+fjZv3uzUB7Mu4kSX5hRhkEQAy0pJEsAU4HLgtaM2irgYQNLrgPcD55TWS3qm\nc82tX2l4+Lp1M9m9exKTJ/8zr3/9IaxYcaFH3ZkVnC/53rwiFChoLEniGOC+9jepedOm9QJ7EiB2\n717EmjUwbVqyPN6E2r6+PpYtW9bGFprZRPmS780rSoGqi6TvANOBOVm3pVIpxaGZ+1QWrOXLl7tA\nmVnX6qoCFRHHSfpD4LNArs45VTsa2ruLb727+My6RLMpClbcAlXr0/tXFOR5lWZZDwwMsGXLFmbM\nOIdZs2a5OJl1AV/yvXmF+CAvN1aShKQvknTv/Qb468q7dbaV9ZPE7NmzmT17dtZNMbMWayZFwYpR\noB4DPihpSkSsjogrgCsqN4qIt1S7c1mSxC1tbaWZmbVU7gtURFxDE5fPSIefX9y6FuVHeRafmVm3\ncdQR+Y4qiYiyPuyZ7sM2q+Coo2IrdNRRvUkS6W1TJP1c0rvK1uU6SaKW0ii/uXMfYO7cB1iy5HLP\nRjcrOKdL1C/3XXzUmSSRegfwo/Q+yZ1znCQxnsqJvP39+NLtZgXmdInGFKFAQR1JEpKeA5xKcl2o\nng62reUmMqnXzPLP6RKNKUqBqsd5wMeB52fdkGaVjpA2b97Mccd9HSgNQR8ALhspYD6SMrNulvtz\nUPWQdBBwbER8gxzPeWrUrFmzmD9/F1OnbmXq1K0sWLCLbdu2MzSU/PT19WXdRDNrQJIucRc9Pavo\n6VmVpkvMzLpZuZX7UXySLgHujIh1ZetOAF4fEReky6cB7wV+DbyI5Mjw7IgYGOsxyuV5FE6tUXy9\nvb0MDQ1l1TSzXCjaKD6PzB2t1utXuC6+akkS6QCKm9LbzwYOKBWn0t0yaGpLeCa6WXfx33T9ilCg\n6kqSKImI68qXnSRhZlZMue/i64Q8d/HV4i4+s+J18dlo+8REXUkrJX1f0qa0m6+0vrATdRvlCYBm\n1k2K0MVX70TdAN4SEb8YtbLAE3XHU57FVzkB8NWv/jc+9KGFvnyHmY1SpEEaRShQUP8l3/P7P90G\npavp7pnYuyd1YtOmhdx++4eZO3eNZ6qbGVC8JIuiFKh6PAlcL2kIeE9E/CzrBrVTPWkTO3deOBKP\nBJ7Ya7avK1qSRdcUqIg4D0DSK4ArgTdm26L2Ki82SbEK4CIqUyfKDyprFTUXLzPLm6IWqFrHo08D\nXXnOaSyl4hJxPgMDA3zwg//CnXeegvQ5Tj757lwfwptZ5yRJFjeycWOynCRZnJZto2ooXIGqccn3\nLwCHknT1vavybp1tZTZKl49fvfqjZSdBT3NxMjMg+YxYseLCwnw+5H4elKR3A28FlkfE6gncvzRR\n970RcUu1bYo6j6Gvr29koITZvsrzoIqt1uuX+wLVCUV9g3mirpkLVNEVeqKumZntm3JfoBpIknhB\nmiJxq6Sry9Y7ScLMrICKMEii3iSJq4CLIuJ7o+7cpUkSpQI0MDAwcj2ZIk3AM7N9QzPJFbk/gkqN\nSpIARp0PpSU2AAAQDklEQVR4kTQZmFFZnLpVaTY4wLx5W1my5HIGBgZYs+YqduxYyI4dC+nvv3Lk\nTWFmloXSZ9W8eVtHPqsa6d0pSoEaz28Bz5Z0o6SNkrp6ku7g4CBr1lwFXDJSjLZs2bLXdnPmHFtX\n4oSZWTuUJ1fs2LGQDRuObuiLcxG6+OqxDXgceDMwGbhN0jciYme2zWq3S0d+W7z4bOBCenqS7j5P\n0DWzoitqgRr1qRsRz0h6CDgkIn4paVdG7eqImTNnsmDB+WXnm5JiBOcXZgKemXW/ZpMrcj8PStIl\nwJ0RsS5dHkmSAL5dliRxBHAtcBDwpYi4puwxLgXuKD1GpSLOYyhSZL5ZO3keVL6N91lV6Im6TpIw\ns1pcoIqt0AWqE/wGMysuF6hic5JEl+rr68u6CWZmbZP7AlVPkoSk56YpEqWfx8tu68okiYhg+fLl\n9Pf3s3nzZidHmBWA014aU4RRfOMmSUTEE8CJAJJeDvxN2W1dlyQRESxa9PcALF4Mkyffz+mn+9Lu\nZnlWtMut50ERChRUJElIOr7GtucB/9D+JmVn2rRe4CPA1cAidu9m5NLuvjKuWT4V7XLreVCUAlUX\nSdOAwyPinqzb0g71pEKUb+NiZWZF1lUFCng78KmsG9Euey7tnnTxrV0LsJLJk4c5/fT73F1glmNF\nu9x6HhS1QO31KSxpP+B0YE7nm9NZkli58iLe974hjj8eZsw4ilmz3uDiZJZjRbvceh4UrkCVJ0lI\nem4pSQJ4A7A2Ioar3a1jDewQSVx99dXjb2hmuSHJ55waUIQC9RjwQUlTImJ1RFwBXFG5UUR8udqd\ny5IkbmlrK83MrKWcJIFngpsVmZMkis1JEmZmVji5L1D1JEmk68+S9ANJt0k6sWx9IZMkIoLNmzc7\nKcKsgJwY0Rq5L1DsSZIoJZmXkiQqnU+SLjEP+PDInSMuBla2uY0tVRpGfuKJ17J4MZx44v1VL5Xs\nLD6z/Gn2Mue2RxEGSUB9SRKDwPEk14m6vVMNa4c9SRGJsZIili9fzrJlyzJooZmNxYkRrVOUAlWP\nbwF/CzwL+ETGbWlYPSkR1bYrLTs1wsy6TVcUKElHAqdHxIJ0+TuS1kfEzoybVrfyAlPq4rvppl+x\ne/epYyZF9Pa6MJnljRMjWqeoBapyWOJk0uei5BN8Ksm5q0IqJUUMDAywZcsWJ0WYFYgTI1qncAWq\nWpJERNwv6XZJN5EM/PhERDxdfrdMGtsEScyePZvZs2dn3RQza5ATI1qjCAWq3iSJD+991+5Okli6\ndGnWTTAzaxsnSeCZ4GZF5iSJYqv1+hXhCGqfEBFlfdYz3WdtZvu83E/UbSBJ4hxJ35d0s6Sjytbn\nPkmi2sS+zZs3exa6mdWtG9Mrct/FJ+kS4I6IuCldPgU4EHhtRFyQrnsOsDEijpE0HfiniDiz4jHu\njIh11faR5SH6ePOfzjjjAl+I0KwGd/Ht+ZK7ceMrgGRoe1E+N7qhi2+8JAkB+0uaQjKo4hBJ+0fE\nM51sZD3qnZBb0t9/5V4JEmZm5bo1vSL3XXz1iIinSPL3vg58BTgYeF6mjRpD9UITwIXA9enPRVRO\n4+rtPXjUDziLz8y6W1GOoMYVEV8hKU5I+nFE/DrjJo2pWpGKOJ/BwUEigiuv/BGbNn0OqH2o7iw+\nM4PuTa8oaoEas89S0mnAXR1sS0uUT+xbuXKWZ6GbWd26Nb2icAWqWpJEuv5fgZcAO4CFlXfrbCub\n41noZtaobvzcKEKBqjdJ4i+r3bmbkyTMzLpZ7oeZd0IeholORG9vL0NDQ1k3wyxTHmZebLVev9yP\n4qsyUfdaSZskfTu9zEZpu1Mk3Zr+nFS2PvcTdSeqMouvGyfqmdm+K/dHUJUTdcvWnwScGRHnSpoE\n3Aqckt78TeD4SJ9cnifqtkppot7NN78cgLlz7y3MRD2zZuTxCMrRZfXrqom6ZZ4EfpP+fhTw09IF\nCiVtAV4M3N+Z5mUvuUz8VSPL1S4Rb2btt3eqQ7+/LE5QUQpUNUuAa9Lfe4HHJH00XX4cmEaXF6h6\nUikqt3HBMmuvbk11yEIhC5Sk+cB/RMRP0lXbSJIj3klytPVJ4NGMmpeRIEmgKF3gcAC4jIKNsDcz\nG1G4AiXpVSTnl84vW70F+L2y5aMi4medbVnnlR8NRQSLFu1i/fqtAJx66i5WrNjubgWzDuvWVIcs\nFK5AATcAD0naBNwbEedFxG5JHwJuTre5tOI+Xfkp3dfXNxJ1JImVKy/yiVmzjHVrqkMWijCK793A\nW4HlEbF6AvcvTdR9b0TcUm2boo7i8zwos3yO4rP61Xr9cl+gOqGobzAXKDMXqKIr9ERdMzPbN+W+\nQDWQJDHWpeALnSThdAgz21cVYZBEAMtKSRIR8Q4YSZK4ADg33W4KcDnw2lF3jrhYUu6urFuP8gl/\nu3cPc8wx/8YNN3yYSZNy/73CzMo4WWJiilCgYPwkibEuBV9olekQmzadxfTpsG3bEJL2yuIzs/xx\nssTEFaVAVVOeJNF1aqVEJIULhoZ8NV2zvHOyxMQVskBVSZLoKnuKU+10iN7egx1dZGZdq3AnM8qS\nJD5W7eZOt6cdhoa2pz+P8eij7+P447/J1KlbOeCA3+eMM5aybdvQyDZmlm9JssRd9PSsoqdnVZos\nMTPrZhVCEY+g9kqSgLEvBZ8qbOGaNGkSq1d/tOwE6+vdd21WIE6WmLjcT9R1koSZ1eKJusVW6Im6\nEXFNRBwzkeKU3v/iiPiDsYpTkfX19WXdBDOztsn9EVQnFPUbkKOOzHwEVXSFPoJqIElirPW5TJIY\nHh6mv7+f/v5+hoeHs26OmbWYU2CaV4RBEnUlSdRYn7skieHhYV72srN45JEzATjssLO5557rnBBh\n1iU8Obc1ilCgoI4kiTrW58b06dOAdSPLDz+8iOnTfTl2s27hybmtUZQCVc1YSRK5S5iolQox1nYu\nVma2rytkgRorSSKvCROVxWbvLr4vT6iLz1l8Zvnky763Ru5H8Um6BLgzItaly68C3hoR51dsV3V9\ntceolMUonOHhYdauXQvA/Pnzff7JbILyOorPCeb1qfX6FfEIqmqSRI31kMMkiUmTJnHGGWdk3Qwz\naxNJPufUpCIUqMeAD0qaEhGrI+LIahuNtb4sSeKW9jXRzMxaLfddfJ3giXZmxZXXLj6rT1MTdSX9\nVWubU3Ufd9SxzYqxJtx2oo1mZtZZ9ZyZf3vbW1GHiFgMrBzj5ly0sdOcxTdxzc7yd0qAtZrfU3ur\nWaAkfRZ4SRohdHG6brGk/jR+6Lyybe+o9nuNx/6ApDskrQR6ytb/iaTvSbpN0pvreJzyNn6wbH3V\ndnaT5cuXZ92EQirN8p83byvz5m1lyZLLG/pAaPb+ZpX8nqpu3HNQku6IiFeXLe8fEc9ImgL8MCKO\nrtyu8j5VHvNQYDXwOuAAklF3L5Q0Cfhxun4XsBE4NSJ2pferOly82v7Gamc1Re1D3tfDYuudAN1t\nPIl7tG44BzUwMMC8eVtHkid6elbx9a+/cJ8YBdjqYebHSTod2AE8Z4JtOpyk0AwDT0r6dbp+OvAC\n4Gvp8vOAw4AHM2qn1WFfLRRZ2Vf+v12IrZ4C9SxJk9JiAkmM0MuAI4A/L9tuEoCk5zB+QXgAeGV6\nxNQLHJqu/zVwH3BGRDxR31Oo2sZa7bQWK9oHyd5Bnnc3FOTZ7P3NKjl5orp6uviuAmYDD0TEOyX9\nM8kH/78Dr4qIY9Lt/gl4iuSI5U8jYvY4j7sUOBPYDLy6tL2k1wF/T5Ji/nBEvK3sPpcAbwaui4iP\nVGnjgxFxbrquajurcRffvqfZWf5OCciPbujig333PVXr9fM8KIpboPr6+li2bFnWzTDLVLcUqH1V\nZgUqjR2qFBFxUtt2OgF+g5kVlwtUsWWWxRcRJ7bz8c3MrHsVIYsvN/bVPmIzsyz4Gg91Ko3cmjv3\nAebOfcAT6cxsTE6FaI2uKFCSFklaJWlA0jsl3SfpdxtNt6hlYGCAdetmsnPnWezceRZf+9pLGRgY\naL7xZtZVnArROl1RoEiGpD8IfAY4EPg88IpW7mDLli3s3r3nv2v37iVs2bKllbtomLP4zPJncHCQ\njRtfwY4dC9mxYyEbNhw9cmrAGtMtBQrgV8DT6b+7aPH5tRkzZjB58s3A9enPRcyYMaOVu2iYs/jM\nrJt1U4GqppF0i5pmzZrFaacdwtSpW5k6dSsLFuzaJ3KyzKwxSSrEXfT0rKKnZ1WaCjEz62YVUldM\n1JV0NnsS0XcAhwA/A06hjnSLeucx5G0Un5MkzPI5DypvnxV55iSJcRR1op0LlFk+C5TVr6kr6pqZ\nmWXBR1DAhg0b/J9gVmAnn3xyZkdR/vxo3livnwuUmZnlkrv4zMwsl1ygzMwsl1ygzMwsl1ygzMws\nl1ygWkDStZI2Sfq2pCM7tM9TJN2a/nT0ApBZPN+K/U+R9HNJ78pg3y9In/utkq7u8L7PkvQDSbdJ\nauu11iTNkfRDSVeWrevIe26MfWf6nmtUteeQh31JWinp++n/5dltak/LXiuP4muh9I/2zIg4t837\nmQTcSpKUAfBN4Pjo8IvZqedbZb/vBo4H1kfEJzu87y8A/xAR3+vkftN93wO8EjgA+GZEvKaN+zqF\nJHj5tRFxQSffc5X7rrgtk/dco2o9hyz3JWkFcElE/KKdbUr31fRr5SOo1noS+E0H9nMU8NOI2BkR\nO4EtwIs7sN9KnXq+I9JcxVOBfqCjc18kTQZmZFGcUoMkhfl04PZ27igi1gPlMSUde89V2Xe5jr/n\nJmKc55D1vjr1dzPqtZJ0QqNHbb6ibgMknQosrVj9voi4J/19CXBNB5rSCzwm6aPp8uPANOD+Duy7\nXKeeb7nzgI8Dz+/wfgF+C3i2pBuB5wL/GBFf7eD+vwX8LfAs4BMd3C/s2++5bvIkcL2kIeA9EfGz\nNu5rCXCNpKOBjwHPI/n7WQT8fVpYa3KBakBE3AzcXO02SfOB/4iIn3SgKdtIXux3knwb+iTwaAf2\nO6LDz7e0z4OAYyOiL32Td9o2kg/mNwOTgdskfSM9omirtC//9IhYkC5/R9L6Tuw7tU++57pNRJwH\nIOkVwJXAG9uxnyqv1YmSjgdeGBHX1fs4LlAtIOlVJP3x53dol1uA3ytbPqrN34RGyeD5lryO5BvY\n54EXAftJ2hQRHbkaXEQ8I+kh4JCI+KWkXZ3Yb2oy6d+rkmjsqSQX6myn8q6gTr/nRnVDZfiea0Yn\nu6Ab3dfTwDNtaUgLXysPkmgBSQ8ADwHDwL2lbylt3udc4O/SxQ+lR3cdkcXzrdKGs4EDMhgkcQRw\nLXAQ8KWI6Fh3k6QLgWNJzh1/ISJWtnFf7wfmkVy65tsRcU6n3nNj7Dvz91wjqj2HTu9L0pnA/4uI\ndWXbfgE4lKSr710R8fM2tKdlr5ULlJmZ5ZJH8ZmZWS65QJmZWS65QJmZWS65QJmZWS65QJmZWS65\nQJmZWS65QJmZWS65QJmZWS79fyvuf+HB62KoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e7bcf98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.Matplot.summary_plot([M_random_treat.m_delta, M_random_treat.delta, M_random_treat.tau_delta, M_random_treat.mu])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def random():\n",
    "    \n",
    "    m_mu = pm.Normal('m_mu', 0, 1e-5, value=0)\n",
    "    tau_mu = pm.Gamma('tau_m', 1, 0.01, value=1)\n",
    "    mu = pm.Normal('mu', m_mu, tau_mu, value=[0]*N)\n",
    "    \n",
    "    m_delta = pm.Normal('m_delta', 0, 1e-5, value=0)\n",
    "    tau_delta = pm.Gamma('tau_delta', 1, 0.01, value=1)\n",
    "    delta = pm.Normal('delta', m_delta, tau_delta, value=[0]*N)\n",
    "    \n",
    "    theta_t = pm.Lambda('theta_t', lambda mu=mu, d=delta: pm.invlogit(mu + d))\n",
    "    theta_c = pm.Lambda('theta_c', lambda mu=mu, d=delta: pm.invlogit(mu))\n",
    "    \n",
    "    treat = pm.Binomial('treat', n_t_obs, theta_t, value=r_t_obs, observed=True)\n",
    "    control = pm.Binomial('control', n_c_obs, theta_c, value=r_c_obs, observed=True)\n",
    "    \n",
    "    treat_pred = pm.Binomial('treat_pred', n_t_obs, theta_t, value=r_t_obs)\n",
    "    control_pred = pm.Binomial('control_pred', n_c_obs, theta_c, value=r_c_obs)\n",
    "    \n",
    "    return locals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 10000 of 10000 complete in 8.7 sec"
     ]
    }
   ],
   "source": [
    "M_random = pm.MCMC(random())\n",
    "M_random.sample(10000, 5000)\n",
    "M_random.sample(10000, 5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEgCAYAAADhUed1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2YXGV9//H3ZxfEQHxgoS1EoUiMtokglp9WxQghkBpg\nDT5Q259cgY2XolixPgARsYT+KEQIVeoT1UoCTRGfkG0AFUJAEVGiFIENKgYiKNBqFhBoQMx+f3+c\nM8uZ2ZndmZ1zdufsfl7XtVd2zpw59z05s/Odc+b+nFsRgZmZWafpmuwOmJmZ1eMCZWZmHckFyszM\nOpILlJmZdSQXKDMz60guUGZm1pFcoCaYpG5J50i6V9JDktZI2rlmneMlPSHp/vTnPkm7pfd1SVqb\nLv++pD3T5Z+Q9KEW+/LHki5O+3G/pE2SPprfs23Y7lpJZzS470JJb6yzfIukQ0fZ5k7pc3hA0pCk\nwl7bktZn2rlf0qfGuZ0Vkv497/7lTVJ/+jyflrRwsvsz1aWvi8rf/68l/UTSW1t4/CGS7s+pL8dL\nujGPbY3HDpPV8DR2EvA6YH9gG/AF4BPACZl1AvhyRCyr8/i/Al4QEXtJWgG8X9LFwELgFS325Urg\nJ8DsiHhC0nOBvVvcxnhE+jPyjoh3j/KYxhuMeArYS9KfAve2173RRcRhmXb+NCKGimxvskXEEgBJ\n9zLGfrBcBPD1iFgKIOlVwHpJd0bETye3axPLR1AT7xDgKxHxWET8AfgA8HZJu2fWUfpTz9NAt6Ru\nYCfgD8AngQ9HxPZmOyHpCOCFwAkR8QRARPwuIu6sWW9NesT36fRI7teS9s/cL0nLJd2dHuVcKOnZ\nNdt4m6RfpJ8ILwFm1j4/ST9O739S0jsadHs/STdI+m9J35M0p95TG+N5L5N0V/pcviKpZ7T1R9vU\nKG3cIOlkSavT/5N7JB1Ss86tJPv+TZmj5Fk16xyVfnq+X9K3Je1dc/8+6VHcfEm3pftmTeb+z0v6\nl5rHfEnSWZnbu0r6pKQBSQ9Kul3S61r+z0j+X+9In8d9rR7NW5Wqv/+IuIXkw9CftbQR6U2SNkr6\njaRLJSlz35j7XdJ/AucDr9QzZ3P+sq1n1qqI8M8E/gCnA98DdiN5Eb4KeARYkFnnb9NlW4DvAgfX\nbOMzwGbgm8C7gCvG0Y8zgG80sd4a4H+Ad6W3n1Nz/weAO4A9ST7wrAY+nbn/z0mOFA9Jb78G2Ar8\nQ4P2rgeW1Vm+JX2+u6b/b+cDP66z3j7AENBV5743AQ8AL0lv/yNw5Tj342jt3ABsAv4i086NDfbB\nJQ22fyDwO+DV6e13kBztqk4fvgXskS57bub+/wP8Btgxvf084HGSo77KOjsBiyrPAzgTGGjQp3uB\nQ+ssn5Pu47mZZbtM9N/WVPkBVgD/nv7eBRwL/ArYtcnHH0LyQfbsdP/uBgwCC1vd78Bx9V67E/Xj\nI6iJdw5wHUnh+TnwHuAp4I8rK0TEl4A5EbEPyYv1Ckn7Zu5/b0TMBo4G3k9SJFr1ApI3q7EEyRHf\n59O2H6u5/wTgrIh4MJJTXaeRvJlW/DWwLiJuSB9/M3A1YxzpNOjHeRHxcCR/Of8AvDw91dasE4AL\nIuLn6e1/BA5R+j1ejgJYFRG3prdvpP6p09GOlN8J/EdE/AAgIr5I8qby6jrrHhsRD6Xr/W64ExE/\nAn4NHJUu+mvgexHxy8w6T0XENfHMacqvAS8d+ylWeZTkNbxI0gvS7T7R4jas2tHpKdVtJPvvtRHx\ncAuPfygiTkv371aSD0x7Ve5sYb+3+neaKxeoCRYR2yPijIiYFxFzSN7MnwP8d816v0n/3UBSzI4a\nsTE4GfgG8B5JN0v6qqRdm+zKo2SK4hhG+8PYCzhfyaCPe4HvA9syb/p/QnL0k6v0DfBh4I9aeNhe\nwAcyfb0beJJivnfL/mE/Tet/a3sBx1T6mva3h+S0bK3BUbbzBeD49Pfj0tvPdDLxAUnflfQ9kqPz\nLrUwyCQi/ofkTMDewAZJt8qDKdr1jYh4EfAV4HkRcV/lDkmvyZxyu1/STU1sr+o1mMd+nwgd1Zlp\n6o0knz5/MMo6zyL5JDVM0guBtwM/BF4ZEa8BBoC/b7LdH5KcW96l5R5X2wwsjYgXZX56IuLB9P5f\nAS+qecwOjO/L9h0rvygZ1dgD3Nd49bp9/WhNX3ePiB+Ooy9jaXcwwWbgczV9/eOI+GqL2/kPkqPE\n1wIvBvpr7n8f8FbgmIh4HUkRa1lE/DwiPhgRLyU56l8naafxbMsInvmAcyIwR9LJw3dG3BwRe2V+\nDhpHG7ns96K5QE0wSbtIen76+77AuSTfxzyZWWdvJYMgkHQ4yYi/K2o2tYrkjeAP6XqV00XPpjn9\nJN/H/KuS0XtI2l1S7ZHaWIf4/wx8UtLw6YHK80t9GThC0qvTT21/DfSOsc16bQr4uKQ/kbQDyf/b\nt9JP781u4xPAGUpGRVX62uwRZyt9HW151iDJwI8dJe1Qc6rxs8AJkv5qeIPj6GtEPEKyr78EXBzJ\nwJysFwIPAf+jZKDOeenyHamv7vPKnoIm+UD1FND0oB2rkh0g8Rjwf4EV2ddtu9ul+f3+MEmBfE76\n9/uCNvvQEheoibcfcJukXwPXAJ+KiNoczYeA+yX9EjgFWFQ55QcgaT7wJxHxZeDbwIMkRxJ/BfwL\nTYiIp0mGpm8H7pb0K5LvSmqHqjccEp5uZw3Ji/uydPTWvcBZmfs3k3xa+0bax0NIvoMatXsNlv0H\nSaH+b2Bfku9pGtki6baavl5PMqjkgvTUyL0kgzpaImk9yanMSNupl4OqfQ71ntOlwP+SHGXeQfKd\nYqWvPyUp5Kdm/l//s85RSTNHal8geUP6tzr3nQ/sTvJh5Srg8yQfehq9EX0p/b8bHkGZfsD5iqRf\nKcnfvBdYXKcYWnOq/uYiGcX3/4BLJT2nhW2MtqzZ/f4t4FbgHuCnJH8/E0bpSA0zM7OO4iMoMzPr\nSC5QZmbWkVygzMysI7lAmZlZR+qYi8VK2hgRrxxjndUko8D+LiKuqrnvnRHxhboPHMN1113nkSJm\nJbZw4cJJu+KB3z/a12j/dcwovmYKVLreGcCP6hSoph5fz8MPP9z0f0JEsGnTJubPn8/WrVvJXH/R\nzCbBrrvuOql/hK28f9hIo+2/XE/xKZk7ZG16hdwTlVw1uuG10iR9RMnVdteQXOG6svwNSuY6uknS\nW5po9xLgpZKul/SxzPI+JXPZ3C7ppDafHhFBX9/ZLFp0DwDLlp1DpxR4M+scEcHAwAADAwN+j2hD\nrkdQko4juZzKIySnD58N3BYRtZdXIU3NXw4cBOwC3BER+6TXgro1Xf4UsAE4PJL5flo6gpK0Y0Q8\nnYYbb4mIl9frd7OfgO68804WLLib7duXAaK7+4tcf/0cXvaylzXzcDMrQKcdQVU+yG7YcAAACxf+\nhIsu+ojPtjQw2v4r4juoh9J/Hwf2GKWNvUgKzRDwmKTKlRJ2J0m9X5nefj4wi/FNQvf69NI9jwM7\nj7XyWDZv3sz27c8cdG7f3sXmzZtdoMxs2KZNm9iw4QAef/xYAK67Llk2b968Se5Z+UzmIIl7gFek\nR0w9JPMJQTJ/zV3AkuzUAU14lqSuqJ7d9AKSSwvtDfxNux2ePXs23d0Xsn37s4A30929ntmzTxjz\ncWZm1rqJKFCNpvb+rZIZG38I3ElyWpCICEnLSa47FsADEfH2moefI+nPIuL8zLJrgasl3RsR70mX\n3ZT+/BfJJHltmTdvHkccsQfr128BDuTwwwf9qcjMqsydO5dDD72CDRuS2wsX/oS5c4+Y3E6VVMeM\n4ptM4xnFB8kL0eeVzSZXp30HBX6faMVo+6/wAiXp+jqLIyIOLbThFniYqFl5dWKBsuZN9CCJKhGx\noOg2zMxs6mk5ByVpYxPrrFYyTfWRde4bbQ6fVvqxIs03vTmzbL6kWySdV7PuWY360yrnG8zMJkYh\n1+KLiD5gTYO785rwKoDlEXF5ZtlOwDl1+nP6KP1pvsFMUHfBgvMd1DWzuvxBNh9NFagOvuJD1bnL\niFhPMo12IQYGBrjqqrls27aUP/zhCq688s8ZGBgoqjkzK6HKB9nFi7ewePEWf5Btw5jfQaVXfHgj\n8JekV3xIl3cBK8lc8UHSlZUrPtQTEUvTKz7Ufi+1NiJWV674QJPTlk+0kUHdZWzevMZBXTMb5qBu\nfpoZJFG6Kz4UpTqoC93dxzqoa2ZWkGYKVNmu+FDYkNNsUHfbNjjyyD38qcjMqjiom58xC1SZrvgg\n6VRgMbCHpOdGRPbwpu3CJYk1az6aTrcBq1ef5gCemVWRxOrVp2WCukf4fWKcSnslCUnvB/4WOLdm\nJF+9dc8CjgA+GBE31N4/nqDdypUrWb58easPM7OcOahbbpN6JYky8AvMrLxcoMptwiYsnEi1QV1J\nF6bD178jad/MermFdMH5BjOziTKZ0220qxLUvRogIt4NIOlQ4GTgPeny0yU9nUuDab7h2mv3B2DR\non5PRGZmI/hisfko7RFUqt5efwz4fRGNZYO627YtdVDXzEZwUDc/ZS9Q9SwDPlfEhhvNqGtmVpEN\n6j7++LFcd93Lh4+mrDVlPsU3gqRe4GcR8dMitl8d1P063d0zHNQ1MyvIlDmCknQgcHBEfLKoNipB\n3RkztgCXO6hrZiMkQd3bmDlzLTNnrk2DunMnu1ulVNph5pLOILkE01Xp7XuA+4Eh4I6IOCmz7gpg\nY2XdWuOZUXf+/Pls3brVX36aTbJOHGbuQRLNm9QJCwv0CPAxSTtFxOURsW+9lTIh3RvyaFTS8FGT\nX3RmVk/2fcLGr7RHUHkaT9Cup6eHwcHCZvYwsyZ14hGUNW+6BHXPkrRB0vqJCOpWfjczs2KUtkBR\nM6NuRJweEYcCZwCnDq+U02y66baGZ9TdYYejnW8ws3HzVWnGVuYCBfWDuq8mmQYkd9Uz6n7DQV0z\nGxeHeZtT5kESI0j6LslEivOL2H6joK5n1DWzVnjW3eZMqQIVEa+X9CrgEiCX75yyRs6ou95BXTOz\ngpT9FF89D1FQ4c0GdWfM2OKgrpmNi8O8zZkyR1CSvkxyeu/3wN/V3p1TG8Mz6oIDeGY2Pp51tzll\nLlC1Qd231VupqKDuypUrffRkZuPmMO/YHNTFQV2zMnNQt9ymZFB3sjioa2Y2MUpboGqvJJEu20nS\nLyW9N7MstytJZIO6gLMLZlaXQ7j5KG2BouZKEql3Az9O70tWyvFKEtmgLuCgrpmN4BBufspcoCAz\nOk/SzsDhQD85jdqr5Rl1zWwsnlE3P2UvUFknAZ8usoEkqHstcCnw5jSoO7vIJs3Mpq0pUaAkPQ94\nXUR8i4KOnqA2qHugg7pmNoJDuPkpcw4Knvmu6SDg2ZK+BLwI2EHSDRGR6xdEDuqa2Vgcws1P2QuU\nACLiauBqAEnHAbvUFKfcXh0O15nZWPw+kY8yF6iqK0lUFkbExdmV8r6ShJmZTQxfSYLWkuAR4VN8\nZh3EV5Iotyl5JYk6U76vkXSzpOvT03yV9QoJ6i5YcL7zDWbWkMO67SvzKb5KUPfqzO23RcR9VStF\nnC7p6TwarAR1t29fChzHlVd+kYGBAU9YaGZVKh9mN2w4AICFC/u56KKP+IxLi0p7BJWq3duF7n0H\ndc2sGQ7r5qPsBSrrMeBSSeskvbiIBqqDujioa2ZWoClToCLipIg4CPgYcF4RbWSDuoCDumZWl8O6\n+Sjzd1CNPAnk8p1TrWxQd/58WL36NJ9TNrMRHNbNx5QpUJIuA/YkOdX33tq7c2yHefPmccopp/gF\nZ2YNOazbvjIXqNop3/+m3kpFBXWXL1+e5+bMzKyGg7o4qGtWZg7qltt0Ceq+MA3p3ijpnzPrFRLU\nXbToHgd1zWxUDuu2p8yn+GqDuquAj0bE96tWKiyoC1deeZGDumZWl8O67SvtEVRKAJK6gdm1xSlv\nDuqaWbMc1m1f2QtUxR+RzAd1haQNkt5URCPVQd23OKhrZlagqVKgtgKPAm8B3gCcJmlG3o1UB3Uv\nd1DXzBpyWLd9Zf4OalhEPC3pfmCPiPi1pKeKaMdBXTNrlsO67ZsSBSp1KvAFSc8DvhIR2zL3FTKj\nrl9sZjYah3XbU+YCVRvUvY8kkFvFM+qamZWTg7qML2jX09PD4OBgEd0xsxY4qFtuUz6oK+m5aUi3\n8vNoZr3cgroAQ0ND9Pf3s2TJEoaGhvLYpJlNMQ7o5qPMp/hqg7oLACTtD7xveKUcg7pDQ0Pst99S\nHnzwGOBoNm48jttvv5iurtLWeTPLmQO6+Sn7O2u9PX4S8KkiGlu3bl1anI4HjueBB97KunXrimjK\nzErKAd38lL1AVZG0G7BXRNw+2X0xM7P2lPkUXz3vAj5f1MZ7e3vZc8+lPPhgcnvWrK/R23txUc2Z\nWQklAd0r2LAhuZ0EdEcMMLYmTJkCJWkH4ChgflFtdHV1cccdlwyf1uvt9fdPZlbNAd38TKV316OB\ndRFRb2hdbq+Orq4ulixZwl133eXiZGZ1VQK68+bNc3FqQ5nfYStB3TcDRMTXImJl7UppUPeNwBN5\nNn7uuefmuTkzM6vhoC4O6pqVmYO65Tblg7rp7aWSfijpJkkLMusVEtSt/G5m1ogDu+0pbYHimaDu\n5entDwOvBRYDZw+vFHE6sCaPBitB3b6+5Pb++x/nImVmdVUCu4sXb2Hx4i0sW3aOi1SLylygoHrw\nwybgYJKRfD8oorHqoC4O6ppZQw7stm/KDDMHrgH+HngW8Jnimzuj+CbMzKaxKVGgJO0LHBURb0xv\nf1fS+po5odpWHdTdx0FdM2vIgd32TYkCBXSTPhcloYMZJN9R5cpBXTNrlgO77ZsSBSoi7pb0A0lX\nk3yv9pmIeDKzSu5BXTOzsXhG3faUuUDVzqh7dr2VPKOumVk5OaiLg3ZmZeagbrlNl6DuCZJulnSt\npDmZ9QoJ6vb19TkDZWZ1OaCbjzKf4hueUVfSzkBfRLxa0u7A54BjoMgZdfvZuPEPnlHXzKp4Rt38\nlP2dVZl/d5S0E8l3U3tI2jHvxhzUNbOxOKCbnzIfQQ2LiCcknQ18E3gM2BV4PvCbSe2YmZmN25Qo\nUAAR8XXg6wCSbo2I3IuTZ9Q1s7E4oJufKVOgKiQdAdxWxLazQd2+Pvz9k5mN4IBufqZMgZL0ReCl\nwOPAsbV359VOJah7yimnuDiZWV0O6OajzAWqNqj7jnorFRXUXb58eZ6bMzOzGg7q4qCdWZk5qFtu\npQ7q1gnkzpd0i6TzatY7TNKN6c+hmeWFBHX7+/sd1DWzhhzWbV8ZTvENB3LT2zsB55DMnguApC7g\nTOCwdNG3JV0fiYKCujBr1nEeKGFmIzism4+yvLMO79WIWA8M1tw/B/h5RGxL54DaDLw4705UB3WP\nd1DXzOpyWDcfZTiCakYP8IikT6S3HwV2A+4urskVwD7Fbd7MbJqbKgVqK8mVI04kOdr6LPDbvBup\nDuqeyaxZRzqoa2YjOKybj7IWqNoTuZuBl2Ruz4mIX+TdqIO6ZtYMh3XzUboCJelUYDHJBWGfGxEn\nRMR2SWcC16arrah9WF7tZ2fUdXEys0Yc1m1fGQpUbSD348DHa1eKiGuAa2qXe0ZdM7NyclCX8QXt\nenp6GBysHUxoZhPNQd1ymy5B3UbLCwnqLlmyxEFdMxuTA7vj1/EFimeCupentytB3Vp1l0fE6cCa\nPDpSCer29UF//9Hsv/9xLlJm1lAlsLt48RYWL97CsmXnuEi1oAwFCsYO6jZcnicHdc2sFQ7stqcs\nBcrMzKaZMozi6xieUdfMWuHAbnvKWqAajfoodDRPNqgL0NvroK6ZNebAbntK9+6aBnVXAL2S/nWs\n5ZW782q/EtS96667XJzMbEyVwO68efNcnFrU8TkoSe8H/hY4NzOSr5XHV4K6H4yIG+qt4xyUWXk5\nB1Vuo+2/ji9QE8EFyqy8XKDKbboEdS+UdL2k70jaN7O8kKBu5Xczs0Yc0m1PxxcomgzqRsS7I2IB\nycy6J2eWFxLUBRzUNbOGHNJtXxkKFDQR1M14DPh9EZ2oDurioK6ZNeSQbvvKUqBasQz4XPHNnEGl\nUJmZWf7KmoOqS1Iv8LOI+GkR268O6u7DrFlHOahrZnU5pNu+shaoEaM+JB0IHBwRHy6qUQd1zaxZ\nDum2r3QFqt6MuuldXwXul3Q9cEdEnJR9WF7tZ2fUNTMbjWfVbU8ZClSzM+ruO/KhnlHXzKysHNTF\nQTuzMnNQt9ymS1D3LEkbJK2fiKBuX1+fM1Bm1pBDuu0rwym+SlD36vR2Jaj72qqVkkAukg4CTgVO\nqCyX9HQeHakEdZMsVD8bN/6B22/3QAkzq1YJ6W7YcAAACxf2c9FFH/EgiRaV5Z21laDuq4G7iuiE\ng7pm1gyHdPNRhiOopkn6LrA7MH+y+2JmZu2ZUgUqIl4v6VXAJUAu3zlleUZdM2uGQ7r5KGuBGu1E\n7kMU9LyyQd2+Pvz9k5nV5ZBuPkpXoBoFdSV9meT03u+Bv6t9WF7tV4K6p5xyiouTmTXkkG77Oj4H\n1akz6ppZZ3AOqtw8o+4Y/AIzKy8XqHKbFkHd9L6dJP1S0nszywoJ6vb39zuoa2Ytc4C3eWX4Dqqp\noG7q3cCP08ckDy4sqAuzZh3ngRJm1jQHeFtTlnfWMYO6knYGDgf6yXFQRFZ1UPd4B3XNrCUO8Lam\nLAWqGScBn5645lbgGXXNzIozJQqUpOcBr4uIb1HQ0RNUgrpfBdYAZ6Yz6vYW1ZyZTTFJgPc2Zs5c\ny8yZa9MA79zJ7lbHKsN3UPXUFqGDgGdL+hLwImAHSTdExECejTqoa2btcIC3NaUrUPWCuukAiqvT\n+48DdqkpToXMqOviZGatcoC3eWUoUE3NqFsREVUXx/OMumZm5eSgLuML2vX09DA4ONqsH2Y2ERzU\nLbdpEdSVtEbSzZKuT0/zVZbnGtSthOz6+vocsjMzK1AZTvE1G9QN4G0RcV/VwhyDupWQ3bXX7g8c\nwuDgOQ7ZmdkIEZEZCDHX7xHj1PFHUKlmZ9Qt9FUwMDDAVVfNZdu2pWzbtpQrr/xzBgZyHShoZiVX\n+SC7ePEWFi/ewrJl5/hsyziVpUA14zHgUknrJL24iAY2b97M9u3P/Jdt397F5s2bi2jKzErKV4vI\nTxlO8TUlIk4CkHQAcB7wprzbmD17Nt3dF7J9+7MA6O5ez+zZJ+TdjJmZUd4jqNFO5T0J5PKdU615\n8+ZxxBF7MGPGFmbM2MKRR+7hPIOZVfHVIvJTuiOoUWbUvQzYk+RU33trH5ZT26xZ81E2bdrERRdd\nxKpVq/zlp5lV8dUi8tPxOahOnVHXOSizzuAcVLl5Rt0xuECZlZcLVLlNl6DuC9OQ7o2S/jmzvJCg\nbuV3MzMrRscXKJ4J6lZO71WCurVWAR+NiPkR8cHhB0ecTjI/RvsdSfMNixbdA+B8g5nlzlPCP6MM\nBQrGCOpK6gZmR8T3i+xENqgLOKhrZrlyyLdaWQrUWP6IZD6oKyRtkJR7Bgpqg7pnOKhrZrlyyLda\n6YaZN7AVeBR4C9AN3CTpWxGxLc9GqoO6L6G7+2oHdc3MClLWI6iqUR8R8TRwP7BHRPweeKqIRh3U\nNbMiOeRbrXRHUI2CusCpwBckPQ/4Ss3RU+5BXfBVis0sXw75Vuv4HFSnBnXNrDM4B1VuDuqOwS8w\ns/JygSq3KR/UlfTcNKRb+Xk0c18hQd0PfehD03r4p5lZ0crwHdSYM+pGxO+ABQCS9gfel7mvkBl1\nt21bzeDgbp5R18wa8sy67en4I6hUszPqApwEfKqITjioa2bNcui2fWUpUE2RtBuwV0TcXsT2PaOu\nmTXLodv2TakCBbwL+HxRG0+CutcClwKVGXVnF9Wcmdm0VtYCNeJErqQdgKOAbxTVaDaoCzioa2YN\nOXTbvjIMkqgySlD3aGBdRAzVe1hObWdm1O1j1arT/KWnmdXl0G37Oj4H5aCumY3GOahyc1B3DH6B\nmZWXC1S5Tfmgbrp8qaQfSrpJ0oLM8kKCup5MzMysWGX4DmrMoG7qw8ArgF2AbwOvgeKCugCLFvU7\nqGtmTXNwtzUdfwSVaiaouwk4mGQk3w+K6EQ2qLtt21IHdc2saQ7utq4sBaoZ1wB/DywFNhTRQHVQ\ndwXbty9zUNfMmuLgbuvKcIpvTJL2BY6KiDemt78raX2xM+qeSXf3LzyjrplZQcp6BFV74rabtNgq\nOak7g+S7q1w5qGtm4+Xgbus6fpi5pDOAH0XEVent4aAu8J1KUFfSacDrSIruZRGxJrONFcDGyjZq\ntTJMtPIl5/z589m6dau/5DSbZGUaZu5BEiOVOgfVqUHdnp4eBgdHu6i6mU2EMhUoG6nUBWoiuECZ\nlZcLVLmVOqjbaSpB3b6+Pg8RNTMrUMcXqBauJHGCpJslXStpTmZ5bleSqOQYFi26h8suO8Q5BjMb\nla88056OL1A8cyWJyvdPlStJDJO0M9AXEa8h+b7q7OEHR5wOrMmjIw7qmlmzHMxtXxkKFIx9JQkB\nO0raCXiEZCqOHfPuhGfUNbNmOZjbvikR1I2IJySdDXwTeAzYFXg+8Js826kO6lZm1HVQ18ysCFOi\nQAFExNeBrwNIujUici1O8ExQd/36LQAcfriDumZWXxLMvYIN6YXXkmDuEZPbqZIpa4FqPG5eOgK4\nrZBGq2bUvYhVq1Y5aGdmdXlG3faVrkA1mvJd0heBlwKPA8fWPizH9pk3bx6rV6/m/PPPz2uzZjYF\nVd4vbHzKUKAeAT4maaeIuDwiPg58vHaliHhHvQdnriRxQ6G9NDOzXPlKEvhafGZl5itJlFupryRR\nJ6h7oaTrJX0nnWajst5hkm5Mfw7NLC8kqAs412BmY3JYd/w6vkBRE9SNiHdHxALgTOBkAEld6e1F\n6c+KdNqNwoK6gIO6ZjYqh3XbU4YCBfUHOTwG/D79fQ7w84jYlk5SuBl4cd6dqA7qnuGgrpmNymHd\n9pRhkEQjy4AL0t97gEckfSK9/SiwG3B3ng1WB3VfQnf31Q7qmpkVpJQFSlIv8LOI+Gm6aCvJlSNO\nJDna+iz7UvnhAAAKpElEQVTw27zbdVDXzFrhsG57SlegJB0IHBwRH84s3gy8JHN7TkT8ooC2h4O6\n4BkxzWx0Duu2p3QFCvgqcL+k64E7IuKkiNgu6Uzg2nSdFTWPyT2oa2bWDL9njF8ZClRtUHffeitF\nxDXANbXLHdQ1MysnB3UZX1DX1+Iz6wwO6pbbdAnqNpppt5Cg7urVq51pMLO6HM7NR8cXKJoI6qZG\nzLSbru+grplNGIdz81OGAgVjB3UbzbSbK8+oa2ZjcTg3P2UpUPUsAz43kQ0mQd1rgUuByoy6syey\nC2Zm00YpC1SdoO6EqAR1Z8zYAsCRRzqoa2bVknDubcycuZaZM9em4dy5k92tUirDMPMqDYK6w3cX\n3HZmRt0+Vq06zaP4zKyKw7n5KV2Bok5QFxrPtJvKPajr2XTNrBGHc/NRhgLVbFC37ky7DuqamZWT\ng7o4aGdWZg7qltt0Ceo2Wp5bUBdgaGiI/v5++vv7GRoaymOTZjYNOLzbujKc4qsEda+GJKgLkE7r\nfjLwnjGWny7p6Tw6MjQ0xH77LeXBB48BYNas47j99ovp6ur4Om9mk6gS3t2w4QAAFi7s56KLPuLB\nE2MoyzvrmEHdJpa3bd26dWlxOh7YwgMPvJV169YV0ZSZTSEO745PWQpUPY2CuhMU4D2z+CbMzKax\nMpziG6FRULfoAG9vby977rmUBx9Mbs+a9TV6ey8uoikzm0I8s+74lK5ANQrqjhHgzUVXVxd33HEJ\n69ato68Pf/9kZk1xeHd8SlegaBDUHWU55BjU7erqYsmSJcO/m5k1w+Hd1pWhQDUb1K273EFdM7Ny\nGvMQQNI7i+6EpI2N7ouICyLi1UBvozzTaH2MiNMj4i8i4oZ8eps45ZRT8tycmZnVaOYc1bsK70UT\nIqKPxhMPTngfly9fPtFNTriigoUOLNp049f8+IxaoCRdArw0vULD6emyPkn96dUdTsqsu7He76Ns\n+yOSNkpaA8zMLH+DpO9LuknSW5rYTraPH8ssr9tPa05Rs4J6tlGbbvyaH78xr8UnaWNEvDJze8eI\neFrSTsAtEfHy2vVqH1Nnm3sClwMHAbuQDGrYR1IXcGu6/ClgA3B4RDyVPu4M4EcRcdVofRytn/WU\n9VpaPT27TnYXprzBwYcnuws2hk6/Ft/AwACLF2/h8cePBWDmzLV885v7eMBEarT9N55BEq+XdBTw\nOLDzOPu0F0mhGQIek/SbdPnuwAuBK9PbzwdmAfdOUj/HzcVjavB+nDz+cGDNFKhnSepKiwnABcB+\nwN7A32TW6wKQtDNjF4R7gFekR0w9wJ7p8t8AdwFLIuJ3zT2Fun0crZ8Toux/XCOvHfaTXK4dVtR2\nzTqVQ7rj18wpvlXAy4B7IuJESf9K8sb/X8CB6Qg7JH0OeILkiOWtEfGyMbZ7CnAMcCfwysr6kg4C\n/onkIrEPRMTbM485A3gLcHFEnJ9ZXunjvRHxnnRZ3X7WM55TfCtXrpzyAyUiIhMsnJtbESlquzY9\ndfopPvBrfjSj7T/PB8X4ClRPTw+Dg4NFdMfMWlCGAmWN5f0dVNPSqzrUiog4tMh2zcys/AotUBGx\noMjtm5nZ1FWGSx11hKGhoeG5n3p7eye5N2bWyfydUz6mRIGSdDxwGPAK4DPA+4A3AF9rNps1mnoz\n6ZqZ1ePZc/MzJQoUyYi/e0lG7D0H+BJwQF4br55JFx54AJYsmSr/dWaWp+zsuQDXXZcsczC3dVPp\nXfah9N/HgT0o+LkdffTRRW7ezGzam0oFqp5WwsMNeSZdM2uWg7n5mRI5KEnH8cwFZytHUL8g+V5q\nzPBwMzmG2kESnqzQrDN0Yg7KgySa56DuGBy0MyuvTixQ1rzR9p8PA8zMrCO5QI3TypUrJ7sLZmZT\nmk/xMb7voHbffXdfi8+sA3TiKT5/B9W80p3ik3S8pLWSBiSdKOkuSX/awrp7p/f9KLPemLP8NlIJ\n6vb1QV8f7L+/g7pmVp9n0M1Ppw4zrw3eXkYSvP1lk+u+Argvva9t9YK6cNUojzCz6cpB3fx0aoGC\n1oK3ExrSNTOz4k31N/LCgrrJUZSZWTUHdfNTpgLVyum6yrq3pLPtPt7i46t0dXVxxx2XZAZJXMy5\n55473s2Z2RQmidWrT8sMkjjCgyTGyaP4cNDOrMw6cRSfNW/SZtTNk2fnNTObXia9QEl6Z0R8Yaz1\nPDuvmdn00gk5qHdNdgfaERHceeed9Pf3c+eddzrvYGYNRQQDAwMMDAz4vaIJk1qgJF0CvFTS9ZJO\nT5f1SeqXdLukkzLrbqz3e4PtNh30bUdEcPzx/8SCBRfS1wcLFtztUJ6Z1eUAb+sm9RRfRCxNp2LP\nnr5bGxGrJe0E3AL8y3g2TfNB33FZuXIlvb29rF+/G9u3nw/A9u3Q3w+bNn3PoTwzq+IAb+s64RRf\nrddL+gRwOm1kl0jCu09m/s21GFeGmW/bdtqI++bPfx09Pbvm2ZyZ2bTTCQXqWZKy/bgA+CDwbzXr\n5RK6zdPcuXPp7f0Q3d1vB9bQ3X0RS5aczNatgwwOPjzZ3TOzDpIEeG9j5sy1zJy5Ng3wzp3sbnW0\nSR/FB1wLXC3pnog4Ebgp/fkvYGtmvXZDt7mf7JXEmjUfZWBggM2bNzN79hzmzTvaoTwzG8EB3tY5\nqMv4gnY9PT2ebsOsAzioW25TIqhbj8O7ZmZTl4+ggOuuu87/CWYltnDhwkk7ivL7R/sa7T8XKDMz\n60idMIrPzMxsBBcoMzPrSC5QZmbWkVygzMysI7lAjYOkwyTdmP6UYki7pDWSbk4vzLs0XVaK5yFp\nvqRbJJ2XWVa37536nBo8h+w+OS6zvFOfw4VpX78jad90Wan2w0Spt787oa1Gr7mc+zPidTJuEeGf\nFn5IivpNwIz057ukoyE7+QdYDexdxucBHAa8CTivUd87/TnVPod6+6TTn0Omj4cCnwNUtv0wmfu7\nE9qq95or+nXSzjZ8BNW6OcDPI2JbRGwDNgMvnuQ+NSubNSjN84iI9UD2sh0j+i5pTr3ldMhzqvMc\nKmrzHx37HDIeA35PCffDRBllf3dCWxOVGau8TpJGpUNaPWor9ZUkJkkP8Eh6xXWAR4HdgLsnr0tN\neQy4VNIg8AHK+zygcd/VYHmnPqeqfRIRv6Ac+2UZyUWdd2Nq7IfppN5rrijLgAskvRz4JPB84NmS\njgf+KS2so3KBat1Wkv/oE0n+ED8L/HZSe9SEiDgJQNIBwHnAqZTweaQa7YOuBss7Up198iY6/PUl\nqRf4WUT8VNJLmAL7YTpp8JrLXfZ1ki5aIOlgYJ+IuLjZ7bhAtW4z8JLM7TkFfwrJ25PA08AvKNfz\nyJ6WqLsPJHXXWz4hvWtOo1MrlX0CHfz6knQgcHBEfDhdVNb9MFEm8vJLrbaVfc3l25GRr5Nxc4Fq\nUURsl3QmyTQhACsmsTtNk3QZsCfJdCUnRsRQWZ6HpFOBxcAekp4bESfU63sn75sGz+HLwB4kp13e\nC539HICvAvenF2m+PSLeX7b9MFHq7e+JbkvSMcD/RsRVmXUr7wPDr7kCZF8nd1SO2iLiO8B3WtmQ\nr8VnZmYdyaP4zMysI7lAmZlZR3KBMjOzjuQCZWZmHckFyszMOpILlJmZdSQXKDMz60guUGZm1pH+\nPy8uIB4zRxUWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f7404a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.Matplot.summary_plot([M_random.m_delta, M_random.delta, M_random.tau_delta, M_random_treat.mu, M_random_mean.m_mu, M_random_mean.tau_mu])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DIC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 527.29959017,  281.75127105,  439.03270863,  285.24870453])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dic = np.array([M_fixed.dic, M_random_mean.dic, M_random_treat.dic, M_random.dic])\n",
    "dic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.   ,  0.852,  0.   ,  0.148])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_dic = dic - dic.min()\n",
    "np.round(np.exp(-0.5*delta_dic) / np.exp(-0.5*delta_dic).sum(), 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, 85% of model weight is on the random mean model, while 15% is on the full random model."
   ]
  }
 ],
 "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.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}