pymc_and_pandas.ipynb

{
 "metadata": {
  "name": "pymc_and_pandas"
 },
 "nbformat": 3,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "source": [
      "PyMC Pandas Example",
      "===================",
      "",
      "This example project shows how to fit a fixed effects Poisson model with PyMC.  It uses pandas Series and DataFrame objects to store data in a classy way."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import pylab as pl",
      "import pymc as mc",
      "import pandas"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "source": [
      "1. Simulate Noisy Data",
      "----------------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# simulate data with known distribution",
      "",
      "N = 100",
      "X = pandas.DataFrame({'constant': pl.ones(N), 'cov_1': pl.randn(N)})",
      "",
      "beta_true = pandas.Series(dict(constant=100., cov_1=20.))",
      "mu_true = pl.dot(X, beta_true)",
      "",
      "Y = mc.rpoisson(mu_true)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# explore the data a little bit graphically",
      "",
      "pl.figure(figsize=(11,4.25))",
      "",
      "pl.subplot(1,2,1)",
      "pl.hist(Y)",
      "pl.xlabel('Observed Count')",
      "",
      "pl.subplot(1,2,2)",
      "pl.plot(X['cov_1'], Y, '.')",
      "pl.xlabel('Covariate 1')",
      "pl.ylabel('Observed Count')"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "source": [
      "2. Model data with PyMC",
      "-----------------------",
      "",
      "The following code creates a fixed effect Poisson model ",
      "where the observed data stored in Y is explained by the",
      "covariate data in X, according to the formula:",
      "",
      "$$",
      "Y_i \\sim \\text{Poisson}(\\mu_i),",
      "$$",
      "$$",
      "\\mu_i = X_i\\cdot \\beta.",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# the simplest approach doesn't work with PyMC 2.1alpha, but it does with 2.2grad",
      "print 'pymc version:', mc.__version__",
      "",
      "beta = mc.Uninformative('beta', value=[Y.mean(), 0.])",
      "mu_pred = mc.Lambda('mu_pred', lambda beta=beta, X=X: pl.dot(X, beta))",
      "Y_obs = mc.Poisson('Y_obs', mu=mu_pred, value=Y, observed=True)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "m = mc.Model([beta, mu_pred, Y_obs])",
      "%time mc.MCMC(m).sample(10000, 5000, 5, progress_bar=False)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "mc.Matplot.plot(beta, common_scale=False)",
      "print '\\ntrue value of beta\\n', beta_true",
      "print '\\npredicted:'",
      "print pandas.DataFrame({'mean':beta.stats()['mean'],",
      "                        'lb':beta.stats()['95% HPD interval'][:,0],",
      "                        'ub':beta.stats()['95% HPD interval'][:,1]},",
      "                       columns=['mean','lb','ub'])"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "source": [
      "2a. TODO: Integrate PyMC and Pandas further",
      "---------------------------------------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# making beta.value a pandas.Series would be slightly cooler than the above",
      "",
      "@mc.stochastic",
      "def beta(value=pandas.Series(dict(constant=Y.mean(), cov_1=0))):",
      "    return 0.",
      "mu_pred = mc.Lambda('mu_pred', lambda beta=beta, X=X: pl.dot(X, beta))",
      "Y_obs = mc.Poisson('Y_obs', mu=mu_pred, value=Y, observed=True)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "beta.value"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# unfortunately the pandas.Series becomes a numpy.array during MCMC",
      "m = mc.Model([beta, mu_pred, Y_obs])",
      "mc.MCMC(m).sample(10000, 5000, 5, progress_bar=False)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "beta.value # in a pandas-centric version of PyMC, this would still be a pandas.Series"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [],
     "language": "python",
     "outputs": [],
     "prompt_number": "&nbsp;"
    }
   ]
  }
 ]
}