Lesson 0 - Coin Toss.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Statistics\n",
    "\n",
    "Bayes rule\n",
    "$$\n",
    "p(\\theta| D) = \\frac{p(D|\\theta)p(\\theta)}{p(D)}\n",
    "$$\n",
    "However, in the samples that we will soon see it is far more convenient to write the above as\n",
    "$$\n",
    "p(\\theta| D) \\propto p(D|\\theta)p(\\theta)\n",
    "$$\n",
    "since $p(D)$ is a constant **that does not depend on $\\theta$**.\n",
    "\n",
    "However, we Bayesian analysis is far more involved in solving the following integral:\n",
    "$$\n",
    "p(D^*|D) = \\int p(D^*|\\theta)p(\\theta|D) d\\theta\n",
    "$$\n",
    "We denote $D^*$ to be test data, $D$ to be the train set and $\\theta$ to be a hidden parameter. This can often be thought of as doing model averaging if $p(\\theta|D)$ is thought of as a weighting function.\n",
    "\n",
    "## Coin Tossing Example\n",
    "\n",
    "Suppose we have a biased coin where the true probability of landing heads is 0.6. We assume however, (for some odd reason that it is biased towards landing tails). A suitable prior for this would be a beta distribution which has the following functional form:\n",
    "$$\n",
    "p(\\theta) = \\frac{\\Gamma(\\alpha)\\Gamma(\\beta)}{\\Gamma(\\alpha + \\beta)}\\theta^{\\alpha - 1}(1-\\theta)^{\\beta - 1}\n",
    "$$\n",
    "and the distribution looks like (https://en.wikipedia.org/wiki/Beta_distribution):\n",
    "![beta](./images/Beta_distribution_pdf.svg.png)\n",
    "\n",
    "Since we believe it is left skewed we will choose $\\alpha < \\beta$. Let $\\alpha = 2$, $\\beta = 3$ and hence, **this forms our prior distribution**. Assume that after 10 coin tosses we observe 7 heads and 3 tails. Each (individual) coin toss is distributed as a bernoulli distribution, conditional on $\\theta$, the probability of achieving a heads. At this point a frequentist method (I don't mean this in a derogative fashion) would state that the probability of receiving a heads is 0.7.\n",
    "\n",
    "The likelihood of the data can be states as follows:\n",
    "$$ p(y_1, \\cdots, y_{10}| \\theta) = \\prod_{i=1}^{10} p(y_i|\\theta)$$ since each draw is independent **given $\\theta$** where, $p(y_i|\\theta) = \\theta^{1(y_i=1)}(1-\\theta)^{1(y_i=0)}$. Hence in our example, the likelihood $ p(y_1, \\cdots, y_{10}| \\theta) = \\theta^7(1-\\theta)^3$. **Keep in mind that $\\theta$ is considered fixed as far as the likelihood is concerned**, hence it is a bernoullii distribution, **not a beta distribution**.\n",
    "\n",
    "### The posterior\n",
    "\\begin{align}\n",
    "p(\\theta | y) \\propto & \\quad p(y|\\theta) p(\\theta) \\\\\n",
    "\\propto & \\quad \\theta^7(1-\\theta)^3 \\theta^2(1-\\theta)^3 = \\theta^{9}(1-\\theta)^6\n",
    "\\end{align}\n",
    "The normalising constant is irrelevant as we can see that the posterior $p(\\theta | y)$ is a beta distribution simply by looking at its functional form. In fact the normalising constant (also known as the partition function) is simply,\n",
    "$$\\int_0^1 \\theta^{9}(1-\\theta)^6 d\\theta$$ since all probability distributions integrate out to one. In this case this turns out to be $\\frac{\\Gamma(9)\\Gamma(6)}{\\Gamma(15)}$. What is far more important to recognise than the normalising constant is that $\\theta \\sim Be(9, 6)$.\n",
    "\n",
    "### The probability of $p(y^*=1|y)$\n",
    "\\begin{align}\n",
    "p(y^*=1|y) &= \\int p(y^*=1|\\theta)p(\\theta|y) d\\theta\n",
    "\\end{align}\n",
    "\n",
    "A frequentist method would have chosen the Maximum-A-Posteriori (MAP) estimate and plugged into $\\theta$ without integrating it out. \n",
    "$$\n",
    "\\begin{align}\n",
    "p(y^*=1|y) &= \\int p(y^*=1|\\theta)p(\\theta|y) d\\theta \\\\\n",
    "&= \\int \\theta \\frac{\\Gamma(15)}{\\Gamma(9)\\Gamma(6)} \\theta^{9}(1-\\theta)^6 d\\theta\\\\\n",
    "&= \\frac{\\Gamma(15)}{\\Gamma(9)\\Gamma(6)} \\int\\theta^{10}(1-\\theta)^6 d\\theta\\\\\n",
    "&= \\frac{\\Gamma(15)}{\\Gamma(9)\\Gamma(6)} \\frac{\\Gamma(10)\\Gamma(6)}{\\Gamma(16)} \\int\\frac{\\Gamma(16)}{\\Gamma(10)\\Gamma(6)}\\theta^{10}(1-\\theta)^6 d\\theta \\\\\n",
    "&= \\frac{\\Gamma(15)}{\\Gamma(9)\\Gamma(6)} \\frac{\\Gamma(10)\\Gamma(6)}{\\Gamma(16)}  \\quad 1 \\\\\n",
    "&= \\frac{9}{15}\n",
    "\\end{align}\n",
    "$$\n",
    "Note the rule $\\Gamma(n) = (n-1)\\Gamma(n-1)$ was used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pymc3 as pm\n",
    "\n",
    "from scipy.special import gamma\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using ADVI...\n",
      "Average Loss = 8.1227:   1%|          | 1292/200000 [00:00<00:15, 12915.73it/s]\n",
      "Convergence archived at 2500\n",
      "Interrupted at 2,500 [1%]: Average Loss = 7.9508\n",
      "100%|██████████| 100500/100500 [00:37<00:00, 2691.01it/s]\n"
     ]
    }
   ],
   "source": [
    "alpha = 2\n",
    "beta = 3\n",
    "y_obs = np.zeros(10)\n",
    "y_obs[:7] = 1\n",
    "niter = 10000\n",
    "\n",
    "with pm.Model() as model:\n",
    "    θ = pm.Beta('θ', alpha=alpha, beta=beta)\n",
    "    y = pm.Bernoulli('y', p=θ, observed=y_obs)\n",
    "    \n",
    "    trace = pm.sample(niter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# norm_const = gamma(9)*gamma(6)/gamma(15)\n",
    "theta = np.linspace(0,1,1000)\n",
    "post_theta = (theta**9)*(1-theta)**6\n",
    "C = np.sum(post_theta*(theta[1]-theta[0]))\n",
    "post_theta = post_theta/C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd8VHW+//HXZyYJSFcJiIBE6SUR\nMKKgwQIqYtu1rOWu5a67yOre9W67F9e9rust692irqs/V1xde8VCL4poolKM9CqRGkENLYC0zDmf\n3x8zcEMIZJLMzHfK5/l4zIPJnJMz70OST775nu/5fkVVMcYYk14CrgMYY4yJPSvuxhiThqy4G2NM\nGrLibowxaciKuzHGpCEr7sYYk4asuBtjTBqy4m6MMWnIirsxxqShLFdv3LZtW83Ly3P19sYYk5I+\n++yzLaqaW9d+zop7Xl4epaWlrt7eGGNSkoisj2Y/65Yxxpg0ZMXdGGPSkBV3Y4xJQ1bcjTEmDVlx\nN8aYNGTF3Rhj0pAVd2OMSUNW3I1JFr4PXsh1CpMmnN3EZIwBNsyBBS/A2hLYsSH8WquOkHcuDLwZ\nupwDIm4zmpRkxd0YF7atgaljYPV0yGkJXS+Agu8BAlvL4PNpsPhVyCuCyx+Btt1cJzYpxoq7MYm2\nfAK8c2f4+UUPwJk/hJzmh+9TtRfmvwCz/gueHApXPwm9r0h8VpOyrM/dmESa8zd4/WbI7QF3zYFz\n7j6ysANkHwdnjYI750C7XvDa92HeU4nPa1KWFXdjEmXOEzDt36HX5fDPU6F1p7o/p9XJcNsU6DkS\npvwSPnsu/jlNWrDibkwirJgE0+4JF/brnoWsJtF/bnbT8Od0Gw6TfgZrPohTSJNOrLgbE2+bF8Nb\nP4KOA+Gav0Mwu/7HyGoSLvBte8Abt8H2qGZ9NRnMirsx8bR/d7gYN20DN7wc7ktvqCYt4YaXwuPh\n3x4NvhezmCb9WHE3Jp6m/zo87PHqsdDypMYf78SucOn/woZPwn34xhyFFXdj4mXVNJj/HJzzUzi1\nKHbHPf0G6HkZzHwAtpTF7rgmrVhxNyYeDuyBKb+C3N5wwW9ie2wRuOIRyGoaHn2jGtvjm7Rgxd2Y\neCj5E1RugMsfgqyc2B+/RTs4fwyUvQefT4/98U3Ks+JuTKxVfA4fPwqn3wRdhsTvfQb9CNr2hGlj\nIHQgfu9jUpIVd2NibcZvILtZeGqBRsobM5m8MZNr3xjMhkv+B7avhQXPN/q9THqx4m5MLK3/JDwZ\nWNHPoEVu/N+v2zA4ZTAU/yk8H40xEVbcjYkVVXjvfmjZAQbdUe9PP2YrvZb98sZMDl9cveBe2LUZ\nSp9pQGiTruqcFVJEmgLFQJPI/uNU9bc19mkCPA+cAWwFrlfVdTFPa0wy+3wabJwbnqI3p1lMD33M\non9qEZx6Hnz0MBTeHp6uwGS8aFru+4ELVfV0oD8wQkTOrrHP7cB2Ve0GPAz8b2xjGpPkVGHWf8MJ\np8GA7yf+/Yf+Er6tgEWvJP69TVKqs7hr2O7Ih9mRR82BtVcBB6erGwcME7HlY0wGKXsPvloCRb84\n6twx0Xa7NEheEXToD7MfC09PYDJeVIt1iEgQ+AzoBjyuqnNr7NIR2AigqiERqQROBLbEMKsxyavk\nIWjVCfK/l9C3rf7L4vJAEY/l/BVWTYHelyc0h0k+UV1QVVVPVfsDnYBBItKvxi61tdKPuG1OREaJ\nSKmIlFZUVNQ/rTHJaP3s8FwvQ/4lPjcsRWmqP4gNfi58/BdnGUzyqNcye6q6Q0Q+AEYAS6ttKgc6\nA+UikgW0BrbV8vljgbEAhYWFds+0SQ8fPQTNToSBt8TkcA3tuvEI8rQ3kt+VPwdfzg9PMWwyVp0t\ndxHJFZE2kefHAcOBlTV2mwDcGnl+LfC+qk14YTLA18th9Qw4+8cxHyHTEG95ReEbqEqfdh3FOBZN\nt0wHYJaILAY+Bd5V1Uki8oCIXBnZ52ngRBEpA34OjIlPXGOSzLwnwxN4Fd7uOgkAu2jGy3vPZu/8\n12HvdtdxjEN1dsuo6mJgQC2v31ft+T7guthGMybJ7d0Oi1+H/Oug2QkNOkQ8Rs+86A3npqz3YeEr\nMPjOmB/fpIZ69bkbY6pZ8CJU7YGz6n83ajwt1zzm+91oPfVRho0/BRDWPXiZ61gmwWz6AWMawvdg\n3lPQ5Rw4Kd91miO8GBpO18Bmzg6scB3FOGLF3ZiG+Hw67FgPg0a5TlKrqf4gdulxXBModh3FOGLd\nMsY0xLyx0Koj9Kr/zUJxu0u1mr00ZbJ3FlcEZ/Pb0G1xfz+TfKzlbkx9bVsLa2bBGbdBMHnbR296\nQ2ku+xkRmOc6inEgeb8zjUlWC18CCUD/f6pz10S00o/mU+3Jer8d1watayYTWcvdmPrwPVjwEnQb\nDq07uk5TB2GcN5QhweWwY4PrMCbBrOVuTH2UzYRdm+DSY89q7bLFXt3bfhG/YBx//uPv+Kt3NYAN\ni8wQ1nI3pj7mPwfNc6HHCNdJolKuucz1e3FlcDa1zOVn0pgVd2Oitfub8GpLp9/gdPbH+proDaZ7\n4Et6ykbXUUwCWXE3JlqLXgU/BANiM/tjokz1BuGpcHlwjusoJoGsuBsTDVVY8AJ0Phtye7hOUy9b\nac0nfl+uCFjXTCax4m5MNDYtgC2fQ/+bXCdpkIn+YPICX9NP1rqOYhLEirsx0Vj8OgRzoM9VrpM0\nyHTvTKo0aF0zGcSGQhpTFy8ES8eFR8gc1+aouyXL8MfaVNKCEj8/XNxVwdavT3vWcjemLms/gG8r\noOB610kaZaI3mE6yBco/BcK/jA4+TPqxlrsxdVn8OjRtA90vqnVzqhTH9/wzOKBBclZMgM6DXMcx\ncWYtd2OOZf9uWDER+n4Hspq4TtMou2jGJ34/WDEp3DVj0poVd2OOZdWU8GpLKd4lc9AMvxC2r4Vv\nlruOYuLMirsxx7L4NWh9Snh8exp41zsDX4U/P/qQ6ygmzqy4G3M0u7+BL2ZBwXUQSI8flQraMF+7\nc0mw1HUUE2fp8R1rTDwsHw/qQf51rpPE1AzvDPoF1tFJKlxHMXFko2WMOZplb0Nub2jX+4hNqTJC\npjbT/TP5Na9wUaCUf3iXuo5j4qTOlruIdBaRWSKyQkSWicjdtexzvohUisjCyOO++MQ1JkF2fQXr\nPwmPkkkz6/UkVvqdrWsmzUXTcg8Bv1DV+SLSEvhMRN5V1ZqX20tUtf6rBRuTjJZPABT6pF9xB5jh\nn8FdwfEcz06208p1HBMHdbbcVXWzqs6PPN8FrACSfX0xYxpn2dvQrg+06+U6SVxM984kKMrw4HzX\nUUyc1OuCqojkAQOAubVsHiwii0Rkqoj0PcrnjxKRUhEpraiwizkmSe3cBBtmQ9/vuk4SN8s0j816\nAhcEFrqOYuIk6uIuIi2AN4F/VdWdNTbPB7qo6unAX4F3ajuGqo5V1UJVLczNzW1oZmPiK827ZMKE\nWV5/igJLyCbkOoyJg6iKu4hkEy7sL6nqWzW3q+pOVd0deT4FyBaRtjFNakyiLHsb2vVNuUU56ut9\nfwAtZS+FgVWuo5g4qPOCqogI8DSwQlVrva1NRE4CvlZVFZFBhH9pbI1pUmMSYecm2DgHLvjNEZtS\nefhjbT72+7JfsxkWmH/o3NY9eJnjVCZWomm5nwPcDFxYbajjSBEZLSKjI/tcCywVkUXAo8ANqjYz\nkUlBy8eH/03DIZA17aUpc/ze1u+epupsuavqR8AxZ/ZX1ceAx2IVyhhnlr0N7fOhbXfXSRLifX8A\nv8t+ji7yFev1JNdxTAzZ9APGHFRZDhvnQt/vZMxCFu/7/QG4MLDAcRITa1bcjTnoUJdM+g6BrGmj\ntme139GKexqy4m7MQSsmQvt+cGJX10kS6n2/P2cFVtCcva6jmBiyicOMgfD0vhvmwHn/7jpJws3y\nB3BH1mTODSwlb8xxh163kTOpzYq7MRBecQmF3kdOj5Tu/e6lfg92ajMuCCxgun+m6zgmRqxbxhgI\nryvapku4WybDhMii2C/gwuBCBN91HBMjVtyN2bcT1n4Iva8AOeao37Q1y+tPO9lBb9ngOoqJESvu\nxqyeAd4B6JW5M1YX+/kAnBdY7DiJiRUr7sasnATNc6HzINdJnKngeFb4pzDUinvasOJuMlvVPlj9\nLvQcCYGg6zROfegXcEZgFc3Y5zqKiQEr7iazrf0QDuwO97dnuGK/gBzxGBxYBpARd+imMyvuJrOt\nmAhNWsGpQ10nca7U78kebWJdM2nCirvJXL4Hq6ZC94shq4nrNM4dIJvZfh8r7mnCbmIymWvDHNiz\npdYblzJVsV/AsOwFnCJfs0HbA4ffxGV3raYOa7mbzLVyEgSbQLeLXCdJGsV+AYC13tOAtdxNZlIN\n35Xa9QJo0uLQy5l+AXGtnsRGP5fzAot50bNfeqnMWu4mM321GCo3ZPSNS7UTiv0CBgeW2cLZKc6K\nu8lMKyaBBKDnpa6TJJ1iv4AWso+Bstp1FNMIVtxNZlo5CU4ZAs3buk6SdD7x+xLSAEODi1xHMY1g\nxd1knq1fwDfLbZTMUeyiGfO1u11UTXFW3E3mWTEx/G8vG9Z3NMVeAfmBdbSl0nUU00BW3E3mWTkJ\nOpwObU5xnSRpfeifDsC5gSWOk5iGqrO4i0hnEZklIitEZJmI3F3LPiIij4pImYgsFpGB8YlrTCPt\n3Azln0Ivm0vmWJZqHlu1JUOD1jWTqqIZ5x4CfqGq80WkJfCZiLyrqsur7XMp0D3yOAt4IvKvMcll\nVXgc+0XTWrF6avi53XV5JCXAx34/igJLAAUycxGTVFZny11VN6vq/MjzXcAKoGON3a4CntewOUAb\nEekQ87TGNNaKSXBCV1ZrzW9hU1OxX0CuVNrqTCmqXn3uIpIHDADm1tjUEdhY7eNyjvwFYIxbe7fD\nupLIKBlridalxAuvzlRko2ZSUtTFXURaAG8C/6qqO2turuVTtJZjjBKRUhEpraioqF9SYxrr8xng\nh6y/PUpfcwKr/E6Rrpmwg3O8Z/o0DakgquIuItmEC/tLqvpWLbuUA52rfdwJ2FRzJ1Udq6qFqlqY\nm5vbkLzGNNzKidDiJOh4huskKaPEz2dQYBVN2e86iqmnaEbLCPA0sEJVHzrKbhOAWyKjZs4GKlV1\ncwxzGtM4VXuhbGZ4bHvARgBHq8QvoIlUMSiw0nUUU0/RjJY5B7gZWCIiCyOv/Ro4BUBV/wZMAUYC\nZcAe4J9jH9WYRvhiFlTtsbtS62mu34v9mkVRYAnFkbHvJjXUWdxV9SPquPqkqgrcFatQxsTcyknQ\ntDXkFR2xyfqPj24fTfjU73lYv7tJDfb3qUl/Xii8nF6PERDMdp0m5RT7BfQKbKQ921xHMfVgxd2k\nvw2fwN5tNpdMA5VEVmcqClrrPZVYcTfpb+VkyGoK3Ya7TpKSVmpnKrS1dc2kGCvuJr2phot71wsh\np7nrNClJCVDi53NuYAmC7zqOiZIVd5PeNi+Eyo2HltOzG3AapsTL50TZRR9Z7zqKiZIVd5PeDi6n\n12OE6yQp7SM/PBXB0FruVjXJKZpx7sakrpWTocs50PxE10lSWgVtWO53oSiwmCe8Kw/bVr3A2wyb\nycOKu0lfW7+AihXcv+kWnrUWZqMV+/n8IDiVZuxjD01dxzF1sG4Zk74iy+m969lcMrFQ4ueTIx5n\nBVa4jmKiYMXdpK/IcnpfYpPUxUKp35N9mm0LZ6cIK+4mPdlyejG3nxzm+r1tvHuKsOJu0tOqKeF/\nbaKwmCr28+kW2MTJbHEdxdTBirtJTysnwQmnQW4v10nSysGpCM61qQiSnhV3k3727oC1xeEbl8SW\n04ulz7UTX+nxh413N8nJirtJP6sjy+n1tv722BNKvHzOCSwlYFMRJDUr7ib9rJwELdpDx0LXSdJS\niV/A8bKbfFnjOoo5BivuJr1U7YXV70HPkbacXpx85PcDsFEzSc7uUDXpZc0HUPUt9L7c5j2Jk220\nYomfR1FwCY9533UdxxyFNW1MelkxEZq0hryhrpOktRK/gIGymhbscR3FHIUVd5M+QgfC/e29RkJW\njus0aa3YLyBbPM62qQiSlhV3kz7WFcO+SuhzleskaW++351vtQlFNaYisGmAk4cVd5M+lo+HnJZw\n2gWuk6S9A2Qzx+9jF1WTWJ3FXUSeEZFvRGTpUbafLyKVIrIw8rgv9jGNqYMXCi/M0XMEZNt0tIlQ\n4udzWuArOsk3rqOYWkQzWuZZ4DHg+WPsU6KqNomHcWf9R7B3G3d81onpn1q3QCKUVFud6WVv2GHb\nbAEP9+psuatqMbAtAVmMabjl4yG7OR/4p7tOkjG+0JP5Uk88ot/dJIdY9bkPFpFFIjJVRPrG6JjG\nRMf3wkMge1zMfmyUTOIcnIpgGUE812FMDbEo7vOBLqp6OvBX4J2j7Sgio0SkVERKKyoqYvDWxgAb\nZsO3FTZKxoFiv4BWsofT5QvXUUwNjS7uqrpTVXdHnk8BskWk7VH2HauqhapamJtrq+OYGFk+HrKO\ng24XuU6ScT72++Gr2KiZJNTo4i4iJ4mE51UVkUGRY25t7HGNiYrvw/IJ0H04NGnhOk3GqaQFi/U0\nimx+96RT52gZEXkFOB9oKyLlwG+BbABV/RtwLfBjEQkBe4EbVFXjltiY6srnwe6voM93XCfJWMV+\nPncGJ9CKb9lJc9dxTESdxV1Vb6xj+2OEh0oak3jLx0OwCXS/2HWSjFXiFfDTrHcYHFjGdH+Q6zgm\nwu5QNanrYJdMt2HQtJXrNBlrgXZjtza11ZmSjBV3k7rKP4Wd5TZKxrEQWcz2+zI0sBiwHtlkYcXd\npK6l4yCraXhhDuPUh34BnQMVdJGvXUcxEVbcTWryQrDs7XBfu3XJOHdwKgIbEpk8rLib1LSuJHzj\nUv61rpMYYL22Z4OfG+maMcnAltkzqWnpuPD0vt0vtvnDk4JQ4hdwZfATsggRstLinLXcTeoJ7Q/P\nJdPrMsg+znUaE1Hs59NS9jJAyg57/eACHvZLOLGsuJvUUzYzvOJSv2tcJzHVzPb7EtIARUHrmkkG\nVtxN6ln6Jhx3AnS1FZeSyU6as1C72Xj3JGHF3aSWA9/Cqinhse3BbNdpTA0lXj4FsobW7HYdJeNZ\ncTep5fNpULUH+l1j/bhJqMTPJyDKOYFaV+U0CWTF3aSWJW9Cyw7QZYjrJKYWi7QrO7UZ59mQSOes\nuJvUsWcblL0Lfb8LgaDrNKYWHkGK/QIuCC5E8F3HyWhW3E3qWPY2eAeg4HrXScwxzPQG0E52kC9r\nXUfJaFbcTepY9Crk9oYOtgh2MvvAPx1PhWHBBa6jZDQr7iY1bP0ivDBH/xshvPCXSVLbacV87c6w\nwHzXUTKaFXeTGha9ChKA/O+5TmKiMNMbSL/AOtqzzXWUjGXF3SQ/34fFr8Jp50OrDq7TmCjM9AcC\nHNE1Y8NXE8dm9zHJb8Ns2LEBLvwPKwwpYrV2ZIOfy4WB+bzsDXMdJyNZy90kv0UvQ04L6HW56yQm\nasJMfyDnBpbSlP2uw2QkK+4muR3YA8vGQ5/vQE4z12lMPcz0B9JUqhgSWOY6Skay4m6S26opcGAX\nnH6D6ySmnub5vditTRluo2acqLO4i8gzIvKNiNQ6WYSEPSoiZSKyWEQGxj6myVgLXoTWp3Dqkzus\nvz3FHCCbYr+AC4MLsIWzEy+alvuzwIhjbL8U6B55jAKeaHwsY4Dt62DNLBjwfdT+yExJM72BnCTb\n6Wd3qyZcnT8xqloMxxysehXwvIbNAdqIiI1XM403/4Xw2PYB/+Q6iWmgWX5/PBUuDpa6jpJxYtEc\n6ghsrPZxeeQ1YxrOC/FV8dO8Hyog7/eLXKcxDbSNVszzezMi8KnrKBknFsW9tnvBa+1gE5FRIlIq\nIqUVFRUxeGuTtlbP4CTZzquerbaU6qb6Z9Ij8CVd5ctDr9m6qvEXi+JeDnSu9nEnYFNtO6rqWFUt\nVNXC3NzcGLy1SVvzn+MbbcP7/gDXSUwjzfAKAbjEWu8JFYviPgG4JTJq5mygUlU3x+C4JlNVfgmr\nZ/CGN5SQ3USd8r7iRBb43RgRtOKeSNEMhXwFmA30FJFyEbldREaLyOjILlOANUAZ8BRwZ9zSmsyw\n8CVQn9esSyZtTPPOpCCwlo5Yd2yi1NksUtUb69iuwF0xS2Qym+9R/v6TrPP7skHbu05jYmSafyb3\n8AqXBEt5xrv0sG0H+93XPXiZi2hpywYPm+Ty+XQ6yRabbCrNrNeTWOGfwiXWNZMwVtxNcpn3JJv1\nBGb4ha6TmBib5p3JmbKKtlS6jpIRrLib5FGxCtZ8wIuh4XYhNQ1N888kIGo3NCWIFXeTPOY9BcEc\nG9ueplZpZ9b4JzEyMMd1lIxgxd0kh307YdEr0O8attLadRoTF8JEfzCDA8vJZbvrMGnPirtJDote\ngQO7YdCPXCcxcTTBG0JQlMuDR7be7a7V2LLibtzzfZg3FjoWQsczXKcxcfSFdmSZ34Urg7NdR0l7\nVtyNe6unw9YyOGt03fualDfBG8KAQBmd5WvXUdKaFXfj3sd/oVzb0u3lHPuTPANM9AYDcGXAWu/x\nZMXdOJM3ZjJX3/MwbJjN30MjbfhjhthEW+b5Pbkq+DG2QlP8WHE3To3Kmsx2bcFr3vmuo5gEmuAN\noUfgS3rJxrp3Ng1ixd04c5ps4uJAKc97F7GXpq7jmASa4p1FSAOR1vuRbORM41lxN878MDiZA2Tx\nfOhi11FMgm2jFcV+Ad8JfkwA33WctGTF3bixczPXBEsY5w21m5Yy1BveeXSQbRQFlriOkpasuBs3\nPn6EAMqT3uWukxhHZvoD2aYtuC74oesoacmKu0m8nZuh9B+86RWx0eZsz1gHyOYd71wuCpTShl2u\n46QdK+4m4f7xh7upCoV43LvKdRTj2BveeTSREFcFP3EdJe1YcTeJtXMzNwXft1a7AWCFdmGJn2dd\nM3Fgxd0k1kcPE8C3Vrs55HXvfPoF1tFX1rmOklasuJvE2b4OPvsH47yh1mo3h0zwhrBPs7kxONN1\nlLRixd0kzvv/BRLkkdA1rpOYJFJJCyZ6g/lu8CNasueI7XYzU8NYcTeJsWkBLHkDBt/J15zgOo1J\nMs97F9Nc9nN1sMR1lLRhxd3Enyq8ex80OxHOudt1GpOEluhpLPS7cktwBjaZWGxEVdxFZISIrBKR\nMhEZU8v220SkQkQWRh4/jH1Uk7LKZsLaYn5beTl593/kOo1JUs+FLqZrYDPnBJa6jpIW6izuIhIE\nHgcuBfoAN4pIn1p2fU1V+0cef49xTpOqQgdg2hjW+u152RvmOo1JYlP8s9iqLbk1OMN1lLQQTct9\nEFCmqmtU9QDwKmDj2Ex05jwOW1fzu9CtVNl87eYY9pPDK96FDA/MJ082u46T8qIp7h2B6pMul0de\nq+kaEVksIuNEpHNtBxKRUSJSKiKlFRUVDYhrUkplOd+++3tmeGfwgd/fdRqTAp4NjaCKLEYFbXRM\nY0VT3KWW12pe8ZgI5KlqAfAe8FxtB1LVsapaqKqFubm59UtqUs/0ewni80DoFtdJTIrYQmvGeUO5\nJlhCLjtcx0lp0RT3cqB6S7wTsKn6Dqq6VVX3Rz58CrAl7DNd2Xuw/B0eD11FudovchO9p7yRZBPi\ntqxprqOktGg6QT8FuovIqcCXwA3ATdV3EJEOqnqwk+xKYEVMU5qUkj/mDaY1+Xf26smMtSl9TT2t\n0w5M9c/k5uB7/L/QVXzLcQCH3ci07sHLXMVLGXW23FU1BPwEmE64aL+uqstE5AERuTKy209FZJmI\nLAJ+CtwWr8Am+d2T9RInsY1fVd3BfnJcxzEp6MnQFbSSPdwanO46SsqKaviCqk4BptR47b5qz+8B\n7oltNJOSymZyU9Ys/ha6nAXa3XUak6IWa1fe8wZwR9YkXvQuYifNXUdKOXaHqomdvTtg4t2U+Sfz\ncOha12lMins4dB2tZQ+3Z009YpvNN1M3K+4mJvLGTGLy/1xH1Y5N/LJqtHXHmEZbpnlM9gbxg+BU\njmen6zgpx4q7iYl/Cs7ksuA8/hT6Hgu1m+s4Jk08HLqW5uxjdNZE11FSjhV303hfLeW+rBf40Ctg\nrGejGEzslGkn3vKLuC04nVPka9dxUooVd9M4e7bBqzexg+b8vOrHqH1LmRj7Q9X1hAhyb9ZLrqOk\nFPtJNA3nVfHJ7y9j//ZN3HHg52yltetEJg19w/E8HvoOlwRLOSew5LBtBy+s2sXVI1lxNw2jCtPu\nYUhwOfdU3W797CaunvYuZb3fjvuznieHKtdxUoIVd9MwHz8Cnz7F2NBlvOUPdZ3GpLn95PDb0K10\nD3zJXVnvuI6TEqy4m3r71b2/gvfuZ7w3hN+HbnQdx2SID/wBvOWdy53BCfSW9Udst+6Zw1lxN/Wz\n7B0ezHqKYi+fX1aNtguoJqEeqLqZHTTnD9lPkk3IdZykZj+ZJnpLxsG4H7BAuzO66me2+IZJuB20\n5DdVPyA/sI5fZr3mOk5Ss+JuorPwZbxxP2SO14NbDoxhD01dJzIZaro/iBdCw7kjazIXBBa4jpO0\nrLibY8obM4k/33s7vPNjPvH7ctuBf7PCbpz7r9D3We534c/ZT9BJvnEdJylZcTdHF9rPH7Oe5BfZ\n43jTK+IHVf/GPpq4TmUM+8nhzqqfEkB5JvuPtOLbQ9ts7HuYFXdTu+3r4JlLuC6rmIeqruUXVaOt\nj90klXXagdFVP+NU+YrHs/9Cll1gPYwVd3Ok5eOpfGQwO79cxagDP+NR72pqX0rXGLfm+H24J/RD\nioJLeTT7MSvw1VhTzBxSOOZl7s9+lsuDc1mnp/GTqn9ho7Z3HcuYYxrnnUcr9nBf9gsoj3F31U8I\nRUpbJi/NZ8XdgBeCBc/zbpP/oBn7+EPV9Yz1Ljv0A2JMsnvGuxRB+Y/sF2nJH7mr6m520cx1LKfs\npzeTqcLqGfDufVCxktXak3uqfsgX2tF1MmPq7WlvJLs4jv/OeoY3c37L7VW/POwvz4Ot+ExpwVtx\nz0S+B8vHh+eH2byItX57Hgz9jOl+Ida3blLZ694FbNR2PJH9CFNyfs39Vbfypl9EJn5fW3HPAIda\nLGMKYNErsOAF2LGBL/wOjPUZQusrAAAGt0lEQVR+xFtekY2EMWljtt+Xy/b/Dw/lPMGfc/7GCG8e\n/xm6mQ0Zdv1IVNXJGxcWFmppaamT984oOzZw3x//xMWBUoYElhMQ5SOvLy94F/GuX4hvA6ZMmgrg\nc3twCndnvUU2IZ73LubvoZF8zQmH9knFLhoR+UxVC+vcL5riLiIjgL8AQeDvqvpgje1NgOeBM4Ct\nwPWquu5Yx7TiHnt5YyaTy3YKA5/zRNEBWPcRfB1e3OALvwOT/LN5wzuPcm3nOKkxidOO7fwq6zWu\nDpbgEWCiP4Rx3lDm+r1Z8+AVruPVW8yKu4gEgc+Bi4By4FPgRlVdXm2fO4ECVR0tIjcA31XV6491\nXCvujaAK31ZAZTlsWc0Tb0ykh5TTM7CRTrIFgL2aw0K/G+/7/ZnpD2SNnuw4tDFudZJv+FFwMlcH\nP6Kl7GWznsAs73Q+8fsx1+9FBW0ASfrWfCyL+2DgflW9JPLxPQCq+vtq+0yP7DNbRLKAr4BcPcbB\nk7K4H4yrCujhrxHjbd4BCO2Dqn3hf0P7/+/fqj2wrxL2bv+/x74d8O0W1q35nA6yjSbyf6vRHNAg\na/RkVmtHFvpd+czvyTLNs350Y2rRhAMMD8zniuBshgSW0kr2ArBFW7HCP4X12p5NeiKb9UQevv1i\naNIKmrQMP3JaQCCr2iPx3ZrRFvdofvo7AhurfVwOnHW0fVQ1JCKVwInAluji1sOKifDWHcS02Cax\nKg1SSXMqtTnbackmPY1p/pmHvvnWaAfWa3sbk25MlPaTw2T/bCb7ZxPEI1/W0j9QRm/ZQO/AekYG\n5nKC7A7v/MITdRxNIBAMF3o5WqGvZaTO4Lvgwnsbcxp1iqYi1DaGqGZVjGYfRGQUMCry4W4RWRXF\n+9emLfH4xZG0tkPGnTNg55wpnJ7zGmB8wt/1N23hNw095y7R7BRNcS8HOlf7uBOw6Sj7lEe6ZVoD\n22oeSFXHAmOjCXYsIlIazZ8l6cTOOTPYOWeGRJxzNB1GnwLdReRUEckBbgAm1NhnAnBr5Pm1wPvH\n6m83xhgTX3W23CN96D8BphMeCvmMqi4TkQeAUlWdADwNvCAiZYRb7DfEM7Qxxphji+oqnKpOAabU\neO2+as/3AdfFNtoxNbprJwXZOWcGO+fMEPdzdnaHqjHGmPixe8+NMSYNJXVxF5ERIrJKRMpEZEwt\n25uIyGuR7XNFJC/xKWMrinP+uYgsF5HFIjJTRKIaFpXM6jrnavtdKyIqIik/siKacxaR70W+1stE\n5OVEZ4y1KL63TxGRWSKyIPL9PdJFzlgRkWdE5BsRWXqU7SIij0b+PxaLyMCYBlDVpHwQvnj7BXAa\nkAMsAvrU2OdO4G+R5zcAr7nOnYBzvgBoFnn+40w458h+LYFiYA5Q6Dp3Ar7O3YEFwPGRj9u5zp2A\ncx4L/DjyvA+wznXuRp7zUGAgsPQo20cCUwnfJ3Q2MDeW75/MLfdBQJmqrlHVA8CrwFU19rkKeC7y\nfBwwTERSeeLmOs9ZVWep6p7Ih3MI33eQyqL5OgP8J/AHYF8iw8VJNOf8I+BxVd0OoKrfJDhjrEVz\nzgq0ijxvzZH306QUVS2mlvt9qrkKeF7D5gBtRKRDrN4/mYt7bdMe1Fwi6LBpD4CD0x6kqmjOubrb\nCf/mT2V1nrOIDAA6q+qkRAaLo2i+zj2AHiLysYjMiczMmsqiOef7ge+LSDnh0Xn/kphoztT3571e\nknlCkphNe5BCoj4fEfk+UAicF9dE8XfMcxaRAPAwcFuiAiVANF/nLMJdM+cT/uusRET6qeqOOGeL\nl2jO+UbgWVX9c2TCwhci5+zHP54Tca1fydxyr8+0Bxxr2oMUEs05IyLDgXuBK1V1f4KyxUtd59wS\n6Ad8ICLrCPdNTkjxi6rRfm+PV9UqVV0LrCJc7FNVNOd8O/A6gKrOBpoSnncmXUX1895QyVzcM3Ha\ngzrPOdJF8SThwp7q/bBQxzmraqWqtlXVPFXNI3yd4UpVTbL5ouslmu/tdwhfPEdE2hLuplmT0JSx\nFc05bwCGAYhIb8LFvSKhKRNrAnBLZNTM2UClqm6O2dFdX1Gu42rzSMILhXwB3Bt57QHCP9wQ/uK/\nAZQB84DTXGdOwDm/B3wNLIw8JrjOHO9zrrHvB6T4aJkov84CPAQsB5YAN7jOnIBz7gN8THgkzULg\nYteZG3m+rwCbgSrCrfTbgdHA6Gpf48cj/x9LYv19bXeoGmNMGkrmbhljjDENZMXdGGPSkBV3Y4xJ\nQ1bcjTEmDVlxN8aYNGTF3Rhj0pAVd2OMSUNW3I0xJg39f+C0Imkii5ElAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x123497e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(trace['θ'][1000:],100, normed='true')\n",
    "plt.plot(theta, post_theta)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\int \\theta p(\\theta | y) d\\theta \\approx \\sum_{\\theta\\sim p(\\theta|y)} \\theta$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.59970305353900688"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(trace['θ'][1000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}