fonnesbeck/pymc_bug.ipynb

## pymc_bug.ipynb
{
 "metadata": {
  "name": "",
  "signature": "sha256:c748160671b52b3e183301ca9b5e70ecbd31ce02e611fca3412fbdd61d17dc3a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "from sklearn.datasets import load_boston\n",
      "import numpy as np\n",
      "import pymc as pm\n",
      "import pandas as pd\n",
      "\n",
      "boston = load_boston()\n",
      "features = ['INDUS', 'NOX', 'RM', 'TAX', 'PTRATIO', 'LSTAT']\n",
      "df = pd.DataFrame(boston.data, columns=boston.feature_names)\n",
      "X = np.array(df.ix[:, features])\n",
      "y = boston.target\n",
      "\n",
      "gamma = pm.Binomial('gamma', 1, 0.5, size=len(features))\n",
      "var = pm.Lambda('var', lambda gamma=gamma: (1-gamma)*0.001 + gamma*10)\n",
      "prec = pm.Lambda('prec', lambda var=var: 1.0/var)\n",
      "b = pm.Normal('b', 0, prec)\n",
      "int_ = pm.Normal('int_', 0, 0.01)\n",
      "taue = pm.Gamma('taue', 0.1, 0.1)\n",
      "mu = int_ + X[:,0]*b[0] + X[:,1]*b[1] + X[:,2]*b[2] + X[:,3]*b[3] + X[:,4]*b[4] + X[:,5]*b[5]\n",
      "observed = pm.Normal('obs', mu, taue, observed=True, value=y)\n",
      "M = pm.MCMC([observed, mu, int_, b, prec, taue, var, gamma])\n",
      "M.sample(10000, 500, 5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " [-----------------100%-----------------] 10000 of 10000 complete in 3.5 sec"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "int_.stats()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "{'95% HPD interval': array([-82.23023941,  34.05807169]),\n",
        " 'mc error': 4.0775156170050106,\n",
        " 'mean': -9.3811298932046085,\n",
        " 'n': 1900,\n",
        " 'quantiles': {2.5: -82.182744891019681,\n",
        "  25: -43.780936763079715,\n",
        "  50: 7.0756300003831951,\n",
        "  75: 24.127800859958935,\n",
        "  97.5: 34.250495843107721},\n",
        " 'standard deviation': 40.788646352591847}"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pm.Matplot.plot(int_)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Plotting int_\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFxCAYAAACiM6b4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYJGV5+P3vvauAoMKCB4jkILIYEk2MpI0nwPQuSYyg\nSYzR14gYIhExAQ+0GNQso8ZTawhRiP4MxlOihvf1kCjKLvQGVkRp/UXRGDlsogKKBlgQ3EVg937/\nqBroHbpndnpq+vj9XFdfXfV0ddVdPd3V91Td/TyRmUiSJGn5rRh2AJIkSdPCxEuSJGlATLwkSZIG\nxMRLkiRpQEy8JEmSBsTES5IkaUBMvCRJkgZkWRKviNg9Ir4bES8r59dGxKbyVl+ObUrSroqID0TE\nZRGxMSJeWLZ1PU55/JJUpfst03pPBL4KZEQEMAOsLR+7ICI2pj23ShqeBJ6bmd8DiIgVzDlOAa1u\n7R6/JC1F5We8ImJP4Cjg00AAq4GrMnNbZm4DNgMHV71dSVqk6Ji+z3EqIlZ3a8fjl6QlWI4zXicD\n7wYeXs7vB9wSEWeW87eWbVcvw7YlaVfcBvxzRNwMvALYl+7HqejR7vFLUl8qTbwiYm/gqZn51oh4\nUdl8E7APcBLFQewc4MZe67jooos8hS9NmTVr1sTCS1UnM08GiIjHAU3gNLofp1b0aO/K45c0nRZz\nDKv6jNdTgD0i4qPAI8v1bwIO6VhmdWZeM99KBn0QljQ8Q05W7gDuAq6hy3EqIlZ2a59vhR6/NM0i\nIgEyc2o+B4s9hlWaeGXm+cD5ABFxHLBXZl4RETPAhnKxM6rcpiQtVkR8DDgAuB04KTN3dDtOZeb2\ncTl+zczMJMC6deum5gtPWg7L/Vlarl81kpkf7JheD6xfrm1J0mJk5vO6tHU9To3L8cuES6rGcn+W\nli3xkjRaas3W6cBfdzT9kOJHMI9rN+pfrzVbDwL2AnYAN7cb9buHEKYkTTQTL2lC1ZqtvYC/AN7S\nY5HZXx5/rdZsdXv+6yl+obymXM9akzFJWhoTL2nC1JqtiykKxtd0efgPgI3As4CfApdT9E0F8H7g\neOCzwDOAN5a3WacDbyi3sSfwaOA7FJ2Rvhj4artR31jt3mhXWeMlVWNsa7wkDVat2QrgT4Aj5jy0\nD/DYdqP+hY62D3ZM33NwqTVbrwF+DBwI/C7whx3rm6k1W5cChwPryrZbyvVDcYly5dL3RP0w4ZKq\nYY2XpF31p8D7gDOBD5S3I9qN+u3AF3o/7V7tRv1/y8nNwLvKG7Vm6wnAl4ELOxa/jXuTrs8Cz6g1\nW09vN+qfW9JeSNIEM/GSRlSt2Xol8E6KcQOP7lZfVWu2HkvR+ec7gdcBr2836m8qH358VbG0G/XL\na83W/sANZdPRwBeBm8v5F1J0lnx+rdnavd2o31nVtiVpkph4SaPrMcDdwG8Dd9Warb3ajfrW2Qdr\nzdahwBXl7B+X919drmDajfoP2Xl8Q2rN1orysaw1Wz8HfAn4OYoOSTVA1nhJ1bDGqw8R8WLgXzPz\nR/Ms8zTglzPz7D7W3/dzpW5qzdbu5eQdwC9SjPzwJ8DvUSQyfwf8pDzD9S1gz/JxgO+VywB8ZVAx\nQ5FwdUxfW2u2frbdqO8YZAwqmHBJ1bDGa5Ei4k+BY4HLgB9FxC8A/wi0gKdS/JorKX699dCIuH9m\n/u0i1r97t+dGxAcovgB/BjgjM6+LiBOAnwV+CXheZt4dEa8BVlH8ouzNmXlHRJwO7F22vyozb1vi\ny6ARVmu2jgAuBl7YbtQ/XGu2jgfO7Vjk34EHU5y92tBu1LfWmq2jgGOAbwDXA48ol50dZ/AA4Aed\nidAwmHRJ0vwmLvHKzHMj4sDOJuDrmfnGiHgD8MjM/K+IOBd4zGLPWmXmT3s8N4H3Z+Z3Oto+TFFn\ncyjwqHIMq1WZedrsAhHxSxS/EtsI3B/4VXaxEFrDU/6CcG/gxwslG7VmazXFrwN3595fAwJ8qNZs\nfQjoPDN7HHAU8ALguNlLi+1G/Zm1ZuurFO+nR3Qs/69lsvX9Je6SJGkAJi7x6uGu8n47sKKczo7p\nxer13NlfhBERDwD+maKn8B9RvNZ3dXne3cB3MvPtfcaiPpSJ08p+OgStNVuPAd5Ecfb007Vm643t\nRr1rbVWt2fo7is5H51oNXEJxpuphwB9RvH8uBrYAz2s36v/Z+YR2o35Yuc4nUvQsf9ViY9fkssZL\nqsbY1XhFxHsoOlZcAfxJZv53RKzl3v/012XmfbvJXl7dLr9cDfxVROwPvDMzb+6yTC+zzz0AeEeP\n526nKET+Q4rLj3tl5uURsSUi3koxOO9ZmXlVRPwoIt5B8fP8dy0yFvXndOB5wGMX86SyxuqKjqZn\nlbeoNVsPo+ig9FntRv3ptWbrWO5Nui6luJz4z0C2G/U7a81WrXzsOopE6t/L+X+jOPvZVbtR/9Ji\nYtZ0MOGSqrHcn6XIXJ6SkIioA88BTqK4dLa2fOgC4MjsseGLLroo16xZ4wFEy6bWbN2Pok+qxwPH\nthv1j5Tt+wLb2o36tjnLr6BIpDt9i6K39i+W84cCnwN+oZxfTZGgzz723726WCh/DXjtsOuzhmWS\nPvNV70tEPBB4eWa+acGFi+VPyMz3LXGb7cys9Xhs4PFovJQlNWTmRHymd8ViP/fLeanxNuBOii+g\nqzJzG0BEbAYO5t4vpZEQEa/t0vyWzLRYeALUmq1HUdRBHQJ8rWy+Gfj7WrO1iuJXg7PLvrrdqDfL\n6d3YuZd3gEe0G/Xvl48/nOIM2H/NWWb2/b1gn1btRv17i98jTYPMvJ3isvau+jOKTnSXxajFI42j\n5Uy8jgfOAvYDbomIM8v2W8u2kUq8MvOvhx2DlqbWbJ1C8WvRjwKHUVza+wLwVuCbcxb/T4ozVpfR\nkXSV3l5rtmZr7r5argtgX2C3sj8rANqN+o/K/rSuovi14R8BpwI/oCh8tyNR9SUijgNeBDyw8wxU\nRFwBbACeArQy83SAX/3VX83ddtuNiNhYtr+xy2p3xW4R8c65658nnpMofgyyHfhAZp5btn8IeHQF\n8UgDNXY1XgARcQxwZWZ+OyIOoRhW5CSKmqdzgBuXY7uaLrVmax+KM6prKIa2me0W5O87FvsFii+F\nWS9pN+r/p2MdHwX+H4r35WUU9XjfAT5eLnIYcH67UX9GrzjajfoW4KEdTX+2+L2RdpaZHwQ+GBHt\nOQ+tAt5OcRz9GkW9Il//+tejvEz4m0vc9L7d1j9PPM8Bfjczb5kT/wsrikcaqLHrxysiDqOo4Tq1\nbNpMcXln1urMtFdrVeFMiv/AAd5S3j8beHp5S4rBnv+F4kzUL7Yb9SvnrON44HXtRv2/57T/S63Z\n+ifg+RRdPEij4obM/CFARNwxAus/FnhJRDwE+FRmXroMMUm75D82f/9Fd+7Iw6ta3wPut+Ljv/LI\nA9ZXtT5YnjNe5wHXlqeXr8jMUyJihuLUOMAZy7BNTahas3UAcEM5JM3eFGegZofSeVHHoscD7283\n6p8APtFlPae3G/XNc9vbjfodwNyka/axP+beoXikcbBbRKwYZG1qZl4HvC0i9qD49e5hHQ8PPB5N\nt5u23vX7p39+8zOrWNduK4O/OeaQO4DRTrwy86AubeupOHBNtlqz9SSKS4dPKOffB5zQZdHHthv1\nb5bLfKzX+rolXdKYmPtr1+w2PTMzk0960pO47LLLzo+I/8nMl1awvW6/tN2prazffRzFaAtzO6Te\nACw1HmmXJd3ftItx9NZNAFzwoCOWHE8309KBqsZIrdk6CHghZdJVOoHioH4GxXiGRwNPnk26AOZ2\nAyFNgsx8Qq/5zul169bFunWdAyMsfXtzt90jnlfMs65Tez0mjarP7FlcqezZmeISmXhpZNSard+g\n6Lrh0WXTcRS/UFwF/N92o/7nHYt/rLxJ6qEs+ZgrM7M+8GAkASZeGhG1Zuu57JxI/QS4uN2o30Ux\n5NKBXZ8oqSd/USiNHhMvDV2t2Xo0O3dSune7Uf/xsOKRxpFjNUrVsMZLE6ccsucjwHMp+sv6A4rL\n6Y8Grmk36v4CSlokEy6pGtZ4aaLUmq25Pzh5bnn/+najftWg45EkaZBMvDRMf0GReH0G+OSQY5Ek\nadmZeGlgas1W55nbfdqN+q3Au4cVjzRJrPGSqmGNlyZCOa7i6eXsSuu4pGqZcEnVsMZLY6nWbD0e\n2BP4LnA9cBHweOB7Jl2SpGll4qXK1ZqtFwPv6/Hw6kHGIknSKBlY4hURa4HZ8SzWZWZrUNvW4NSa\nrRdQJF0XAF8GjgA+DVzRbtT9m0vLxBovqRoTUeMVESuAGWBt2XRBRGzMzKWOZanRM9s9xMfajfoH\nhhmINE1MuKRqLHeN14plWu9cq4GrMnNbZm4DNgMHD2jbWga1Zitrzdan5rR9mGLw6quAfxpKYJIk\njbBBXWrcF7glIs4s528F9gOuHtD2VaFas3VqOfmsWrPVBJ5KkWi9oGw/vBxjUZIkdRhU4nUTsA9w\nEhDAOcCNA9q2KlRrtt4F/DlwIvAGYDYJeyKwBTiy3aj/aEjhSVPLGi+pGhNR40VxafGQjvnVmXnN\ngLatJag1W3sA0W7Ut9WarY8Afwwc0W7UN5WDW78CeDJwLXBAu1H/xhDDlaaWCZdUjYnoxyszt0fE\nDLChbDpjENtVJT4OPLPWvOcHiSe2G/VN5XQDeG27Ud9Wzl836OAkSRonA+tOIjPXA+sHtT1V5tc6\npr/dbtTfOzvTbtS3A9vu+xRJktTNoH7VqDFUa7ZWAA8DasAe7Ub90CGHJKmHmZmZnK3zktS/o7du\nuqfOaznYc73m85cUZ7S+2m7UPaBLI8waL6kaE1HjpfFSa7YeBvywnD3CpEuSpGp4qVE7qTVbAbx0\ndr6jkF6SJC2RZ7ymXK3Z+nngz4CPAq8Efh14LPAV4PlDDE3SItiPl1SNSenHSyOoPLv1nXL29I6H\nPtJu1I8dfESS+mXCJVVjUsZq1GhaU95fWt4/D1hp0iVJ0vLwjNd0eyrwmXajfsywA5EkaRqYeE23\nBwIXDzsISUtnjZdUDWu8VLlas3UXcCXwy8DLhhyOpAqYcEnVsMZLlao1W6+kSLh/Gfg34J+HG5Ek\nSdOj0sQrIt4TERsj4uKIOKijfW1EbCpv9Sq3qd5qzdbKWrP16DnNK4DrgROAP2436rcMPjJJkqZT\npZcaM/NEgDK5agAvjYgVwAywtlzsgojYmJn2hl6hWrN1MLA/cBlwN/AG4BXAg2rN1hOADwG3Ag8C\nzm036v8wrFglVc8aL6ka41rjdRtwZzm9GrgqM7cBRMRm4GDg6mXa9lSpNVvHAR/o8tBfdUxfPucx\ne6OXJowJl1SNkRyrMSKOAl49p/lVmXlFOX08cFY5vS9wS0ScWc7fCuyHiVdfas3WY4B9KC4XHktx\nNrGbLwO/QZF0PaFsOw94fbtRv3K545QkSffVV+KVmRuADd0ei4hjgCsz89tl000UicJJQADnADf2\ns91pV2u2HgL8B93/bj8CDqToBPV24NPAAe1G/frBRShJ0q7bsmXLoV/7/m1n/fTuHXdVsb6EQ6pY\nz3Kq9FJjRBwGHJmZp3Y0b2bnF2J1Zl5T5XanQa3Zej7wT10e+hjFsD9/227U7wI+3PGYSZc0Jazx\n0ph6yEe/dsOa//j+7SPTy8K41XidB1wbERuBb2TmyZm5PSJmuPcM2RkVb3NavLW8v4piMOtbgV8A\nPtFu1LcOKyhJo8GES6rGSNZ49ZKZB/VoXw+sr3JbU+gK4GfbjXpn9xBfGFYwkiRp8Ubm1J56qzVb\nfwY8A3j+sGORJEn9c8ig8fAn5f15Q41C0siyxkuqxrjVeKlitWbrAcATgZe3G/W7hx2PNEkiYneK\nusm3Z+bZEbEWWFc+vC4zW+VyXdtHiQmXVI2xqvHSsrgeoN2on7XQgpIW7UTgq0BGRDBnlA2g5egb\nkqpkjdcIqzVbBwCrhh2HNIkiYk/gKIo+74KOUTbKkTY2R8Tqbu0Uo29I0qJ5xmtE1ZqtvYC/pBh3\n8fFDDkeaRCcD7wYeXs7vR/dRNqJH+0iNvmGNl1QNa7ym19nAccC32436N4YdjDRJImJv4KmZ+daI\neFHZ3GuUjRU92keKCZdUDWu8ptdu5f1LhxqFNJmeAuwRER8FHklxLNxEl1E2ImJlt/bBhSppkph4\nja5HAGvajfq/DzsQadJk5vnA+QARcRywV2Ze0W2UDUffkFQlE68RVGu2fg44AvjfYcciTbrM/GDH\ndNdRNsZh9A1rvKRqWOM1nb5Y3pt4SdolJlxSNZa7xqvy7iQiYveI+G5EvKyjbW1EbCpv9aq3OUlq\nzdaDKS4zwggW8EqSpP4txxmv2Q4JAbDzwV1Xa7b2Ba4DvgI82Z7qJUmaLJWe8ZrTIeEsOx/cdU8G\nHgCc1G7U7xp2MJLGx8zMTM7WeUnq39FbN91T57Uc+jrjFRFHAa+e0/wq4Ons3CEhwL6MQeeDI2J3\n4JPtRr097EAkjRdrvKRqjGQ/Xpm5gXt/Wg3c0yHh4Zn5to4OCaF3p4Qq1ZqtFcAHgF8Arh1qMJIk\nadlUWeN1nw4JI2IjcCV2PriQXwKOLadbwwxEkiQtn8oSrx4dEn6rnLfzwS5qzVYC/wY8Dvgp8DTg\nW8OMSdJ4sh8vqRpj2Y9XZ4eE5fzIdz44aLVma2U5ecxsW7tR/9KQwpE05ky4pGqMXT9e2mVPKO+/\nXt5fNKxAJEnSYJh4Dc/vA//abtQfB6wEfmvI8UiSpGXmkEHD80TgXIB2o75jyLFIGnPWeEnVGMsa\nL82v1mw9ENgT8NedkiphwiVVYyT78dKS3QDsBfxk2IFIkqTBscZrOPYq7zcPNQpJkjRQJl4DUmu2\notZsdXYkS7tRv21Y8UiaLI7VKFVjJMdq1OLUmq0HAF8CfqXWbH0duJNiXEtJqoQ1XlI17MdrMjwG\n+JVy+leBDwL/PrRoJEnSUHjGa5nVmq3fAZ4PXAY8G7h/u1H/3nCjkiRJw2DitYxqzdYvAp8rZ89u\nN+o/GGY8kiaX/XhJ1Rirfrwi4kDgw+V625n5yrJ9LbCuXGxdZraq3O6oqTVbKyheg/eXTX8BvGd4\nEUmadCZcUjXGrR+vdwCvzcwvzjZExApgBlhbNl0QERszc+J+fVNrtu7fbtTvAl4LvIFiv1vtRv3d\nw41MkiSNgsqK6yNiJfCozqSrtBq4KjO3ZeY2ir6rDq5qu6Oi1mytAe6sNVtJkXRBcZbvluFFJUmS\nRklfZ7wi4ijg1XOa3wjsERGfAh4MvCszPwnsC9wSEWeWy90K7Adc3V/II2v9nPlrgW8AHxpCLJKm\njDVeUjVGssYrMzcAGzrbIuJ+FEnVs4GVwKUR8XngJmAf4CQggHOAG5cQ86i6E9gDeBiQ7UZ9EvdR\n0ogy4ZKqMTb9eGXm3RRnefbPzDuBn5YPbQY6e2xfnZmTODj018r7G026JElSN1UX158GvC8i9gb+\npazpIiJmuPcM2RkVb3NUXAu8t92oT9yPBiRJUjUqTbwy83vA73ZpX899a6AmRq3ZOh54DvA3w45F\n0nSyxkuqxkjWeOletWbrQODcctZLjJKGwoRLqsbY1HhNsc4zeZuHFoUkSRp5nvFauh3AZ4BPWd8l\nSZLmY+LVh1qz9UfAKcArgL2Bk9uN+v8MNypJ08waL6ka1niNpo+X918u728eViCSBCZcUlWs8Rpt\n/wA8qt2o3zrsQCRJ0ujzjNci1ZqtnysnX95u1M8aajCSJGmsjGTiVWu2dms36nfWmq0TgCPbjfoL\nOh57DfDLwN3A8UMoaH9mef++AW9XknqyxkuqxrTWeP201mzdQTH2IbVm65x2o/7FWrP1JOAtHctd\nBHxkwLHdCZzbbtS3Dni7ktSTCZdUjeWu8RrVxAvKpKt0aa3Z6rbMoQOKhVqztQ9wJPBeYCvw4kFt\nW5IkTYZRLq6/YZ7H/qq8H2QXDluAT5XTlw1wu5IkaUJUmnhFxAsj4ssRcWlE/GZH+9qI2FTe6ruw\nqgPbjfoBwNry9mzgNcA7gP/TbtTfSDFMz7BOrf/1kLYrSV3NzMzkbJ2XpP4dvXXTPXVey6HqS42n\nAr8G7AVcADwpIlYAMxQJFMAFEbExM3seINqN+vXl/UXzbGs7sLKSqBdQa7YOKic/3W7Uf28Q25Sk\nxbDGS6rGuNV4fYuiDmp/4Etl22rgqszcBhARm4GDgauXuK2BJV4Uv6LcAvzhgLYnSZImUF+JV0Qc\nBbx6TvOrKAaMfjmwG3B22b4vcEtEnFnO3wrsxwgmXrVm6zeAPduN+sY5D+0FXNhu1O+ucnuSJGm6\n9JV4ZeYGYENnW0QcBBydmc8s5y+JiAuBm4B9gJMoarLOAW5cStCl7cBDyw5NP0mR8P0B8Lx2o35A\nrdk6DnhMu1Fv1Jqth7cb9R8C1JqthwE/BW4DDgMuB+7XbtS3AxcCDyzjpNZsvQf4AXB9ubwkjST7\n8ZKqMU79eK2cXV9EBPAAIIHNwCEdy63OzGsq2N4vAr8NvK6cv2T2gVqz9QzgA+X0qeX9DPAvwH92\nWdfdnd1V1JqtdwKvnLPMByuIWZKWhQmXVI2xqfHKzKsj4ksRcT7FryXPzsw7ACJihnvPkJ1R0SZX\nz/PYZ7q0rStvu2Ju0gXwsF18riRJUleVFtdn5pt7tK+nqP+q0v4d01uAVeX0ayjqza6j+IXlXwAn\nUhTI/1+Kff4E8MZ2o76j1mz9KXAQ8BKK2jOALwLXAF8u1/s14JsVxy9JkqbMKPdcv5B/oKjzehXF\nZcdDKbp72A68rdZs7QE8t92of7Cs1VoN3NRu1H/UuZJ2o35uOfnaWrP1SeCb7Ub99QPbC0mqgDVe\nUjXGqcZroNqN+ikds/9V3jofv4OyLqv8NeJOj/dY5+9XGaMkDYoJl1SN5a7xGuUhgyRJkiaKiZck\nSdKAjO2lRknSvazxkqphjZckaUEmXFI1rPGSJEmaECZekiRJA+KlRkmaANZ4SdWwxkuStCATLqka\n1nhJkiRNCBMvSZKkAekr8YqIwyPi8ohozmlfGxGbylt9oXZJGoaIeFNEtCLiwog4qGwb6+PXzMxM\nztZ5Serf0Vs33VPntRz6rfHaHXgL8OTZhohYAcwAa8umC4BWt/aI2JiZHiAkDUVmvg4gIp4CnBYR\nJzLmxy9rvKRqjGSNV2ZeCNw8p3k1cFVmbsvMbcDmiFjdrR04eClBS1JFngj8Fx6/JA3IvGe8IuIo\n4NVzml+VmVd0WXxf4JaIOLOcvxXYD4ge7Vf3HbUkLVFEXAI8BDgcOIQxPX595ZrrL/jGDT95dBXr\n2n1lrHjs/g8897GPPGCmivVJuq95E6/M3ABs2MV13QTsA5xEcbA6B7iR4qxat/aeLrroopE7jS9p\nsmTmERHxBOBDwCuo6Pg1aP97+107zrnsup+frUmZvUzSjwfvvpK3Pv3g3aqKTRpHo9yP19x6gs0U\n/zXOWp2Z10TEym7tvVa6Zs0a6xQkDcoNFMfBa6jg+DVMS0m4JN1ruWu8+kq8IuI04OnA/hHx4Mx8\nSWZuj4gZ7j1DdgZAr3ZJGpaI+DjFZcY7gT/PzB0evyQNQl+JV2a+DXhbl/b1wPpdbZekYcjM53Zp\n8/gladk5ZJAkTYAqarwkjXaNlyRpRJhwSdUYyX68lsu49BDdr4j4QERcFhEbI+KFZdtY95a9kG6j\nHCx2n8f1teix753vgeM62idt399T7uPF/fYMP677LknzGZkzXuPUQ/QSJPDczPweTE1v/zuNcrCY\nfe7VPkavxX1GeGDOewAm832QmScClAlTIyJOYnr+7pLU08gkXnT0EA0QEbM9RI9MR4UV6ewu4z77\nXPaWvWJuO2P6WmTmhRFxZEfTLu/zuL8WXfZ91twuUyb5fXAbxS8Hp+bvPizWeEnVmKYar14930/S\ngfY24J8j4maKDhunsbf/xe7zpL0WO70Hyj6hJvl9cDxwFkXc0/x3X3YmXFI1RrIfr2XSq+f7iZGZ\nJwNExOOAJnAaY9pb9hIsdoSDiXoturwHfp8KR30YJRFxDHBlZn47Ig5hiv/ukjRrlBKvrj3fDyuY\nZXYHcBcT0Fv2Luq8tLaoEQ4m4LXoNRLD7HsAKhr1YZRExGHAkZl5atk0bX93SepqZBKvaeghOiI+\nBhwA3A6cNA29ZUeXUQ4Ws8/j/Fr02PePA/tTXHJ8GUzmvgPnAddGxEbgisw8ZVr+7sNijZdUjWmq\n8Zr4HqIz83ld2ia6t+xuoxwsdp/H9bXose/36TG9bJ+0fT+oS9tU/N2HxYRLqsZU9eMlSZI0yUy8\nJEmSBsTES9LYi4jXRMShw45jmI7euume2hRJ/Vvuz9JI1XhJUp/+EXh6RBwLbAUuyMz2kGMaKGu8\npGpY4yVJCzsAeDjFMe0nwMM7x8iUpFHhGS9Jk+BxwLmZeU8nqxFx5xDjkaSuPOMlaRJ8ejbpioiH\nwj3dUUwNa7ykaljjJUkLOxmYKadfArxpiLEMhTVeUjWmaaxGSerXQyJiT2Al8DPDDkaSejHxkjQJ\n3gG8rpx++zADkaT5mHhJGnuZ+V3g9GHHMUyO1ShVY6rGapSkfkTEPwH/CyRAZr5iuBENngmXVA1r\nvCRpYd/MzLcMOwhJWoiJl6RJ8MsR8XfAdpjOM16SxoOJl6Sxl5kvGHYMw2aNl1QNa7wkaQERcTzw\n6Mw8LSKem5kfH3ZMg2bCJVXDsRolaWGHADeV078yzEAkaT6e8ZI0CW4G1kTEI4Abhh2MJPVi4iVp\n7GXm2yPiPcCOzLx92PEMgzVeUjWs8ZKkBUTEuo7pzMw3DDOeYTDhkqphP16StLC/Le/vD7xsmIFI\n0nxMvCRNgqeV9yuAxw0xDkmal4mXpElwS3m/A3jJMAMZFmu8pGpY4yVJu24FcGhEHAqQmRcPOZ6B\nMeGSqmE/XpK0sBOAfYGHl9P7lDdJGime8ZI0CW7OzE8CRMRTMvPTww5Ikrox8ZI0CT5XDpK9O/DJ\nYQczDNZ4SdWwxkuSFpCZn4uIL2bmrcOOZVhMuKRqWOMlSQuIiDcDbymnTxtyOJLUk4mXpEmwAvhO\nOX3gEON3UfQyAAATdUlEQVSQpHmZeEmaBJcDvx4R/wh8ftjBDMPRWzfdU5siqX/L/VmyxkvS2MvM\nTwCfGHYcw2SNl1QNa7wkaQER0Rh2DJK0K0y8JE2CfSJiz2EHIUkL8VKjpEnwIKAZEXcCZOYrhhzP\nwNmPl1QN+/GSpHlExBMz8+SIeGRm/s+w4xkWEy6pGtZ4SdL8nlXenzDUKCRpF5h4SRp3vx4RpwC1\niDilnJakkeSlRknj7gQggU8NO5BhssZLqoY1XpI0j8z8zrBjGAUmXFI1rPGSJEmaECZekiRJA+Kl\nRkmaANZ4SdWwxkuStCATLqka1nhJkiRNCBMvSZKkAfFSoyRNAGu8pGpY4yVJWpAJl1QNa7wkSZIm\nhImXpKkTEe+JiI0RcXFEHFS2rY2ITeWt3rFs13ZJ6oeXGiVNncw8EaBMpBoRcRIwA6wtF7kAaEXE\nirntEbExM3PQMS/EGi+pGtZ4SdLyuQ24E1gNXJWZ2wAiYnNErKa4KrBTO3AwcPWQ4u3JhEuqxnLX\neJl4SZpmxwNnAfsBt0TEmWX7rWVb9GgfucRL0ngYucTroosuGrlT+JKW15o1a2LQ24yIY4ArM/Pb\nEXEIsA9wEkWydQ5wI8UZr27tktSXkUu8YDgHYUnDMYx/tiLiMODIzDy1bNoMHNKxyOrMvCYiVnZr\nH1Sci2GNl1QNa7wkqXrnAddGxEbgisw8JSJmgA3l42cAZOb2bu2jyIRLqoY1XpJUscw8qEvbemD9\nrrZLUj/sx0uSJGlAPOMlSRPAGi+pGtZ4SZIWZMIlVWNkx2qMiMMj4vKIaO7Csg65IUmSpt5Sznjt\nDrwFePJ8C43TkBuSJEnLqe8zXpl5IXDzLix6z1Ac5bAbs0NuSJIqcvTWTffUpkjq33J/lgZR47Uv\nDrkhScvKGi+pGiNb47UIN1EMuXE68NpyesEhNyLijLnT3dp6PPfFEfGwBdb/tIh42UJxSJIkVWWp\nideuDO3TdSiOXXjeui7T3dp2Sswi4k+BY2cfj4i/jYiNEfH6iLgmIvaIiDdSDI57dER8vnMdc6e7\niYgTIuINEfH/RsT9yrbXRMTbyvY9ImJlRJwdEW+KiPdFxF7lch8ol/mHiHhEl7YDd+G1kSRJY2gp\nv2o8jWL4jGMi4r0d7c+JiGfMzmfmdori+g0UvT+f0e8253FPYpaZ5wItikFtAU4Bvp6ZbwQeBTwS\neB1wLvAZ4LfnrOOe6XnOun0Y+DxwKPCoiHgXsCozTwN2ZOYdwG8BB2bm64A9KJJBgATen5kvzszr\nu7Rdt7SXQtI0ssZLqsZyf5aWUlz/tsx8Wmb+Yma+pKP9vMz87Jxl12fmU8vbhvuubdnd1TE9u8/J\nwvvf66zbR4E7gV+iqJP78451zS6XwDPL6ReU8wCPA/53zlm1+7TNd3l1MWfnJE2Hz+x5uHVeUgWW\n+7M0LUMGdeu64mrgWQARse8i1xfAH5bTe5X3WyLireX6HkQ5tls5wC7AR8r7x5X3nUldt7b5Lq/u\n0tm5zul+2iRJUrUmMvHKzJk586/umP7P8v4Hmbm2nN6VbjE61/d7mfmacvry8v7NHW23ZeaOcnpd\nef+Tjuf/pMs679O2SF3r31hkAkeP2rkq2iRJmnYTmXj1KyJeGxGvnZ0edjwjoO+krVubCZy0fKzx\nkqoxsjVekygz/zoz/3p2etjxTKCRORMnTRprvKRqWOMl3VcVZ+K8lCpJGjgTL02zZbmUupizbiZw\nkjRdTLyk6iz6rFu3x72Uqn5Y4yVVwxovafos+UycSdv0scZLqoY1XpL6MbT6N0lSbyZekroZaK2b\nJE0LEy9JVeun1k1LZI2XVI3l/izdb9nWLEkaGOu7pGrMfpbuv0zr94yXJEnSgJh4SZIkDUjfiVdE\nrI2ITeWtvsCyL4yIL0fEpRHxm/1uU5LUnTVeUjVGssYrIlYAM8DasumCiNiYmdnjKacCvwbsBVwA\nPKmf7UqSurPGS6rGctd49Vtcvxq4KjO3AUTEZuBg4Ooey38LOBLYH/hSn9uUJEkaa/0mXvsCt0TE\nmeX8rcB+9E681gMvB3YDzu5zm5IkSWOt38TrJmAf4CQggHOAG7stGBEHAUdn5jPL+Usi4sLZs2WS\npKWbrUnxkqO0NLOfpQsedMSyrL/fxGszcEjH/OrMvKbHsitntxMRATwA6FULJknqgwmXVI2RrPHK\nzO0RMQNsKJvOmH0sIp4DbM3Mz5bLXh0RX4qI8yl+RXl2Zt6xtLAlSZpMW7ZsWQkcVvFqv7pq1art\nFa9Tfei75/rMXE9RuzW3/bwubW/udzuSJE2ZvS757y3/tuHqm/euYmVHrd73liMOWrUauK2K9Wlp\nHDJIkiaANV6T5Ttb7rjry9f+ePcq1nXIQ/e8q4r1TItRrfGSJI0QEy6pGo7VKEmSNCFMvCRJkgbE\nS42SNAGs8ZKqYY2XJGlBJlxSNazxkiRJmhAmXpIkSQPipUZJmgDWeEnVsMZLkrQgEy6pGtZ4SZIk\nTQgTL0mSpAHxUqMkTQBrvKRqWOMlSVqQCZdUjeWu8eo78YqItcC6cnZdZrbmWfZA4MPl9tqZ+cp+\ntytJkoZny5YtP88CecmWLVsO3sXVHbj0iMZLX4lXRKwAZoC1ZdMFEbExM7PHU94BvDYzv9jP9iRJ\n0mj4xg9u/8z6q296xHzLvPOS716+q+v77pY7pqrevN8zXquBqzJzG0BEbAYOBq6eu2BErAQeZdIl\nScvHGi/18t0tdzxo4+Yt63fkzb1OjizK97bcse8FV928ar5lFnp8lI1qjde+wC0RcWY5fyuwH10S\nL+ChwB4R8SngwcC7MvOTfW5XktSFCZd6ueR/btn7kv+55YnDjmNcjGqN103APsBJQADnADfOs+yt\nwLOBlcClEfH52bNlkiRJ06Lf66qbgUM65ldn5jXdFszMu4Brgf0z807gp31uU5Ikaaz1dcYrM7dH\nxAywoWw6Y/axiHgOsDUzP9vxlNOA90XE3sC/eLZLkqpljZdUjVGt8SIz1wPru7Sf16Xte8Dv9rst\nSdL8TLikajhWoyRJ0oQw8ZIkSRoQhwySpAlgjZdUjZGt8ZIkjQ4TLqka1nhJUsUi4vCIuDwimh1t\nayNiU3mrL9QuSf3wjJekabQ78BbgydB9/Fmg1ce4tJI0L894SZo6mXkhcHNH0z3jz5b9DG6OiNXd\n2inGpR05R2/ddE9tiqT+LfdnyTNektR7/Nno0d5tXNqhssZLqsaojtUoSZOk1/izK3q0S1JfTLwk\nTavomO46/mxErOzWPpDoJE0kEy9JUyciTgOeDuwfEQ/OzJd0G392vnFpR439eEnVsB8vSapYZr4N\neNuctl7jz3ZtHzUmXFI17MdLkiRpQph4SZIkDYiXGiVpAljjNVxbtmzZA1hZ0er2qmg96sPI1nhF\nxFpgXTm7LjNbCyy/O3AV8PbMPLvf7UqS7suEa7iuvnHrJy6/9sePrWp97et+/JCq1qXFGcl+vPoc\nRuNE4KuAQ21IkibKlm137fjgV39w4LDj0Ojrt8ZrUcNoRMSewFHAp9m57xxJkqSp0e+lxl7Da/Qa\nRuNk4N3Aw/vcniRpHtZ4SdUY1RqvXsNr3EdE7A08NTPfGhEv6nN7kqR5mHBJ1RjJGi96DK/RY9mn\nAHtExEeBRwL3K+vBvtXntiVJY2DLli0rqbbboly1atXdFa5PGri+Eq/5htGIiOcAWzPzs+Wy5wPn\nl48dB+xl0iVJk+97W+5451eu//EfbN8x7w+vdtnjH/Ggr/36qlXPqmJd0rD03Z3EPMNrnDfPcz7Y\n7/YkSb2NYo3XHdt37H7u5d//2bt2VPNj9jf/zqO+VsmKpHmMao2XJGmEjFLCJY0zx2qUJEmaECZe\nkiRJA+KlRkmaAKNY4yWNI2u8JEkLqiLhunN7csUNt7/omz+86vkVhMSOHbn79mp+0CgNzKj24yVJ\nmjB33L2D93zp+p8ZdhzSJLPGS5IkaUA84yVJE8AaL6ka1nhJkhZkwiVVw368JEmSJoSJlyRJ0oB4\nqVGSJoA1XlI1rPGSJC3IhEuqhjVekiRJE2JJiVdErI2ITeWtvsCy74mIjRFxcUQctJTtSpIkjaO+\nLzVGxApgBlhbNl0QERszu48PkZknls+rAw3gpf1uW5K0M2u8pGqMco3XauCqzNwGEBGbgYOBqxd4\n3m3AnUvYriRpjmlIuB60+8pV1//wxt+pan177rbyy6tWrdpS1fo0GUZ5rMZ9gVsi4sxy/lZgPxZO\nvI4HzlrCdiVJU+jvLr328PuvXPG5Ktb16IfsmS990oFHApuqWJ+0q5aSeN0E7AOcBARwDnDjfE+I\niGOAKzPz20vYriRpCl1947bK1rX7yuhaFiMtt6UU128GDumYX52Z1/RaOCIOA47MzL9dwjYlSV0c\nvXXTPbUpkvq33J+lvs94Zeb2iJgBNpRNZ8w+FhHPAbZm5mc7nnIecG1EbAS+kZkn97ttSdLOpqHG\nSxqEUa7xIjPXA+u7tJ/Xpc0uJCRJ0lSzA1VJkqQBccggSZoA9uMlVWOU+/GSJI0IEy6pGo7VKEmS\nNCFMvCRJkgbES42SNAGs8Vqc7Ul8d8u213/zhttvqmJ9e9xvxSELL6VxYI2XJGlBJlyLc8UPbo8T\n/r9vHzXsODR6rPGSJEmaECZekiRJA+KlRkmaANZ4SdWwxkuStCATLqka1nhJkiRNCBMvSZKkAfFS\noyRNAGu8pGqMbI1XRKwF1pWz6zKzVcWykqTFM+GSqrHcNV59JV4RsQKYAdaWTRdExMbMzKUsK0mS\nNMn6rfFaDVyVmdsycxuwGTi4gmUlSZImVr+XGvcFbomIM8v5W4H9gKuXuKwkqQ/WeEnVGNUar5uA\nfYCTgADOAW6sYFlJUh9MuKRqjGSNF8Xlws6R2Fdn5jUVLAtABAlZ3sO9093aFnp8WG3GMNpxGcOo\nxHXhhUjS1Ogr8crM7RExA2wom86YfSwingNszczPLrRs7/UT/cQlafxcdNFsIiZJk6/v7iQycz2w\nvkv7ebu6rCSpGtZ4SdUY1RovSdIIMeGSquFYjZIkSRPCxEuSJGlAvNQoSRPAGi+pGtZ4SdIIGPUx\nZ024pGqMaj9ekjQ1HHNWUlWs8ZKkhTnmrKRKeMZLkhY2smPO7rvX/Vce+/j9r9vyhfMOBFj11Odc\nN+yYNL2+AgcCHPv4/cf2fTj7WXroEX90/W4r486q1x+jdqb8oosuGq2AJC27NWvWjPRoFRFxCPCX\n7Dzm7JvmDn/m8UuaTos5ho1c4iVJoyYiVgKXUNR4BbAhM58y3KgkjSMvNUrSAvoZc1aSuvGMlyRJ\n0oD4q0ZJkqQBMfGSJEkakJFKvCJibURsKm/1YcezqyLiPRGxMSIujoiDyrax3JdZEbF7RHw3Il5W\nzo/l/kTEgeXfZlNEvLNsG7t9iYgXRsSXI+LSiPjNsm0s9iMiDo+IyyOi2dHWNfZx2Sfovl9l+9D3\nrdv7ZdAxzImn83P4N8OOp2P7Ox3nhhlTt++RYcYz7G2X2x/Z79YlfUdm5kjcKJLAS4EHlLdLKGvQ\nxuUG1IG/p/jV07jvyynAJ7j35/NjuT/Ax4And8yP5fsMuAJYCTwYuGyc/iYUvwT8faDZ628wjn+b\nufs1Svs29/0y7Nd37udw2PF0xDB7nHvZCMVUB/5+2PGMwmsx9zUZpePeUr4jR+mM1yT0DH0bcCdj\nvi8RsSdwFPBpijfUWO5P2QXAozLzix3NY7kvwLeAI4GjgS8xRvuRmRcCN3c03Sf2iFjdrZ0R3Sfo\nul8wOvs29/3SNbZljgHo+TkcWjwdcXUe50YiptLs98iw4xmF12LWSH23LvU7cpS6kxjZnqEX4Xjg\nLIq4x3lfTgbeDTy8nB/X/XkosEdEfIriP/93ATcwnvuyHng5xbit5zC+fxPo/VmPHu1D3aeIOAp4\n9ZzmV2XmFV0WH+i+9YqNe98vuwFnLxBbZa9vj3jeyJzPYWZ+chDxzBPTq4Cns/NxjkHEtAvvp9nv\nkYHEM49R+k4ete/WJX1HjlLidROwDzv3DH3jUCNahIg4BrgyM78dRS/XY7kvEbE38NTMfGtEvKhs\nHte/zU0UH4BnU1x2uRT4U8ZsX8rahqMz85nl/CXAnzNm+9Gh1/tpRY/2ocrMDdzbf9dCBrpv3WLr\n9n6JiAvnia0yPeK5H3M+hxHx+UHEM09MewOHZ+bbOo5zDCKm+d5Pnd8jg4pnHiNx3B+179YqviNH\nKfHaDBzSMb865wzHMaoi4jDgyMw8tWwa230BnkLx3+lHgUdSvEc2MYb7k5l3RcS1wP6ZeX1E/BS4\nhvHbl5WUn9WICIoagnHbj87hNLp+PspLUuO0T7DzfsFo7Fu390v2im2ZYrhHZt7d5XPIsOIp3ec4\nFxEbgSuHFVOX7xEY7ms09O+xEf1uXfp35DCK0uYpVvst4Avl7ahhx7OIuP8buBjYCJw1zvsyZ7+O\nA04a5/0Bfg44n+Js1ynjui/A6eV+fB540TjtB3Aa8O/At4H3zhf7uOxTr/0alX3r9n4Z5uvb7XM4\nKn/vzuPckF+jzu+Rvxt2PMPedpfXZOS+W/v9jrTnekmSpAEZpV81SpIkTTQTL0mSpAEx8ZIkSRoQ\nEy9JkqQBMfGSJEkaEBMvSZKkATHxkiRJGhATL0mSpAH5/wHXsyU4Qo8EtQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106cab6d0>"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gamma.stats()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "{'95% HPD interval': array([[ 0.,  1.,  1.,  0.,  0.,  0.],\n",
        "        [ 1.,  1.,  1.,  1.,  1.,  1.]]),\n",
        " 'mc error': array([ 0.02982943,  0.        ,  0.01438836,  0.01606565,  0.03076783,\n",
        "         0.02588083]),\n",
        " 'mean': array([ 0.38315789,  1.        ,  0.95157895,  0.06947368,  0.88631579,\n",
        "         0.79684211]),\n",
        " 'n': 1900,\n",
        " 'quantiles': {2.5: array([ 0.,  1.,  0.,  0.,  0.,  0.]),\n",
        "  25: array([ 0.,  1.,  1.,  0.,  1.,  1.]),\n",
        "  50: array([ 0.,  1.,  1.,  0.,  1.,  1.]),\n",
        "  75: array([ 1.,  1.,  1.,  0.,  1.,  1.]),\n",
        "  97.5: array([ 1.,  1.,  1.,  1.,  1.,  1.])},\n",
        " 'standard deviation': array([ 0.48615627,  0.        ,  0.21465427,  0.25425792,  0.31742733,\n",
        "         0.40234906])}"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pm.Matplot.plot(gamma)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Plotting gamma_0\n",
        "Plotting gamma_1\n",
        "Plotting gamma_2\n",
        "Plotting gamma_3\n",
        "Plotting"
       ]
      }
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pm.Matplot.plot(var)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pm.Matplot.plot(prec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pm.Matplot.plot(b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pm.Matplot.plot(taue)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}