2014_12_31a_pymc_initial_values_examples.ipynb

{"nbformat": 3, "worksheets": [{"cells": [{"cell_type": "code", "language": "python", "outputs": [{"output_type": "stream", "stream": "stdout", "text": "Wed Dec 31 09:55:52 PST 2014\r\n"}], "collapsed": false, "prompt_number": 1, "input": "!date\nimport matplotlib.pyplot as plt, numpy as np, seaborn as sns, pandas as pd\n%matplotlib inline\nsns.set_context(\"poster\")\nsns.set_style('whitegrid')", "metadata": {"trusted": true}}, {"cell_type": "code", "language": "python", "outputs": [], "collapsed": true, "prompt_number": 2, "input": "# set random seed for reproducibility\nnp.random.seed(12345)", "metadata": {"trusted": true}}, {"source": "From http://stackoverflow.com/questions/27713254/fitting-logistic-regression-with-pymc-zeroprobability-error:\n\n> To teach myself PyMC I am trying to define a simple logistic regression. But I get a ZeroProbability error, and does not understand exactly why this happens or how to avoid it.\n\n> Here is my code:\n\n", "cell_type": "markdown", "metadata": {}}, {"cell_type": "code", "language": "python", "outputs": [{"ename": "ZeroProbability", "evalue": "Stochastic observed's value is outside its support,\n or it forbids its parents' current values.", "traceback": ["\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mZeroProbability\u001b[0m                           Traceback (most recent call last)", "\u001b[1;32m<ipython-input-3-e2203c482542>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     12\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[1;36m1.0\u001b[0m \u001b[1;33m/\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mw0\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mw1\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     13\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m \u001b[0mobserved\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpymc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mBernoulli\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'observed'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlogistic\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mobserved\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/pymc/distributions.pyc\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, *args, **kwds)\u001b[0m\n\u001b[0;32m    269\u001b[0m                 \u001b[0mrandom\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdebug_wrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    270\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 271\u001b[1;33m                 \u001b[0mStochastic\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlogp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mlogp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrandom\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlogp_partial_gradients\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlogp_partial_gradients\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0marg_dict_out\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    272\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    273\u001b[0m     \u001b[0mnew_class\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.pyc\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, logp, doc, name, parents, random, trace, value, dtype, rseed, observed, cache_depth, plot, verbose, isdata, check_logp, logp_partial_gradients)\u001b[0m\n\u001b[0;32m    757\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mcheck_logp\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    758\u001b[0m             \u001b[1;31m# Check initial value\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 759\u001b[1;33m             \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlogp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfloat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    760\u001b[0m                 raise ValueError(\n\u001b[0;32m    761\u001b[0m                     \u001b[1;34m\"Stochastic \"\u001b[0m \u001b[1;33m+\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.pyc\u001b[0m in \u001b[0;36mget_logp\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    914\u001b[0m                     (self._value, self._parents.value))\n\u001b[0;32m    915\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 916\u001b[1;33m                 \u001b[1;32mraise\u001b[0m \u001b[0mZeroProbability\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0merrmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    917\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    918\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mlogp\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mZeroProbability\u001b[0m: Stochastic observed's value is outside its support,\n or it forbids its parents' current values."], "output_type": "pyerr"}], "collapsed": false, "prompt_number": 3, "input": "import pymc\nimport numpy as np\n\nx = np.array([85, 95, 70, 65, 70, 90, 75, 85, 80, 85])\ny = np.array([1., 1., 0., 0., 0., 1., 1., 0., 0., 1.])\n\nw0 = pymc.Normal('w0', 0, 0.000001)    # uninformative prior (any real number)\nw1 = pymc.Normal('w1', 0, 0.000001)    # uninformative prior (any real number)\n\n@pymc.deterministic\ndef logistic(w0=w0, w1=w1, x=x):\n    return 1.0 / (1. + np.exp(w0 + w1 * x))\n\nobserved = pymc.Bernoulli('observed', logistic, value=y, observed=True)", "metadata": {"trusted": true}}, {"source": "What has gone wrong?", "cell_type": "markdown", "metadata": {}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "pyout", "prompt_number": 4, "metadata": {}, "text": "array(-519.4387150567381)"}], "collapsed": false, "prompt_number": 4, "input": "# when you initialize the Normal Stochastics\n# their values are drawn from the prior\nw0 = pymc.Normal('w0', 0, 0.000001)    # uninformative prior (any real number)\nw0.value", "metadata": {"trusted": true}}, {"source": "Since you have a diffuse prior, this can lead to a value that has such small posterior probability that it is effectively zero.  The solution is simple: start from a better value:", "cell_type": "markdown", "metadata": {}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "pyout", "prompt_number": 5, "metadata": {}, "text": "array(0.0)"}], "collapsed": false, "prompt_number": 5, "input": "w0 = pymc.Normal('w0', 0, 0.000001, value=0)    # uninformative prior (any real number)\nw0.value", "metadata": {"trusted": true}}, {"source": "This can also make MCMC convergence faster, although any starting point that does not raise a ZeroProbability error should yield the same posterior distribution if your burnin period is long enough.", "cell_type": "markdown", "metadata": {}}, {"cell_type": "code", "language": "python", "outputs": [], "collapsed": true, "prompt_number": 6, "input": "w0 = pymc.Normal('w0', 0, 0.000001, value=0)    # uninformative prior (any real number)\nw1 = pymc.Normal('w1', 0, 0.000001, value=0)    # uninformative prior (any real number)\n\n@pymc.deterministic\ndef logistic(w0=w0, w1=w1, x=x):\n    return 1.0 / (1. + np.exp(w0 + w1 * x))\n\nobserved = pymc.Bernoulli('observed', logistic, value=y, observed=True)", "metadata": {"trusted": true}}], "metadata": {}}], "metadata": {"name": "", "signature": "sha256:efc4a659c9b37c524bbda8c39015a003db4fc61780a046b217d0b94def2cb08a"}, "nbformat_minor": 0}