BayesWithDad_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### July 13, 2015\n",
    "\n",
    "####A brief attempt at a bayesian tutorial for my dad."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pymc\n",
    "\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "For this we are going to use a simple example from this website: http://sabermetricinsights.blogspot.com/2014/05/bayesian-linear-regression-with-pymc.html\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to try to model the mpg of a car as a function of the car's weight, i.e., \n",
    "$$mpg = b0 + b1 * weight + error$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cars_data =  ({\"weight\": [ 4.36, 4.054, 3.605, 3.94, 2.155, 2.56, 2.3, 2.23, 2.83, 3.14, 2.795,\n",
    "               3.41, 3.38, 3.07, 3.62, 3.41, 3.84, 3.725, 3.955, 3.83, 2.585, 2.91, \n",
    "               1.975, 1.915, 2.67, 1.99, 2.135, 2.67, 2.595, 2.7, 2.556, 2.2, 2.02,\n",
    "               2.13, 2.19, 2.815, 2.6],\n",
    "              \"mpg\": [16.9, 15.5, 19.2, 18.5, 30, 27.5, 27.2, 30.9, 20.3, 19.5, 21.6, 16.2, \n",
    "               20.6, 20.8, 18.6, 18.1, 17, 17.6, 16.5, 18.2, 26.5, 21.9, 34.1, 35.1, \n",
    "               27.4, 31.5, 29.5, 28.4, 28.8, 26.8, 33.5, 34.2, 31.8, 37.3, 30.5, 22, 21.5]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# NOTE: the linear regression model we're trying to solve for is \n",
    "# given by:\n",
    "# y = b0 + b1(x) + error\n",
    "# where b0 is the intercept term, b1 is the slope, and error is \n",
    "# the error\n",
    "\n",
    "# model the intercept/slope terms of our model as \n",
    "# normal random variables with comically large variances\n",
    "b0 = pymc.Normal(\"b0\", 0, 0.0003)\n",
    "b1 = pymc.Normal(\"b1\", 0, 0.0003)\n",
    "\n",
    "# model our error term as a uniform random variable\n",
    "err = pymc.Uniform(\"err\", 0, 500)\n",
    "\n",
    "\n",
    "# \"model\" the observed x values as a normal random variable\n",
    "# in reality, because x is observed, it doesn't actually matter\n",
    "# how we choose to model x -- PyMC isn't going to change x's values\n",
    "x_weight = pymc.Normal(\"weight\", 0, 1, value=np.array(cars_data[\"weight\"]), observed=True)\n",
    "\n",
    "# this is the heart of our model: given our b0, b1 and our x observations, we want\n",
    "# to predict y\n",
    "@pymc.deterministic\n",
    "def pred(b0=b0, b1=b1, x=x_weight):\n",
    "    return b0 + b1*x\n",
    "\n",
    "# \"model\" the observed y values: again, I reiterate that PyMC treats y as \n",
    "# evidence -- as fixed; it's going to use this as evidence in updating our belief\n",
    "# about the \"unobserved\" parameters (b0, b1, and err), which are the \n",
    "# things we're interested in inferring after all\n",
    "y = pymc.Normal(\"y\", pred, err, value=np.array(cars_data[\"mpg\"]), observed=True)\n",
    "\n",
    "# put everything we've modeled into a PyMC model\n",
    "model = pymc.Model([pred, b0, b1, y, err, x_weight])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#This requires installing some things so we won't worry about it for now.\n",
    "#import pymc.graph\n",
    "#graph = pymc.graph.graph(model)\n",
    "#graph.write_png('univariate.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 50000 of 50000 complete in 4.6 sec"
     ]
    }
   ],
   "source": [
    "# prepare for MCMC\n",
    "mcmc = pymc.MCMC(model)\n",
    "\n",
    "# sample from our posterior distribution 50,000 times, but\n",
    "# throw the first 20,000 samples out to ensure that we're only\n",
    "# sampling from our steady-state posterior distribution\n",
    "mcmc.sample(50000, 20000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-8.48565608161\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEFdJREFUeJzt3X2sZVddxvHvA6WmSmOtkOlrYjHTyOAbVDtNiDIq1CGa\nvoTQl8RatVHDCBhNDB0VOlGDgEIETRsTWlpUhkxUmiG8dSDcpBpwgBQZmA5tDRM6F2ZKDArGiFP5\n+cfdlx6Gc1/Pyz7nru8nOZm919l7n7Xm3vPcddbea59UFZKktjyt7wpIkqbP8JekBhn+ktQgw1+S\nGmT4S1KDDH9JatCq4Z/k0iQfTfK5JJ9N8uqu/Pwkh5I8kuSBJOcN7LM3yaNJjiW5eqD8iiRHuufe\nOrkmSZLWslbP/zTwO1X1POAq4LeSPBe4HThUVZcDH+nWSbIDuBHYAewG7kyS7lh3AbdV1XZge5Ld\nY2+NJGldVg3/qjpZVZ/ulv8LeBi4GLgGuK/b7D7gum75WmB/VZ2uquPAY8DOJBcC51bV4W67dw7s\nI0masnWP+Sf5AeD5wL8A26rqVPfUKWBbt3wRcGJgtxMs/bE4s3yxK5ck9WBd4Z/kmcA/AL9dVV8f\nfK6W7g/hPSIkaY6ctdYGSZ7BUvD/TVXd3xWfSnJBVZ3shnSe6MoXgUsHdr+EpR7/Yrc8WL445LX8\nIyJJm1BVWXurp6wa/t3J2ruBo1X1FwNPHQRuBd7Y/Xv/QPm7kryFpWGd7cDhqqokX0uyEzgM3AK8\nbRwNmCdJ9lXVvr7rMSm2b35t5bZBE+3bcMd5rZ7/C4FfAj6T5KGubC/wBuBAktuA48ANAFV1NMkB\n4CjwJLCnnrpt6B7gXuAc4P1V9cGNVlaSNB6rhn9V/RMrnxd48Qr7vB54/ZDyTwE/stEKSpLGzxm+\n07XQdwUmbKHvCkzYQt8VmKCFviswYQt9V2DWZJa+zCVJbeUxf0mahM1kpz1/SWqQ4S9JDTL8JalB\na07yktSPla7d9ryYxsHwl2bamflv7ms8HPaRpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalB\nhr8kNchJXtKcW+1bnJwNrJUY/tIUjSOohx9j2GHNfa3M8JembtSg9pYPGp1j/pLUIMNfkhpk+EtS\ngwx/SWqQ4S9JDfJqH2lGrHYZqDRuhr80M7yEU9PjsI8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lq\nkOEvSQ0y/CWpQYa/JDXIGb7ShHi7Bs0yw1+aKG/ZoNnksI8kNcjwl6QGOewjbWHDzjtUlWNPMvyl\nrc1zDhrOYR9JapDhL0kNMvwlqUGGvyQ1yBO+0hg4m1fzxvCXxsYrazQ/HPaRpAatGf5J7klyKsmR\ngbJ9SU4keah7vHTgub1JHk1yLMnVA+VXJDnSPffW8TdFkrRe6+n5vwPYfUZZAW+pqud3jw8AJNkB\n3Ajs6Pa5M8nyZ9+7gNuqajuwPcmZx5QkTcma4V9VDwJfHfLUsAHNa4H9VXW6qo4DjwE7k1wInFtV\nh7vt3glct7kqS5JGNcqY/6uS/GuSu5Oc15VdBJwY2OYEcPGQ8sWuXJLUg81e7XMX8Efd8h8DbwZu\nG0eFkuwbWF2oqoVxHFeStooku4BdoxxjU+FfVU8MVOLtwHu71UXg0oFNL2Gpx7/YLQ+WL65w7H2b\nqZMktaLrFC8srye5Y6PH2NSwTzeGv+x6YPlKoIPATUnOTnIZsB04XFUnga8l2dmdAL4FuH8zry1J\nGt2aPf8k+4EXAc9K8jhwB7AryY+zdNXPF4DfBKiqo0kOAEeBJ4E9VbU882UPcC9wDvD+qvrgmNsi\nSVqnPJXN/UtSftGE5tHS7R2GzfBdT9lGth39mL7Htp7NZKczfCWpQYa/JDXI8JekBhn+ktQgw1+S\nGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yC9wl1awdMsGaWsy/KVVbeQ+PNL8cNhHkhpk+EtSgwx/\nSWqQY/5SY4adyPYe/+0x/KXmeLJaDvtIUpMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+S\nGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalB\nhr8kNcjwl6QGndV3BST1L0mdWVZV6aMumg7DXxJwZvab+1udwz6S1CDDX5Ia5LCPxPAx79at9H/i\nuYCtwfCXvsVx7283LPtb/z/ZOhz2kaQGGf6S1KA1wz/JPUlOJTkyUHZ+kkNJHknyQJLzBp7bm+TR\nJMeSXD1QfkWSI91zbx1/UyRJ67Wenv87gN1nlN0OHKqqy4GPdOsk2QHcCOzo9rkzyfIg4V3AbVW1\nHdie5MxjSpKmZM3wr6oHga+eUXwNcF+3fB9wXbd8LbC/qk5X1XHgMWBnkguBc6vqcLfdOwf2kSRN\n2Wav9tlWVae65VPAtm75IuDjA9udAC4GTnfLyxa7cklzxltBbA0jX+pZVTXOa6ST7BtYXaiqhXEd\nW9I4eEls35LsAnaNcozNhv+pJBdU1cluSOeJrnwRuHRgu0tY6vEvdsuD5YvDDlxV+zZZJ0lqQtcp\nXlheT3LHRo+x2Us9DwK3dsu3AvcPlN+U5OwklwHbgcNVdRL4WpKd3QngWwb2kSRN2Zo9/yT7gRcB\nz0ryOPA64A3AgSS3AceBGwCq6miSA8BR4ElgT1Utf0bcA9wLnAO8v6o+ON6mSJLWK09lc/+SlCeO\n1Iel81bDxrJXusXBerZdb9n8H9P3bb82k53O8JWkBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kN8pu8\n1By/slEy/NUs70+jtjnsI0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSg7zOX9LI/FL3+WP4\na8tyJu80OWlu3hj+2uJW+jYqqW2O+UtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5Ia\nZPhLUoMMf0lqkOEvSQ0y/CWpQd7YTdJEeJvn2Wb4S5oQb/M8yxz2kaQGGf6S1CDDX5IaZPhLUoMM\nf0lqkOEvSQ0y/CWpQV7nry1h2IQiSSsz/LWFOKlIWi+HfSSpQYa/JDXI8JekBhn+ktQgw1+SGmT4\nS1KDRgr/JMeTfCbJQ0kOd2XnJzmU5JEkDyQ5b2D7vUkeTXIsydWjVl6StDmj9vwL2FVVz6+qK7uy\n24FDVXU58JFunSQ7gBuBHcBu4M4kfvKQpB6MI3zPnElzDXBft3wfcF23fC2wv6pOV9Vx4DHgSqQN\nSlJnPvqukzRvxtHz/3CSTyb59a5sW1Wd6pZPAdu65YuAEwP7ngAuHvH11aw64yFpI0a9vcMLq+rL\nSZ4NHEpybPDJqlqrV/YdzyXZN7C6UFULI9ZR0oxYKQ/8YveNSbIL2DXKMUYK/6r6cvfvV5K8h6Vh\nnFNJLqiqk0kuBJ7oNl8ELh3Y/ZKu7Mxj7hulTpJm2bDsN/c3qusULyyvJ7ljo8fY9LBPku9Ocm63\n/D3A1cAR4CBwa7fZrcD93fJB4KYkZye5DNgOHN7s60uSNm+Unv824D1Jlo/zd1X1QJJPAgeS3AYc\nB24AqKqjSQ4AR4EngT1V5WCtJPUgs5S/ScqxP61ladx42O2b11O2kW0nsb/HHHZM3/ej2Ux2ep29\nJDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JatCot3eQpJENu+2Dl39OluEvaQYMmxOgSTL8Jc0kPw1M\nluEvaUb5aWCSPOErSQ0y/CWpQQ77aKb5FY3SZBj+mgOO/Urj5rCPJDXI8JekBhn+ktQgw1+SGmT4\nS1KDDH9JapCXemomeD2/NF2Gv2bIsPz3mn49xZu9jY/hL2mOOOFvXBzzl6QG2fPX1Dm+L/XP8FdP\n/Pgu9clhH0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDnOSliXI2ryZtpd8xb/i2\nOsNfU+BsXk2Sd4PdDId9JKlBhr8kNchhH0lbkl/8sjrDX2PjyV3NFs81rcbw15j5hpPmgWP+ktQg\nw1+SGuSwj6RmeBL4KYa/pIZ85zmpVv8gGP7aMK/q0dbS5kUKUx3zT7I7ybEkjyZ5zTRfW+NWQx6S\n5sXUwj/J04G/AnYDO4Cbkzx3Wq8/C5Ls6rsOq0lS63n0Xc/+LPRdgQla6LsCE7awoa3X+16Y5/fD\nNHv+VwKPVdXxqjoNvBu4doqvPwt29V2BtQ3rzdvDX7LQdwUmaKHvCkzYwga3X+mT7dZ5L0xzzP9i\n4PGB9RPAzim+fi+SXAg8s1s9P8n2bvk/q+qJEY67rtvYrtYzaeGklqThphn+8/1nctO+6274xksH\nCl61vJCMmr3ru3JhpVvezvNHVmlWrPd9NGudrVRN5/2f5CpgX1Xt7tb3At+sqjcObGMYSdImbPSP\nyzTD/yzg88DPAV8CDgM3V9XDU6mAJOlbpjbsU1VPJnkl8CHg6cDdBr8k9WNqPX9J0uyYiRu7JXl5\nks8l+b8kVwyUvyTJJ5N8pvv3Z/qs52ac0bYXDJSfn+SjSb6e5C/7rOMoVmpf99zebkLfsSRX91XH\ncUnyY0k+1v0+Hkxybt91GqckVyY5nOShJJ9I8pN912mckry7a9tDSb6Q5KG+6zROSV6V5OEkn03y\nxrW2n5XbOxwBrgf+mm+/NOUrwC9W1ckkz2NpyOiSHuo3isG2Dfof4A+BH+4e82po+5LsAG5kaULf\nxcCHk1xeVd+cfhXH5u3A71bVg0l+Ffg94HU912mc3gS8tqo+lOSl3frcdbhWUlU3LS8n+XPgP3qs\nzlh1HeNrgB+tqtNJnr3WPjPR86+qY1X1yJDyT1fVyW71KHBOkmdMt3ajWaVt/11V/wx8o4dqjc1K\n7WNpAt/+qjpdVceBx1ia6DfPtlfVg93yh4GX9VmZCfgy8L3d8nnAYo91mZgsXWN9A7C/77qM0SuA\nP+0m0FJVX1lrh5kI/3V6GfCp5cZtIVv1pMtFLE3kW3aCpU8A8+xzSZZnpb8cuLTPykzA7cCbk3wR\n+DNgb8/1mZSfAk5V1b/1XZEx2g78dJKPJ1lI8hNr7TC1YZ8kh4ALhjz1+1X13jX2fR7wBuAlk6jb\nqEZp2zwYY/tm/g/dam0Ffg14W5LXAgeB/51m3cZhlfb9AfBq4NVV9Z4kLwfuYUbfcytZ5+/qzcC7\npler8VjjZ3cW8H1VdVV3ruYA8JzVjjfNSz039UuU5BLgH4FbquoL463VeGy2bfNik+1b5Nt7xpcw\nB8MI62jrzwMkuRz4hcnXaLxWa1+Sv62qF3erf8/SOY65stbPr5tvdD3wgtW2m0Vr/OxewVJOUlWf\nSPLNJN9fVf++0j6zOOzzrVlqSc4D3ge8pqo+1l+VxmbYDLyZmvI9osG2HARuSnJ2kstY+lh6uJ9q\njcfySbQkT2PpZP1d/dZo7B5L8qJu+WeBYedy5t2LgYer6kt9V2TM7mfpZ7bcMTl7teCHGQn/JNcn\neRy4Cnhfkg90T70S+EHgjoFLtJ7VW0U3YZW2keQ48GbgV5J8MckP9VTNTVupfVV1lKWPnkeBDwB7\nav4nldyc5PPAw8CJqrq35/qM228Ab0ryaeBPuvWt5ka21oneZfcAz0lyhKX2/fJaOzjJS5IaNBM9\nf0nSdBn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ16P8B2JuZ1wcRkAAAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1130cdfd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print np.mean(mcmc.trace('b1')[:])\n",
    "plt.hist(mcmc.trace('b1')[:], bins=50);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}