PyMC + Plotly.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Using Plotly to Display MCMC Output"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[Plotly](https://plot.ly) is a collaborative data analysis and graphing platform. Plotly's [Scientific Graphing Libraries](https://plot.ly/api) interface Plotly's online graphing tools with the following scientific computing languages:\n",
      "\n",
      "* Python\n",
      "* R\n",
      "* Matlab\n",
      "* Julia\n",
      "\n",
      "You can think of Plotly as \"Graphics as a Service\". It generates interactive, publication-quality plots that can be embedded in the locaiton of your choice. You can style them locally with code or via the online interface; plots can be shared publicly or privately with a url, and your graphs are accessible from anywhere.\n",
      "\n",
      "You can install Plotly on Python via pip:\n",
      "\n",
      "    pip install plotly"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import plotly"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In order to use Plotly, you need an account. You can sign up using the API, without visiting the website:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "reg = plotly.signup(\"_pymc_demo\", \"_pymc_demo@vanderbilt.edu\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Thanks for signing up to plotly!\n",
        "\n",
        "Your username is: _pymc_demo\n",
        "\n",
        "Your temporary password is: p143f. You use this to log into your plotly account at https://plot.ly/plot.\n",
        "\n",
        "Your API key is: vqg7cjezgb. You use this to access your plotly account through the API.\n",
        "\n",
        "To get started, initialize a plotly object with your username and api_key, e.g. \n",
        ">>> py = plotly.plotly('_pymc_demo', 'vqg7cjezgb')\n",
        "Then, make a graph!\n",
        ">>> res = py.plot([1,2,3],[4,2,1])\n",
        "\n",
        ">>> print(res['url'])\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ply = plotly.plotly(username=reg['un'], key=reg['api_key'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Consider the following dataset, which is a time series of recorded coal mining \n",
      "disasters in the UK from 1851 to 1962:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc.examples import disaster_model\n",
      "import numpy as np\n",
      "\n",
      "data = {'y': disaster_model.disasters_array,\n",
      "        'x': np.arange(1851, 1963),\n",
      "        \"type\":\"scatter\",\"mode\":\"markers\"}\n",
      "layout = { \n",
      "'xaxis':{'showgrid':False,'zeroline':False, 'title': 'Year'},\n",
      "'yaxis':{'showgrid':False,'zeroline':False, 'title': 'Disaster count'} }\n",
      "\n",
      "ply.iplot(data, layout=layout, width=850, height=400)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "High five! You successfuly sent some data to your account on plotly. View your plot in your browser at https://plot.ly/~_pymc_demo/0 or inside your plot.ly account where it is named 'plot from API'\n"
       ]
      },
      {
       "html": [
        "<iframe height=\"450\" id=\"igraph\" scrolling=\"no\" seamless=\"seamless\" src=\"https://plot.ly/~_pymc_demo/0/850/400\" width=\"900\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<IPython.core.display.HTML at 0x111c36a50>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Occurrences of disasters in the time series is thought to be derived from a \n",
      "Poisson process with a large rate parameter in the early part of the time \n",
      "series, and from one with a smaller rate in the later part. We are interested \n",
      "in locating the change point in the series, which perhaps is related to changes \n",
      "in mining safety regulations.\n",
      "\n",
      "We represent our conceptual model formally as a statistical model:\n",
      "\n",
      "$$\n",
      "\\begin{aligned}  (D_t | s, e, l) \\sim\\text{Poisson}\\left(r_t\\right), &\\text{where } r_t=\\left\\{\\begin{array}{lll}             e &\\text{if}& t< s\\\\ l &\\text{if}& t\\ge s             \\end{array}\\right.,&t\\in[t_l,t_h]\\\\         s\\sim \\text{Discrete Uniform}(t_l, t_h)\\\\         e\\sim \\text{Exponential}(r_e)\\\\         l\\sim \\text{Exponential}(r_l)     \\end{aligned}\n",
      "$$\n",
      "\n",
      "The symbols are defined as:\n",
      "    \n",
      "* $D_t$: The number of disasters in year $t$.\n",
      "* $r_t$: The rate parameter of the Poisson distribution of disasters in year $t$.\n",
      "* $s$: The year in which the rate parameter changes (the switchpoint).\n",
      "* $e$: The rate parameter before the switchpoint $s$.\n",
      "* $l$: The rate parameter after the switchpoint $s$.\n",
      "* $t_l$, $t_h$: The lower and upper boundaries of year $t$.\n",
      "* $r_e$, $r_l$: The rate parameters of the priors of the early and late rates, respectively.\n",
      "\n",
      "We can fit this model using [PyMC](http://github.com/pymc-devs/pymc)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc import MCMC, graph\n",
      "\n",
      "M = MCMC(disaster_model)\n",
      "\n",
      "graph.dag(M, format='png')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<pydot.Dot at 0x10ee4bd10>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "![DAG](MCMC.png)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nchains = 3\n",
      "for i in range(nchains):\n",
      "    M.sample(5000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---------        24%                  ] 1204 of 5000 complete in 0.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------54%                  ] 2708 of 5000 complete in 1.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------82%-----------       ] 4118 of 5000 complete in 1.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------100%-----------------] 5000 of 5000 complete in 1.7 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-------------    35%                  ] 1792 of 5000 complete in 0.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------65%----              ] 3273 of 5000 complete in 1.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------98%----------------- ] 4921 of 5000 complete in 1.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------100%-----------------] 5000 of 5000 complete in 1.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [----------       26%                  ] 1333 of 5000 complete in 0.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------52%                  ] 2613 of 5000 complete in 1.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------87%-------------     ] 4358 of 5000 complete in 1.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------100%-----------------] 5000 of 5000 complete in 1.7 sec"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The MCMC sampler generates a *trace* for each unknown parameter, which can be visualized either as a time series or a histogram:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "trace = M.early_mean.trace()\n",
      "\n",
      "data = [{'y': trace}]\n",
      "data.append({'y': trace,\n",
      "             'xaxis': 'x2',\n",
      "             'yaxis': 'y2',\n",
      "             'type': 'histogramy'})\n",
      "\n",
      "layout = {\n",
      "    \"xaxis\":{\n",
      "        \"domain\":[0,0.5],\n",
      "        \"title\": \"Iteration\"\n",
      "    },\n",
      "    \"yaxis\":{\n",
      "        \"title\": \"Value\"\n",
      "    },\n",
      "    \"xaxis2\":{\n",
      "        \"domain\":[0.55,1],\n",
      "        \"title\": \"Value\"\n",
      "    },\n",
      "    \"yaxis2\":{\n",
      "        \"anchor\":\"x2\",\n",
      "        \"side\": \"right\",\n",
      "        \"title\": \"Frequency\"\n",
      "    },\n",
      "    \"showlegend\": False,\n",
      "    \"title\": \"Posterior samples of early Poisson mean\"\n",
      "}\n",
      "\n",
      "ply.iplot(data, layout=layout, width=850,height=400)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "\n"
       ]
      },
      {
       "html": [
        "<iframe height=\"450\" id=\"igraph\" scrolling=\"no\" seamless=\"seamless\" src=\"https://plot.ly/~_pymc_demo/1/850/400\" width=\"900\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<IPython.core.display.HTML at 0x111cf9050>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We ran the model three times. One easy way to check for lack of convergence is to examine the samples from each chain together. If they don't look similar, then they have not yet converged."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nparams = len(M.stochastics)\n",
      "mcmc_data = {'chain {0}'.format(i): {s.__name__:s.trace(chain=i)[-100:] for s in M.stochastics} for i in range(nchains)}\n",
      "\n",
      "attr = mcmc_data['chain 1'].keys()\n",
      "colors = {'chain 0': 'rgb(31, 119, 180)', \n",
      "          'chain 1': 'rgb(255, 127, 14)',\n",
      "          'chain 2': 'rgb(44, 160, 44)'}\n",
      "\n",
      "data = []\n",
      "for i in range(nparams):\n",
      "    for j in range(nparams):\n",
      "        for chain in mcmc_data.keys():\n",
      "            x = mcmc_data[chain][attr[i]]\n",
      "            y = mcmc_data[chain][attr[j]]\n",
      "            if i==j:\n",
      "                data.append({\"name\": chain, \n",
      "                     \"x\": x + np.random.randn(len(x))*np.mean(x)/100., \n",
      "                     \"y\": y + np.random.randn(len(y))*np.mean(y)/100.,\n",
      "                     \"type\":\"scatter\",\"mode\":\"markers\",\n",
      "                     'marker': {'color': colors[chain], 'opacity':0.2},\n",
      "                     \"xaxis\": \"x\"+(str(i) if i!=0 else ''), \n",
      "                     \"yaxis\": \"y\"+(str(j) if j!=0 else '')})\n",
      "            else:\n",
      "                data.append({\"name\": chain, \n",
      "                     \"x\": x, \"y\": y,\n",
      "                     \"type\":\"scatter\",\"mode\":\"markers\",\n",
      "                     'marker': {'color': colors[chain], 'opacity':0.2},\n",
      "                     \"xaxis\": \"x\"+(str(i) if i!=0 else ''), \n",
      "                     \"yaxis\": \"y\"+(str(j) if j!=0 else '')})\n",
      "padding = 0.04;\n",
      "domains = [[i*padding + i*(1-3*padding)/nparams, i*padding + ((i+1)*(1-3*padding)/nparams)] for i in range(nparams)]\n",
      "\n",
      "layout = {\n",
      "    \"xaxis\":{\"domain\":domains[0], \"title\":attr[0], \n",
      "             'zeroline':False,'showline':False,'ticks':'', \n",
      "             'titlefont':{'color': \"rgb(67, 67, 67)\"},'tickfont':{'color': 'rgb(102,102,102)'}},\n",
      "    \"yaxis\":{\"domain\":domains[0], \"title\":attr[0], \n",
      "             'zeroline':False,'showline':False,'ticks':'', \n",
      "             'titlefont':{'color': \"rgb(67, 67, 67)\"},'tickfont':{'color': 'rgb(102,102,102)'}},\n",
      "    \"xaxis1\":{\"domain\":domains[1], \"title\":attr[1], \n",
      "              'zeroline':False,'showline':False,'ticks':'', \n",
      "              'titlefont':{'color': \"rgb(67, 67, 67)\"},'tickfont':{'color': 'rgb(102,102,102)'}},\n",
      "    \"yaxis1\":{\"domain\":domains[1], \"title\":attr[1], \n",
      "              'zeroline':False,'showline':False,'ticks':'', \n",
      "              'titlefont':{'color': \"rgb(67, 67, 67)\"},'tickfont':{'color': 'rgb(102,102,102)'}},\n",
      "    \"xaxis2\":{\"domain\":domains[2], \"title\":attr[2], \n",
      "              'zeroline':False,'showline':False,'ticks':'', \n",
      "              'titlefont':{'color': \"rgb(67, 67, 67)\"},'tickfont':{'color': 'rgb(102,102,102)'}},\n",
      "    \"yaxis2\":{\"domain\":domains[2], \"title\":attr[2], \n",
      "              'zeroline':False,'showline':False,'ticks':'', \n",
      "              'titlefont':{'color': \"rgb(67, 67, 67)\"},'tickfont':{'color': 'rgb(102,102,102)'}},\n",
      "    \n",
      "    \"showlegend\":False,\n",
      "    \"title\":\"Posterior samples for coal mining disasters model\",\n",
      "    \"titlefont\":{'color':'rgb(67,67,67)', 'size': 20}\n",
      "    }\n",
      "\n",
      "ply.iplot(data,layout=layout)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "\n"
       ]
      },
      {
       "html": [
        "<iframe height=\"650\" id=\"igraph\" scrolling=\"no\" seamless=\"seamless\" src=\"https://plot.ly/~_pymc_demo/2/600/600\" width=\"650\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<IPython.core.display.HTML at 0x10edccad0>"
       ]
      }
     ],
     "prompt_number": 8
    }
   ],
   "metadata": {}
  }
 ]
}