{
 "metadata": {
  "name": "SpaceShuttleBayesFactor"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Bayes Factor for Log Regression/"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#####Is there really a relationship between failure and temperature?\n",
      "\n",
      "An critism of our above analysis is that *assumed* that the relationship followed a logistic model, this we implictly assumed that the probabilities change over temperature. Let's look at the data again. (Top figure)\n",
      "\n",
      "Could it be that infact this particular sequence of defects occured by chance?  After all, I can produce similar plots. (Bottom figure) This might explain the large overlap in temperatures."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figsize( 12.5, 6)\n",
      "subplot(211)\n",
      "plt.scatter( challenger_data[:,0], challenger_data[:,1], s = 75, \\\n",
      "    color=\"k\", alpha = 0.5) \n",
      "plt.yticks([0,1])\n",
      "plt.ylabel(\"Damage Incident?\")\n",
      "plt.title(\"(Real) Defects of the Space Shuttle O-Rings vs temperature\")\n",
      "\n",
      "subplot(212)\n",
      "n = challenger_data.shape[0]\n",
      "plt.scatter( challenger_data[:,0], stats.bernoulli.rvs(0.6, size = n) ,\n",
      "     s = 75, color=\"k\", alpha = 0.5) \n",
      "plt.yticks([0,1])\n",
      "plt.ylabel(\"Damage Incident?\")\n",
      "plt.xlabel(\"Outside temperature (Farhenhit)\" )\n",
      "plt.title(\"(Artificial) Defects of the Space Shuttle O-Rings vs \\\n",
      "temperature\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'figsize' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-1-a4aeded293db>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mfigsize\u001b[0m\u001b[1;33m(\u001b[0m \u001b[1;36m12.5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m6\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0msubplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m211\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m plt.scatter( challenger_data[:,0], challenger_data[:,1], s = 75, \\\n\u001b[0;32m      4\u001b[0m     color=\"k\", alpha = 0.5) \n\u001b[0;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0myticks\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNameError\u001b[0m: name 'figsize' is not defined"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Introducing the Bayes factor \n",
      "\n",
      "The *Bayes factor* is a measure comparing two models. In our example, on the one hand we believe that temperature plays an important role in the probability of O-ring failures. On the other hand, we believe that there is no significant connection, and the pattern appears by coincidence. We can compare these model using the ratio of the *probabilities of observing the data, given the model*:\n",
      "\n",
      "$$ \\frac{ P( \\text{observations} | \\text{Model}_1 ) }{  P( \\text{observations} | \\text{Model}_2 ) }$$\n",
      "\n",
      "\n",
      "Calculating this is not at all obvious. For simplicity, let's select a set of parameters from the posterior distribution (which is tantamount to selecting a particular model)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alpha_hat = alpha_samples[0,0] #select the first alpha sample\n",
      "beta_hat = beta_samples[0,0] #select the first beta sample\n",
      "\n",
      "p_hat = logistic( temperature, beta_hat, alpha_hat)\n",
      "print \"estimates of probability at observed temperature, defects: \"\n",
      "print np.array(zip(p_hat, temperature, D) )[:3, :]\n",
      "print \"...\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'alpha_samples' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-2-bf332447957b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0malpha_hat\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0malpha_samples\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;31m#select the first alpha sample\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mbeta_hat\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbeta_samples\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;31m#select the first beta sample\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mp_hat\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlogistic\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mtemperature\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbeta_hat\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha_hat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"estimates of probability at observed temperature, defects: \"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNameError\u001b[0m: name 'alpha_samples' is not defined"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To calculate the numerator in the Bayes factor, we start by **assuming a model**, in our case assume `alpha_hat`, `beta_hat` are true, and calculate the probability of the defects we observe, equal to the product of:\n",
      "\n",
      "$$ P(\\; \\text{Ber}(\\; p(t_i \\; | \\; \\hat{\\beta}, \\hat{\\alpha} )\\; ) = D_i \\; ),\\;\\; i=1..N $$\n",
      "\n",
      "For example, using the output above, the first computation, $i=0$, would look like \n",
      "\n",
      "\\begin{align}\n",
      "& p = p(t_1 \\; | \\; \\hat{\\beta}, \\hat{\\alpha} ) = 0.667 \\\\\\\\\n",
      "& d = D_0 = 0 \\\\\\\\\n",
      "& \\Rightarrow \\; P( \\; \\text{Ber}(p) = d ) = ?? \n",
      "\\end{align}\n",
      "\n",
      "The probability of this ocurring is $(1-0.667) = 0.333\\; $ (Recall the definition of a Bernoulli random variable $\\text{Ber}(p)$: it is equal to $1$ with probability $p$ and equal to 0 with probability $1-p$). As each observation is independent, we multiply all the observations together. A clever way to do this for a specific $i$ is, recalling the $D_i$ can only be 1 or 0:\n",
      "\n",
      "$$\\left( p(t_1 \\; | \\; \\hat{\\beta}, \\hat{\\alpha} )\\right)^{D_i} \\left( 1 - p(t_1 \\; | \\; \\hat{\\beta}, \\hat{\\alpha} ) \\right)^{(1-D_i)}$$\n",
      "\n",
      "\n",
      "(it is possible that the quantity can overflow, so it is recommended to take the $\\log$ of the above::\n",
      "\n",
      "$$ D_i\\log( p(t_1 \\; | \\; \\hat{\\beta}, \\hat{\\alpha} ) ) + (1-D_i)\\log( 1- p(t_1 \\; | \\; \\hat{\\beta}, \\hat{\\alpha} ) ) $$\n",
      "\n",
      "and sum the terms instead of multiplying. Be sure to use this transform for both models you are comparing.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print p_hat\n",
      "_product = p_hat**( D )*(1-p_hat)**(1-D)\n",
      "prob_given_model_1 = _product.prod()\n",
      "print \"probability of observations, model 1: %.10f\"%prob_given_model_1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'p_hat' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-3-be149a40cb07>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mprint\u001b[0m \u001b[0mp_hat\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0m_product\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp_hat\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mD\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mp_hat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mprob_given_model_1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_product\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"probability of observations, model 1: %.10f\"\u001b[0m\u001b[1;33m%\u001b[0m\u001b[0mprob_given_model_1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNameError\u001b[0m: name 'p_hat' is not defined"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The second model we are comparing against is that $\\beta=0$, i.e. there is no relationship between probabilites and temperature. We perform the same calculations as above. Notice that `beta=0` here in the below PyMC code:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "beta = 0\n",
      "alpha = pm.Normal( \"alpha\", 0, 0.001, value = 0 )\n",
      "\n",
      "@pm.deterministic\n",
      "def p( temp = temperature, alpha = alpha, beta = beta):\n",
      "    return 1.0/( 1. + np.exp( beta*temperature + alpha) ) \n",
      "\n",
      "\n",
      "observed = pm.Bernoulli( \"bernoulli_obs\", p, \\\n",
      "                        value = D, observed=True)\n",
      "\n",
      "model = pm.Model( {\"observed\":observed, \"beta\":beta, \"alpha\":alpha} )\n",
      "\n",
      "#mysterious code to be explained in Chapter 3\n",
      "map_ = pm.MAP(model)\n",
      "map_.fit()\n",
      "mcmc = pm.MCMC( model )\n",
      "mcmc.sample( 260000, 220000, 3 )\n",
      "######\n",
      "\n",
      "alpha_samples_model_2 = mcmc.trace(\"alpha\")[:, None]\n",
      "alpha_hat = alpha_samples_model_2[0] #use the first 'model'\n",
      "beta_hat = 0\n",
      "\n",
      "p_hat = logistic( temperature, beta_hat, alpha_hat )\n",
      "print \"estimates of probability at observed temperature, defects: \"\n",
      "print np.array(zip(p_hat, temperature, D ) )[:3, :]\n",
      "print\n",
      "print \"Notice the probability is constant for all temperatures?\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'mc' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-4-2fd082360270>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mbeta\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0malpha\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mNormal\u001b[0m\u001b[1;33m(\u001b[0m \u001b[1;34m\"alpha\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0.001\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvalue\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;33m@\u001b[0m\u001b[0mmc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdeterministic\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mtemp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtemperature\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0malpha\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbeta\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbeta\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNameError\u001b[0m: name 'mc' is not defined"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#compute the probability of observations given the model_2\n",
      "\n",
      "_product = p_hat**( D )*(1-p_hat)**(1-D)\n",
      "prob_given_model_2 = _product.prod()\n",
      "print \"probability of observations, model 2: %.10f\"%prob_given_model_2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'p_hat' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-5-de7da8dccc70>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m#compute the probability of observations given the model_2\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0m_product\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp_hat\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mD\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mp_hat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[0mprob_given_model_2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_product\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"probability of observations, model 2: %.10f\"\u001b[0m\u001b[1;33m%\u001b[0m\u001b[0mprob_given_model_2\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNameError\u001b[0m: name 'p_hat' is not defined"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Bayes factor = %.3f\"%(prob_given_model_1/prob_given_model_2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'prob_given_model_1' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-7-c40e6307b245>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"Bayes factor = %.3f\"\u001b[0m\u001b[1;33m%\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mprob_given_model_1\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mprob_given_model_2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[1;31mNameError\u001b[0m: name 'prob_given_model_1' is not defined"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Is this good? Below is a chart that can be used to compare the computed Bayes factors to a degree of confidence in M1. \n",
      "\n",
      "<table><tbody><tr><th>Bayes factor</th><th>Supporting M1</th></tr><tr><td>&lt;1:1  </td><td>Negative (supports M2)</td></tr><tr><td>  1:1 to 3:1 </td><td>Barely worth mentioning</td></tr><tr><td>  3:1 to 10:1 </td><td>Substantial</td></tr><tr><td>  10:1 to 30:1 </td><td>Strong</td></tr><tr><td>  30:1 to 100:1</td><td> Very strong</td></tr><tr><td>  &gt; 100:1</td><td> Decisive</td></tr></tbody></table>\n",
      "\n",
      "We are not done yet. If you recall, we selected only one model, but recall we actually have a distribution of possible models (sequential pairs of ($\\beta, \\alpha$) from the posterior distributions). So to be correct we need to average over samples from the posterior:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p = logistic( temperature[None,:], beta_samples, alpha_samples )\n",
      "_product = p**( D )*(1-p)**(1-D)\n",
      "prob_given_model_1 = _product.prod(axis=1).mean()\n",
      "print \"expected prob. of obs., given model 1: %.10f\"%prob_given_model_1\n",
      "\n",
      "    \n",
      "p_model_2 = logistic( temperature[:,None], \n",
      "                np.zeros_like(temperature), alpha_samples_model_2)\n",
      "\n",
      "_product = p_model_2**( D )*(1-p_model_2)**(1-D)\n",
      "prob_given_model_2 = _product.prod(axis=1).mean()\n",
      "print \"expected prob. of obs., given model 2: %.10f\"%prob_given_model_2\n",
      "print \n",
      "print \"Bayes factor: %.3f\"%(prob_given_model_1/prob_given_model_2)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'logistic' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-9-5e4bf3a1b066>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlogistic\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mtemperature\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbeta_samples\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha_samples\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0m_product\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mD\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mprob_given_model_1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_product\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"expected prob. of obs., given model 1: %.10f\"\u001b[0m\u001b[1;33m%\u001b[0m\u001b[0mprob_given_model_1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNameError\u001b[0m: name 'logistic' is not defined"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It looks we have a pretty good model, and temperature *is* significant. Look at you, you're a rocket scientist now."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}