/file1.py

## file1.py
{
 "metadata": {
  "name": "2013_04_12d_naive_bayes_in_pymc_revisited"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!date"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Sat Apr 13 13:33:11 PDT 2013\r\n"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "    \n",
      "    \n",
      "    From: pymc@googlegroups.com [mailto:pymc@googlegroups.com] On Behalf Of Tommy Engstr\u00f6m\n",
      "    Sent: Wednesday, April 10, 2013 2:40 PM\n",
      "    To: pymc@googlegroups.com\n",
      "    Subject: [pymc] Naive bayes in PyMC\n",
      "    \n",
      "    I'm completely new to PyMC so this is most likely a trivial question but I haven't managed to figure it out. \n",
      "    \n",
      "    I thought I'd start my pymc journey by building a naive bayes classifier using pymc but I'm having trouble understanding how to do that.\n",
      "    \n",
      "    The problem:\n",
      "    \n",
      "    1000 observations of 10 attributes + 1 target attribute, all booleans.\n",
      "    How can fit the network once I have it? I need to be able to set all 11 attributes at once. Or am I just thinking completely wrong here?\n",
      "    \n",
      "    Thankful for any guiding words\n",
      "    \n",
      "    \n",
      "    \n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pymc as mc"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 42
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Naive Bayes in its simplest form is a prediction model for categorial data.\n",
      "\n",
      "I happen to have code from a recent project which simulates data of this form from an interesting distribution:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import random\n",
      "def simulate(n, p, seed):\n",
      "    random.seed(123456+seed)\n",
      "    mc.np.random.seed(123456+seed)\n",
      "    \n",
      "    # make A clusters, beta distributed\n",
      "    A = 5\n",
      "    X_true = mc.rbeta(.5, .5, size=(A,p))\n",
      "    y_true = mc.rbeta(1, 1, size=A)\n",
      "    \n",
      "    \n",
      "    X = zeros((n,p))\n",
      "    p_true = zeros(n)\n",
      "    for i in range(n):\n",
      "        a_i = random.randrange(A)\n",
      "        X[i] = mc.rbernoulli(X_true[a_i])\n",
      "        p_true[i] = y_true[a_i]\n",
      "    \n",
      "    y = mc.rbinomial(1, p_true)\n",
      "    \n",
      "    test = random.sample(range(n), n/4)\n",
      "    train = list(set(range(n)) - set(test))\n",
      "    \n",
      "    X_train = X[train]\n",
      "    y_train = y[train]\n",
      "    \n",
      "    X_test = X[test]\n",
      "    y_test = y[test]\n",
      "\n",
      "    return locals()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 1000\n",
      "p = 10\n",
      "data = simulate(n=n, p=p, seed=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 44
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As mentioned on the PyMC list, it is not necessary to use PyMC for naive bayes prediction, and scikits-learn has a very simple way to handle this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sklearn.naive_bayes"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nb = sklearn.naive_bayes.BernoulliNB(alpha=0, fit_prior=False)\n",
      "nb.fit(data['X_train'], data['y_train'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "BernoulliNB(alpha=0, binarize=0.0, class_prior=None, fit_prior=False)"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y_pred = nb.predict_proba(data['X_test'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# accuracy\n",
      "mean((y_pred[:,1] >= .5) == data['y_test'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "0.66000000000000003"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "But there is nothing wrong with implementing the naive bayes model in PyMC; it could be a good project for getting started.\n",
      "\n",
      "Formulating the naive bayes predictor in PyMC requires writing this machine learning standard in the language of Bayesian statistics, specifically, the \"data likelihood\" of $X_j | y$:\n",
      "$$X_j \\sim Be(\\alpha_j y + \\beta_j (1-y)),$$\n",
      "where $\\alpha_j$ is the probability of $X_j$ being 1 when $y$ is 1, and $\\beta_j$ is the probability of $X_j$ being 1 when $y$ is 0.\n",
      "\n",
      "The tradition in ML is not to think much about the prior distribution of $y$, and so I won't:\n",
      "$$y_i \\sim Be(.5).$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n_test = len(data['X_test'])\n",
      "y = [mc.Bernoulli('y_%d'%i, .5) for i in range(n_test)]\n",
      "\n",
      "alpha = empty(p)\n",
      "beta = empty(p)\n",
      "for j in range(p):\n",
      "    # alpha[j] is Pr[X_j = 1 | y = 1] in training data\n",
      "    alpha[j] = (data['X_train'][:,j] * data['y_train']).sum() / data['y_train'].sum()\n",
      "    \n",
      "    # beta[j] is Pr[X_j = 1 | y = 0] in training data\n",
      "    beta[j] = (data['X_train'][:,j] * (1-data['y_train'])).sum() / (1-data['y_train']).sum()\n",
      "    \n",
      "X = [mc.Bernoulli('X_%d_%d'%(i,j), alpha[j]*y[i]+beta[j]*(1-y[i]), value=data['X_test'][i,j], observed=True) for i in range(n_test) for j in range(p)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 49
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mcmc = mc.MCMC([y, X])\n",
      "%time mcmc.sample(10000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r",
        "[****************100%******************]  10000 of 10000 complete"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " CPU times: user 19min 19s, sys: 168 ms, total: 19min 19s\n",
        "Wall time: 19min 23s\n"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot([y_i.stats()['mean'] for y_i in y], y_pred[:,1], 's', mec='k', color='none')\n",
      "plot([0,1], [0,1], 'k--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 51,
       "text": [
        "[<matplotlib.lines.Line2D at 0x1b076a10>]"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD9CAYAAACoXlzKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH9hJREFUeJzt3X10U+UdB/Bv2mAiqA0vgo60zXTI7BultqxVCnFnG7a8\nKahMlNMCnmkZKtYjA88wwcEG08NwigU23naGOlYVKWB3hjF4cFKOFKRliHO2kHSnoJVCUk1Km2d/\nYGLSlzRNk9zk5vs55x5J7yX5+dh+ffrc5z6PQgghQEREspEgdQFERBRaDHYiIplhsBMRyQyDnYhI\nZhjsREQyw2AnIpIZv8G+YMECjBo1CpmZmb1e8/jjjyM9PR05OTk4duxYyAskIqL+8Rvs8+fPR3V1\nda/n33jjDZw9exYnT57Eli1bMH/+/JAXSERE/eM32AsLCzF06NBez+/fvx/z5s0DAIwfPx4dHR2w\nWq2hrZCIiPpFOZC/bLVakZyc7Hmt1WphtVqh1Wp9rlMoFAP5GCKiuBXM4gADvnna9UN7C3EhBA8h\nYDAYJK8hWg62BduCbdH9uHTpEh599NFuHeT+GFCwa7VaWCwWz+ueeutERBQYk8mErKwsOJ1O1NXV\nBf0+Awr24uJi7Ny5EwBQW1uLxMREjB49eiBvSUQUl5566imUlJRgw4YN2Lp1KzQaTdDv5XeM/YEH\nHsDBgwfx5ZdfIjk5GStXrsTly5cBAI888ghmz56N9957D+np6VCpVNi2bVvQhcQLvV4vdQlRg23x\nHbbFd+K1LaZMmYIVK1YMKNDdFEKIsC/bq1AoEIGPISKSlWCzk0+eEhHJDIOdiChEjEYjdDod1Gp1\nt0Oj0WDixIkoKyvD3//+97DWMaB57EREdIVGo8HFixc9r72nfl933XVITU3F0aNHccstt+CnP/1p\nWGvhGDsRUQi4b3pmZ2cDAMxmMwAgKSkJDocDCoUCY8aM6dc0xmCzkz12IqIgJSQkdAvegwcPArgS\nyklJSfj666+hUCiQl5cHpTIykcsxdiKiAUpKSvIckydP9jk3ePBgDB48OGKhDjDYiYjCSoq1shjs\nRER90Ov10Ol00Gg0UCqVnsM9DHPp0iXYbDYIIdDY2Oj5excvXvQcx48fj1i9vHlKROSHXq/H4cOH\n0d7eDuDKgoaJiYkAgM7OTgDAoEGDPE/ld807hUKBlJQUn8APVLDZyWAnIuqFTqeD1Wr1BLg/GRkZ\nnhkv7puqSUlJaG1tDfrzGexERCGm0+kAAGfOnAEAz43RxsZGOBwOnDt3zu/flyrYOd2RiCgI33zz\njefP7sB3z12XGm+eEhEFwd+2oVJjj52IKECHDx9GR0dHtzH348ePe544jQbssRNR3NNoNNBoNEhI\nSPA5zpw547NLnHtmDOA7P33JkiVRMwwDMNiJiAAAra2tSElJgcvlgsvl8izo5X3zUgjh6a27v+6e\n+hhNOBRDRNSFyWTCwoULMWjQIJw/fx7Z2dlBzUOXCnvsREReli1b5tl7dPDgwSHZqi7SGOxERF6m\nTp2Kuro6FBcXS11K0DgUQ0TkpbCwUOoSBozBTkSEKzNjLl26hIQE34EM99OnsYRLChBR3LHZbFi6\ndClyc3OxcOFCqcvpVbDZyTF2IoorJpMJWVlZcDqdmD17ttTlhAWHYogoLrh76Xv37sWmTZti+uZo\nXxjsRBQXFixYgGuvvRZ1dXUxOYWxPzjGTkRx4ZtvvsHVV18tdRn9wvXYiYhkhjdPiYhwZSz9q6++\nkroMSTHYiUg23DNeKisrpS5FUrx5SkQxL55mvASCPXYiimne89JjfY2XUGGPnYhiWk1NDTZs2MBA\n98JZMUQUVYxGY1Dn5CjY7GSPnYiiTk8BHm+hPhAcYyeimNDQ0IAPPvhA6jJiQp/BXl1djczMTKSl\npWHt2rXdzn/yySf40Y9+hIyMDKSlpeHtt98OS6FEFJ9sNhvKysqwe/duOJ1OqcuJCX6D3el0oqys\nDNXV1Thx4gQqKytx7Ngxn2tWrVqFBQsWoL6+Hm+88QYWL14c1oKJKH54z3gpKyvDj3/8Y6lLigl+\nx9hramqQnp6O0aNHAwDmzJmDffv2Yfz48Z5rkpOTPbt5t7a2IjU1NYzlElGsMxqNMJvNPW4O7b2p\nxYoVK7B9+3bPvHSOsQfOb7BbrVYkJyd7Xmu1WpjNZp9rli9fjoKCArz00ktoa2vDu+++2+N7ef9H\n0ev10Ov1QRdNRLHN/fPvnQvuP5vNZhiNRpw7dw4PPvggjhw5giNHjkS+SAmYzeZuGRsMv8GuUCj6\nfIPy8nI8/PDDePLJJ3H48GE89NBDOHnyZLfr+H9bIgqEXq+P27zo2ulduXJlUO/jN9i1Wi0sFovn\ntcVi8enBA8ChQ4c8H56fnw+Hw4Hz589j5MiRQRVERPLjHdTePVKj0Qij0cjnXELMb7Dn5eWhvr4e\nTU1NGDlyJHbt2oVNmzb5XHPzzTfjwIEDKCkpwalTp9DW1obhw4eHtWgiij3ucPcOefdN0TFjxkhT\nlEz5nRWjVqtRUVGBKVOmYNy4cZg1axZycnJgMBhQVVUFAFi3bh02btyI9PR0zJ49G3/+85+RmJgY\nkeKJKHY1NDRg48aNcDqdWLBggdTlyEqfT54WFRWhqKjI52ve4z5jx47Fhx9+GPrKiEiWnE4nduzY\ngfPnz0OlUsFkMiE7OxvAlVkxnFgxcFxSgIgiqrm5GXfddRfWrVuH9evXx+2N0nBisBNRRG3cuBEq\nlUrqMmSNwU5EEcGeeeRw2V4iCgubzYavv/4ao0aNkrqUmMXNrIkoYoxGIzQaDdRqdbdDo9GgpKQE\nWVlZeP3116UuNS6xx05EAdHpdLBarXC5XD4/zwqFAtdddx2WLFkCp9OJl19+GQkJCXjttde4q9EA\ncaMNIgo7rVaL0tJSAMD27dtRWlqK7du3A7gyL33Pnj0QQqCsrIyhLiEOxRBRSHzxxRcoLi7GiBEj\noFarpS4nrjHYiSgkJkyYwKUBogSDnYhIZjjGTkT90tDQAJfLhdbWVqxZs8azXd2qVaugVCo5DBMF\nOCuGiAKSkpICq9Xq+Vn23q/BPSuGDyGFFmfFEFHYmEwmJCYmorS0FOvWrYNGo5G6JPKDwU5Efv3m\nN7/B5s2bPXuPUvTjUAwR+VVfXw+tVsteugSCzU4GOxFRlOJaMUQ0YOyAyQODnYhgs9lQVlaG5557\nTupSKAQY7ERxzmQyISsrC06nE0888YTU5VAIcFYMUZyy2WxYunQp9u7dyxkvMsNgJ5Ih94NCZrO5\n2zn3ZtEtLS1wOp2oq6vjjBeZ4VAMkUwZjUbo9XqYzWbPodfrPaH/hz/8AVu3bmWoyxCDnShOKZX8\nhV2uGOxEMmez2XD27Fmpy6AIYrATyVhDQwP3Ho1D/F2MSIacTifKysqwe/duvP766ygqKpK6JIog\nLilAJDMmkwmzZs3CTTfdhMGDB3cbS3fPiuESu9GPy/YSEQDAarXitddeYy89jrHHThQj3D3tnvQ0\nX51iH3vsRHHA3wNHRG4MdqIY9e677+LixYtSl0FRiNMdiWKMeyXG0tJSXHPNNVKXQ1GIwU4UQ959\n911kZmaivb0ddXV1+NnPfiZ1SRSFOBRDFCPOnj2L0tJSbN68mTNeyC/OiiGKERMmTMBVV13V4xov\nnBUjT9zzlIhIZsK252l1dTUyMzORlpaGtWvX9njNrl27MH78eGRlZWHu3Ln9LoKIfLlcLqlLoBjm\nt8fudDrxwx/+EIcOHcKoUaNQUFCAzZs3Y/z48Z5rPv74Y/ziF7+AyWTCkCFD8NVXX2HYsGG+H8Ie\nO1FA3LsaDRkyBC+88ILU5ZDEwtJjr6mpQXp6OkaPHg2lUok5c+Zg3759Ptds27YNixcvxpAhQwCg\nW6gTUWDcM16cTid+/etfS10OxTC/s2KsViuSk5M9r7VabbebNKdPn0ZiYiJefPFFCCFgMBgwY8aM\nbu/lveCQXq/n03JE3+Leo+Tm3ulqoPwGu0Kh6PMNXC4XGhsbUVNTA4vFgttvvx0TJ07s1nPnSnJE\nPVu9ejX3HiUA3Tu9K1euDOp9/Aa7VquFxWLxvLZYLD49eABITk7GxIkTkZiYCJ1Oh7S0NHz66afI\nz88PqiCieLN69WokJiZKXQbJiN8x9ry8PNTX16OpqQmXL1/Grl27uj0YMXXqVM+vDl9++SVOnTqF\nm2++OWwFE8kNQ51CzW+wq9VqVFRUYMqUKRg3bhxmzZqFnJwcGAwGVFVVAQDuueceDB8+HOnp6Zg4\ncSLWrFmD66+/PiLFE8USu92Ozz77TOoyKA7wASWiEPB3D8loNMJkMmHhwoWYP38+nn322cgVRjGN\n67ETSayncH/mmWewaNEiVFVVccYLRQxXdyQKE7PZjIqKCjgcDtTV1THUKWLYYycKk6+++grFxcXY\nunWr1KVQnGGPnShMZs2ahTFjxkhdBsUhBjsRkcxwKIZogEwmEz7++GM+XU1Rg8FOFCS73Y6lS5ei\nqqoKf/rTn3DXXXdJXRIRAA7FEAXFZDIhMzPTM+OFoU7RhD12oi76etjoxRdfxAsvvMB56RS1GOxE\nPegp3N1fu+eee1BSUsKVGClqMdiJ+iklJUXqEoj84hg7kR+dnZ1Sl0DUbwx2oh7Y7XYsWrQIixYt\nkroUon5jsBN10dDQ4JnxsnbtWqnLIeo3jrETfcs9L3337t2YNm0aUlJSsH79eqnLIuo3BjvRt/74\nxz/C4XDgzJkznPFCMY0bbRB9SwgR0AbuRJESbHZyjJ3oWwx1kgsGO8Udu92OkydPSl0GUdgw2Cmu\nuNd4efXVV6UuhShsePOUZEGn0+HMmTO9nl++fDlaW1u59yjFBfbYSVYMBgOEEJ4jNTUVALj3KMUV\n9tgpLnDvUYon7LFTXODeoxRPGOxERDLDYCdZOXLkiNQlEEmOY+wU8+x2O1paWgAA77zzDh80orjH\nHjvFNPe89Pvuuw8XLlzwmRHjffjb7o5Ibthjp5i1adMmrFq1ivPSibrgImAUs86dOweVSsWVGEm2\nuAgYxTS9Xg+NRgOFQtHtSEhIgFqt7vZ3Ro0axVAn6gGDnaJGdnY2VCqVz9j45MmTMWnSJP7GR9QP\nHGOnqNbZ2YnPP/8cHR0dUpdCFDPYY6eoZTKZ8NFHH8HlckGpZB+EKFAMdoo6drsdixYtQklJCX7w\ngx9g7NixnJtO1A/sBlFEaTQaXLp0qdvXhRBQqVQAgG3btnlWYrz77rsjXSJRzOtzumN1dTWefvpp\ndHZ2oqSkBL/61a96vO6NN97Afffdh48++gg5OTm+H8LpjvQtjUYDjUaDxsZGn68rlUq4XK4ev08U\nCgWuuuoqOByOCFVJFB3CMt3R6XSirKwM1dXVOHHiBCorK3Hs2LFu19lsNrz44ovIz8/vdwFEAHDN\nNdcgJSWlx6dGXS4XQ52oH/wGe01NDdLT0zF69GgolUrMmTMH+/bt63bdihUrsGzZMs9UNaJA2O12\nHD16VOoyiGTH7xi71WpFcnKy57VWq4XZbPa5pra2Fk1NTSguLsbzzz/f600u77U69Ho99Hp90EVT\n7DOZTFi4cCHuvfde3HbbbVKXQxQVzGZzt4wNht9g72smgsvlQnl5OXbs2OH5Wm89di7CRMCV74+W\nlhaUlJRwjReiLrp2eleuXBnU+/gNdq1WC4vF4nltsVh8evA2mw0nT570FNLc3IwZM2agqqqq2w1U\nokOHDsFms0EIAbvdjmnTpvmc5/IARKHhd4w9Ly8P9fX1aGpqwuXLl7Fr1y4UFRV5ziclJeGLL75A\nQ0MDGhoakJ+fz1CnXiUmJmLv3r0+N0W9j64zZYgoOH577Gq1GhUVFZgyZQpcLhfmzZuHnJwcGAwG\n5ObmYvr06ZGqk2SgoKBA6hKI4gKX7SUiilLBZiefPKV+0+l0aG1t7fHckiVLMGnSJBw+fBjPPPNM\nhCsjIoBrxVCQlixZgtbWVp/juuuuw759+1BSUoLs7GypSySKW+yxU0iYTCb873//w/Dhw1FXV8cZ\nLkQSYo+dBmzHjh0oKSnB8OHDMXPmTIY6kcQY7DRgM2bMQF1dHa6++mqpSyEicCiGQmDo0KFSl0BE\nXthjp35xr7K4fv16zxK87qO3mTJEFFnssVNA7HY7li5diqamJj4hShTl2GOnPplMJmRmZsLhcPgs\n+EZE0Yk99jhlNBqxatUqdHZ29nhepVKhvLwcra2tqKqq4kqMRDGEPfY4plQqkZSUBIPB4FmYy2Aw\nQKVSQa1W49SpU569RxnqRLGDa8XEmb566omJiQCubFW3ZMkSrqNPJCGuFUMBUyqVUCqV6Ojo8AQ4\nAKxZswYA0NHRIWV5RDRADHZCe3s7mpubpS6DiEKEwR7nOjo6UFFRgbFjx0pdChGFCIM9Tgkh0NnZ\niba2NgBATU2Nz3n3g0hEFHs4KyYOuVwuXL58GdnZ2bhw4YLPbBj3LBmHw8Ebp0Qxij32ONTZ2Qkh\nBI4fP851XohkiNMdiYiiVLDZyaEYIiKZYbDLmMlkwrJly6Qug4gijMEuQ3a7HYsWLUJJSQkmTZok\ndTlEFGEM9hhhNBqhVquhUCh6PdRqNcaNG+dZiZFrvBDFJwZ7DFGr1VCpVEhNTcXkyZM90xSTkpKg\nUqlw7bXX4vTp09iwYQO2bt3KvUeJ4hSnO8rI1VdfjbFjx7KXThTn2GOXkYSEBCiV/H81UbxjsMeo\n3pbdJSJi9y7KGI1GbN++vdvG0A6HA+3t7VAoFGhpacHFixclqpCIoh2DPQrpdDro9XqftVr0ej0O\nHjwIIQTsdjuAK0+leWtuboZOp4tgpUQUjTgUEwPsdjva2tqgUqkwd+5czzZ2XQ+HwwGz2Sx1uUQk\nMQZ7DDhw4AA6OjqQm5uLMWPGSF0OEUU5BnsMuPvuuzFz5kzOeCGigDDYiYhkhl3AKGK32/Hf//4X\njY2NWL9+PdavX+9znk+SElEgGOxRwmQyYeHChSgqKkJjY6PU5RBRDOtzKKa6uhqZmZlIS0vD2rVr\nu51//vnnkZ6ejoyMDEyaNAkNDQ1hKVSuvFdi3LBhA1555RWpSyKiGOc32J1OJ8rKylBdXY0TJ06g\nsrISx44d87kmPz8ftbW1qK+vx9y5c1FeXh7WguXkyJEjXImRiELOb7DX1NQgPT0do0ePhlKpxJw5\nc7Bv3z6fawoLC6FSqQAAd9xxB5qamsJXrcyMGDGCKzESUcj5HWO3Wq1ITk72vNZqtX4fgNm0aRNm\nzpzZ47muT1Hq9fp+FSpHN910E2666SapyyCiKGE2m0PykKHfYO/6yLo/O3fuRG1tLQ4ePNjjee9g\nlyP3Gi890el0fCKUiPrUtdO7cuXKoN7H71CMVquFxWLxvLZYLD49eLcDBw5g9erV2LNnDwYNGhRU\nIXKg0+nQ2Njoc5SWluLChQsoKysLardxIqL+8ttjz8vLQ319PZqamjBy5Ejs2rULmzZt8rnm2LFj\nePTRR/GPf/wDI0aMCGux0cTdQ29ubkZHR4dnGd2EhAQkJCTgmmuuQWZmJj7//HO0tLTgd7/7Xb9+\nAyIiCpbfYFer1aioqMCUKVPgcrkwb9485OTkwGAwIC8vD9OmTcPSpUvR1taGe++9FwCQmpqK3bt3\nR6R4qel0Os9KjO6hFvefL1y4gPr6eqjVauTm5nLGCxFFTJ8PKBUVFaGoqMjna97jPv/85z9DX1WM\na2lpwX/+8x/Mnj2bs4SIKOK4VkwYDB06lCsxEpFkuKRACLnXeHE4HACA+vp6AEB2draUZRFRnGGw\nD1BHRwcAeHY8MhqNPuPtgPynehJRdGGwB8npdOLo0aNoa2vDhx9+CCEEnnvuOQDwzIoBwAexiCji\nFCICk6sVCoWs5nC7V2K88847sW7dOi4HQERhEWx2ssfeDzabDUuXLsXevXuxefPmbrOFiIiiAYO9\nH44cOQKn04m6ujr20okoanEohogoSnEopp/cNzV72q2otLQUAGezEFFsiusHlMxmM0pLS30W7dLp\ndHjqqadw+vRpqcsjIgpK3PbYe3PhwgVkZWVh2LBhEEJw4S4iijkM9m+5Z7ycPn0ab731Fmpqahjq\nRBSTZD0UYzQaodFokJCQAIVC4XMcPHgQarUa27dvR21tLbKystDe3o7c3FxOYySimCbrYAcAjUaD\nlJQUTJ48GUIICCFgMBigUqlwww03AABuvPFGvPLKK9iyZQuUSv4SQ0SxLa5TrLm5GQBQUFDg83XO\nhiGiWCa7YPfee7S1tRV2ux3AlRD33hDD3VvnfqREJDeyHIpx7z26ZMkSaLVaaLVa3HDDDbhw4QIe\neughPixFRLImy2DvyuVyoaWlBadPn8aDDz7I2S5EJGuyG4rp6ssvv0RbW5vntffeoyqVCjqdToKq\niIjCR9Y99okTJ2Lo0KHYv38/DAaDz8wYIQQcDgfH14lIdmTdY7/zzjtRX1+PpKQk1NTUSF0OEVFE\nyDLY3Wu+dMVhFyKKBzGzbK/3XqJd6fV6PPbYYxg+fPiAPoOIKJoEm50xFeze/3Rbvnw5Dhw4gI6O\nDhw9ehQJCbK+bUBEcSTY7IzpFDSZTNi4cSM6Ozvx3nvvMdSJiBCjY+zee48WFxdjzJgx3KqOiOhb\nMdnF/eSTT9De3o66ujqMGTNG6nKIiKJKTPbY8/LykJeXJ3UZRERRSfJg1+l0sFqt6Ozs7PWayZMn\nexbwcu9V6q2nrxERxStJg91oNOLMmTMBX0tERH2TbIxdoVBg5cqVPZ4zGAw+r7kaIxFR4CS9eeoO\ncPdqi5MnTwYA7Nu3z+c6rsZIRBS4qJoV494Uwz3ezkAnIuq/qAp2tVoNAJgxY4bElRARxa6oCHb3\nGPoHH3wAAJ6xd46tExH1n+TB7m+4RaVSwWAwyGrNdDn9uwwU2+I7bIvvsC0Grs9gr66uRmZmJtLS\n0rB27dpu551OJ+bMmYPMzEzccccdAU1fNJlMAK70zHvqlXtvhCG3aY78pv0O2+I7bIvvsC0Gzm+w\nO51OlJWVobq6GidOnEBlZSWOHTvmc83LL7+MG2+8EXV1dXj66afx+OOP9/p+NpsNZWVlKCkpwf79\n+312M/I+iIgoeH6DvaamBunp6Rg9ejSUSiXmzJnTbSri/v37MW/ePABXbnr+61//6jGcTSYTsrKy\n4HQ6UVdXh6KiohD+axARkYfwY+fOneLRRx/1vH7ttdfEI4884nPNLbfcIs6dO+d5PXbsWNHc3Oxz\nDQAePHjw4BHEEQy/SwqEah45h1eIiCLH71CMVquFxWLxvLZYLEhOTu52zdmzZwEALpcLLS0tuP76\n68NQKhERBcJvsOfl5aG+vh5NTU24fPkydu3a1W1svLi4GH/9618BAG+//TYKCgq4kxERkYT8DsWo\n1WpUVFRgypQpcLlcmDdvHnJycmAwGJCbm4vp06dj8eLFmDdvHjIzM3Httdfi1VdfjVTtRETUk6BG\n5nvxzjvviIyMDHHrrbeKNWvWdDvvcDjE/fffLzIyMsTtt98uGhsbQ/nxUaWvtvj9738v0tLSRHp6\nuigsLBSff/65BFVGRl9t4VZZWSkUCoU4evRoBKuLrEDa4m9/+5vIzs4WmZmZ4oEHHohwhZHTV1uc\nOnVKTJgwQaSnp4tbb71V7N69W4Iqw2/+/Pli5MiRIiMjo9drHnvsMZGWlibGjx8vamtr+3zPkAW7\nw+EQOp1OWK1WcfnyZZGbm9utgBdeeEE88cQTQggh3nrrLTFjxoxQfXxUCaQt3n//feFwOIQQQlRU\nVIi7775bilLDLpC2EEKIS5cuicLCQlFQUCDbYA+kLY4fPy4mTJgg7Ha7EEKIlpYWKUoNu0Da4sEH\nHxQbN24UQgjx73//W2i1WilKDbv3339f1NbW9hrslZWVYubMmUIIIWpra8W4ceP6fM+QDYaHcs57\nrAukLQoLC6FSqQAAd9xxB5qamqQoNewCaQsAWLFiBZYtWwaVSiXL7wkgsLbYtm0bFi9ejCFDhgAA\nhg0bJkWpYRdIWyQnJ+PixYsAgNbWVqSmpkpRatgVFhZi6NChvZ73zs3x48ejo6MDVqvV73uGLNit\nVqvPjBmtVtvtw72vSUhIwPDhw3H+/PlQlRA1AmkLb5s2bcLMmTMjUVrEBdIWtbW1aGpqQnFxMQD5\nLtccSFucPn0ax48fR25uLm677Tbs2bMn0mVGRCBtsXz5cuzYsQPJycmYOnUqXnrppUiXGRX6mydA\nCLfGk+sPYzD60xY7d+5EbW0tDh48GMaKpNNXW7hcLpSXl2PHjh2er8m1xx7I94XL5UJjYyNqampg\nsVhw++23Y+LEibLruQfSFuXl5Xj44Yfx5JNP4vDhw3jooYdw8uTJCFQXfbr+TPTVfiHrsXPO+3cC\naQsAOHDgAFavXo09e/Zg0KBBkSwxYvpqC5vNhpMnT0Kv1+P73/8+Dh8+jBkzZqC2tlaKcsMqkO+L\n5ORkTJ8+HYmJidDpdEhLS8Onn34a6VLDLpC2OHToEO6//34AQH5+PhwOhyx/w+9L17ayWq3QarX+\n/1KobgB88803IjU1VVitVtHe3i5yc3O73QTzvnn65ptviunTp4fq46NKIG1RW1srbr75ZvHZZ59J\nVGVkBNIW3vR6vWxvngbSFm+++aYoKSkRQgjxxRdfiO9973vi/PnzElQbXoG0RXFxsdi+fbsQ4srN\n01GjRomOjg4pyg27hoYGvzdP3ZMrjh49KrKysvp8v5BOd9y/f79natJvf/tbIYQQzz77rNizZ48Q\n4sqd8Pvuu09kZGSIgoIC0dDQEMqPjyq9tUVVVZUQQoif/OQn4oYbbhDZ2dkiOzvbc9dbjvr6vvAm\n52AXIrC2KC8vF2lpaWLs2LHiL3/5i1Slhl1fbfHJJ5+I/Px8kZaWJm699VbPz47c/PznPxc33nij\nGDRokNBqtWLLli1i48aNnhlBQgjxy1/+0jPdMZCfD4UQMh3QJCKKU3z2n4hIZhjsREQyw2AnIpIZ\nBjsRkcww2ImIZIbBTkQkM/8HIVRTE0kfTMAAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10d0b0d0>"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    }
   ],
   "metadata": {}
  }
 ]
}
	{
	"metadata": {
	"name": "2013_04_12d_naive_bayes_in_pymc_revisited"
	},
	"nbformat": 3,
	"nbformat_minor": 0,
	"worksheets": [
	{
	"cells": [
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"!date"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "stream",
	"stream": "stdout",
	"text": [
	"Sat Apr 13 13:33:11 PDT 2013\r\n"
	]
	}
	],
	"prompt_number": 41
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	" \n",
	" \n",
	" From: pymc@googlegroups.com [mailto:pymc@googlegroups.com] On Behalf Of Tommy Engstr\u00f6m\n",
	" Sent: Wednesday, April 10, 2013 2:40 PM\n",
	" To: pymc@googlegroups.com\n",
	" Subject: [pymc] Naive bayes in PyMC\n",
	" \n",
	" I'm completely new to PyMC so this is most likely a trivial question but I haven't managed to figure it out. \n",
	" \n",
	" I thought I'd start my pymc journey by building a naive bayes classifier using pymc but I'm having trouble understanding how to do that.\n",
	" \n",
	" The problem:\n",
	" \n",
	" 1000 observations of 10 attributes + 1 target attribute, all booleans.\n",
	" How can fit the network once I have it? I need to be able to set all 11 attributes at once. Or am I just thinking completely wrong here?\n",
	" \n",
	" Thankful for any guiding words\n",
	" \n",
	" \n",
	" \n"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"import pymc as mc"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 42
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"Naive Bayes in its simplest form is a prediction model for categorial data.\n",
	"\n",
	"I happen to have code from a recent project which simulates data of this form from an interesting distribution:"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"import random\n",
	"def simulate(n, p, seed):\n",
	" random.seed(123456+seed)\n",
	" mc.np.random.seed(123456+seed)\n",
	" \n",
	" # make A clusters, beta distributed\n",
	" A = 5\n",
	" X_true = mc.rbeta(.5, .5, size=(A,p))\n",
	" y_true = mc.rbeta(1, 1, size=A)\n",
	" \n",
	" \n",
	" X = zeros((n,p))\n",
	" p_true = zeros(n)\n",
	" for i in range(n):\n",
	" a_i = random.randrange(A)\n",
	" X[i] = mc.rbernoulli(X_true[a_i])\n",
	" p_true[i] = y_true[a_i]\n",
	" \n",
	" y = mc.rbinomial(1, p_true)\n",
	" \n",
	" test = random.sample(range(n), n/4)\n",
	" train = list(set(range(n)) - set(test))\n",
	" \n",
	" X_train = X[train]\n",
	" y_train = y[train]\n",
	" \n",
	" X_test = X[test]\n",
	" y_test = y[test]\n",
	"\n",
	" return locals()"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 43
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"n = 1000\n",
	"p = 10\n",
	"data = simulate(n=n, p=p, seed=0)"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 44
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"As mentioned on the PyMC list, it is not necessary to use PyMC for naive bayes prediction, and scikits-learn has a very simple way to handle this:"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"import sklearn.naive_bayes"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 45
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"nb = sklearn.naive_bayes.BernoulliNB(alpha=0, fit_prior=False)\n",
	"nb.fit(data['X_train'], data['y_train'])"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "pyout",
	"prompt_number": 46,
	"text": [
	"BernoulliNB(alpha=0, binarize=0.0, class_prior=None, fit_prior=False)"
	]
	}
	],
	"prompt_number": 46
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"y_pred = nb.predict_proba(data['X_test'])"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 47
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"# accuracy\n",
	"mean((y_pred[:,1] >= .5) == data['y_test'])"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "pyout",
	"prompt_number": 48,
	"text": [
	"0.66000000000000003"
	]
	}
	],
	"prompt_number": 48
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"But there is nothing wrong with implementing the naive bayes model in PyMC; it could be a good project for getting started.\n",
	"\n",
	"Formulating the naive bayes predictor in PyMC requires writing this machine learning standard in the language of Bayesian statistics, specifically, the \"data likelihood\" of $X_j \| y$:\n",
	"$$X_j \\sim Be(\\alpha_j y + \\beta_j (1-y)),$$\n",
	"where $\\alpha_j$ is the probability of $X_j$ being 1 when $y$ is 1, and $\\beta_j$ is the probability of $X_j$ being 1 when $y$ is 0.\n",
	"\n",
	"The tradition in ML is not to think much about the prior distribution of $y$, and so I won't:\n",
	"$$y_i \\sim Be(.5).$$"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"n_test = len(data['X_test'])\n",
	"y = [mc.Bernoulli('y_%d'%i, .5) for i in range(n_test)]\n",
	"\n",
	"alpha = empty(p)\n",
	"beta = empty(p)\n",
	"for j in range(p):\n",
	" # alpha[j] is Pr[X_j = 1 \| y = 1] in training data\n",
	" alpha[j] = (data['X_train'][:,j] * data['y_train']).sum() / data['y_train'].sum()\n",
	" \n",
	" # beta[j] is Pr[X_j = 1 \| y = 0] in training data\n",
	" beta[j] = (data['X_train'][:,j] * (1-data['y_train'])).sum() / (1-data['y_train']).sum()\n",
	" \n",
	"X = [mc.Bernoulli('X_%d_%d'%(i,j), alpha[j]y[i]+beta[j](1-y[i]), value=data['X_test'][i,j], observed=True) for i in range(n_test) for j in range(p)]"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 49
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"mcmc = mc.MCMC([y, X])\n",
	"%time mcmc.sample(10000)"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "stream",
	"stream": "stdout",
	"text": [
	" \r",
	"[**************100%****************] 10000 of 10000 complete"
	]
	},
	{
	"output_type": "stream",
	"stream": "stdout",
	"text": [
	" CPU times: user 19min 19s, sys: 168 ms, total: 19min 19s\n",
	"Wall time: 19min 23s\n"
	]
	}
	],
	"prompt_number": 50
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"plot([y_i.stats()['mean'] for y_i in y], y_pred[:,1], 's', mec='k', color='none')\n",
	"plot([0,1], [0,1], 'k--')"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "pyout",
	"prompt_number": 51,
	"text": [
	"[<matplotlib.lines.Line2D at 0x1b076a10>]"
	]
	},
	{
	"output_type": "display_data",
	"png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD9CAYAAACoXlzKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH9hJREFUeJzt3X10U+UdB/Bv2mAiqA0vgo60zXTI7BultqxVCnFnG7a8\nKahMlNMCnmkZKtYjA88wwcEG08NwigU23naGOlYVKWB3hjF4cFKOFKRliHO2kHSnoJVCUk1Km2d/\nYGLSlzRNk9zk5vs55x5J7yX5+dh+ffrc5z6PQgghQEREspEgdQFERBRaDHYiIplhsBMRyQyDnYhI\nZhjsREQyw2AnIpIZv8G+YMECjBo1CpmZmb1e8/jjjyM9PR05OTk4duxYyAskIqL+8Rvs8+fPR3V1\nda/n33jjDZw9exYnT57Eli1bMH/+/JAXSERE/eM32AsLCzF06NBez+/fvx/z5s0DAIwfPx4dHR2w\nWq2hrZCIiPpFOZC/bLVakZyc7Hmt1WphtVqh1Wp9rlMoFAP5GCKiuBXM4gADvnna9UN7C3EhBA8h\nYDAYJK8hWg62BduCbdH9uHTpEh599NFuHeT+GFCwa7VaWCwWz+ueeutERBQYk8mErKwsOJ1O1NXV\nBf0+Awr24uJi7Ny5EwBQW1uLxMREjB49eiBvSUQUl5566imUlJRgw4YN2Lp1KzQaTdDv5XeM/YEH\nHsDBgwfx5ZdfIjk5GStXrsTly5cBAI888ghmz56N9957D+np6VCpVNi2bVvQhcQLvV4vdQlRg23x\nHbbFd+K1LaZMmYIVK1YMKNDdFEKIsC/bq1AoEIGPISKSlWCzk0+eEhHJDIOdiChEjEYjdDod1Gp1\nt0Oj0WDixIkoKyvD3//+97DWMaB57EREdIVGo8HFixc9r72nfl933XVITU3F0aNHccstt+CnP/1p\nWGvhGDsRUQi4b3pmZ2cDAMxmMwAgKSkJDocDCoUCY8aM6dc0xmCzkz12IqIgJSQkdAvegwcPArgS\nyklJSfj666+hUCiQl5cHpTIykcsxdiKiAUpKSvIckydP9jk3ePBgDB48OGKhDjDYiYjCSoq1shjs\nRER90Ov10Ol00Gg0UCqVnsM9DHPp0iXYbDYIIdDY2Oj5excvXvQcx48fj1i9vHlKROSHXq/H4cOH\n0d7eDuDKgoaJiYkAgM7OTgDAoEGDPE/ld807hUKBlJQUn8APVLDZyWAnIuqFTqeD1Wr1BLg/GRkZ\nnhkv7puqSUlJaG1tDfrzGexERCGm0+kAAGfOnAEAz43RxsZGOBwOnDt3zu/flyrYOd2RiCgI33zz\njefP7sB3z12XGm+eEhEFwd+2oVJjj52IKECHDx9GR0dHtzH348ePe544jQbssRNR3NNoNNBoNEhI\nSPA5zpw547NLnHtmDOA7P33JkiVRMwwDMNiJiAAAra2tSElJgcvlgsvl8izo5X3zUgjh6a27v+6e\n+hhNOBRDRNSFyWTCwoULMWjQIJw/fx7Z2dlBzUOXCnvsREReli1b5tl7dPDgwSHZqi7SGOxERF6m\nTp2Kuro6FBcXS11K0DgUQ0TkpbCwUOoSBozBTkSEKzNjLl26hIQE34EM99OnsYRLChBR3LHZbFi6\ndClyc3OxcOFCqcvpVbDZyTF2IoorJpMJWVlZcDqdmD17ttTlhAWHYogoLrh76Xv37sWmTZti+uZo\nXxjsRBQXFixYgGuvvRZ1dXUxOYWxPzjGTkRx4ZtvvsHVV18tdRn9wvXYiYhkhjdPiYhwZSz9q6++\nkroMSTHYiUg23DNeKisrpS5FUrx5SkQxL55mvASCPXYiimne89JjfY2XUGGPnYhiWk1NDTZs2MBA\n98JZMUQUVYxGY1Dn5CjY7GSPnYiiTk8BHm+hPhAcYyeimNDQ0IAPPvhA6jJiQp/BXl1djczMTKSl\npWHt2rXdzn/yySf40Y9+hIyMDKSlpeHtt98OS6FEFJ9sNhvKysqwe/duOJ1OqcuJCX6D3el0oqys\nDNXV1Thx4gQqKytx7Ngxn2tWrVqFBQsWoL6+Hm+88QYWL14c1oKJKH54z3gpKyvDj3/8Y6lLigl+\nx9hramqQnp6O0aNHAwDmzJmDffv2Yfz48Z5rkpOTPbt5t7a2IjU1NYzlElGsMxqNMJvNPW4O7b2p\nxYoVK7B9+3bPvHSOsQfOb7BbrVYkJyd7Xmu1WpjNZp9rli9fjoKCArz00ktoa2vDu+++2+N7ef9H\n0ev10Ov1QRdNRLHN/fPvnQvuP5vNZhiNRpw7dw4PPvggjhw5giNHjkS+SAmYzeZuGRsMv8GuUCj6\nfIPy8nI8/PDDePLJJ3H48GE89NBDOHnyZLfr+H9bIgqEXq+P27zo2ulduXJlUO/jN9i1Wi0sFovn\ntcVi8enBA8ChQ4c8H56fnw+Hw4Hz589j5MiRQRVERPLjHdTePVKj0Qij0cjnXELMb7Dn5eWhvr4e\nTU1NGDlyJHbt2oVNmzb5XHPzzTfjwIEDKCkpwalTp9DW1obhw4eHtWgiij3ucPcOefdN0TFjxkhT\nlEz5nRWjVqtRUVGBKVOmYNy4cZg1axZycnJgMBhQVVUFAFi3bh02btyI9PR0zJ49G3/+85+RmJgY\nkeKJKHY1NDRg48aNcDqdWLBggdTlyEqfT54WFRWhqKjI52ve4z5jx47Fhx9+GPrKiEiWnE4nduzY\ngfPnz0OlUsFkMiE7OxvAlVkxnFgxcFxSgIgiqrm5GXfddRfWrVuH9evXx+2N0nBisBNRRG3cuBEq\nlUrqMmSNwU5EEcGeeeRw2V4iCgubzYavv/4ao0aNkrqUmMXNrIkoYoxGIzQaDdRqdbdDo9GgpKQE\nWVlZeP3116UuNS6xx05EAdHpdLBarXC5XD4/zwqFAtdddx2WLFkCp9OJl19+GQkJCXjttde4q9EA\ncaMNIgo7rVaL0tJSAMD27dtRWlqK7du3A7gyL33Pnj0QQqCsrIyhLiEOxRBRSHzxxRcoLi7GiBEj\noFarpS4nrjHYiSgkJkyYwKUBogSDnYhIZjjGTkT90tDQAJfLhdbWVqxZs8azXd2qVaugVCo5DBMF\nOCuGiAKSkpICq9Xq+Vn23q/BPSuGDyGFFmfFEFHYmEwmJCYmorS0FOvWrYNGo5G6JPKDwU5Efv3m\nN7/B5s2bPXuPUvTjUAwR+VVfXw+tVsteugSCzU4GOxFRlOJaMUQ0YOyAyQODnYhgs9lQVlaG5557\nTupSKAQY7ERxzmQyISsrC06nE0888YTU5VAIcFYMUZyy2WxYunQp9u7dyxkvMsNgJ5Ih94NCZrO5\n2zn3ZtEtLS1wOp2oq6vjjBeZ4VAMkUwZjUbo9XqYzWbPodfrPaH/hz/8AVu3bmWoyxCDnShOKZX8\nhV2uGOxEMmez2XD27Fmpy6AIYrATyVhDQwP3Ho1D/F2MSIacTifKysqwe/duvP766ygqKpK6JIog\nLilAJDMmkwmzZs3CTTfdhMGDB3cbS3fPiuESu9GPy/YSEQDAarXitddeYy89jrHHThQj3D3tnvQ0\nX51iH3vsRHHA3wNHRG4MdqIY9e677+LixYtSl0FRiNMdiWKMeyXG0tJSXHPNNVKXQ1GIwU4UQ959\n911kZmaivb0ddXV1+NnPfiZ1SRSFOBRDFCPOnj2L0tJSbN68mTNeyC/OiiGKERMmTMBVV13V4xov\nnBUjT9zzlIhIZsK252l1dTUyMzORlpaGtWvX9njNrl27MH78eGRlZWHu3Ln9LoKIfLlcLqlLoBjm\nt8fudDrxwx/+EIcOHcKoUaNQUFCAzZs3Y/z48Z5rPv74Y/ziF7+AyWTCkCFD8NVXX2HYsGG+H8Ie\nO1FA3LsaDRkyBC+88ILU5ZDEwtJjr6mpQXp6OkaPHg2lUok5c+Zg3759Ptds27YNixcvxpAhQwCg\nW6gTUWDcM16cTid+/etfS10OxTC/s2KsViuSk5M9r7VabbebNKdPn0ZiYiJefPFFCCFgMBgwY8aM\nbu/lveCQXq/n03JE3+Leo+Tm3ulqoPwGu0Kh6PMNXC4XGhsbUVNTA4vFgttvvx0TJ07s1nPnSnJE\nPVu9ejX3HiUA3Tu9K1euDOp9/Aa7VquFxWLxvLZYLD49eABITk7GxIkTkZiYCJ1Oh7S0NHz66afI\nz88PqiCieLN69WokJiZKXQbJiN8x9ry8PNTX16OpqQmXL1/Grl27uj0YMXXqVM+vDl9++SVOnTqF\nm2++OWwFE8kNQ51CzW+wq9VqVFRUYMqUKRg3bhxmzZqFnJwcGAwGVFVVAQDuueceDB8+HOnp6Zg4\ncSLWrFmD66+/PiLFE8USu92Ozz77TOoyKA7wASWiEPB3D8loNMJkMmHhwoWYP38+nn322cgVRjGN\n67ETSayncH/mmWewaNEiVFVVccYLRQxXdyQKE7PZjIqKCjgcDtTV1THUKWLYYycKk6+++grFxcXY\nunWr1KVQnGGPnShMZs2ahTFjxkhdBsUhBjsRkcxwKIZogEwmEz7++GM+XU1Rg8FOFCS73Y6lS5ei\nqqoKf/rTn3DXXXdJXRIRAA7FEAXFZDIhMzPTM+OFoU7RhD12oi76etjoxRdfxAsvvMB56RS1GOxE\nPegp3N1fu+eee1BSUsKVGClqMdiJ+iklJUXqEoj84hg7kR+dnZ1Sl0DUbwx2oh7Y7XYsWrQIixYt\nkroUon5jsBN10dDQ4JnxsnbtWqnLIeo3jrETfcs9L3337t2YNm0aUlJSsH79eqnLIuo3BjvRt/74\nxz/C4XDgzJkznPFCMY0bbRB9SwgR0AbuRJESbHZyjJ3oWwx1kgsGO8Udu92OkydPSl0GUdgw2Cmu\nuNd4efXVV6UuhShsePOUZEGn0+HMmTO9nl++fDlaW1u59yjFBfbYSVYMBgOEEJ4jNTUVALj3KMUV\n9tgpLnDvUYon7LFTXODeoxRPGOxERDLDYCdZOXLkiNQlEEmOY+wU8+x2O1paWgAA77zzDh80orjH\nHjvFNPe89Pvuuw8XLlzwmRHjffjb7o5Ibthjp5i1adMmrFq1ivPSibrgImAUs86dOweVSsWVGEm2\nuAgYxTS9Xg+NRgOFQtHtSEhIgFqt7vZ3Ro0axVAn6gGDnaJGdnY2VCqVz9j45MmTMWnSJP7GR9QP\nHGOnqNbZ2YnPP/8cHR0dUpdCFDPYY6eoZTKZ8NFHH8HlckGpZB+EKFAMdoo6drsdixYtQklJCX7w\ngx9g7NixnJtO1A/sBlFEaTQaXLp0qdvXhRBQqVQAgG3btnlWYrz77rsjXSJRzOtzumN1dTWefvpp\ndHZ2oqSkBL/61a96vO6NN97Afffdh48++gg5OTm+H8LpjvQtjUYDjUaDxsZGn68rlUq4XK4ev08U\nCgWuuuoqOByOCFVJFB3CMt3R6XSirKwM1dXVOHHiBCorK3Hs2LFu19lsNrz44ovIz8/vdwFEAHDN\nNdcgJSWlx6dGXS4XQ52oH/wGe01NDdLT0zF69GgolUrMmTMH+/bt63bdihUrsGzZMs9UNaJA2O12\nHD16VOoyiGTH7xi71WpFcnKy57VWq4XZbPa5pra2Fk1NTSguLsbzzz/f600u77U69Ho99Hp90EVT\n7DOZTFi4cCHuvfde3HbbbVKXQxQVzGZzt4wNht9g72smgsvlQnl5OXbs2OH5Wm89di7CRMCV74+W\nlhaUlJRwjReiLrp2eleuXBnU+/gNdq1WC4vF4nltsVh8evA2mw0nT570FNLc3IwZM2agqqqq2w1U\nokOHDsFms0EIAbvdjmnTpvmc5/IARKHhd4w9Ly8P9fX1aGpqwuXLl7Fr1y4UFRV5ziclJeGLL75A\nQ0MDGhoakJ+fz1CnXiUmJmLv3r0+N0W9j64zZYgoOH577Gq1GhUVFZgyZQpcLhfmzZuHnJwcGAwG\n5ObmYvr06ZGqk2SgoKBA6hKI4gKX7SUiilLBZiefPKV+0+l0aG1t7fHckiVLMGnSJBw+fBjPPPNM\nhCsjIoBrxVCQlixZgtbWVp/juuuuw759+1BSUoLs7GypSySKW+yxU0iYTCb873//w/Dhw1FXV8cZ\nLkQSYo+dBmzHjh0oKSnB8OHDMXPmTIY6kcQY7DRgM2bMQF1dHa6++mqpSyEicCiGQmDo0KFSl0BE\nXthjp35xr7K4fv16zxK87qO3mTJEFFnssVNA7HY7li5diqamJj4hShTl2GOnPplMJmRmZsLhcPgs\n+EZE0Yk99jhlNBqxatUqdHZ29nhepVKhvLwcra2tqKqq4kqMRDGEPfY4plQqkZSUBIPB4FmYy2Aw\nQKVSQa1W49SpU569RxnqRLGDa8XEmb566omJiQCubFW3ZMkSrqNPJCGuFUMBUyqVUCqV6Ojo8AQ4\nAKxZswYA0NHRIWV5RDRADHZCe3s7mpubpS6DiEKEwR7nOjo6UFFRgbFjx0pdChGFCIM9Tgkh0NnZ\niba2NgBATU2Nz3n3g0hEFHs4KyYOuVwuXL58GdnZ2bhw4YLPbBj3LBmHw8Ebp0Qxij32ONTZ2Qkh\nBI4fP851XohkiNMdiYiiVLDZyaEYIiKZYbDLmMlkwrJly6Qug4gijMEuQ3a7HYsWLUJJSQkmTZok\ndTlEFGEM9hhhNBqhVquhUCh6PdRqNcaNG+dZiZFrvBDFJwZ7DFGr1VCpVEhNTcXkyZM90xSTkpKg\nUqlw7bXX4vTp09iwYQO2bt3KvUeJ4hSnO8rI1VdfjbFjx7KXThTn2GOXkYSEBCiV/H81UbxjsMeo\n3pbdJSJi9y7KGI1GbN++vdvG0A6HA+3t7VAoFGhpacHFixclqpCIoh2DPQrpdDro9XqftVr0ej0O\nHjwIIQTsdjuAK0+leWtuboZOp4tgpUQUjTgUEwPsdjva2tqgUqkwd+5czzZ2XQ+HwwGz2Sx1uUQk\nMQZ7DDhw4AA6OjqQm5uLMWPGSF0OEUU5BnsMuPvuuzFz5kzOeCGigDDYiYhkhl3AKGK32/Hf//4X\njY2NWL9+PdavX+9znk+SElEgGOxRwmQyYeHChSgqKkJjY6PU5RBRDOtzKKa6uhqZmZlIS0vD2rVr\nu51//vnnkZ6ejoyMDEyaNAkNDQ1hKVSuvFdi3LBhA1555RWpSyKiGOc32J1OJ8rKylBdXY0TJ06g\nsrISx44d87kmPz8ftbW1qK+vx9y5c1FeXh7WguXkyJEjXImRiELOb7DX1NQgPT0do0ePhlKpxJw5\nc7Bv3z6fawoLC6FSqQAAd9xxB5qamsJXrcyMGDGCKzESUcj5HWO3Wq1ITk72vNZqtX4fgNm0aRNm\nzpzZ47muT1Hq9fp+FSpHN910E2666SapyyCiKGE2m0PykKHfYO/6yLo/O3fuRG1tLQ4ePNjjee9g\nlyP3Gi890el0fCKUiPrUtdO7cuXKoN7H71CMVquFxWLxvLZYLD49eLcDBw5g9erV2LNnDwYNGhRU\nIXKg0+nQ2Njoc5SWluLChQsoKysLardxIqL+8ttjz8vLQ319PZqamjBy5Ejs2rULmzZt8rnm2LFj\nePTRR/GPf/wDI0aMCGux0cTdQ29ubkZHR4dnGd2EhAQkJCTgmmuuQWZmJj7//HO0tLTgd7/7Xb9+\nAyIiCpbfYFer1aioqMCUKVPgcrkwb9485OTkwGAwIC8vD9OmTcPSpUvR1taGe++9FwCQmpqK3bt3\nR6R4qel0Os9KjO6hFvefL1y4gPr6eqjVauTm5nLGCxFFTJ8PKBUVFaGoqMjna97jPv/85z9DX1WM\na2lpwX/+8x/Mnj2bs4SIKOK4VkwYDB06lCsxEpFkuKRACLnXeHE4HACA+vp6AEB2draUZRFRnGGw\nD1BHRwcAeHY8MhqNPuPtgPynehJRdGGwB8npdOLo0aNoa2vDhx9+CCEEnnvuOQDwzIoBwAexiCji\nFCICk6sVCoWs5nC7V2K88847sW7dOi4HQERhEWx2ssfeDzabDUuXLsXevXuxefPmbrOFiIiiAYO9\nH44cOQKn04m6ujr20okoanEohogoSnEopp/cNzV72q2otLQUAGezEFFsiusHlMxmM0pLS30W7dLp\ndHjqqadw+vRpqcsjIgpK3PbYe3PhwgVkZWVh2LBhEEJw4S4iijkM9m+5Z7ycPn0ab731Fmpqahjq\nRBSTZD0UYzQaodFokJCQAIVC4XMcPHgQarUa27dvR21tLbKystDe3o7c3FxOYySimCbrYAcAjUaD\nlJQUTJ48GUIICCFgMBigUqlwww03AABuvPFGvPLKK9iyZQuUSv4SQ0SxLa5TrLm5GQBQUFDg83XO\nhiGiWCa7YPfee7S1tRV2ux3AlRD33hDD3VvnfqREJDeyHIpx7z26ZMkSaLVaaLVa3HDDDbhw4QIe\neughPixFRLImy2DvyuVyoaWlBadPn8aDDz7I2S5EJGuyG4rp6ssvv0RbW5vntffeoyqVCjqdToKq\niIjCR9Y99okTJ2Lo0KHYv38/DAaDz8wYIQQcDgfH14lIdmTdY7/zzjtRX1+PpKQk1NTUSF0OEVFE\nyDLY3Wu+dMVhFyKKBzGzbK/3XqJd6fV6PPbYYxg+fPiAPoOIKJoEm50xFeze/3Rbvnw5Dhw4gI6O\nDhw9ehQJCbK+bUBEcSTY7IzpFDSZTNi4cSM6Ozvx3nvvMdSJiBCjY+zee48WFxdjzJgx3KqOiOhb\nMdnF/eSTT9De3o66ujqMGTNG6nKIiKJKTPbY8/LykJeXJ3UZRERRSfJg1+l0sFqt6Ozs7PWayZMn\nexbwcu9V6q2nrxERxStJg91oNOLMmTMBX0tERH2TbIxdoVBg5cqVPZ4zGAw+r7kaIxFR4CS9eeoO\ncPdqi5MnTwYA7Nu3z+c6rsZIRBS4qJoV494Uwz3ezkAnIuq/qAp2tVoNAJgxY4bElRARxa6oCHb3\nGPoHH3wAAJ6xd46tExH1n+TB7m+4RaVSwWAwyGrNdDn9uwwU2+I7bIvvsC0Grs9gr66uRmZmJtLS\n0rB27dpu551OJ+bMmYPMzEzccccdAU1fNJlMAK70zHvqlXtvhCG3aY78pv0O2+I7bIvvsC0Gzm+w\nO51OlJWVobq6GidOnEBlZSWOHTvmc83LL7+MG2+8EXV1dXj66afx+OOP9/p+NpsNZWVlKCkpwf79\n+312M/I+iIgoeH6DvaamBunp6Rg9ejSUSiXmzJnTbSri/v37MW/ePABXbnr+61//6jGcTSYTsrKy\n4HQ6UVdXh6KiohD+axARkYfwY+fOneLRRx/1vH7ttdfEI4884nPNLbfcIs6dO+d5PXbsWNHc3Oxz\nDQAePHjw4BHEEQy/SwqEah45h1eIiCLH71CMVquFxWLxvLZYLEhOTu52zdmzZwEALpcLLS0tuP76\n68NQKhERBcJvsOfl5aG+vh5NTU24fPkydu3a1W1svLi4GH/9618BAG+//TYKCgq4kxERkYT8DsWo\n1WpUVFRgypQpcLlcmDdvHnJycmAwGJCbm4vp06dj8eLFmDdvHjIzM3Httdfi1VdfjVTtRETUk6BG\n5nvxzjvviIyMDHHrrbeKNWvWdDvvcDjE/fffLzIyMsTtt98uGhsbQ/nxUaWvtvj9738v0tLSRHp6\nuigsLBSff/65BFVGRl9t4VZZWSkUCoU4evRoBKuLrEDa4m9/+5vIzs4WmZmZ4oEHHohwhZHTV1uc\nOnVKTJgwQaSnp4tbb71V7N69W4Iqw2/+/Pli5MiRIiMjo9drHnvsMZGWlibGjx8vamtr+3zPkAW7\nw+EQOp1OWK1WcfnyZZGbm9utgBdeeEE88cQTQggh3nrrLTFjxoxQfXxUCaQt3n//feFwOIQQQlRU\nVIi7775bilLDLpC2EEKIS5cuicLCQlFQUCDbYA+kLY4fPy4mTJgg7Ha7EEKIlpYWKUoNu0Da4sEH\nHxQbN24UQgjx73//W2i1WilKDbv3339f1NbW9hrslZWVYubMmUIIIWpra8W4ceP6fM+QDYaHcs57\nrAukLQoLC6FSqQAAd9xxB5qamqQoNewCaQsAWLFiBZYtWwaVSiXL7wkgsLbYtm0bFi9ejCFDhgAA\nhg0bJkWpYRdIWyQnJ+PixYsAgNbWVqSmpkpRatgVFhZi6NChvZ73zs3x48ejo6MDVqvV73uGLNit\nVqvPjBmtVtvtw72vSUhIwPDhw3H+/PlQlRA1AmkLb5s2bcLMmTMjUVrEBdIWtbW1aGpqQnFxMQD5\nLtccSFucPn0ax48fR25uLm677Tbs2bMn0mVGRCBtsXz5cuzYsQPJycmYOnUqXnrppUiXGRX6mydA\nCLfGk+sPYzD60xY7d+5EbW0tDh48GMaKpNNXW7hcLpSXl2PHjh2er8m1xx7I94XL5UJjYyNqampg\nsVhw++23Y+LEibLruQfSFuXl5Xj44Yfx5JNP4vDhw3jooYdw8uTJCFQXfbr+TPTVfiHrsXPO+3cC\naQsAOHDgAFavXo09e/Zg0KBBkSwxYvpqC5vNhpMnT0Kv1+P73/8+Dh8+jBkzZqC2tlaKcsMqkO+L\n5ORkTJ8+HYmJidDpdEhLS8Onn34a6VLDLpC2OHToEO6//34AQH5+PhwOhyx/w+9L17ayWq3QarX+\n/1KobgB88803IjU1VVitVtHe3i5yc3O73QTzvnn65ptviunTp4fq46NKIG1RW1srbr75ZvHZZ59J\nVGVkBNIW3vR6vWxvngbSFm+++aYoKSkRQgjxxRdfiO9973vi/PnzElQbXoG0RXFxsdi+fbsQ4srN\n01GjRomOjg4pyg27hoYGvzdP3ZMrjh49KrKysvp8v5BOd9y/f79natJvf/tbIYQQzz77rNizZ48Q\n4sqd8Pvuu09kZGSIgoIC0dDQEMqPjyq9tUVVVZUQQoif/OQn4oYbbhDZ2dkiOzvbc9dbjvr6vvAm\n52AXIrC2KC8vF2lpaWLs2LHiL3/5i1Slhl1fbfHJJ5+I/Px8kZaWJm699VbPz47c/PznPxc33nij\nGDRokNBqtWLLli1i48aNnhlBQgjxy1/+0jPdMZCfD4UQMh3QJCKKU3z2n4hIZhjsREQyw2AnIpIZ\nBjsRkcww2ImIZIbBTkQkM/8HIVRTE0kfTMAAAAAASUVORK5CYII=\n",
	"text": [
	"<matplotlib.figure.Figure at 0x10d0b0d0>"
	]
	}
	],
	"prompt_number": 51
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 51
	}
	],
	"metadata": {}
	}
	]
	}