HW2Solns.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Homework 2 Solutions"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Question 1 (a)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The function we are looking to integrate is:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$I=\\int_{0}^{\\pi/2} \\theta^2 {\\rm sin}\\theta d\\theta$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We approximate this integral using the standard Monte Carlo method and using $N = 1000$ samples. We repeat the approximation $N_{MC} = 100$ times to determine the variance of the approximation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "## Function to integrate\n",
      "fun = lambda x: x**2*np.sin(x)\n",
      "thetamin=0.\n",
      "thetamax=np.pi/2.\n",
      "v=thetamax-thetamin\n",
      "N=1e3\n",
      "Nmc=1e2\n",
      "\n",
      "Imc = np.zeros(Nmc)\n",
      "for i in np.arange(0,Nmc):\n",
      "    theta=v*np.random.uniform(low = 0, high=1,size=N)+thetamin\n",
      "    f=fun(theta)\n",
      "    Imc[i]=v*np.mean(f)\n",
      "\n",
      "outmean=np.mean(Imc)\n",
      "outvar=np.var(Imc)\n",
      "print outmean, outvar"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.14011725903 0.00133217082033\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We actually know how to integrate this function, so we can check our result here."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$I=\\int_{0}^{\\pi/2} \\theta^2 {\\rm sin}\\theta d\\theta = 2 \\theta {\\rm sin}\\theta - (\\theta^2 -2){\\rm cos}\\theta |_{0}^{\\pi/2}$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fint = lambda x: 2*x*np.sin(x)-(x**2-2)*np.cos(x)\n",
      "Iexact = fint(thetamax) - fint(thetamin)\n",
      "print Iexact"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.14159265359\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Question 1 (b)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next, we use one control variate to reduce the variance of the estimation (and thus increase the quality of the approximation). We use $g(\\theta) = \\theta^2$. Indeed, it seems that the variance of $f(\\theta)\u2212g(\\theta)$ is reduced with respect to $f(\\theta)$. We now apply the method of control variates and find:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gun= lambda x: x**2 #Control variate\n",
      "gint= lambda x: x**3/3 #integral of CV\n",
      "tau=(1/v)*(gint(thetamax)-gint(thetamin))\n",
      "\n",
      "I2 = np.zeros(Nmc);\n",
      "for i in np.arange(0,Nmc):\n",
      "    theta=(thetamax-thetamin)*np.random.uniform(low = 0, high=1,size=N)+thetamin;\n",
      "    f=fun(theta)\n",
      "    #Control variate\n",
      "    g=gun(theta)\n",
      "    q=np.cov(f,g)\n",
      "    b= - q[0,1]/np.var(g)\n",
      "    I2[i]=v*(np.mean(f)+b*(np.mean(g)-tau))\n",
      "\n",
      "mean2=np.mean(I2) #Mean of MC with control variate\n",
      "var2=np.var(I2) #Variance of MC with control variate\n",
      "print mean2, var2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.14139100704 9.70330598283e-06\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(Imc,10, histtype='step', label=u'Basic MC')\n",
      "plt.hist(I2,5, histtype='step', label = u'Control Variate')\n",
      "plt.plot([Iexact, Iexact], [0 ,50], 'y.-') \n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<matplotlib.legend.Legend at 0x106706b10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHkJJREFUeJzt3XtUVOXeB/Dv5hYqCAPIIIhOB60kL7CiPFrIJJFliph5\nQQXM6uhZJ3vVTmKZItapscxT65xlr3mo8J7dFBXMKMZeMQ+to2YnU8vyAsHIHVGR2/P+QU4ilxmG\nucDj97PWLJ299+znN8/a853NM/uiCCEEiIhIGk6OLoCIiKyLwU5EJBkGOxGRZBjsRESSYbATEUmG\nwU5EJBkXcxbSaDTo3bs3nJ2d4erqiry8PJSVlWHatGk4e/YsNBoNtm/fDm9vb1vXS0REJpi1x64o\nCvR6PY4cOYK8vDwAgE6nQ0xMDE6dOoXo6GjodDqbFkpEROYxeyjmxvOYMjIykJSUBABISkrCjh07\nrFsZERFZxOw99gceeAARERFYv349AMBgMECtVgMA1Go1DAaD7aokIiKzmTXGnpubi759+6K4uBgx\nMTG44447ms1XFAWKotikQCIi6hizgr1v374AgD59+mDSpEnIy8uDWq1GUVERAgICUFhYCH9//xav\nGzhwIE6fPm3diomIJBcSEoKffvrJ4tebHIq5fPkyLl68CAC4dOkS9u3bh6FDhyI2Nhbp6ekAgPT0\ndMTFxbV47enTpyGE4MMKj5SUFIfX0NUeR48+jJwc4JtvIlBbW87+dOCD/WndR2d3iE3usRsMBkya\nNAkAUF9fj5kzZ+LBBx9EREQEpk6dirS0NOPhjkT2FBq6Bbm5Kgwf/jlcXXmoLdE1JoP91ltvxdGj\nR1tM9/HxQXZ2tk2KIjLHtTBnqBM1xzNPuwmtVuvoEqTC/rQu9mfXogghbHajDUVRYMPVE0GvV6DV\nchsjuXQ2O806KoaIOs/Hxwfl5eWOLoO6EJVKhbKyMquvl8FOZCfl5eX8C5aasdX5PxxjJyKSDIOd\niEgyDHYiIskw2InIZsaNG4eNGzc6uoybDoOdyIF8fABFsd3Dx8e8OjQaDXr27AlPT0/4+Phg/Pjx\nyM/P7/T7y8zMREJCQodf5+TkBLVajYaGBuO0uro6+Pv7w8mpeWx99tlnGD16NHr37g1/f39otVrs\n2rWr07V3Zwx2IgcqLweEsN3D3KMrFUXB7t27cfHiRRQWFkKtVmP+/Pm2ffMm+Pj4ICsry/g8KysL\nPj4+zY4k+eijjzB16lTMnj0bBQUFuHDhAlauXMlgd3QBRNS13HLLLZg8eTKOHz9unLZnzx6Eh4fD\ny8sL/fv3R2pqqnFeTU0NZs2aBT8/P6hUKtxzzz0oLi4G0HRGalpamnHZ9evXIzQ0FL1798add96J\nI0eOtFlHQkICNmzYYHy+YcMGJCYmGg8ZFUJg0aJFWL58OebMmQNPT08AwOjRo/HOO+9YpzO6KQY7\nEQH4/S5ply9fxgcffICRI0ca53l4eGDTpk2orKzEnj178Pbbb2Pnzp0Amq7uWlVVhfz8fJSVlWHd\nunVwd3cH0PxeDR9++CFSU1OxceNGVFVVYdeuXfD19W2znokTJ+Krr75CVVUVysvLceDAAUycONE4\n/+TJk8jPz8djjz1m9b7o7niCEhFBCIG4uDi4uLjg0qVL8Pf3x969e43zo6KijP8fOnQopk+fjv37\n92PixIlwc3NDaWkpfvzxRwwdOhTh4eGttvGvf/0LycnJuOuuuwAAf/jDH9qtyd3dHRMmTMC2bdvQ\n2NiIiRMnGr8wAKC0tBTA7/eLoN9xj52IoCgKdu7cifLycly9ehX/+Mc/EBUVZbzl5b///W/cf//9\n8Pf3h7e3N9atW2cM1oSEBIwdOxbTp09HUFAQkpOTUV9f36KN/Px8hISEdKimxMREpKenY+PGjc2G\nYQAY9/YLCws789alxGAnomYURcGkSZPg7OyM3NxcAMCMGTMQFxeH/Px8VFRUYN68eWhsbAQAuLi4\nYPny5fj+++9x8OBB7N69u9nY+DXBwcEdvitQZGQkioqKcOHCBdx7773N5t1+++0IDg7GRx99ZOE7\nlReDnYgAoNmPktf23gcPHgwAqK6uhkqlgpubG/Ly8rBlyxbj2Ller8d3332HhoYGeHp6wtXVFc7O\nzi3W/+STT2L16tU4fPgwhBD46aefcO7cOZN17dq1CxkZGS2mK4qCNWvW4KWXXsL777+PqqoqNDY2\n4sCBA5g7d25nuqLb4xg7kQOpVE3Hm9ty/eaaMGECnJ2doSgKNBoNNmzYYAz2tWvX4tlnn8XTTz+N\nqKgoTJs2DRUVFQCAoqIizJs3D/n5+fDw8MD06dNbPXb9scceQ2lpKWbMmIGCggLceuut2LhxI/r3\n799i2esPaQwNDW1z3uTJk+Hh4YG//e1vmD9/Pnr06IEhQ4bgueeeM/+NS4jXY6durTtdj52fB7pR\nW9tEZ7cVDsUQEUmGwU5EJBkGOxGRZBjsRESSYbATEUmGwU5EJBkGOxGRZBjsRESSYbBTt1LfWI+n\nMp5C/MfxiP84HgAQ/3E8Dpw74ODKyJY0Gg2++OILu7Qlw+38eEkB6laqa6ux+bvNSIv97eYNpdtw\nqfYSDpw7gPv63+fY4izgs8oH5TVm3ubIAip3FcqSy8xadsuWLVizZg1OnjwJT09PhIWFYenSpS0u\nvtVRs2fPRnBwMF566SWL13H9dd2vp9PpkJWVhf379zebXlJSgsDAQBw9erTFJQlMyczMNHtZjUaD\nd999F2PGjOlQG7bGYKdux83ZDfFDm/bW9foZCO3TsQ9uV1JeUw6RYrvLDCip5l2IZs2aNVi1ahXW\nrVuHsWPHws3NDXv37kVGRkang92U+vp6uLhYFkUJCQlYtmwZzpw5A41GY5y+bds2DB8+vEOhfu0U\n/ta+QNrSVS8TwaEYoptcZWUlUlJSsHbtWsTFxaFHjx5wdnbGI488glWrVgEArl69igULFiAoKAhB\nQUFYuHAhamtrATRd3bFfv35Ys2YN1Go1AgMD8f777wMA3nnnHWzZsgWvvfYaPD09jXdA0mg0eO21\n1zBs2DB4enqioaEBGRkZuPPOO6FSqXD//ffjxIkTJmsPCgrCmDFjWgydXLuNXkVFBcaPHw9/f3/4\n+PhgwoQJKCgoMC6n1Wrx4osv4t5774WHhwd+/vnnZrfzO336NMaMGQM/Pz/06dMHs2bNQmVlJYCm\nL5Vz585hwoQJ8PT0xOrVqwEAhw4dwqhRo6BSqRAWFtbirwm7EDZk49XTTaj8SrnwetXL+DwnByL5\n82Tx6v+96sCqzNPa5wErbPsZMWf9WVlZwsXFRTQ0NLS5zLJly8TIkSNFcXGxKC4uFqNGjRLLli0T\nQgiRk5MjXFxcREpKiqivrxeZmZmiZ8+eoqKiQgghxOzZs43LXjNgwAARHh4u8vPzRU1NjTh58qTo\n1auXyM7OFvX19eK1114TAwcOFHV1dUIIITQajfjiiy9arW3z5s1i0KBBxucnTpwQbm5uoqSkRJSW\nlopPPvlEXLlyRVy8eFFMmTJFxMXFGZeNiooSAwYMEMePHxcNDQ2irq5OaLVakZaWJoQQ4qeffhLZ\n2dmitrZWFBcXi9GjR4sFCxYYX39jXfn5+cLX11dkZWUJIYT4/PPPha+vryguLm619rYysrPZyT12\noptcaWkp/Pz84OTUdhxs2bIFy5cvh5+fH/z8/JCSktJsL9nV1RXLly+Hs7MzHn74YXh4eODkyZPG\n+eKG4QpFUfDMM88gKCgIt9xyCz744AOMHz8e0dHRcHZ2xl//+ldcuXIFBw8eNFl/XFwcDAYDvv76\nawBNe+vjxo2Dr68vfHx8MGnSJLi7u8PDwwMvvPBCsz1oRVEwe/ZsDB48GE5OTi2GhEJCQhAdHQ1X\nV1f4+flh4cKF7e6Bb9q0CePGjcNDDz0EAHjggQcQERHRoXF7a2CwE93kfH19UVJSYrwjUmt+/fVX\nDBgwwPi8f//++PXXX5ut4/ovhp49e6K6urrddoODg43/LywsbHZddkVREBwc3GzYpC09e/bElClT\njHdt2rx5MxITEwE03Zh77ty50Gg08PLyQlRUFCorK5t90Vxfx40MBgOmT5+Ofv36wcvLCwkJCcZb\nArbm7Nmz+PDDD6FSqYyP3NxcFBUVmXwf1sRgJ7rJjRw5Erfccgs+/fTTNpcJDAzEmTNnjM/PnTuH\nwMBAs9bf1o+R108PDAzE2bNnjc+FEDh//jyCgoLMaiMpKQnbt2/Hvn37UF1djQkTJgAA3njjDZw6\ndQp5eXmorKzE/v37IYRoFuzt/Vj6wgsvwNnZGf/9739RWVmJjRs3NvsCvPG1/fv3R0JCAsrLy42P\nixcvYvHixWa9D2thsBPd5Ly8vLBy5Ur85S9/wc6dO3H58mXU1dUhKysLycnJAID4+Hi8/PLLKCkp\nQUlJCVauXNnqXZJao1ar8fPPP7e7zNSpU7Fnzx58+eWXqKurwxtvvAF3d3eMGjXKrDYiIyPh7e2N\nuXPnIj4+3jikUl1djR49esDLywtlZWVITU1t8dobh4muV11djV69eqF3794oKCjA66+/3uK9nT59\n2vh81qxZ2LVrF/bt24eGhgbU1NRAr9eb9ZeHNfFwRyIHUrmrzD4k0dL1m2PRokUICAjAyy+/jJkz\nZ8LT0xMRERFYunQpAODFF19EVVUVhg0bBqApiF988UXj69vb633iiScwZcoU49Eun3zySYtlbrvt\nNmzatAnz589HQUEBwsPDsWvXrg4dBpmYmIiVK1cah2EAYMGCBZgxYwb8/PwQFBSERYsWtbh/anu1\np6SkIDExEV5eXhg0aBBmzZqFN9980zj/+eefx/z587F48WIsW7YMixYtws6dO7F48WLEx8fD2dkZ\nI0aMwNq1a81+H9Zg1q3xGhoaEBERgX79+mHXrl0oKyvDtGnTcPbsWWg0Gmzfvh3e3t4tV95Fj/Gk\n7quipgKaNzWoWNJ0v029XsHe+mR4u3tjyX1LHFxd+/h5oBs59NZ4b731FkJDQ43fbDqdDjExMTh1\n6hSio6Oh0+ksLoCIiKzLZLDn5+cjMzMTTz75pPEbJCMjA0lJSQCafrTYsWOHbaskIiKzmQz2hQsX\n4vXXX292KJPBYIBarQbQ9OOBwWCwXYVERNQh7f4ysXv3bvj7+yM8PBx6vb7VZdq6OM81K1asMP5f\nq9VCq9VaUicRkbT0en2bGWuJdoP94MGDyMjIQGZmJmpqalBVVYWEhASo1WoUFRUhICAAhYWF8Pf3\nb3Md1wc7ERG1dONOb2uHZXZEu0Mxr7zyCs6fP49ffvkF27ZtM15sJzY2Funp6QCA9PR0xMXFdaoI\nIiKyng4dx35tyGXJkiWYOnUq0tLSjIc7ElH7VCpVhy4JS/JTqcw7z6CjzA72qKgoREVFAQB8fHyQ\nnZ1tk4KIZFVWZt4NL4g6i5cUICKSDIOdiEgyDHYiIskw2ImIJMNgJyKSDIOdiEgyDHYiIskw2ImI\nJMNgJyKSDIOdiEgyDHYiIskw2ImIJMNgJyKSDIOdiEgyDHYiIskw2ImIJMNgJyKSDIOdiEgyDHYi\nIskw2ImIJMNgJyKSDIOdiEgyDHYiIskw2ImIJMNgJyKSDIOdiEgyDHYiIskw2ImIJMNgJyKSDIOd\niEgyDHYiIskw2ImIJMNgJyKSDIOdiEgyDHYiIsm0G+w1NTUYMWIEwsLCEBoaiueffx4AUFZWhpiY\nGNx222148MEHUVFRYZdiiYjItHaD3d3dHTk5OTh69CiOHTuGnJwcHDhwADqdDjExMTh16hSio6Oh\n0+nsVS8REZlgciimZ8+eAIDa2lo0NDRApVIhIyMDSUlJAICkpCTs2LHDtlUSEZHZTAZ7Y2MjwsLC\noFarcf/99+POO++EwWCAWq0GAKjVahgMBpsXSkRE5nExtYCTkxOOHj2KyspKjB07Fjk5Oc3mK4oC\nRVFsViAREXWMyWC/xsvLC4888gj+85//QK1Wo6ioCAEBASgsLIS/v3+br1uxYoXx/1qtFlqttjP1\nEhFJR6/XQ6/XW219ihBCtDWzpKQELi4u8Pb2xpUrVzB27FikpKTgs88+g6+vL5KTk6HT6VBRUdHq\nD6iKoqCd1RN1WEVNBTRvalCxpOlILL1ewd76ZHi7e2PJfUscXB2RdXQ2O9vdYy8sLERSUhIaGxvR\n2NiIhIQEREdHIzw8HFOnTkVaWho0Gg22b99ucQFERGRd7Qb70KFDcfjw4RbTfXx8kJ2dbbOiiIjI\ncjzzlIhIMgx2IiLJMNiJiCTDYCcikgyDnYhIMgx2IiLJMNiJiCTDYCcikgyDnYhIMgx2IiLJMNiJ\niCTDYCcikgyDnYhIMgx2IiLJMNiJiCTDYCcikgyDnYhIMgx2IiLJMNiJiCTDYCcikgyDnYhIMgx2\nIiLJMNiJiCTDYCcikgyDnYhIMgx2IiLJMNiJiCTDYCcikgyDnYhIMgx2IiLJMNiJiCTDYCcikgyD\nnYhIMgx2IiLJMNiJiCTj4ugCiNpzpPAIcs7kGJ9fqbviwGoAvR44fNi+bXp5AXPmAIpi33ap+zIZ\n7OfPn0diYiIuXLgARVHwpz/9Cc888wzKysowbdo0nD17FhqNBtu3b4e3t7c9aqabyLr/rMOPZT9i\nuHq4cdqy0cscVo9OB/ToAdx6q/3afO45YOZMwN3dfm1S92Yy2F1dXfH3v/8dYWFhqK6uxl133YWY\nmBi89957iImJweLFi7Fq1SrodDrodDp71Ew3mSmhUzAvYp6jyzCaOxd46CH7tbd2rf3aIjmYHGMP\nCAhAWFgYAMDDwwODBw9GQUEBMjIykJSUBABISkrCjh07bFspERGZpUM/np45cwZHjhzBiBEjYDAY\noFarAQBqtRoGg8EmBRIRUceYHezV1dWYPHky3nrrLXh6ejabpygKFP6yQ0TUJZh1VExdXR0mT56M\nhIQExMXFAWjaSy8qKkJAQAAKCwvh7+/f6mtXrFhh/L9Wq4VWq+100UREMtHr9dDr9VZbn8lgF0Lg\niSeeQGhoKBYsWGCcHhsbi/T0dCQnJyM9Pd0Y+De6PtiJiKilG3d6U1NTO7U+k8Gem5uLTZs2Ydiw\nYQgPDwcAvPrqq1iyZAmmTp2KtLQ04+GORETkeCaD/b777kNjY2Or87Kzs61eEJElPjv9Gaprq43P\nfXr4YNHIRQ6siMhxeEkB6vamD5mO6Fuj4e7ibnw8u+9ZR5dF5DC8pAB1e2EBYQgLCGs2bVmO485O\nJXI07rETEUmGwU5EJBkGOxGRZBjsRESSYbATEUmGwU5EJBkGOxGRZBjsRESSYbATEUmGwU5EJBle\nUoCsQgigoACor7fuei9eBEpdgTNnrLteIpkx2Mkqjh4F/vhHoG9f6663dCSwtwxYf7LlvMpK4NNP\nrdsekQwY7GQVV68C4eHAoUPWXe+83UBYADAvouU8nc66bRHJgmPsRESSYbATEUmGwU5EJBkGOxGR\nZBjsRESSYbATEUmGwU5EJBkGOxGRZBjsRESSYbATEUmGwU5EJBleK0ZSCQlAUZH92quoAFxd7dce\nADj9tlsSE9PKzPvamN5Jhw//3i5RV8Vgl9SHHwIffwy4udmvzQED7NcWADz1FPDtt8DixS3nZR9s\nfXpnOTkBkZHWXy+RNTHYJRYdDbi7O7oK21Gpmv5tdc/8oG322Im6A/5RSUQkGQY7EZFkGOxERJJh\nsBMRSYbBTkQkGQY7EZFkeLgjURfn4gL88Y/2PTFKqwXWrLFfe2RdJoN9zpw52LNnD/z9/fHdd98B\nAMrKyjBt2jScPXsWGo0G27dvh7e3t82LJboZffstUFlpv/aOHQP+93/t1x5Zn8lgf/zxxzF//nwk\nJiYap+l0OsTExGDx4sVYtWoVdDoddDqdTQslulmFhNi3vdpa+7ZH1mfyj7vIyEiorp3i95uMjAwk\nJSUBAJKSkrBjxw7bVEdERB1m0Ri7wWCAWq0GAKjVahgMBqsWJZtjx4AffrBvmw0N9m2PiLqOTv94\nqigKFEVpc/6KFSuM/9dqtdBqtZ1tsttJTgYuXQL69rVfm48/bt8LgBGR5fR6PfR6vdXWZ1Gwq9Vq\nFBUVISAgAIWFhfD3929z2euD/WYlBPDCC8BDDzm6EiLqim7c6U1NTe3U+iw6gCo2Nhbp6ekAgPT0\ndMTFxXWqCCIish6TwR4fH49Ro0bh5MmTCA4OxnvvvYclS5bg888/x2233YYvv/wSS5YssUetRERk\nBpNDMVu3bm11enZ2ttWLIbIWN2c3aN7UNJu2NHIpnrrrKccURGRHPPOUpJS/MB+X6i4Zn7/9zds4\nXX7agRUR2Q+DnaTUp1cf9EEf43NVDxUqaiocWBGR/fAiYEREkuEeO3UZDY0NSP82HVfrrxqnfV/8\nPcICwhxYFVH3w2CnLuN81Xn8z97/wayhs4zThvQZgsj+kQ6siqj7YbBTl+Lbwxdvj3/b0WUQdWsc\nYycikgyDnYhIMgx2IiLJMNiJiCTDYCcikgyDnYhIMjfd4Y4nTgAffWTfNk/zEiVEZEc3XbDv2AFk\nZgJjxtivzfh44K677NceEd3cbrpgB4DRo4GXX3Z0FUREtsExdiIiyTh8j72+vumeoPbS0GC/toiI\nHMGhwf7998CwYYCTnf9u0Ons2x4RkT05NNgvXgTuvhs4dMiRVRARyYVj7EREkmGwExFJhsFORCQZ\nBjsRkWQY7EREkmGwExFJxuEnKBFR11NbCxQV2bdNPz/AhYlkFexGImrG1xe4cAEIC7Nfm9XVwHPP\nASkp9mtTZgx2Impm0CAgP9++bep0QEWFfduUGcfYiYgkw2AnIpIMh2KIqEvYvRsoLLRfe05OwMqV\nQHCw/dq0FwY7ETlcQgLQt69921y9uukKswx2IiIbCAoCkpLs2+bWrfZtz544xk5EJBkGOxGRZDoV\n7Hv37sUdd9yBQYMGYdWqVdaqiYiIOsHiYG9oaMDTTz+NvXv34vjx49i6dSt++OEHa9ZG19Hr9Y4u\nQSrsT+tif3YtFgd7Xl4eBg4cCI1GA1dXV0yfPh07d+60Zm10HX5wrIv9aV3sz67F4mAvKChA8HXH\nCfXr1w8FBQVWKYqIiCxncbAritLpxnv1Am6/vdOrIUm4ObthiP8Qs5c/ceJPAIBvvx2Hurr2LzSi\n7qVGkGdQp+ojuQwcCPTu7egqbERY6OuvvxZjx441Pn/llVeETqdrtkxISIgAwAcffPDBRwceISEh\nlkazEEIIRQghYIH6+nrcfvvt+OKLLxAYGIh77rkHW7duxeDBgy1ZHRERWYnFZ566uLjgn//8J8aO\nHYuGhgY88cQTDHUioi7A4j12IiLqmiz68XTOnDlQq9UYOnRom8s888wzGDRoEIYPH44jR44Yp2s0\nGgwbNgzh4eG45557LGleKqb68sSJExg5ciTc3d3xxhtvNJvHE8Ra6kx/cttsyVR/bt68GcOHD8ew\nYcNw77334tixY8Z53D6b60xfdnjbtGRg/quvvhKHDx8WQ4YMaXX+nj17xMMPPyyEEOLQoUNixIgR\nxnkajUaUlpZa0qyUTPXlhQsXxDfffCOWLl0qVq9ebZxeX18vQkJCxC+//CJqa2vF8OHDxfHjx+1V\ndpdlaX8KwW2zNab68+DBg6KiokIIIURWVpbxs87tsyVL+1KIjm+bFu2xR0ZGQqVStTk/IyMDSb9d\nqm3EiBGoqKiAwWC4/svEkmalZKov+/Tpg4iICLi6ujabzhPEWmdpf17DbbM5U/05cuRIeHl5AWj6\nrOf/dk89bp8tWdqX13Rk27TJRcDaO3lJURQ88MADiIiIwPr1623R/E2BJ4hZH7fNzklLS8O4ceMA\ncPvsrOv7Euj4tmmz67G39e1y4MABBAYGori4GDExMbjjjjsQGRlpqzKkZY0TxKi53Nxc9O3bl9um\nBXJycvDuu+8iNzcXALfPzrixL4GOb5s22WMPCgrC+fPnjc/z8/MRFNR01l9gYCCApj+JJ02ahLy8\nPFuUIL0b+/j8+fPo16+fAyvq/vr+dgsfbpsdc+zYMTz11FPIyMgwDjVw+7RMa30JdHzbtEmwx8bG\nYsOGDQCAQ4cOwdvbG2q1GpcvX8bFixcBAJcuXcK+ffvaPbKGfnfjX0ARERH48ccfcebMGdTW1uKD\nDz5AbGysg6rrfm7sT26bljl37hweffRRbNq0CQMHDjRO5/bZcW31pSXbpkVDMfHx8di/fz9KSkoQ\nHByM1NRU1NXVAQDmzp2LcePGITMzEwMHDkSvXr3w3nvvAQCKiorw6KOPAmg6c3XmzJl48MEHLSlB\nGqb6sqioCHfffTeqqqrg5OSEt956C8ePH4eHhwdPEGuFpf154cIFbputMNWfK1euRHl5Of785z8D\nAFxdXZGXl8cTGFthaV9akps8QYmISDK8NR4RkWQY7EREkmGwExFJhsFORCQZBjsRkWQY7EREkmGw\nExFJhsFORCSZ/wdBGLOJwBogtgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106706a10>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Question 3"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bayes theorem states that $P (\\theta | D) = P (D | \\theta) \\frac{P(D)}{P(\\theta)} $ or expressed more clearly in the form of this problem:\n",
      "\n",
      "$P({\\rm healthy~|~PosTest}) = P({\\rm PosTest~|~healthy}) \\frac{P({\\rm healthy})}{P({\\rm PosTest})}$\n",
      "\n",
      "where PosTest is the situtation where a patient has tested positive for the disease. From the question we know that the probability of someone not having the disease (i.e, that they are healthy) is $P({\\rm healthy}) = 1 - 10^{-6} = 0.999999$. Also the probability that someone is healthy inspite of receiving a positive test result (i.e., a false positive) is $P({\\rm PosTest~|~healthy}) = 0.0001$. Also, we are missing from the Bayes equation is probabiliy that someone will test positive, $P({\\rm PosTest})$.\n",
      "\n",
      "There are only two possible outcomes, either a person is healthy or sick. So the marginal reduces to:\n",
      "\n",
      "$P({\\rm PosTest}) = P({\\rm PosTest~|~healthy})P({\\rm healthy}) + P({\\rm PosTest~|~sick})P({\\rm sick})$\n",
      "\n",
      "but we know form the question that $P({\\rm PosTest~|~sick}) =1 $ and $P({\\rm sick})=10^{-6}$. Therefore, \n",
      "\n",
      "$P({\\rm PosTest}) = (0.0001 \\times 0.999999) + (1 \\times 10^{-6}) = 0.000101$\n",
      "\n",
      "Finally,\n",
      "\n",
      "$P({\\rm healthy~|~PosTest}) = \\frac{0.0001 \\times 0.999999}{0.000101} = 0.990098$\n",
      "\n",
      "So we can clearly see that if you were to go and take the test, the chance of you still being healthy inspite of the test informing you that you may be sick is $>99\\%$. In other words, it's probably a waste of your valuable time that could be better spent in class!\n",
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Question 4"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Part (a): The priors should be based on our previous knowledge, but not specifically take into account the data. So, from the information given in the question, we know that Challenger failed at a relatively low temperature (31F). We also know from common knowledge that 31F is a relatively low temperature, which is unusual for the typical temperatures seen in Florida and that this was not the first shuttle launch, there were 50 before the STS-51 mission - these launches did not fail. With this intuition in mind, the regression gradient should probably be negative, i.e., it is less likely that the shuttle will fail at more typical Florida temperatures. In turn, this means that the intercept should probably be positive. So the parameters should have mean values that can be normally distributed. So, we choose extremely simple and loose priors for $\\Theta = (\\theta_1,\\theta_2)$ of $\\theta_1 \\sim N(31,30)$ and $\\theta_2 \\sim N(-0.5,1)$. Any sensible positive intercept for $\\theta_1$ would be fine, although less prior information could also be assumed if we were to uniformly distribute $\\theta_1$ above zero, and any mean negative gradient for $\\theta_2$ that allows for the more unlikely possibility of a positive gradient would also be acceptable."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Part (b) We use the data to define the likelihood (L) as:\n",
      "\n",
      "$L(D|\\Theta) = \\prod\\limits_{i=1}^{n=23} f(fail_i | \\theta_1,\\theta_2, T_i) \\times f(success_i|\\theta_1,\\theta_2,T_i)$\n",
      "    \n",
      "And implement a multivariate regection alogrithm to sample from the posterior given by:\n",
      "$P (\\Theta|D) = L(D|\\Theta)P (\\Theta)$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from scipy.stats import norm\n",
      "import pandas as pd\n",
      "\n",
      "# PART A: GENERATE THE PRIORS\n",
      "def prior(theta1, theta2):\n",
      "    a=norm.pdf(theta1, loc= 31, scale=30)\n",
      "    b=norm.pdf(theta2, loc=-0.5, scale=1)\n",
      "    c=a*b\n",
      "    return c\n",
      "\n",
      "# PART B: CALCULATE THE POSTERIOR\n",
      "data=np.loadtxt('./files/chall.txt')\n",
      "T=data[:,0]\n",
      "p_fail=data[:,1]\n",
      "\n",
      "\n",
      "def mylikelihood(pfail, temp, thet1, thet2):\n",
      "    logistic_fun= 1./(1 + np.exp(-1*(thet1 + thet2*temp)))\n",
      "    p=np.prod(logistic_fun**pfail*(1-logistic_fun)**(1-pfail))\n",
      "    return p\n",
      "\n",
      "#def mylikelihood(pfail, t, a, b):\n",
      "#    logistic_fun= lambda (temp,thet1,thet2): 1./(1 + np.exp(-1*(thet1 + thet2*temp)))\n",
      "#    p=np.prod(logistic_fun(t, a, b)**pfail*(1-logistic_fun(t, a, b))**(1-pfail))\n",
      "#    return p\n",
      "\n",
      "# Now lets use rejection to draw random samples, t1 and t2 from the posterior\n",
      "\n",
      "foo=200\n",
      "x = np.linspace(0,50,foo)\n",
      "y = np.linspace(-1,0.3,foo)\n",
      "xx,yy = np.meshgrid(x,y)\n",
      "zz=np.zeros((foo,foo))\n",
      "for i in np.arange(0,foo):\n",
      "    for j in np.arange(0,foo):\n",
      "        zz[i,j] = mylikelihood(p_fail, T, xx[i,j], yy[i,j])\n",
      "\n",
      "plt.contourf(xx,yy,zz)\n",
      "plt.xlabel('t_1')\n",
      "plt.ylabel('t_2')\n",
      "plt.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD7CAYAAACG50QgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFilJREFUeJzt3X9QVNfdx/HPWpjpTH+AJmFRlg6tgvwIAoohbWaMhqxU\nTQjVRE06iU0t+mTGpnY6Q9K0fziZUaE/0qZjZ2odkyFpHqOdNpFG5QlWwUaHMhUNbaKjSSQisLQp\nkMSY1Ejv84eAu+wP2LvI7nLerxkme++eyz2emXy+h7P33nVYlmUJAGCEKdHuAABg4hD6AGAQQh8A\nDELoA4BBCH0AMAihDwAGSYh2B4Y4HJmS3op2NwAgrhQUFOjkyZNjbu+Ilev0HQ6HpOZodyNG7JBU\nGe1OxAjG4hrG4hrG4ppbFU6Ms7wDAAYh9AHAIIR+TJob7Q7EEMbiGsbiGsbCLtb0ASCusaYPAAiC\n0AcAgxD6AGAQQh8ADELoA4BBCH0AMEjEoV9fX6/s7GxlZmaqpqbG7/0XXnhBBQUFmjNnjm677Ta1\ntbVFekoAgE0RPXBtYGBAGzZs0MGDB5WWlqb58+ervLxcOTk5w22+8pWv6MiRI0pKSlJ9fb3WrVun\n5mauxweAaIhopt/S0qJZs2YpIyNDiYmJWr16tfbu3evT5qtf/aqSkpIkSSUlJbpw4UIkpwQARCCi\n0O/s7FR6evrwtsvlUmdnZ9D2O3fu1NKlSyM5JQAgAhEt71x9dMLYHD58WM8884yOHj0aySkBABGI\nKPTT0tLU0dExvN3R0SGXy+XXrq2tTZWVlaqvr9fUqVND/MYdXq/nSpoXSfcAYBI6LqnV9tERPXDt\nypUrmj17tv785z9rxowZuuWWW7Rr1y6fD3LPnz+vO+64Q7/73e906623Bu8ID1wDABvCe+BaRDP9\nhIQEbdu2TWVlZRoYGNDatWuVk5Oj7du3S5LWr1+vJ598Un19fXrkkUckSYmJiWppaYnktAAAm3i0\nMgDENR6tDAAIgtAHAIMQ+gBgEEIfAAxC6AOAQQh9ADAIoQ8ABiH0AcAghD4AGITQBwCDEPoAYBBC\nHwAMQugDgEEIfQAwCKEPAAYh9AHAIIQ+ABiE0AcAgxD6AGAQQh8ADELoA4BBCH0AMAihDwAGIfQB\nwCCEPgAYhNAHAINEHPr19fXKzs5WZmamampqArZ59NFHlZmZqYKCAp04cSLSUwIAbIoo9AcGBrRh\nwwbV19frzTff1K5du3Tq1CmfNvv379dbb72ls2fP6re//a0eeeSRiDoMALAvotBvaWnRrFmzlJGR\nocTERK1evVp79+71aVNXV6c1a9ZIkkpKStTf36+enp5ITgsAsCmi0O/s7FR6evrwtsvlUmdn56ht\nLly4EMlpAQA2JURysMPhGFM7y7LGeNwOr9dzJc2z1S8AmLyOS2q1fXREoZ+WlqaOjo7h7Y6ODrlc\nrpBtLly4oLS0tCC/sTKS7gCAAebJd0K8M6yjI1reKS4u1tmzZ9Xe3q7Lly9r9+7dKi8v92lTXl6u\n5557TpLU3Nys5ORkOZ3OSE4LALApopl+QkKCtm3bprKyMg0MDGjt2rXKycnR9u3bJUnr16/X0qVL\ntX//fs2aNUuf+9zn9Oyzz45LxwEA4XNYIxfco+TqOn9ztLsBAHHmVr/PTUPhjlwAMAihDwAGIfQB\nwCCEPgAYhNAHAINEdMlmTFtUEnj/4b9ObD8AIIZM/ks2g4W/RAEAMAmEd8nm5A99KXTwD6EAAIhL\nhH5wYwl/iQIAII4Q+mMTqABsGvwZiSIAIGYR+uEZGf6bArQZ2kf4A4g5hL493uG/KUS7ofcoAABi\nAqEfuaECsGmUdkPvUwAARA2hP35GCf8lC/6oA0eW+7ahAACYUIT++AsQ/ksW/NGvGQUAwMQj9K+v\nEQUgUPgPOXBkOQUAwHVG6E+MEOH/P9qu32i9T3MKAIDrg9CfWKOEv7ehQkABADB+CP3o8SoAQ+E/\nMviHUAAAjA9CP/oCzP5Hhn/566+qrmCxpAAFgPAHMGaEfmwZMfsPFP5D6goW6zda73sVEAUAQEiE\nfuxaVOIX/t6h782nAGwa3EkBAOCH0I99AWb/wcJfogAACIXQjx9Bln7CKgCEP2A4Qj8+DS79SAFm\n/zVe7R679pICAGBCQ7+3t1erVq3Su+++q4yMDO3Zs0fJyck+bTo6OvTQQw/pn//8pxwOh9atW6dH\nH33UvyOmh/6QUEs/NQHaDxYBCgBgqgkM/aqqKt14442qqqpSTU2N+vr6VF1d7dPG4/HI4/GosLBQ\nFy9e1Lx58/Tyyy8rJyfHtyOEvr9gBSBQ+EvSY6z/A+aZwNDPzs5WU1OTnE6nPB6PFi5cqNOnT4c8\npqKiQt/97ndVWlrq2xFCP7hww1+iAADGmMDQnzp1qvr6+iRJlmVp2rRpw9uBtLe36/bbb9cbb7yh\nz3/+874dIfTHZgwF4Oiuq/+97f7BHYEKAOEPTBLjHPput1sej8dv/+bNm7VmzRqfkJ82bZp6e3sD\n/p6LFy9q4cKF+vGPf6yKigr/jjgcktZ67Zkrad6Y/hHGGnHdf6DZ/1ABkAaLwGABuOfI/zH7B+LS\ncUmtXts7J3Z5p7GxUampqeru7taiRYsCLu98+umnuuuuu7RkyRJt3LgxcEeY6dsXYvbvHfpDvP8C\nWFrwB2b/QFwLb6Y/JZJTlZeXq7a2VpJUW1sbcAZvWZbWrl2r3NzcoIGPCB3+69Wf2/+qA3cs1z1H\n/k9LC/6guv9drNtOeoX8oKO7Bn8Kpf0PrJCV5NDeQ2VaYnVKTSX+XxYPYNKI+JLNlStX6vz58z6X\nbHZ1damyslL79u3Ta6+9pgULFmjOnDmDs3lp69at+vrXv+7bEWb64yvM2b/ku/zD+j8QL7g5CyMF\nu/ErQAHw+avAe/1fogAAMYnQRzDjOfuXKABATCD0MRYRfPjL8g8QSwh9hCvEpZ9jKQBc/glEE6EP\nu/jwF4hDhD7GAwUAiBOEPsZbBMs/Pjd/SRQAYNwR+rhewpj9B7r0k9k/cD0Q+pgIYywAI+8GpgAA\n443Qx0QK9ZWPI5aARp39SxQAIGyEPqJljI99Hhn+EgUAsI/QRywI8eiHoEbO/iWWf4BREfqIJSNm\n/5LG9LWPErN/YGwIfcSqcL/0XeKL34FREfqIBzYLALN/YCRCH/HEK/ylEAXgscHXj107lAIASIQ+\n4leo9X/JtwB44cNfmI3Qx2QQ6vr/IOoKFksSs38YhtDHZBJk+UcavQD4zf4lCgAmIUIfk1WYBWBo\n5i8x+8dkRujDBEHW/4d4F4GR4S+JAoBJhNCHSYbCXwpaAKSrRSBo+A8eS/gjPhH6MFWA5R/JvwCM\nxOwf8Y3QB2wVAGb/iE+EPuAr0gIweBwFALGJ0AeCC7MABJz9SxQAxJDwQj/B7ml6e3u1atUqvfvu\nu8rIyNCePXuUnJwcsO3AwICKi4vlcrn0pz/9ye4pgcgNh/XV8D+gwZn8gmtNAhWAoeJwYNNyn+MJ\nf8Qb2zP9qqoq3XjjjaqqqlJNTY36+vpUXV0dsO1TTz2l48eP68MPP1RdXV3gjjDTR7QEmf2Phtk/\nYsMELe9kZ2erqalJTqdTHo9HCxcu1OnTp/3aXbhwQd/61rf0ox/9SE899VTQmT6hj5hgowAQ/oiu\nCQr9qVOnqq+vT5JkWZamTZs2vO3tvvvu0xNPPKEPPvhAP/vZzwh9xIcR1/9LoxeA4fD3OoYCgOtv\nHNf03W63PB6P3/7Nmzf7bDscjsHQ9vXKK68oJSVFRUVFamxsHHOngKjzCesR6/8KXAC897H2j1gV\n0fJOY2OjUlNT1d3drUWLFvkt7zzxxBN6/vnnlZCQoE8++UQffPCBVqxYoeeee86/Iw6HpLVee+ZK\nmmena8D1E+AvAClwEfCZ+Xu3pwAgIscltXpt75yY5Z2qqirdcMMNeuyxx1RdXa3+/v6gH+RKUlNT\nE8s7mFxGKQAHjiz3eR2oLQUAkZugNf3e3l6tXLlS58+f97lks6urS5WVldq3b59P+6amJv385z/n\n6h1MTiM+AB4T77aEP2zj5iwgeoLM/kfl3ZYCgLAQ+kBssFMAvNsR/hgTQh+IPeEWgJFtKAAIitAH\nYlskBYDwhx9CH4gf4RSAke9TACCJ0AfikXf4D9kUoF2gfYS/4Qh9IL4FKgDeNin4XwUUAAMR+sDk\nMVoBCIbwNwihD0xOdgoA4W8AQh+Y/CgAGEboA2YJtwAQ/pMMoQ+YiwJgIEIfgBReASD84xihD2Ck\nsRYAwj8OEfoAQhlLASD84wihD2CsKACTAKEPwI7RCgDhH6MIfQCRClUACP8YE17oJ1zHngCIVyOD\nfagIEPhxj9AHMDrCftKYEu0OAAAmDqEPAAYh9AHAIIQ+ABiE0AcAgxD6AGAQQh8ADGI79Ht7e+V2\nu5WVlaXFixerv78/YLv+/n7de++9ysnJUW5urpqbuesWAKLFduhXV1fL7XbrzJkzKi0tVXV1dcB2\n3/ve97R06VKdOnVKbW1tysnJsd1ZAEBkbD97Jzs7W01NTXI6nfJ4PFq4cKFOnz7t0+b9999XUVGR\n3nnnndE7wrN3AMCG8J69Y3um39PTI6fTKUlyOp3q6enxa3Pu3DnddNNNevjhhzV37lxVVlbq0qVL\ndk8JAIhQyGfvuN1ueTwev/2bN2/22XY4HIMzdV9XrlxRa2urtm3bpvnz52vjxo2qrq7Wk08+GeSM\nO7xez5U0b7T+A4BhjktqtX10yNBvaGgI+t7Qsk5qaqq6u7uVkpLi18blcsnlcmn+/PmSpHvvvTfo\n2v9VlWPrNQAYa558J8Q7wzra9vJOeXm5amtrJUm1tbWqqKjwa5Oamqr09HSdOXNGknTw4EHl5eXZ\nPSUAIEK2P8jt7e3VypUrdf78eWVkZGjPnj1KTk5WV1eXKisrtW/fPknS66+/ru985zu6fPmyZs6c\nqWeffVZJSUn+HeGDXACwgW/OAgCDTNDVOwCA+EPoA4BBCH0AMAihDwAGIfQBwCCEPgAYhNAHAIMQ\n+gBgEEIfAAxC6AOAQQh9ADAIoQ8ABiH0AcAghD4AGITQBwCDEPoAYBBCHwAMQugDgEEIfQAwCKEP\nAAYh9AHAIIQ+ABiE0AcAgxD6AGAQQh8ADGI79Ht7e+V2u5WVlaXFixerv78/YLutW7cqLy9P+fn5\neuCBB/Sf//zHdmcBAJGxHfrV1dVyu906c+aMSktLVV1d7demvb1dO3bsUGtrq/7+979rYGBAL774\nYkQdBgDYZzv06+rqtGbNGknSmjVr9PLLL/u1+eIXv6jExERdunRJV65c0aVLl5SWlma/twCAiNgO\n/Z6eHjmdTkmS0+lUT0+PX5tp06bpBz/4gb70pS9pxowZSk5O1p133mm/twCAiCSEetPtdsvj8fjt\n37x5s8+2w+GQw+Hwa/f222/rl7/8pdrb25WUlKT77rtPL7zwgr75zW8GOeMOr9dzJc0brf8AYJjj\nklptHx0y9BsaGoK+53Q65fF4lJqaqu7ubqWkpPi1+dvf/qavfe1ruuGGGyRJy5cv17Fjx0KEfuXY\new4ARpon3wnxzrCOtr28U15ertraWklSbW2tKioq/NpkZ2erublZH3/8sSzL0sGDB5Wbm2v3lACA\nCNkO/ccff1wNDQ3KysrSoUOH9Pjjj0uSurq6tGzZMklSQUGBHnroIRUXF2vOnDmSpHXr1o1DtwEA\ndjgsy7Ki3QlJg58JNEe7GwAQZ25VODHOHbkAYBBCHwAMQugDgEEIfQAwCKEPAAYh9AHAIIQ+ABiE\n0AcAgxD6AGAQQh8ADELoA4BBCH0AMAihDwAGIfQBwCCEPgAYhNAHAIMQ+gBgEEIfAAxC6AOAQQh9\nADAIoQ8ABiH0AcAghD4AGITQBwCDEPoAYBDbof/73/9eeXl5+sxnPqPW1tag7err65Wdna3MzEzV\n1NTYPR0AYBzYDv38/Hy99NJLWrBgQdA2AwMD2rBhg+rr6/Xmm29q165dOnXqlN1TGuR4tDsQQxiL\naxiLaxgLu2yHfnZ2trKyskK2aWlp0axZs5SRkaHExEStXr1ae/futXtKgwT/y8k8jMU1jMU1jIVd\n13VNv7OzU+np6cPbLpdLnZ2d1/OUAIAQEkK96Xa75fF4/PZv2bJFd99996i/3OFw2O8ZAGDchQz9\nhoaGiH55WlqaOjo6hrc7OjrkcrkCtp05c6befvvWiM43ueyMdgdiCGNxDWNxDWMhXc3OcIQM/bGy\nLCvg/uLiYp09e1bt7e2aMWOGdu/erV27dgVs+9Zbb41HVwAAIdhe03/ppZeUnp6u5uZmLVu2TEuW\nLJEkdXV1admyZZKkhIQEbdu2TWVlZcrNzdWqVauUk5MzPj0HAITNYQWbpgMAJp2o35Fr8s1b3/72\nt+V0OpWfnz+8r7e3V263W1lZWVq8eLH6+/uj2MOJ09HRoUWLFikvL08333yzfvWrX0kyczw++eQT\nlZSUqLCwULm5ufrhD38oycyxGDIwMKCioqLhC0hMHYuMjAzNmTNHRUVFuuWWWySFPxZRDX3Tb956\n+OGHVV9f77OvurpabrdbZ86cUWlpqaqrq6PUu4mVmJioX/ziF3rjjTfU3NysX//61zp16pSR4/HZ\nz35Whw8f1smTJ9XW1qbDhw/rtddeM3Ishjz99NPKzc0dviLQ1LFwOBxqbGzUiRMn1NLSIsnGWFhR\ndOzYMausrGx4e+vWrdbWrVuj2KOJd+7cOevmm28e3p49e7bl8Xgsy7Ks7u5ua/bs2dHqWlTdc889\nVkNDg/Hj8dFHH1nFxcXWP/7xD2PHoqOjwyotLbUOHTpk3XXXXZZlmfv/SUZGhvXee+/57At3LKI6\n0+fmLX89PT1yOp2SJKfTqZ6enij3aOK1t7frxIkTKikpMXY8/vvf/6qwsFBOp3N42cvUsfj+97+v\nn/70p5oy5VpcmToWDodDd955p4qLi7Vjxw5J4Y/FuFyyaRc3b4XmcDiMG6OLFy9qxYoVevrpp/WF\nL3zB5z2TxmPKlCk6efKk3n//fZWVlenw4cM+75syFq+88opSUlJUVFSkxsbGgG1MGQtJOnr0qKZP\nn65//etfcrvdys7O9nl/LGMR1Zl+ODdvmcLpdA7fBd3d3a2UlJQo92jifPrpp1qxYoUefPBBVVRU\nSDJ7PCQpKSlJy5Yt0/Hjx40ci2PHjqmurk5f/vKXdf/99+vQoUN68MEHjRwLSZo+fbok6aabbtI3\nvvENtbS0hD0WUQ1975u3Ll++rN27d6u8vDyaXYq68vJy1dbWSpJqa2uHw2+ysyxLa9euVW5urjZu\n3Di838TxeO+994avwPj444/V0NCgoqIiI8diy5Yt6ujo0Llz5/Tiiy/qjjvu0PPPP2/kWFy6dEkf\nfvihJOmjjz7Sq6++qvz8/PDH4np94DBW+/fvt7KysqyZM2daW7ZsiXZ3JtTq1aut6dOnW4mJiZbL\n5bKeeeYZ69///rdVWlpqZWZmWm632+rr64t2NyfEX/7yF8vhcFgFBQVWYWGhVVhYaB04cMDI8Whr\na7OKioqsgoICKz8/3/rJT35iWZZl5Fh4a2xstO6++27Lsswci3feeccqKCiwCgoKrLy8vOG8DHcs\nuDkLAAwS9ZuzAAATh9AHAIMQ+gBgEEIfAAxC6AOAQQh9ADAIoQ8ABiH0AcAg/w8IPFPQev2qwQAA\nAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x108692a10>"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "t1min = 0 # domain limits\n",
      "t1max = 40\n",
      "t2min = 0.1\n",
      "t2max = -1\n",
      "zmax = 1e-5\n",
      "\n",
      "N = 50000\n",
      "accepted = 0\n",
      "samples = np.zeros((2,N))\n",
      "\n",
      "\n",
      "# generation loop\n",
      "while (accepted < N):\n",
      "    \n",
      "    # pick uniform numbers\n",
      "    t1 = np.random.uniform(t1min, t1max)\n",
      "    t2 = np.random.uniform(t2min, t2max)\n",
      "    z = np.random.uniform(0,zmax)\n",
      "\n",
      "    # Do the accept/reject comparison\n",
      "    if z < mylikelihood(p_fail, T, t1, t2):\n",
      "        samples[:,accepted] = t1,t2\n",
      "        accepted += 1\n",
      "\n",
      "plt.hist2d(samples[0,:], samples[1,:], normed=True, bins=200)\n",
      "plt.xlabel('t_1')\n",
      "plt.ylabel('t_2')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD7CAYAAACG50QgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX90VeWZ7z8H6mWuRAw/zMFAew84cUWTQII4uu4cTfiR\nKW0BScPFdiACJYz1LilQbaCLeAkV2oSCQGWu42oYoMS1ZlikIAIyRSGRtJdqDBRJJ6vUcEQJCRNC\nisEKmpz7x/u+e+9zcgL5RXKS83zWytr77P2es1/3kmc/+3mf5/u4/H6/H0EQBCEiGNDbExAEQRB6\nDjH6giAIEYQYfUEQhAhCjL4gCEIEIUZfEAQhghCjLwiCEEF8pbcnYHC5/g54r7enIQiC0Mf4H/j9\nvnaPdoVLnr7L5YKJfih/BUgKOvsW4Nb7dT07sVaUAGm9PIf2UEL4z7OE8J8jyDy7mxJknt3JGjpi\nxsPG0wegvBh4BriqD7ymt3OB00GDe9v4C4Ig9D3Cy+jjBvKwvPqoZ9S26RVgqh7zVtB4kAeAIAhC\n+5CF3A7j6e0JtBNPb0+gHXh6ewLtxNPbE2gnnt6eQDvx9PYE2omntydwWwivmD7/V3+a23rAtCFq\ne9gcyMP29GP09oPbMzlBEISwpS/H9KlDLeIOCTwcBRwu0x/G6e1U7LDOpR6YmyAIQt8nzIw+KO/9\nFbU7Ucf0y69iGfso/UBoqrOPWR7/a0h8XxAEoW3C0Og7FmrLz7Y+na2N/r5M8JnzJrMnBtvom7RP\nCfkIgiAYwszouwn01I0xv4SVvbNZH/K09X3j/Zt0T2P8LyFvAYIgRDphZvSDjbLx0pOwjXie2viu\nwr/Eqf19cfZXTOi/yYR8zG+4EQRBiHTCzOgH4zTcxmM/a5+rDxp++Kqd5VOl1wN8wYVebV1D3gIE\nQej/hLnRdxri4Nj8OMjVu2eMoXdk/WTrba45luT4PafXL8ZeEITIQYqzBEEQIogw9/SdBMfoP8Aq\n4npUe/PxQBX2PsA8vS3zgs8s8h7R20tIlo8gCJFEHzL6wWEYN1acvkkb/3IgKi5w2Cy1GbDhGi35\njnRP0CmfbcX6BUEQ+h9hJsOwuoPfCpZhcNvHjPE/4RhuFn736W0VDlkHXRAmMX5BEPoUfVqGobM4\nQzM6hNOks3b2ae++DDvkY6hyfjB6P0eQUI8gCP2VPm70g7NxLmEXdHnVpt7xURv9zMwiAIqZZ3v9\nG/TDYUEmNBktH/H6BUHoX/Tx8E4ogkM+ZvH2NMxTsfwBG64BkOHea33rnC7xrbjLCyP0QV+x3hHP\nXxCEcCUiwztO6oK22uh7Mq1q3ZZ9gwHYOyuDNW71oEmjBICKHV7I11+t1wu+TWAreYr3LwhC36Uf\nevptkYT9FqBCPxP8ZdbZhewAYHXzGpIHngTg6MPT1clye5zt9YvxFwQhHOiYpy/FWYIgCBFEBHj6\nzhi/CdFoXZ55MHZXJQAezgEwiRK2sxAAX50H0OEg4+wXmd81KZ4gXr8gCL1HxMf0g3EaZLO4qy14\n0TiqqxIAqN6gtkNTG9nCDwCYFF0CwNeyz9NQP0p9J0r/RNNc4KWgawVLQwuCIIQXEeDpOwlW1Jxr\n70/TKZ7ZMCFTPRSWW+L9sIllgHoTANj4cK6d598kWT6CIPQW4unfhGAv3CHBYFXmemnMHApANFcA\nWM0aMnRC/zHSABj7XiXVWertgKJMx29Klo8gCOFLhBn9UNr5SYHHyqD6VWXMZ/A2AC8+/bwV5/9w\nRSIAFwqGMTr7svqOKQA77JCBsBDjLwhC+BBh4Z32kAdoCQfdkGXymwf4BUsBuEI0AGPw0aj3E0s/\nVOP3AZuDm7aI0RcE4XbSsfBOl4x+Q0MDTz75JB999BEej4fdu3cTHR3datz3vvc9Dh48SExMDB98\nEDru3TtG37nwqnvw8gHWgm+UDtvEY6l1+t93AbD218/hwQfAJI4BMLrysl3YZeL95cVIyEcQhNtH\nDxr9nJwcRowYQU5ODgUFBVy5coX8/PxW444fP05UVBRPPfVUmBn9UCRhG2mH8ddqnW8kTAFgeulR\n+L069nzOiwBsGPUCrl/o22luQzkEpneCGH9BELqPHjT68fHxlJaW4na7qa2tJS0tjaqqqpBjfT4f\nM2bM6ANGPxRJME15/ZPfPABACqes8E7hmSVq2GVYm/ocAHv1q0HFcq+d419ujH8MkukjCEL30IMV\nuXV1dbjdauHS7XZTVycerCAIQjhzy+yd9PR0amtrWx1ft25dwGeXy6W99f5AiCyfw8pLP3qXquY9\nmj2dFzc9r86dVZuyjAmcIhmA92seA2DrpkUsWV6oBtTrSmBfGW17+lLgJQjC7eOWRv/IkSNtnjNh\nnZEjR3Lx4kViYmLaHNs+Shz7Hv3XGwTr9Dto0ha+MI5Tm5SBv5AxDADvmQq8l7PUed2469kz29i+\naQEAFV5dALbAC03G6Jt7Jou9giC0B5/+6xxdytOfOXMmO3fuZMWKFezcuZNZs2Z15edAFz6FD8YA\nOx9mb6lNk5ti1yMAFC9T+fqTNx2wKnbTdEbP0NhGq7J36awxADTMGgUjtNe/2VTzOh80YvgFQWgL\nD4EOcWmHvt3llM05c+Zw/vz5gJTNmpoaFi9ezMGDBwH47ne/S2lpKZcvXyYmJoaf/OQnLFy4MHAi\nYb2Q2xY/VJu1uuvWLHgxQYV8cgs3th5uwkAFE0iteweAlgVK25/D+qQ0ahcEoUP0YPZOd9L3jL4z\n9KO99rVYmTomy2cNq/GhPPx5Ncqr3xq7iGSUZv8OXem7bfmz6ouby7DeJgKuId6/IAihEKPfgxgJ\nB8fCryno2qE2xzMfYh8ZAMxCtWf01lRQGTsWgGgaAXiKnQAcLZ0OacbrN/1+Jb1TEIS2EMG1HsQY\n40v2oSa91Q3XH6t/nwlPK/ffxPmLYjOtKt7PXA0AvP3LGQAMnNFEyx69Cjzb6d2bB4t07hIEofOI\n0e8WnFIOemG2KF1ty4ZQUaaydnbsUqGcPWeyYJf+ygW9PaQ2r7tnciozBYAXSjaog2lnad37VxZ8\nBUHoONIuURAEIYIQT79b+QDb+9YhnxHPWIu7x5onAXAtbgCDH20BYErsGwC8jQrv7CODwhVK1sFT\n4AMg68weKNQhn83O+L5ZU5CYvyAI7UOMfrcSojVjeR5GwbNB3+2oM818kjEcwJJsRn1Uyp26vitF\nZ/g8l7CWjfG56uA0nSl0uAwJ7wiC0FHE6N82PgixP1dtHoWnPlXZOm+fUR7+1oxFACxku/WtY6g3\ng//zNy/g+dwHwJIRWtIBrzb8Ti4hDwJBEG6GpGzedtzYFb3G+OdhYj6T/Spl03j8e8kgQ6d2ntNV\nd9E0UueqACBTL/zOjt1F8avz1IfvO1M8Rc5BECILSdkMM+qwjb5Juyy28vmP6grqRKYD4D/rYna2\nSu3Zk6bjPGvggH8yAGU6rx9g7NOVAFR//6vqQFSmo0m78/qCIAgKMfo9wk287zS93aM22dkvs2ev\nMvZbS1TI59mabappC1CZqoq6trDUKuhCP+Sr70oAxukfNHIOktopCIKNhHfChYl5alsPE86p0M/7\nZ5Q889rE58gtVVo+F1KHWV+Jvq68ft8gDwCPN79Dw4JR6mSR6dX70m2dtiAIvY3IMPQRgmWbHa0Z\ntVjpol1bAdWZqyhRhYNK9KtB4ZkloMQ9Le9/LxnWQvDjKEG3atc7jmuIxy8I/Y8e7JwlCIIg9C0k\npt9rtOF1N5VBkZJt2FallTffwyrYOlWgGreUJU7Au15l9OxNVYJuyZxk1F6l5TMpQ2n7VM97ForM\nj5sUT2c6qXj/ghBJiNEPG5wNW3QD9XIV8tmW9ayl1ZOmm7R4ayqozFFhHZ9O7cw9s5Fr3xwQMK6w\nYAmuWfrVb4Hu3NX0AXZGkRh9QYgkxOiHHSFkFoqusm2E8vqTN51Sxxog8cqHAPjrVW/iytSxzGQ/\nAB+WJqpxcbArczYAWSN0itCCZ8BnUjtD9AMWBKHfIgu5fYI8QGfjLNNdurLB//+Usc/OfhmABWzH\n+20V8tn6a13he307JwepkFDKdfXAiNrRDKawt9wZ8hHDLwh9DynO6ocUA1qqebM2/mVDGHhAifc3\n10SpY/9o5/ab5iyDL7fgi1Wdu7xb1APheM5DPBb/vvrO8zrkUy6ibYIQCYjR7xNcws631+GY8rm0\nFCqv/6FVxwFoLBlqSTiYhi383l7UXZvzHKD0ffxn1VvClPeUyufRbzwDh01uv+nYZdo2CoLQXxCj\n3ycIod5JHeQq41wRr731elj49MyAb07JeIO3C5Womy/bA8Co9Q1cW6oWfBfqvo4pb55ioytNf8vt\n2ErIRxD6E2L0+xzOFo36ATBbL8ruyWS7brRu+vG+XTPDkm02aZ88CqsHrQFgzXW1jjLvRDEevwoN\nLblLB/ybThPYqQvkISAIfRtZyO2zOCt6TTOVcfAvKuTjn6HCN2WxExiDD4BRNSqHvyg2k3mF6kEx\nJVs3cVk/gws5SuJhr27kvmR5oZ3aX553W/4rBEHoKlKRKwiCILSBePr9gql6+wGW1z9PxfnH7qq0\nOnAtYxMAqXXv0PxfKuPnQKKSbG5kKBnXVUho8CHVyvFAxmRmrHtb/V6uuZYuHJMwjyCECZKyGUGY\nEE8IWQUt2la9PIGMTcqYP9GsCrdOu8dhVJnHJPoA8FwvYfAWbexz1INg+vqjZK5SGg7F8bphy2xz\nnSSkN68g9D3E0++3mAfCXFir4vyLVmnVzsIlPJ/9IgAb1r8AQFnOBHYELQL7GEMj0QAc0+qeR9ep\nZi8U4qjqFeMvCL2HePoRTBvetw7NGCkHnra1efim2qxmjZXaOXzhJwCsGbjayvfPXaH0/GcXKBGg\nYu88SDMNWwwfIFk+ghDeiNHvV4RSz3wNK86vq2+38Sw8rQ7NQ3nrjUTDIXXs8vDRegfuyz4DwKQC\nZfy36F6+e4ZnMbBWVQS3jDTXEikHQQh3uhTeaWho4Mknn+Sjjz7C4/Gwe/duoqOjA8Z8/PHHPPXU\nU1y6dAmXy8U//dM/8YMf/KD1RCS8cxsx3rc2/hO9EK92z+y6D4CEM9VMSVTpmymcsr65oUaFfypj\nlaJnQk01ANeGD2D+ILUwUPywjveXg73QG9wMXhCE20MPds7KyclhxIgR5OTkUFBQwJUrV8jPzw8Y\nU1tbS21tLcnJyTQ1NfHQQw+xb98+HnjggcCJiNHvJoI7ctU5jjkatHt0Fa9PbV72Z+PhHGAb/VGj\nGuA/An/tWpzK8h18ooWHUpX8w36eAGB06WW756+V4F+HGH5BuJ30oNGPj4+ntLQUt9tNbW0taWlp\nVFVV3fQ7s2bNYsmSJUyZMiVwImL0e5Ak7GbtzwAwwV9mibT9QodwAK7ohdzNLAcg2fEWYPr2omq+\nmJ2xi+Li4CyfMkTDRxBuJz1YnFVXV4fbrbxIt9tNXd3N47k+n4+TJ0/yyCOPdOWygiAIQie55UJu\neno6tbW1rY6vW7cu4LPL5dLeemiampqYPXs2W7ZsISoqqo1RJY59j/4Tuh9HEZfW6a/I8oJ+SfvB\ne/bI5bqga09pFgBbUx3SzVrTx3j6e2qyIE+Ncy3Tnke919GuUTx+Qeg6Pqy4bCfocninpKSEkSNH\ncvHiRSZNmhQyvPPFF18wffp0vvGNb7Bs2bLQE5HwTi9hjP8liFKhHrPIy0p4MfP5gNGnUA1Z9qzP\noignM+DcvMJizi5W+3EX1NZ1xQ/1ekCaVPMKQvfTg3n6M2fOZOfOnaxYsYKdO3cya9asVmP8fj+L\nFi3iwQcfbNPgC72Jie3XQZPpzTtXbeuH8ELxBsBuuWjF8edBNFcAmJ52VB2Lgzv9SrTt2nW1PuDH\nxZQErdm/TD9UNr+CGH5B6B26nLI5Z84czp8/H5CyWVNTw+LFizl48CBlZWU8/vjjjBs3zgr//Oxn\nP2PatGmBExFPv5cIzvYBy/uP8sJKtfvyqmzA9vSTOWU1bDmnw3BDaWy18OvBZ4WIZvI6oENJRWf1\ntV7rxv8WQYhEejB7pzsRo99bhKqgTWo9bI8K5UzOPADA26UzeCVNnfqqf7I1bPpe2+sHuC/xDO/w\nOADR2vs/NiiNGcVayM30ArDeOILnIgjCzRGjL3QrQQ8AbfwHeK/RvFMvyOsszbLYCQzVaZ/HmATA\nszXbQDv1zrz+pWwBoNilv0wZtrGXvH5BaD+ivSN0G25aGeB9yui3FA5m7Zuq564RZTPVuwDbY5V4\nG2dhSqqK6ZtmLkvZYuX7R/vVQ2Jb1rOwT3+5yXj9Ma2vLwhClxBPX7gFIap5QcX7d6hdE/JZw2p8\njAFgXqkK25SlTsBbWKEGPqo2RYmZlnZ/1J+b1W8kHLAVPHNNvP8tJNQjCLdCOmcJgiAIbSCevtBB\nHHn9nmcCzgw4cU01aAE8132A6sJ17ZvKtxg36DQA+5lpxfxNBtCoMw3MTtSyzd/Qcf7DZ4HT9vUA\n8fwFIRhZyBVuC6GyfMwxZfzH+itZyHYAcs+ofP5rcQMYfEJ15Lr2qDL+vkEeEs4otc4LiSqvfzsL\nLe0fw5J1hZB7NWgeL3X9P0UQ+hWykCvcFkJ52CbOrxQ1q8d4eSFeFXMlv6n68k4/c9TK3jEx/StE\nwwm1P2qY0nDI/f1GijICK3z9D7psOYfNxvhLm0ZB6Ari6QtdIEQ+v0endJ64BsAa92py9yqvvzJj\nrDUs4dvK0y/79QQATpFiefrzzujc/WHwfKxq67ipThV7tYxspnVBl4R8hEhGPH2hx7jU+pBPSTm0\njJwKwAtsYK9fyXMsZzOA0u2Pa/3VSbo1o07hh2Ewq0DF/KPd6oHwwp4NMNusJeTprbOqWB4AgnAz\nxNMXukBwnN+N7f3rJi3z7NGmS9c5h3qq1bBlcQNbf6kUPLezALBbMwJ4v63SPrf+ehFLsgrVQSuv\nP69L/xWC0LeRhVwhLDAPhLlYGThr1YNg8qoDvF2jmrAb2eVrSwfwta+cB2DSQOXx7ynNojJVt2nc\nq8JBUzLeYJKW4H5hnVo/oB7YbDp1iXyzEGmI0RfCAudbgN43KZ7ZsGjVVsDW6z+Hx/L6TYWv57qP\nvYMyALvY60LqMPaijpn2jjOK3wbTpbPcLPhKlo8QKYjRF8KKEPH2qDwYoXbNgm+Gey87r88HVEon\nqCwf75mKgF+rTLQXg61agBMt3Jd6BoBqV4I+m9dN8xeEcEcqcgVBEIQ2EE9f6EHMIq8bUJW7LBui\ntvHg/3vdblPn9VdmjLXSOM3ib2rdOzQ3K3XPsliV7rmZ5da4bcXPqi+HbMweqsBMEPo6Et4Rwo7g\nRi0x2OmeukvXtCFWwo9p2OLhHNNLlT5/Waoy8N6kCvgPNa4yVoV6TpJi5fabCt/Rr15Wg77vrOg9\nordS3CX0J8ToC32CEF732jwAxq6qBOD09XHWQq5Pe/qnSGaN/v/EeP/RNFqLv6dICRi/rfJZqxKY\nJmeGj3j9Qn9BjL7QJwj2/ucCOtSje7NM/vSAZcxNYde89cUcyFGduoyBX359E4PPKn0f50Iv2M1c\nAJbcpfP7mwBeIRAx/kJfRSpyhT6Bs6ALAlIsvXkAHB0z3dLs35fqA6AxJ5qF15WoW+OgoQAM3tJC\nUY6Sf0hBaf5s0j16C2uWcF+syuwZ+6l6g6h+OAHKTVWv0/iL9y/0f8TTF3qJm6l2mmNz7WMl6i3g\njdQploefm6Y0fYpKMq0q3pI65dm/7p4JwKTrJdZbgIn3L2WLLd9sIj5S1Sv0WcTTF/oEobzpYO//\nLUBp+JCvjP4M3sYfp7J8skteBmDL9aV4BqlCrX1u3YZRZ/P4Bnm4kqhCRJu197+FpRR7HfoQAIfB\nzi6ShV6h/yJGXwhDnMZfPwAOa+XNqkxcO5RXsyt2dqtvGmOfcl1X9w6KtuL6e85kqUFbYNcv1XdP\nkQzAxrvyoMlk+ojRF/ovUpwlCIIQQUhMXwhDgjN7nMyFiSrUM+zEBQBeHzgTb41W4YxVSp1WG8YV\nDXbKpi76WpvzHJuaVajn/JdfAyDqRLMd388tdszDeP2yuCuEK5KyKfRr5mJV2E5zCLhlKgG3wsIl\n6tg31aYoNtMq3CpLnNDq13yMAVSYx2rUMn2wOll+ldYFXdK5Swg3ZCFX6Ne8hrXgeli75lVetpUp\n+QXPJh9gx/aTOWm1ZixJVLF9Dz4yrqs3Ae9lLegWC825qkBg7XvPAfDCqxvg++n6us7FZUHou4jR\nF/ogRsJBG2LfVahSIZ8XspTG/uRdBwBl/L0ow56sc/hTONUqjXMSx6wmLmbc2Kcrqc7Xqp0+c23x\n9IW+TaeNfkNDA08++SQfffQRHo+H3bt3Ex0dHTDm888/JzU1levXr3Pjxg2eeOIJfvazn3V50kKk\nUxe0/UCnXAJRKuRztHQ6AI2p0aCkfOzeu/8H0Ik8o1CN2dcmPmf9+vRRSu+HCz9gUlWJ+tk/N6tj\njwJN5qHjzDKSmL/QN+h0TD8nJ4cRI0aQk5NDQUEBV65cIT8/v9W4zz77jDvvvJMvv/wSr9fLhg0b\n8Hq9rSciMX2h3QQXcSVhe/86r9+jm/DugF2pKj3TVOsm1FRbC74mDDSJY4w6E/gASOMYJTrdM7dG\nFYK5VvihyLnQC8rzF6Mv9BY9FNPfv38/paWlAMyfP5+0tLSQRv/OO+8E4MaNGzQ3NzNs2LDOXlIQ\nNMFFXM5wiz5Wrz+mXWW7fwFgd9rKjn2ZBSgpB7OQe5JkTiaix/kAWM0afqH79F6IVf/fTthVRoVX\nST5YPXoP1yFGX+grdNrTHzp0KFeuXAHA7/czbNgw67OTlpYWJkyYwIcffsgzzzzD+vXrQ09EPH2h\nS5hq2iCPHyyvf+w5rb1TmcCweJXuefk/RwPKu1+oHwRGvdNbWMG1rNalLIMvq/UA19/rfzo+AOP9\nS7xf6Gm60dNPT0+ntra21fF169YFfHa5XNpot2bAgAGcOnWKv/zlL3z961+npKSEtLS0dk9QEARB\n6D5uavSPHDnS5jm3201tbS0jR47k4sWLxMTE3PRCd999N9/61rcoLy+/idEvcex79J8gtIdgD9uk\nVk61Mm+qXf9N7SyD6E0qlr81UcX2F7LdknE2cfwx2T5LwiEa+y22MVapey46p2oDtmU9C0U65EPw\nIq8gdDc+HOlkHaZLC7nDhw9nxYoV5Ofn09jY2CqmX19fz1e+8hWio6P561//yte//nVWr17NlClT\nWk9EwjvCbSNYSC0PJqq9F997HoDcvRupzLA7cQE0Em2tA0y6XqLODUq2ftXo+3ztK+etYw1fGaX3\nirGRkI9wO+mhityGhgbmzJnD+fPnA1I2a2pqWLx4MQcPHuT06dMsWLCAlpYWWlpayMrK4kc/+lHo\niYjRF7pEW1r4TkkHc24qVm/GtWozdlUlH65QK7nPF7wIwIYVL1BZENiUJWF9NVtz1NvBszXbAJgd\nu4tlbALgsYffVwOrcMg1m/UFp/GXNwGhuxAZBkEIIlSKpz7m0cZ/Jbz4tPL6TRrn6uY1XP6+WujF\nKDHHwdpYO6UT1BtBCsrrN+GgrKw9dnaPpd7paBQjCN2GyDAIQhDBXrXD4zapnc/DC2WqmpcRarNo\n01ZYo88X6W2c3X/XW1phnXv+l+rtwDwwju96yNLvL3Y9or+cBIzT+691/j9HELqAePpCBOEM9QQn\nHqRDlO7Rq43+5HMHmKSTC0w6J2At+CbUVKsDDXAtTqV2Ng5S5zax3Ho4mBTQCtc4bIKF3AShs0h4\nRxDaQSj5Zv1G4MlT2x0wIP4aAOfdSoJ5VE0D14YrA+8b5AFUSMdIOW9nIQC5pRu59mjguKfYScUY\nHU7ySchH6C4kvCMI7cAZ3wc71RLwaeH9NDcta5X3v3qVivMUsoTVg9T+hjMvAJBwtpq1GSrObxn9\n32/kWGoaANP3Ki2fjIx9VKxURn/s0x8DUO2a65iPSTOVto3C7UM8fUFoheMtYKLW7I/H2j63SqX8\nmIXcMY6caZPuOW9xMdm/VD18C2u0xv/vYUrGG4AtCEea8fjrsOP8JttHZJyF9tAxT1/aJQqCIEQQ\n4ukLwk0xoRbt/e/xMsB7LWDETvd85pXqYiwl1ElZxgS8K1R2z9oCO8XTtHU0Am6jH76svlAP+Ey/\nRsnnFzqCxPQF4TagjW8VtFSpdoovrlJ5/RnX98IhPUz34/XWVKAFOll+XRVuDf5uC1t/HdioxfT5\n9Qw8R4XLSI4bQ3+JtovOBKFziKcvCB1iql3QZWz0SpiQoLz0LdrSe1dUkF2gY/rrdUz/BBz49WQA\nq5hrVKl6NShLncB8fgVA9RjTrasMiesLt0ZSNgWhGwnOpJlKK+/b44Udave5VLXI68PDnjOqPVdl\nopJyiKbRyvE3uftG3M0sCgPMKH5b7ZQBm3UmkbXIK126hGDE6AvCbSYp6HMmLNO7VWqz683ZloGv\ndymNnqn+Cdax6Wd0S0Zt0w9kTLaqeU3c/0DsZGZ8Qz8ATDtIyrAfQGL8BZCYviDcdkLIJ2/W0sra\n+GdV7mFXgmrTOMlv1DvHWKqdRYl6vO7W5cNje/t6MTgl9hRj39SNX+7SIZ8mMfRC1xBPXxA6THA1\n71TgtN7XUgueOIb9WS3Svj5wJgDeMxWWsZ+3Xmf7aFmeh1KPB6wHAJQVTOCUzvtfsrxQDazH1gHi\nFccc5GEQuYinLwi3mWADG8Lg+s7S8LeqTeMTf94PQOM90TSXRgFQljMBgH1kAPB+6WMcSFWLvCbk\nMwYfQ3XI58VNKlNoOwup3me8finiEjqOFGcJgiBEEBLeEYQu48yoMaEfh2b/NOXxsxKaHh0I2CJs\nhmNMYuF1peQ5+KxqvH4hcRgnUZ26fIwBYBLHSKz8UH0p0Ug4DMEO9ZhFZvH+IwcJ7whCD1MXYr8O\n+KHaPWyfjfI2AzBspYr3Xy5TTVoShldTlqhCPilxKof/JMk0MjTgStE08mKCCvXs9c8CoOJVL3w/\nOLVUUjuF0IjRF4TbhtHM19k+h4HDSsCtIV710s3OVAVcy9lkCbcNvqw8/ekctVM6dbx/VGEDZKtj\nZuF3/tOUvz5qAAAcKklEQVS/ovp5XSnWZN40zMIyiPEXnEh4RxBuG6HaNBrS1WaZkm72J7hYm600\neoxEw2rW8P7ix9Q4067x91CZo1JAPdd9gAoHTUnU6p3FWr1z9llswy8Szf0bCe8IQpgQHOf/gFbp\nnpvVxoWfzGyVi5lbsxGAxtihoJcDslO1pEPcElbrHo47mQ/A1sRFLNeN2Y/Ga6M/Lw6K4oLmcwnx\n+gXx9AWhVzBef6Z96F/0Vlf1TthUxjvXUwF74TeaRhXiAQ5kq5DPpOslHBuUpsbpBd8S0ih+WL8e\nlJsF39ew20RK7L//IDIMghAmtEch0xj/cfb4aSrkM/nNA5YH72zUYnR7DD7G8GyNknpwtnJMLNVZ\nPvohwvedIR+DhH76PmL0BaEP4gaP7tJl+ueWDLHOHk99CFDFXOZBYLVm3LuRygwV5zcPhDH4rJh/\n1AmVMTQg/hotI5v1L0pj9v6DdM4SBEEQ2kAWcgWhV3GEgHx5el9v085ahV1PeJWUw6SBxyxvPrdQ\nLfgWZWda6p0m82dc3WlK3Y8D4G9wATB8xCc0RKlUUTu1Eyy9IEu+WejPSHhHEHqFYNE2aL3I+kMr\npROtt8YB+CR1OKCqeAGy1u3h+CoV/jFaPSdJwacfDkayOZmTLGULoAu6APIBX3HQPCTk07eQlE1B\n6AMEL+66CUztBHgJNuep3Yn6UD6MHqH66h5PUIY+c1UR3lKtwa+LuP67q5iFftWHd1SNyvYpis1k\nEiUAVJi2X76zWDUD4ulHBOLpC0JY4PT8nQ8Ek92T7jimvX/TuCUb/H9SIZyiDJUCOolj1pvAPpRc\ngwcf25vV4u/5L78GwLFBacwYoxu1jNC/V15M654B0qs3fOkhT7+hoYEnn3ySjz76CI/Hw+7du4mO\njg45trm5mYkTJzJ69GjeeOONzl5SEPoxbRlTY3xNto3b7tGrC7sog4feOw7Y0gwA8wpV2KYxW/27\nXHh9OxsOvaBOah1/X+wYBpy4BkBLvmr4zqxMyC3rwn+LEM502tPPyclhxIgR5OTkUFBQwJUrV8jP\nzw859qWXXuL999/n008/Zf/+/aEnIp6+INwEp6ethdyitMfvhQE7lOE+7VaLsidJIZorAJZo2yaW\nsVw/Kcy56XuPgooCMSVVSzmUToc0Y/TNdSX0E770UJ5+fHw8paWluN1uamtrSUtLo6qqqtW4Tz75\nhAULFrBq1SpeeumlNj19MfqC0B7c2Au+ppr3FTvHf6XajH26kkm6/eIax78rk9u/ECXjPKqmga2x\niwKusLp5DdED1eJvtav1v2lZ6A03eii8U1dXh9utvAC3201dXejX0+XLl/Pzn/+cq1evhjwvCEJH\ncP47y9PbqXYGzvfVg6Dam0B1veqwFZ2qDHg0jVZK56j1anF34Pwmmterbl4XcpTL7xvoYVPdcvV7\nUaZLl5b7bFXRK/Q1bmr009PTqa2tbXV83bp1AZ9dLpf21AM5cOAAMTExpKSkUFJS0o7pOMd49J8g\nCDahGra4sb1v3Uyl8BkmbzoAwMZXcwHY9fRsJl0vUedV5ibNcVHW/uocJeS2gO0ku5Wmf9aJPerk\ndC3eVh8HTWYuIRrECz2AT/91jpsa/SNHjrR5zoR1Ro4cycWLF4mJiWk15ne/+x379+/n0KFDfP75\n51y9epWnnnqKX/3qV238alpH5i4IghCBeAh0iEs79O0uLeQOHz6cFStWkJ+fT2NjY5sLuQClpaVs\n2LBBYvqC0CVCqWImYXv6c/W2DqJUls+AP6tF3jT3MTLYB2C3ZtzVQmW20u1JWF8NwLWlA5g56HXA\nLvbyaM9y4zdybQG3er1tyuv6f5bQBXoopr9y5UrmzJnDtm3brJRNgJqaGhYvXszBgwdbfSdUCEgQ\nhI4QKpRyybFvsmySoEnF+VsKVZz/aPx0GjN1WvUgtdmevYD31+tGLTqNc+ag163FXyPvcIoUdXIl\ndnWwoeiHiIBb30GKswShX+Es8poaeMwzxA4F68KuNzZNIQUVvx+1Quv0F0xmHxkAbLmu8v6Nnv/j\nze/QWK8eBC0LdF7/YQAj5SA6/T2PSCsLggDY1bz6TSDqGfvUBr2th7GrKgHYwg8AmF56FJS8D0WJ\n6i3BwzkASphkpXuOXnfZ+g3L+28KNv7C7UeMviAIAZ628fjfwkrz9OhDI7By+xdlbtWnfCy/rjT7\nlw5SqT2z2Auohi1GwM0Iuh0jjaMP6zaNs/Tv5pY5ri8PgNuLGH1BEEKmdkLgAwDgGZimd7X3n5lQ\nRJpOn362UHXkui/7DKCKunL3aknnDPst4LHi99WXjX7P80B5cG6/9Oi9PYjKpiAIAcbVafzfCjr2\nClTpsE+iKqAsXjuPk6vUwu3CLBXK+bAwEYC12c9Zv1qiU6y3XN/L8Uyl+PlYqTL+Y9+rpHqMLuzy\nmW+IlEM4IJ6+IEQUxus3dTVO71tr+vzLENDSO8/tWgvAhhot1HYWKlMDWzM+/jdHSflcvQmY2P/R\nyumgIz52aqdTxM2EfMTz7zrSLlEQBEFoA/H0BSFiCBXnr6NVlg9zYa1W8NRe+oRNykvfzxOWbk9R\njhF8s0M9pojrGGn4GANA9cMJ9m9ZLSGd1xet/q4hC7mCILSbJMe+MfoxEKUNerbaDFipqnp3uueT\ncX1vwC88PqiU/TwBwF6d399ItKXoaYV8sqajC4IdqZ3ORWYJ+XQOMfqCIHQKY4CTsBd8dZx/mvb8\n47GE3Axv751BZYaWcqhRUg5FsfZbgDH6qXXv0PKoLujymcyet7BxGnvx/tuPZO8IgtAlPsBO7dTy\nCoe1Ea8CVAq/WqwFBv7PJl5npjoYqzbRXGGMDvUYj7+lfrDdpatQq3bmxmFX84bKOBK6G/H0BUHQ\nhIqzJzmOaabZXjzA5DcPUFKn+vE2v6S0+Z8veNGq3E04o7x/LkN26ssAnCQZgIosryPkY3puvOa4\nroR8bo14+oIgdIpQhjW4mjbGTsEsV5ujWdOtqt4LBaoRyzHSSNaaPtGJqoJ3VGkDC/SDIFkLuI3Z\n5aN4zDz15XgdQiqfi53TL8a+uxGjLwhCCELJNwO4bWlldFy+CChSbwajl2k9njLgPdWPd1SNyvah\nAbw1FQA0xirRtr11W6zFYlMbwMQhUN66P0fozCOho0h4RxCEDuDs0evAZPvEq82wExesPrumV6/H\n0e0pTR9byhZO1SmvvyVfL/KWYcs5HM5zXNcgxj4QKc4SBEEQ2kA8fUEQQtBWCMXZpSt4sRdYprp1\nUQbsULsvJ6j4zXYW8Cvm632V0bOQ7da+Ue0sfnWe/XvP621TGa3XF8TjV0ieviAItxV3iGPG8I/T\n29PWA2DsJqXX39gczZqB6t/4kkolwO+vd3EgdTKA1bhlb3MGnoEqt7/CpR8iXAX0Qm+rhi2Rjhh9\nQRB6nGAht3TVqQvA2O1sGOa9AMDl/zUagMpfj7VSOi8kqsyfc3h47FWl1vny0+otYYmrEGX4nbwU\n4vqR6P1LyqYgCD2OMbbG43/J1tkx1bcj4qwHgOtFZaQWsZXkRJXaaZq1Dx3UiH+g6qf9EMcBmOAv\no6JYf3m2qeadSmupaIO0a2wL8fQFQegg7fWqzThHr949tpwDAPtgwiqVq9nIUAA+LE0EXbB7IVZ5\n/0vZQnGxjvXPNh7/EOw8z0ju0iXhHUEQepzgB0EoRc8kLFc/Sh8agSXNkOFWQm4L2M4mlgOwXGs+\nmAcCqAVhgKOl0yHNWcXrvH6oOfVXJLwjCEKPE2xYnVIOhg+wFnqbsLYthcr7T16lwjwz1r1tPRtK\n4pW8w3n31yzphqPFSvNnQmYZFSZbaPNc/YNaKyigS1ikGP/2IUZfEIRuoD2ePtiLr3l6exZy1d4L\n+bpJ70qsqE1GqvL+R51pYGmiSu0c4FVvBmPwUeE1Rl+HjYwu0OGrRGao59aI0RcEoRsI5ekH77ux\nHwCvOM7rxV9jwMGqyC0uVXH8KanRvF0zAwDXWb/9c2ZtYK3eGl2g7CEwW/f+tVI8Q80z8pCKXEEQ\nhAhCFnIFQegFTAz+tH3I6Pc0XbWbtjjF2PT+sHiV679l4FKyKvcAkJlQBEDxN3SGTz2WCqjNK/RP\nT1+ydwRBCFtCLaoGH5sLUTpn04RvvHZlr6/OA4DH7ePDmkQAXG9oM2aE2vY59jeb65j0Tgjs2NXX\n6aHsnYaGBp588kk++ugjPB4Pu3fvJjo6utU4j8fDkCFDGDhwIHfccQfvvvtuZy8pCEKfJ5SnHUIy\nuUkbfSPj7IXqdQkBx6pXJjAl4Q3rvPNc5q4iil3a6/focz7nNZwVxJG14Ntpo5+fn096ejo5OTkU\nFBSQn59Pfn5+q3Eul4uSkhKGDRvWpYkKgtCfCOXxO7tl6eNN6frYENtRr7e3jWhH0zwctHe/jE0U\nl8wLPFfmhSKT168LxqLioMk0hDf0xxCQTafDO/Hx8ZSWluJ2u6mtrSUtLY2qqqpW48aMGUN5eTnD\nhw+/+UQkvCMIAqCMv3koGC98Llb8f5p2670wdpUK+VQXq7eARZlbAaXf42NMwDmqaP3gKDeSDmAX\nePU1eiimP3ToUK5cuQKA3+9n2LBh1mcnY8eO5e6772bgwIE8/fTTLF68OPRExOgLQgRyK40ch2xz\ncLFXlNda3B2wUjdc141YXt6UzermNQBWMxcP5zjqmq6/q3+j6WyI69fRt0I+3RjTT09Pp7a2ttXx\ndevWBXx2uVzaaLfmt7/9Lffeey//9V//RXp6OvHx8Tz22GPtnqAgCP2Z9ur3OHP89XfisZqqt4zQ\nXbdmqc0mltOwbxQAazKVM7m6eQ3s0T9hItFVcdBkMohMKMnp8fe/at6bGv0jR460ec6EdUaOHMnF\nixeJiQnV0xLuvfdeAO655x4yMjJ49913b2L0Sxz7HuwVGEEQIgtjbN8KcUxT7rblm8sCT5EKw2ap\n1E4T92+oGsWETDXwlFe3aKwfDNN1qqjPKd1sUkrDMcvHp/86R6eLs2bOnMnOnTsB2LlzJ7NmzWo1\n5rPPPuPTTz8F4Nq1a/zmN78hKSmp1TibNMefp7NTEwRB6Md4CLSVHaPTMf2GhgbmzJnD+fPnA1I2\na2pqWLx4MQcPHqS6uppvf/vbAHz55ZfMnTuXH//4x6EnIjF9QRA6RBJWSCZKe/xGyG0edhqn2ZZh\n5f0fT30IgNWs4eg6Heevd4wrDyXZHK6hHinOEgSh3xFqwdcZNRgXeGptnBXvNymbkz89wFDUom5x\npd2Hd3LCAQCOLtfGfx+OUI9TslmMfrciRl8QhPbhzOc3DVpM7H2uY5w20s62jSv1qTJ7+9wupda2\nvVmpeDZMHwWHzW8Y438E0LF/SywuXIy/6OkLgtAvMZ62s5jqg6BzxvgnYRlln47v1NvfmrBLWf0K\nvOzVDdkbqlS2D9nY8g9GsplxBKp1mmsEzyP8EaMvCEIfoZ0SDtY2SHKhyQ1Kl42KEfpBkA3R6Poi\nZ1WvWQewdHvcDkE4Y/wvYTeC7zuI0RcEoQ8TqnlLMKcd53Ts35JmgIpHlYUf1qhSPKMHNlJdqat4\nz+hxidiKn5tNPv8R7LeOcI33t0Zi+oIg9EOcC7/OnHvjmTs0fUx17g693YcdCjJhnhEENHMHoMhZ\nzetszRhs+G9VddxVJKYvCELEEsrjDlXgpbNyJj5jnzKLu/HYBl5X7g44cI2WQlX1O2CDlnwoiqP1\nm0Wo2H54ef/SOUsQBCGCEE9fEIR+xK28auOJ68ybcof0Qr3O1FmJrc2zQ208bh+NK7Wcw3Sd5cNV\nmGiye/TKb/k4AnP7QaWVho+cgxh9QRAihFDa/aexFnd9Or5T5bUXegv1qZUeWqq0qJtZ0GWII7XT\nSDQ72j/yQ70NLwE3MfqCIPRzQrVmNB5/jOO89tY3n4WJunOXfg605A+22jWaDl5j36yk+lWd5ePR\n43112EVc5i0ivKp5JXtHEIQIwbnoGirF02T5DMEy2E5NH48+fUBvHcVeAZIPVjWv8f5fI7CKmBCf\nu4Jk7wiCIITgZsVdoPLuQeXea6NsBNwmeu1hxsAXYss6GKqAaXq/THv/I/LA9wrhghh9QRAEwC60\nqsPO59ehmvI8IE/tGxX5lVgxf8vrX4m9HnBYx4ZWeiHX0ZMXHFW9bse1nF7/7QsHidEXBCHCCTaw\nSdgG+JLjmA75FOqQj8+ZvaPZhx3eWea1j1kxfxPnN9dMJ3Rv3tsX+xejLwhChBNsYNuKs5vwj/b+\noxzZO1WOYR69dfaVMoVfvuDfPI2tFNozjdnF6AuCILTiJuEVk9oJloCblfbpfBDUO7Y79P6jQXn9\nTWBn+XTn4m7bSEWuIAhCBCGeviAIQocITrsEleaJEmYLpugqjDCpn0GpoB7sSuAmIwLnJlDADboz\n11+MviAIQituld4JKuvGPAB0uMZ3FesBsECfihpip3lOc0g9gMr+KXJ25wIVKgruCBZKMrpziNEX\nBEG4JaEeAo7Yu0dvfUPAZ4qy9HeivHac36Db8Srjr8et1QvEuWXYPX+DVEED6Jxks1TkCoIgdJpQ\n6Z5BHbs4bXfdMh6+WeTdh92ly2QAzQJyTR5/cF9e57UMUpErCILQA4TytC9hvQFE5eljcXasP1dn\n/pgK3xE4wjs69FMPdlqoPtXk6PkrMX1BEIText16vylPf54KTdpQe7QxN56+zyHuVq6PzcKh2qnP\nebwOKQfn24VzMbl9iNEXBEHoFHVt7CcFHavDquz1mTEOHf/4oO3mMpim3wSsoq6rEKW7fDWFknFu\nP2L0BUEQupVLQZ/d2N550CIvXkd4R+NxLPwe1ufWDrHfDjYH/3bHkOIsQRCECEKydwRBEG4rzgVf\n003rSNB5VGonKC/fZPIYj78c7LcE490PAYqB2R3K3hGjLwiC0GPcLMVznOOcXsA1Qf1pDj3/w9r4\nR8Vp7R5Xh4x+p8M7DQ0NpKenc//99/MP//APNDY2hhzX2NjI7NmzeeCBB3jwwQc5ceJEZy8pCILQ\nB3DG8Ntiqv67hDLydSjv37wBXCVgsRdUeucIYG2c+mu6iu39t59OG/38/HzS09P505/+xJQpU8jP\nzw85bunSpXzzm9/kP//zPzl9+jQPPPBAZy8pCILQBzBGPBQx+u8tbImFcfrvkv5zYz8A9APkMLbR\nL8Ru3mKasnSATod34uPjKS0txe12U1tbS1paGlVVVQFj/vKXv5CSkkJ1dfWtJyLhHUEQIoZbvQnE\nBH0ehx3ycbwBTBwC5R0L73Q6ZbOurg63W03c7XZTV9f6yXbu3DnuueceFi5cyB/+8AceeughtmzZ\nwp133tnZywqCIPQDblVVa4x+ettDJg5REg7lbQ8JxU2Nfnp6OrW1ta2Or1u3LuCzy+XSnnogX375\nJRUVFWzdupWHH36YZcuWkZ+fz09+8pM2rlji2PdgqxgJgiBECknYuf5O5U3NwxXQUAJ/AK50/Ndv\navSPHDnS5jkT1hk5ciQXL14kJib4dQRGjx7N6NGjefjhhwGYPXt2m7F/RVq7Ji0IgtC3uZk2/ge0\nlla+hOX9vzcOmAvEQTPAmg5dudMLuTNnzmTnzp0A7Ny5k1mzZrUaM3LkSL761a/ypz/9CYC33nqL\nhISEzl5SEARB6CKdXshtaGhgzpw5nD9/Ho/Hw+7du4mOjqampobFixdz8OBBAP7whz+QnZ3NjRs3\nuO+++9i+fTt3331364n0mYVcH30j7OQj/OfpI/znCDLP7saHzLO9OBd8Y4KOuZVYWwcXcqU4q8OU\n0DfCUCWE/zxLCP85gsyzuylB5hmKm7VGdMb5naH0S8D/Fj19QRCEvsfNMnoutbHfei31VojgmiAI\nQtgRXOAV6oHg9P7bT9iEd9LS0igtLe3taQiCIPQpUlNTKSkpaff4sDH6giAIwu1HwjuCIAgRhBh9\nQRCECKLXjf7hw4eJj48nLi6OgoKC3p5Om3g8HsaNG0dKSgp/93d/19vTsfje976H2+0mKclukNxe\n2eueJNQ88/LyGD16NCkpKaSkpHD48OFenKHi448/ZtKkSSQkJJCYmMgvfvELIPzuaVvzDKd7+vnn\nn/PII4+QnJzMgw8+yI9//GMg/O5lW/MMp3vppLm5mZSUFGbMmAF04n76e5Evv/zSf9999/nPnTvn\nv3Hjhn/8+PH+P/7xj705pTbxeDz+y5cv9/Y0WvHOO+/4Kyoq/ImJidaxH/3oR/6CggK/3+/35+fn\n+1esWNFb07MINc+8vDz/xo0be3FWrbl48aL/5MmTfr/f7//000/9999/v/+Pf/xj2N3TtuYZbvf0\n2rVrfr/f7//iiy/8jzzyiP/48eNhdy/9/tDzDLd7adi4caP/H//xH/0zZszw+/0d//feq57+u+++\ny9/+7d/i8Xi44447+M53vsPrr7/em1O6Kf4wXPN+7LHHGDp0aMCx/fv3M3/+fADmz5/Pvn37emNq\nAYSaJ4TfPR05ciTJyckAREVF8cADD3DhwoWwu6dtzRPC654aRd0bN27Q3NzM0KFDw+5eQuh5Qnjd\nS4BPPvmEQ4cOkZ2dbc2to/ezV43+hQsX+OpXv2p9Hj16tPU/brjhcrmYOnUqEydO5Je//GVvT+em\ntEf2Olx4+eWXGT9+PIsWLer11/xgfD4fJ0+e5JFHHgnre2rm+eijjwLhdU9bWlpITk7G7XZb4ahw\nvJeh5gnhdS8Bli9fzs9//nMGDLBNd0fvZ68a/VByzOHKb3/7W06ePMmbb77JP//zP3P8+PHenlK7\naEv2Ohx45plnOHfuHKdOneLee+/lueee6+0pWTQ1NZGZmcmWLVu46667As6F0z1tampi9uzZbNmy\nhaioqLC7pwMGDODUqVN88sknvPPOOxw7dizgfLjcy+B5lpSUhN29PHDgADExMaSkpLT5BtKe+9mr\nRn/UqFF8/PHH1uePP/6Y0aNH9+KM2ubee+8F4J577iEjI4N33323l2fUNkb2GmhT9jociImJsf4n\nzc7ODpt7+sUXX5CZmUlWVpalHhuO99TMc968edY8w/We3n333XzrW9/i/fffD8t7aTDzLC8vD7t7\n+bvf/Y79+/czZswYvvvd73L06FGysrI6fD971ehPnDiRs2fP4vP5uHHjBv/+7//OzJkze3NKIfns\ns8/49NNPAbh27Rq/+c1vArJQwo32yF6HAxcvXrT29+7dGxb31O/3s2jRIh588EGWLVtmHQ+3e9rW\nPMPpntbX11shkb/+9a8cOXKElJSUsLuXbc3T2UCqt+8lwE9/+lM+/vhjzp07x7/9278xefJkdu3a\n1fH7eduWmNvJoUOH/Pfff7//vvvu8//0pz/t7emEpLq62j9+/Hj/+PHj/QkJCWE1z+985zv+e++9\n13/HHXf4R48e7f/Xf/1X/+XLl/1Tpkzxx8XF+dPT0/1Xrlzp7Wm2mue2bdv8WVlZ/qSkJP+4ceP8\nTzzxhL+2tra3p+k/fvy43+Vy+cePH+9PTk72Jycn+998882wu6eh5nno0KGwuqenT5/2p6Sk+MeP\nH+9PSkryr1+/3u/3+8PuXrY1z3C6l8GUlJRY2TsdvZ8iwyAIghBB9HpxliAIgtBziNEXBEGIIMTo\nC4IgRBBi9AVBECIIMfqCIAgRhBh9QRCECEKMviAIQgQhRl8QBCGC+P8uMYgQoUtOmwAAAABJRU5E\nrkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1086b3350>"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Part (c) Now we use the sample draws from part (b) to calculate the failure probability (i.e., the logistic function) at $T=31$. The probability of failure is $\\sim99\\%$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#PART C: CALCULATE THE PROBABILITY OF FAILURE AT 31F\n",
      "t1=samples[0,:]\n",
      "t2=samples[1,:]\n",
      "logistic_fun= 1./(1 + np.exp(-1*(t1 + t2*31.0))) # Likelihood at 31F\n",
      "plt.hist(logistic_fun,50, histtype='step', label='Failure Probability at 31F')\n",
      "plt.legend()\n",
      "print np.mean(logistic_fun) # Mean probability of failure at 31F"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.98637231124\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtQXOd9//H3ussvvgobYkDapUESa+O1EGziImVau6gU\nXXAsy0MKZlpANvLEqNPK2G3luJMYp62gTScZ2TUdt103stoYeZQESGpt5Kpat74I1wTJHa1rbWws\n2OVSywiELARCPL8/NpxwhCwuQgjFn9fMzi7Pnufs91nY89lz4ziMMQYREZFfuOpyFyAiIvOLgkFE\nRGwUDCIiYqNgEBERGwWDiIjYKBhERMRmSsFw9uxZfD4f99xzDwDV1dW43W58Ph8+n489e/ZY09bU\n1ODxeMjIyGDv3r1We0tLC5mZmXg8HrZs2WK1Dw0NUVxcjMfjYeXKlRw9enS2xiYiIjMwpWDYvn07\nXq8Xh8MBgMPh4NFHH6W1tZXW1lbWrVsHQCgUYteuXYRCIQKBAJs3b2bsNInKykr8fj/hcJhwOEwg\nEADA7/eTmJhIOBymqqqKrVu3XopxiojIFE0aDJFIhJdffplNmzZZC3ljDOc7L66xsZGSkhLi4uJI\nS0sjPT2d5uZmurq6GBgYICcnB4CysjIaGhoAaGpqory8HIDCwkL27ds3a4MTEZHpmzQYqqqq+Pa3\nv81VV/1yUofDwTPPPENWVhYVFRX09fUB0NnZidvttqZzu91Eo9EJ7S6Xi2g0CkA0GiU1NRUAp9NJ\nfHw8vb29szM6ERGZtgsGw09+8hOSkpLw+Xy2NYTKykra2to4ePAgCxcu5LHHHrvkhYqIyNxwXujJ\nN954g6amJl5++WVOnz7NiRMnKCsr44UXXrCm2bRpk7VT2uVy0dHRYT0XiURwu924XC4ikciE9rE+\n7e3tLFq0iJGREfr7+0lISJhQS3p6Ou+///7FjVZE5DNm6dKl/PznP59eJzNFwWDQfOUrXzHGGNPZ\n2Wm1f+c73zElJSXGGGMOHz5ssrKyzNDQkPnggw/MkiVLzOjoqDHGmJycHHPgwAEzOjpq1q1bZ/bs\n2WOMMebZZ581Dz/8sDHGmBdffNEUFxef9/WnUepl9eSTT17uEqbkSqjzSqjRGNU521Tn7JrJsvOC\nawznBIh1VNKf/dmfcejQIRwOB4sXL+a5554DwOv1UlRUhNfrxel0UldXZ/Wpq6tj48aNDA4OUlBQ\nwNq1awGoqKigtLQUj8dDYmIi9fX100s2ERGZVVMOhtzcXHJzcwHYuXPnp073xBNP8MQTT0xo/9KX\nvsT//M//TGj/3Oc+x0svvTTVMkRE5BLTmc+zbCw857sroc4roUZQnbNNdV5+jl9sg5r3HA7Hec+d\nEBGRTzeTZeeUNyWJiMj8dMcd8OGHE9vj4mY2PwWDiMgVLhyG1lZYsMDenpsL3d3Tn5+CQUTkV0Bi\nIsTH29ucM1zCa+eziIjYKBhERMRGwSAiIjYKBhERsVEwiIiIjYJBRERsFAwiImKjYBARERsFg4iI\n2CgYRETERsEgIiI2UwqGs2fP4vP5rGs79/b2kp+fzy233MLq1avp6+uzpq2pqcHj8ZCRkcHevXut\n9paWFjIzM/F4PGzZssVqHxoaori4GI/Hw8qVKzl69OhsjU1ERGZgSsGwfft2vF6vdZnO2tpa8vPz\nOXLkCHl5edTW1gIQCoXYtWsXoVCIQCDA5s2brf8DXllZid/vJxwOEw6HCQQCAPj9fhITEwmHw1RV\nVbF169ZLMU4REZmiSYMhEonw8ssvs2nTJmsh39TURHl5OQDl5eU0NDQA0NjYSElJCXFxcaSlpZGe\nnk5zczNdXV0MDAyQk5MDQFlZmdVn/LwKCwvZt2/f7I9SRESmbNJgqKqq4tvf/jZXXfXLSXt6ekhO\nTgYgOTmZnp4eADo7O3G73dZ0brebaDQ6od3lchGNRgGIRqOkpqYC4HQ6iY+Pp7e3dxaGJiIiM3HB\n/9b9k5/8hKSkJHw+H8Fg8LzTOBwOaxPTpVZdXW09zs3N/ZW+5qqIyEwEg0FreT2Ti/TAJMHwxhtv\n0NTUxMsvv8zp06c5ceIEpaWlJCcn093dTUpKCl1dXSQlJQGxNYGOjg6rfyQSwe1243K5iEQiE9rH\n+rS3t7No0SJGRkbo7+8nISHhvPWMDwYREZlo/Jfmhgbo6Xlq2vO44Kakbdu20dHRQVtbG/X19fzO\n7/wOO3fuZP369ezYsQOAHTt2sGHDBgDWr19PfX09w8PDtLW1EQ6HycnJISUlhQULFtDc3Iwxhp07\nd3LvvfdafcbmtXv3bvLy8qY9CBERmT3TuvDb2Cajxx9/nKKiIvx+P2lpabz00ksAeL1eioqK8Hq9\nOJ1O6urqrD51dXVs3LiRwcFBCgoKWLt2LQAVFRWUlpbi8XhITEykvr5+NscnIiLT5DBjhxrNcw6H\ngyukVBGRORUfD+3tE6/5nJ0Nhw5Nf9mpM59FRMRGwSAiIjYKBhERsVEwiIiIjYJBRERsFAwiImKj\nYBARERsFg4iI2CgYRETERsEgIiI2CgYREbFRMIiIiI2CQUREbBQMIiJio2AQEREbBYOIiNhcMBhO\nnz7NihUryM7Oxuv18vWvfx2IXXvZ7Xbj8/nw+Xzs2bPH6lNTU4PH4yEjI4O9e/da7S0tLWRmZuLx\neNiyZYvVPjQ0RHFxMR6Ph5UrV3L06NHZHqOIiEzDBYPh6quvZv/+/Rw8eJB33nmH/fv389prr+Fw\nOHj00UdpbW2ltbWVdevWARAKhdi1axehUIhAIMDmzZutKwdVVlbi9/sJh8OEw2ECgQAAfr+fxMRE\nwuEwVVVVbN269RIPWURELmTSTUnXXnstAMPDw5w9e5abbroJ4LyXimtsbKSkpIS4uDjS0tJIT0+n\nubmZrq4uBgYGyMnJAaCsrIyGhgYAmpqaKC8vB6CwsJB9+/bNzshERGRGJg2G0dFRsrOzSU5OZtWq\nVdx+++0APPPMM2RlZVFRUUFfXx8AnZ2duN1uq6/b7SYajU5od7lcRKNRAKLRKKmpqQA4nU7i4+Pp\n7e2dvRGKiMi0OCeb4KqrruLgwYP09/ezZs0agsEglZWVfPOb3wTgG9/4Bo899hh+v/+SF1tdXW09\nzs3NJTc395K/pojIlSQYDBIMBgHo7p7ZPCYNhjHx8fHcfffdvP3227YF8qZNm7jnnnuA2JpAR0eH\n9VwkEsHtduNyuYhEIhPax/q0t7ezaNEiRkZG6O/vJyEh4bw1jA8GERGZaPyX5oYG6Ol5atrzuOCm\npGPHjlmbiQYHB3nllVfw+Xx0j4uhH/3oR2RmZgKwfv166uvrGR4epq2tjXA4TE5ODikpKSxYsIDm\n5maMMezcuZN7773X6rNjxw4Adu/eTV5e3rQHISIis+eCawxdXV2Ul5czOjrK6OgopaWl5OXlUVZW\nxsGDB3E4HCxevJjnnnsOAK/XS1FREV6vF6fTSV1dHQ6HA4C6ujo2btzI4OAgBQUFrF27FoCKigpK\nS0vxeDwkJiZSX19/iYcsIiIX4jDnO7xoHnI4HOc9EkpE5LMuPh7a22P342Vnw6FD01926sxnERGx\nUTCIiIiNgkFERGwUDCIiYqNgEBERGwWDiIjYKBhERMRGwSAiIjYKBhERsVEwiIiIjYJBRERsFAwi\nImKjYBARERsFg4iI2CgYRETERsEgIiI2FwyG06dPs2LFCrKzs/F6vXz9618HoLe3l/z8fG655RZW\nr15tXf4ToKamBo/HQ0ZGBnv37rXaW1payMzMxOPxsGXLFqt9aGiI4uJiPB4PK1eu5OjRo7M9RhER\nmYYLBsPVV1/N/v37OXjwIO+88w779+/ntddeo7a2lvz8fI4cOUJeXh61tbUAhEIhdu3aRSgUIhAI\nsHnzZuvKQZWVlfj9fsLhMOFwmEAgAIDf7ycxMZFwOExVVRVbt269xEMWEZELmXRT0rXXXgvA8PAw\nZ8+e5aabbqKpqYny8nIAysvLaWhoAKCxsZGSkhLi4uJIS0sjPT2d5uZmurq6GBgYICcnB4CysjKr\nz/h5FRYWsm/fvtkfpYiITNmkwTA6Okp2djbJycmsWrWK22+/nZ6eHpKTkwFITk6mp6cHgM7OTtxu\nt9XX7XYTjUYntLtcLqLRKADRaJTU1FQAnE4n8fHx9Pb2zt4IRURkWpyTTXDVVVdx8OBB+vv7WbNm\nDfv377c973A4cDgcl6zA8aqrq63Hubm55ObmzsnriohcKYLBIMFgEIDu7pnNY9JgGBMfH8/dd99N\nS0sLycnJdHd3k5KSQldXF0lJSUBsTaCjo8PqE4lEcLvduFwuIpHIhPaxPu3t7SxatIiRkRH6+/tJ\nSEg4bw3jg0FERCYa/6W5oQF6ep6a9jwuuCnp2LFj1hFHg4ODvPLKK/h8PtavX8+OHTsA2LFjBxs2\nbABg/fr11NfXMzw8TFtbG+FwmJycHFJSUliwYAHNzc0YY9i5cyf33nuv1WdsXrt37yYvL2/agxAR\nkdlzwTWGrq4uysvLGR0dZXR0lNLSUvLy8vD5fBQVFeH3+0lLS+Oll14CwOv1UlRUhNfrxel0UldX\nZ21mqqurY+PGjQwODlJQUMDatWsBqKiooLS0FI/HQ2JiIvX19Zd4yCIiciEOM3Y86TzncDi4QkoV\nEZlT8fHQ3h67Hy87Gw4dmv6yU2c+i4iIjYJBRERsFAwiImKjYBARERsFg4iI2CgYRETERsEgIiI2\nCgYREbFRMIiIiI2CQUREbBQMIiJio2AQEREbBYOIiNgoGERExEbBICIiNpMGQ0dHB6tWreL2229n\n2bJlPP3000DsMptutxufz4fP52PPnj1Wn5qaGjweDxkZGezdu9dqb2lpITMzE4/Hw5YtW6z2oaEh\niouL8Xg8rFy5kqNHj87mGEVEZBomDYa4uDi++93vcvjwYQ4cOMCzzz7Lu+++i8Ph4NFHH6W1tZXW\n1lbWrVsHQCgUYteuXYRCIQKBAJs3b7YuElFZWYnf7yccDhMOhwkEAgD4/X4SExMJh8NUVVWxdevW\nSzhkERG5kEmDISUlhezsbACuv/56brvtNqLRKMB5rwrU2NhISUkJcXFxpKWlkZ6eTnNzM11dXQwM\nDJCTkwNAWVkZDQ0NADQ1NVFeXg5AYWEh+/btm53RiYjItE1rH8OHH35Ia2srK1euBOCZZ54hKyuL\niooK+vr6AOjs7MTtdlt93G430Wh0QrvL5bICJhqNkpqaCoDT6SQ+Pp7e3t6LG5mIiMzIlIPh5MmT\nfPWrX2X79u1cf/31VFZW0tbWxsGDB1m4cCGPPfbYpaxTRETmiHMqE505c4bCwkL+4A/+gA0bNgCQ\nlJRkPb9p0ybuueceILYm0NHRYT0XiURwu924XC4ikciE9rE+7e3tLFq0iJGREfr7+0lISJhQR3V1\ntfU4NzeX3NzcqY9UROQzIBgMEgwGAejunuFMzCRGR0dNaWmpeeSRR2ztnZ2d1uPvfOc7pqSkxBhj\nzOHDh01WVpYZGhoyH3zwgVmyZIkZHR01xhiTk5NjDhw4YEZHR826devMnj17jDHGPPvss+bhhx82\nxhjz4osvmuLi4gl1TKFUEZHPpAULjOnrm9ielTWzZeekawyvv/46//Iv/8Ly5cvx+XwAbNu2jRdf\nfJGDBw/icDhYvHgxzz33HABer5eioiK8Xi9Op5O6ujocDgcAdXV1bNy4kcHBQQoKCli7di0AFRUV\nlJaW4vF4SExMpL6+foYxJyIiF8thzHkOLZqHHA7HeY+CEhH5rIuPh/b22P142dlw6ND0l50681lE\nRGwUDCIiYqNgEBERGwWDiIjYKBhERMRGwSAiIjYKBhERsVEwiIiIjYJBRERsFAwiImKjYBARERsF\ng4iI2CgYRETERsEgIiI2CgYREbFRMIiIiM2kwdDR0cGqVau4/fbbWbZsGU8//TQAvb295Ofnc8st\nt7B69Wr6+vqsPjU1NXg8HjIyMti7d6/V3tLSQmZmJh6Phy1btljtQ0NDFBcX4/F4WLlyJUePHp3N\nMYqIyDRMGgxxcXF897vf5fDhwxw4cIBnn32Wd999l9raWvLz8zly5Ah5eXnU1tYCEAqF2LVrF6FQ\niEAgwObNm62rB1VWVuL3+wmHw4TDYQKBAAB+v5/ExETC4TBVVVVs3br1Eg5ZREQuZNJgSElJITs7\nG4Drr7+e2267jWg0SlNTE+Xl5QCUl5fT0NAAQGNjIyUlJcTFxZGWlkZ6ejrNzc10dXUxMDBATk4O\nAGVlZVaf8fMqLCxk3759sz9SERGZkmntY/jwww9pbW1lxYoV9PT0kJycDEBycjI9PT0AdHZ24na7\nrT5ut5toNDqh3eVyEY1GAYhGo6SmpgLgdDqJj4+nt7f34kYmIiIzMuVgOHnyJIWFhWzfvp0bbrjB\n9pzD4cDhcMx6cSIiMvecU5nozJkzFBYWUlpayoYNG4DYWkJ3dzcpKSl0dXWRlJQExNYEOjo6rL6R\nSAS3243L5SISiUxoH+vT3t7OokWLGBkZob+/n4SEhAl1VFdXW49zc3PJzc2d9oBFRH6VBYNBgsEg\nAN3dM5yJmcTo6KgpLS01jzzyiK39T//0T01tba0xxpiamhqzdetWY4wxhw8fNllZWWZoaMh88MEH\nZsmSJWZ0dNQYY0xOTo45cOCAGR0dNevWrTN79uwxxhjz7LPPmocfftgYY8yLL75oiouLJ9QxhVJF\nRD6TFiwwpq9vYntW1syWnQ5jfnHI0Kd47bXXuOuuu1i+fLm1uaimpoacnByKiopob28nLS2Nl156\niRtvvBGAbdu28fzzz+N0Otm+fTtr1qwBYoerbty4kcHBQQoKCqxDX4eGhigtLaW1tZXExETq6+tJ\nS0uz1eFwOJikVBGRz6T4eGhvj92Pl50Nhw5Nf9k5aTDMFwoGEZHzm+1g0JnPIiJio2AQEREbBYOI\niNgoGERExEbBICIiNgoGERGxUTCIiIiNgkFERGwUDCIiYqNgEBERGwWDiIjYKBhERMRGwSAiIjYK\nBhERsVEwiIiIjYJBRERsJg2GBx98kOTkZDIzM6226upq3G43Pp8Pn8/Hnj17rOdqamrweDxkZGSw\nd+9eq72lpYXMzEw8Hg9btmyx2oeGhiguLsbj8bBy5UqOHj06W2MTEZEZmDQYHnjgAQKBgK3N4XDw\n6KOP0traSmtrK+vWrQMgFAqxa9cuQqEQgUCAzZs3W1cOqqysxO/3Ew6HCYfD1jz9fj+JiYmEw2Gq\nqqrYunXrbI9RRESmYdJguPPOO7npppsmtJ/vUnGNjY2UlJQQFxdHWloa6enpNDc309XVxcDAADk5\nOQCUlZXR0NAAQFNTE+Xl5QAUFhayb9++ixqQiIhcnBnvY3jmmWfIysqioqKCvr4+ADo7O3G73dY0\nbrebaDQ6od3lchGNRgGIRqOkpqYC4HQ6iY+Pp7e3d6ZliYjIRXLOpFNlZSXf/OY3AfjGN77BY489\nht/vn9XCzqe6utp6nJubS25u7iV/TRGRK0kwGCQYDALQ3T2zecwoGJKSkqzHmzZt4p577gFiawId\nHR3Wc5FIBLfbjcvlIhKJTGgf69Pe3s6iRYsYGRmhv7+fhISE877u+GAQEZGJxn9pbmiAnp6npj2P\nGW1K6urqsh7/6Ec/so5YWr9+PfX19QwPD9PW1kY4HCYnJ4eUlBQWLFhAc3Mzxhh27tzJvffea/XZ\nsWMHALt37yYvL28mJYmIyCyZdI2hpKSEV199lWPHjpGamspTTz1FMBjk4MGDOBwOFi9ezHPPPQeA\n1+ulqKgIr9eL0+mkrq4Oh8MBQF1dHRs3bmRwcJCCggLWrl0LQEVFBaWlpXg8HhITE6mvr7+EwxUR\nkck4zPkOL5qHHA7HeY+EEhH5rIuPh/b22P142dlw6ND0l50681lERGwUDCIiYqNgEBERGwWDiIjY\nKBhERMRGwSAiIjYKBhERsVEwiIiIjYJBRERsFAwiImKjYBARERsFg4iI2CgYRETERsEgIiI2CgYR\nEbFRMIiIiM2kwfDggw+SnJxsXb4ToLe3l/z8fG655RZWr15NX1+f9VxNTQ0ej4eMjAz27t1rtbe0\ntJCZmYnH42HLli1W+9DQEMXFxXg8HlauXMnRo0dna2wiIjIDkwbDAw88QCAQsLXV1taSn5/PkSNH\nyMvLo7a2FoBQKMSuXbsIhUIEAgE2b95sXTmosrISv99POBwmHA5b8/T7/SQmJhIOh6mqqmLr1q2z\nPUYREZmGSYPhzjvv5KabbrK1NTU1UV5eDkB5eTkNDQ0ANDY2UlJSQlxcHGlpaaSnp9Pc3ExXVxcD\nAwPk5OQAUFZWZvUZP6/CwkL27ds3e6MTEZFpm9E+hp6eHpKTkwFITk6mp6cHgM7OTtxutzWd2+0m\nGo1OaHe5XESjUQCi0SipqakAOJ1O4uPj6e3tndloRETkojkvdgYOhwOHwzEbtUyqurraepybm0tu\nbu6cvK6IyJUiGAwSDAYB6O6e2TxmFAzJycl0d3eTkpJCV1cXSUlJQGxNoKOjw5ouEongdrtxuVxE\nIpEJ7WN92tvbWbRoESMjI/T395OQkHDe1x0fDCIiMtH4L80NDdDT89S05zGjTUnr169nx44dAOzY\nsYMNGzZY7fX19QwPD9PW1kY4HCYnJ4eUlBQWLFhAc3Mzxhh27tzJvffeO2Feu3fvJi8vbyYliYjI\nbDGTuP/++83ChQtNXFyccbvd5vnnnzcff/yxycvLMx6Px+Tn55vjx49b0//VX/2VWbp0qbn11ltN\nIBCw2t9++22zbNkys3TpUvNHf/RHVvvp06fN7/3e75n09HSzYsUK09bWdt46plCqiMhn0oIFxvT1\nTWzPyprZstNhzC+OJ53nHA4HV0ipIiJzKj4e2ttj9+NlZ8OhQ9NfdurMZxERsVEwiIiIjYJBRERs\nFAwiImKjYBARERsFg4iI2CgYRETERsEgIiI2CgYREbFRMIiIiI2CQUREbBQMIiJio2AQEREbBYOI\niNgoGERExEbBICIiNhcVDGlpaSxfvhyfz0dOTg4Avb295Ofnc8stt7B69Wr6+vqs6WtqavB4PGRk\nZLB3716rvaWlhczMTDweD1u2bLmYkkRE5CJdVDA4HA6CwSCtra289dZbANTW1pKfn8+RI0fIy8uj\ntrYWgFAoxK5duwiFQgQCATZv3mxdVaiyshK/3084HCYcDhMIBC5yWCIiMlMXvSnp3EvGNTU1UV5e\nDkB5eTkNDQ0ANDY2UlJSQlxcHGlpaaSnp9Pc3ExXVxcDAwPWGkdZWZnVR0RE5t5FrzH87u/+Lnfc\ncQf/+I//CEBPTw/JyckAJCcn09PTA0BnZydut9vq63a7iUajE9pdLhfRaPRiyhIRkYvgvJjOr7/+\nOgsXLuSjjz4iPz+fjIwM2/MOhwOHw3FRBY5XXV1tPc7NzSU3N3fW5i0i8qsgGAwSDAYB6O6e2Twu\nKhgWLlwIwM0338x9993HW2+9RXJyMt3d3aSkpNDV1UVSUhIQWxPo6Oiw+kYiEdxuNy6Xi0gkYmt3\nuVznfb3xwSAiIhON/9Lc0AA9PU9Nex4z3pR06tQpBgYGAPjkk0/Yu3cvmZmZrF+/nh07dgCwY8cO\nNmzYAMD69eupr69neHiYtrY2wuEwOTk5pKSksGDBApqbmzHGsHPnTquPiIjMvRmvMfT09HDfffcB\nMDIywu///u+zevVq7rjjDoqKivD7/aSlpfHSSy8B4PV6KSoqwuv14nQ6qaurszYz1dXVsXHjRgYH\nBykoKGDt2rWzMDQREZkJhzn3sKJ5yuFwTDgCSkREID4e2ttj9+NlZ8OhQ9NfdurMZxERsVEwiIiI\njYJBRERsFAwiImKjYBARERsFg4iI2CgYRETERsEgIiI2CgYREbFRMIiIiM1F/XdVERGZG8aA3w+f\nfDLxuaGh2X0tBYOIyBXg+HH4wz+Ehx+e+Nxjj8H118/eaykYRESuENddB9u3X/rX0T4GERGx0RqD\niMg80tcHTz4JIyP29tOn566GebPGEAgEyMjIwOPx8Nd//deXuxwRkcvigw/ghz8Er9d+++IX4Z/+\naW5qmBcX6jl79iy33nor//7v/47L5eI3fuM3ePHFF7ntttusaa6UC/UEg0Hreqvz2ZVQ55VQI6jO\n2farVufoKHz88cT2UAj+4i/gc5+zt/f3x9YO3n774muc6YV65sWmpLfeeov09HTS0tIAuP/++2ls\nbLQFw5XiV+2P+nK6EmoE1Tnb5nudJ0/CsWPQ0BAkLS3Xat+/Hw4cgKvO2Q7zyivw/vvw+c/b248d\ng7vugqqqia/xhS/Mft3TMS+CIRqNkpqaav3sdrtpbm6+jBWJyMUYHo7dznXmTOx2rpMnYzeAnh54\n553Y49OnYwvaX1we3nLqFAwMTFwIDw3BRx9N/BZ+9iyEw+c/pPOnP4Ubbpg4r6am2P0119jbBwfh\n6qtjr9HQ8Mv2jz+Gr3wltrAfb/lyWLUKMjImvvZ8NS+CwXHub/0yGxiA6mp47z17++gonDgBN930\n6X3few9aWiZ/jRMnYn9wznN+A/3952+frtOnYyfC3HCDvd2Y2DeV/n549VX7c11dcOONEz8gH38c\n+yD8v/9nbz91KlbnudeZ/TT9/RAXN/FD+2na2uB735v4fo+N4dd/fWKfSCQ2hl/7tYmvvWABXHut\nvf3//i82v3PbP82JE7H3Z/z0//d/sROPenshIcE+/chI7PjzpKSJY+jqgvT0ia/x0Udw880Tx/DR\nR7H379wF1fHjsftzF3pjC9qx38/x4/DCC7E6BwYmvq8DA7GF+bm//9HR2H1i4sRaxzaRnPsRNiZW\n67l/M598Eqv/3Fp7e2PfqJOSYsHw+uux9sHB2OvfeKN9+rG/PZfL3j48HAuBc/82zp6NfSZuvXXi\nGLKyYgvuc8fw0EPw27898ffgcMTGUF0du81nt90Ghw7NoKOZB958802zZs0a6+dt27aZ2tpa2zRL\nly41gG666aabbtO4LV26dNrL5Hmx83lkZIRbb72Vffv2sWjRInJycibsfBYRkbkxLzYlOZ1O/u7v\n/o41a9ZGlbeoAAAGeklEQVRw9uxZKioqFAoiIpfJvFhjEBGR+WPenOA2ZrIT3f71X/+VrKwsli9f\nzm/+5m/yztjhC/OoxsbGRrKysvD5fHzpS1/iP/7jP+a8Rpj6SYP//d//jdPp5Ic//OEcVvdLk9UZ\nDAaJj4/H5/Ph8/n4y7/8y8tQ5dTez2AwiM/nY9myZZftkMvJ6vzbv/1b673MzMzE6XTS19c37+o8\nduwYa9euJTs7m2XLlvG9731vzmuEyes8fvw49913H1lZWaxYsYLDhw/PeY0PPvggycnJZGZmfuo0\nf/zHf4zH4yErK4vW1tYLz/Ci9hrPspGREbN06VLT1tZmhoeHTVZWlgmFQrZp3njjDdPX12eMMWbP\nnj1mxYoV867GkydPWo/feeedGe38uVhTqXNsulWrVpm7777b7N69e17WuX//fnPPPffMeW3jTaXO\n48ePG6/Xazo6Oowxxnz00Ufzss7xfvzjH5u8vLw5rDBmKnU++eST5vHHHzfGxN7LhIQEc+bMmXlX\n55/8yZ+Yb33rW8YYY/73f//3sryf//mf/2l+9rOfmWXLlp33+X/7t38z69atM8YYc+DAgUmXm/Nq\njWH8iW5xcXHWiW7jffnLXyb+F8ffrVixgkgkMu9qvO6666zHJ0+e5PPnntkyB6ZSJ8AzzzzDV7/6\nVW6++eY5rxGmXqe5zFs8p1Ln97//fQoLC3G73QDz+vc+5vvf/z4lJSVzWGHMVOpcuHAhJ06cAODE\niRMkJibivNjjuC9Bne+++y6rVq0C4NZbb+XDDz/ko48+mtM677zzTm66wHH0TU1NlJeXA7HlZl9f\nHz09PZ86/bwKhvOd6BaNRj91er/fT0FBwVyUZplqjQ0NDdx2222sW7eOp59+ei5LBKZWZzQapbGx\nkcrKSuDynE8ylTodDgdvvPEGWVlZFBQUEAqF5rrMKdUZDofp7e1l1apV3HHHHezcuXOuy5zWZ+jU\nqVP89Kc/pbCwcK7Ks0ylzoceeojDhw+zaNEisrKy2D4X/2/6HFOpMysry9oM+9Zbb3H06NE5/8I6\nmfON40I1zoujksZMZ8G0f/9+nn/+eV4fOxNmjky1xg0bNrBhwwb+67/+i9LSUt4792y5S2wqdT7y\nyCPU1tZa/4fqcnwrn0qdX/ziF+no6ODaa69lz549bNiwgSNHjsxBdb80lTrPnDnDz372M/bt28ep\nU6f48pe/zMqVK/F4PHNQYcx0PkM//vGP+a3f+i1uPPfssTkwlTq3bdtGdnY2wWCQ999/n/z8fA4d\nOsQN5561eQlNpc7HH3+cLVu2WPtsfD4fv3buWXHzwLmf7wuNbV4Fg8vloqOjw/q5o6PDWi0f7513\n3uGhhx4iEAhccPXpUphqjWPuvPNORkZG+Pjjj0k836mjl8hU6mxpaeH+++8HYjv69uzZQ1xcHOvX\nr59XdY5fEKxbt47NmzfT29tLwrmnGV/mOlNTU/n85z/PNddcwzXXXMNdd93FoUOH5jQYpvP3WV9f\nf1k2I8HU6nzjjTf48z//cwCWLl3K4sWLee+997jjjjvmVZ033HADzz//vPXz4sWLWbJkyZzVOBXn\njiMSieA697Tx8WZ1D8hFOnPmjFmyZIlpa2szQ0ND593Rc/ToUbN06VLz5ptvztsaf/7zn5vR0VFj\njDEtLS1myZIl87LO8TZu3Gh+8IMfzGGFMVOps7u723o/m5ubzRe+8IV5Wee7775r8vLyzMjIiPnk\nk0/MsmXLzOHDh+ddncYY09fXZxISEsypU6fmtL4xU6mzqqrKVFdXG2NifwMul8t8/PHH867Ovr4+\nMzQ0ZIwx5h/+4R9MeXn5nNY4pq2tbUo7n998881Jdz7PqzWGTzvR7bnnngPga1/7Gt/61rc4fvy4\ntV08Li6Ot956a17V+IMf/IAXXniBuLg4rr/+eurr6+esvunUOR9Mpc7du3fz93//9zidTq699tp5\n+35mZGSwdu1ali9fzlVXXcVDDz2E1+udd3VCbB/YmjVruObcf7w0j+p84okneOCBB8jKymJ0dJS/\n+Zu/mdO1xKnWGQqF2LhxIw6Hg2XLluH3++e0RoCSkhJeffVVjh07RmpqKk899RRnfvHfCr/2ta9R\nUFDAyy+/THp6Otdddx3//M//fMH56QQ3ERGxmVdHJYmIyOWnYBARERsFg4iI2CgYRETERsEgIiI2\nCgYREbFRMIiIiI2CQUREbP4/bzuoe4wLHa8AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10868b9d0>"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Part (d): We repeat (c) using 5000 temperatures randomly drawn from a normal distribution, and use a standard vanilla MC. Here for speed we use the posterior draws directly from part (b). The mean failure probability, based on the available data is $\\sim 31\\%$, with 68% of the failures within $\\pm 4.5\\%$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#PART D: Failure probability for T~N(68 ,1)\n",
      "Nmc=5000\n",
      "temp = np.random.normal(68,1,Nmc)\n",
      "fail_prob = np.zeros(Nmc)\n",
      "for n in np.arange(0,Nmc): #Vanilla MC\n",
      "    logistic_fun= 1./(1 + np.exp(-1*(t1 + t2*temp[n])))\n",
      "    fail_prob[n] = np.mean(logistic_fun)\n",
      "    #Technically , we should redraw t1,t2 from the Posterior for each loop using MCMC, \n",
      "    #but this would take an eternity, and the variance and mean are dominated by the \n",
      "    #distribution of temperatures, so the difference between results is negligible.\n",
      "\n",
      "plt.hist(fail_prob,30, histtype='step', label=u'Failure Probability for T~N(68,1)')\n",
      "print np.mean(fail_prob), np.var(fail_prob)\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.324475172322 0.00341114118785\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 41,
       "text": [
        "<matplotlib.legend.Legend at 0x10a06b950>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xtc1HW+x/HXIKioqCgyKJik4OJYKGmJmruYopkbBy/L\neo1crVa33dPaeWTp9vBydnO6bJ7ccncrU7QsLCvYMvJUQhdNWttybXQV1wsXQQ1Q8Aoy5w+Ok6jM\nDJcZ8Mf7+XjMQ5jf9/f7fX5f4T0/vr+byW632xEREcPyaeoCRETEsxT0IiIGp6AXETE4Bb2IiMEp\n6EVEDE5BLyJicG4FfXh4ONHR0cTExHDbbbcBUFxcTHx8PH379mXMmDGUlpY62i9fvpzIyEiioqLY\nsmWLZyoXERG3uBX0JpOJzMxM/vGPf5CdnQ2A1WolPj6effv2MWrUKKxWKwA2m43U1FRsNhsZGRnM\nmzePqqoqz22BiIg45fbQzZXXVaWnp5OcnAxAcnIy7777LgBpaWlMnToVPz8/wsPDiYiIcHw4iIiI\n97m9Rz969GgGDx7MSy+9BEBRURFmsxkAs9lMUVERAAUFBYSFhTnmDQsLIz8/v7HrFhERN/m60+iL\nL76ge/fuHD9+nPj4eKKiompMN5lMmEymWud3Nk1ERDzLraDv3r07AN26dWPChAlkZ2djNpspLCwk\nJCSEo0ePEhwcDEBoaCi5ubmOefPy8ggNDa2xvIiICA4cONBY2yAi0iL06dOHnJycOs/ncujmzJkz\nlJWVAXD69Gm2bNnCzTffTEJCAikpKQCkpKSQmJgIQEJCAm+88QYXLlzg4MGD7N+/33GmziUHDhzA\nbrfrZbezePHiJq+hubzUF+oL9YXzV313kF3u0RcVFTFhwgQAKisrmT59OmPGjGHw4MEkJSWxevVq\nwsPD2bhxIwAWi4WkpCQsFgu+vr6sWrVKQzciIk3IZdDfeOONfPPNN1e936VLFz766KNrzrNw4UIW\nLlzY8OpERKTBdGVsE4uLi2vqEpoN9cUP1Bc/UF80nMlut3v9wSMmk4kmWK2IyHWtvtnp1lk3cv3r\n0qULJSUlTV2GiLghMDCQ4uLiRlue9uhbCPW5yPWjtt/X+v4ea4xeRMTgFPQiIganoBcRMTgFvVy3\n7rrrLtavXw/A2rVrGTFiRBNXVHdLlixh5syZ9ZrX1TY765+AgAAOHTpUr/Ve6V//+hcDBw6kY8eO\nPP/8842yTG85f/48/fv3d9yU0Zvr7devHydOnPDK+hT0LViXLmAyee7VpYv7tYSHh9OuXTsCAgII\nCAigY8eOFBYWOp1n8+bN9Q7Jhli7di2tWrUiICCATp06ERMTw/vvv1+vZXnyqnFn/VNWVkZ4eDgA\n9957L48//ni91/PUU08xatQoTp06xYMPPljv5QD88pe/dPwMtGnThtatWzu+Hz9+vNvLOXToED4+\nPlfNM2PGDJYuXer4/sUXX+QnP/mJ4068AF9//TU//vGPCQgIICQkhJUrVzqmfffdd8TFxdG5c2d6\n9uzJ73//+1pr2Lp1KyNHjqRz587ceOONNaa1adOGX/ziF47neHiagr4FKykBu91zr7qczWkymXjv\nvfcoKyujrKyMU6dOERIS4rmNv8zFixfrPM/w4cMpKyujtLSU2bNnk5SUxMmTJ+u8bCOcCXX48GEs\nFku95r2yf/7yl784fgYWLlzIlClTHN/X58M0Ozub7du3O76/8k67f/3rX2t8GJ44cYJx48Yxd+5c\niouLOXDgAGPGjHFMnzlzJiNGjKCkpISsrCxWrVrF3/72t2uuu0OHDsyZM4enn376mtOnTp1KSkoK\nFRUVdd6uulLQS7NVWlrKT3/6U4KDg+nSpQt33313jWcbxMXFsXr16qvmu7Q3d/mTzS5vu3btWoYP\nH878+fMJCgpi6dKlXLhwgf/6r/+iV69ehISEMHfuXM6dO1drbZcC2mQyMWvWLM6ePUtOTg5Llixh\n8uTJzJw5k06dOpGSkkJBQQEJCQl07dqVyMhIXn75ZcdyTCYT586dY8qUKXTs2JFBgwaxa9cux3Sr\n1UpERAQdO3akf//+jgf8XF7Hr3/9azp37ky/fv345JNPXPYPgI+PDwcOHODFF19kw4YNPPXUUwQE\nBJCQkMAzzzzD5MmTa7T/zW9+w0MPPXTVcu644w4yMzN58MEH6dixIzk5OZw8eZJ77rmH4OBgwsPD\n+cMf/uDor2v1vbM+bugH4SOPPMKiRYuuWi7AkSNH+Pe//82QIUMc05599lnuvPNOx8OT2rdvX+O2\n7DabjenTp2Mymejduze33347Npvtmuu+9dZbmT59+lV785eEhYURGBhY44PIUxT00mxc+UtdVVXF\n7NmzOXLkCEeOHMHf37/G0ICr5yBc7sq22dnZ9OnTh2PHjrFw4UIWLFhATk4O3377LTk5OeTn57Ns\n2TKXy62srOTll18mICCAvn37AtVPX/vZz37GyZMnmTZtGlOmTOGGG27g6NGjvPXWWyxcuJCtW7c6\ntjktLY2kpCRKSkqYNm0aiYmJjj3diIgIPv/8c06dOsXixYuZMWNGjfHkHTt2EBERwffff8/SpUuZ\nOHGi4/nN7jwn4v7772f69OksWLCAsrIy0tPTmTFjBhkZGY6/UCorK0lNTXU8Ue5yn3zyCSNGjOCF\nF17g1KlTRERE8Otf/5qysjIOHjxIVlYW69atY82aNbX2vSfNnTuXffv28fHHH9fYboB//vOf9O7d\nGx+fH2Jwx44dBAYGMnz4cMxmMwkJCTVuuz5mzBhSUlKorKxk7969bN++ndGjR9e7vn79+vHtt9/W\ne353KeilWbDb7SQmJhIYGEhgYCATJ06kS5cuTJgwgbZt29KhQwcWLlxIVlZWo6yvR48e/OpXv8LH\nx4c2bdrw0ksv8eyzz9K5c2c6dOjAY489xhtvvFHr/F9++SWBgYF0796d1NRU3nnnHQICAgAYNmwY\nCQkJABw/fpxt27bx5JNP0rp1awYMGMCcOXNYt26dY1mDBw9m4sSJtGrVivnz53Pu3DnHXt7kyZMd\nQ1hJSUlERkayY8cOx7zBwcH853/+J61atSIpKYkf/ehHvPfee3Xuj8s/ZENCQhgxYgRvvvkmABkZ\nGXTr1o2YmBiX81+8eJHU1FSWL19O+/bt6dWrFw8//LDjoDDU7Pu2bdvWudY1a9YwcOBAQkJCuOee\ne/j0008pKSlh3bp1PPfcczXatmvXjkWLFvG73/3uqu0sLS11/J9dkpubS0pKCitXruTIkSPceOON\nTJ061TF9xYoVpKam4u/vj8ViYc6cOQwaNKjO23BJQECA44PZkxT00iyYTCbS0tIoKSmhpKSEt99+\nmzNnzvDAAw8QHh5Op06d+MlPfsLJkycbZVy7Z8+ejq+PHz/OmTNnGDRokOODZty4cU7PiIiNjaWk\npMQR5HfccYdj2uWP0iwoKKBLly60b9/e8d4NN9xQYwjq8vYmk4mwsDCOHj0KwLp164iJiXHUtXv3\nbr7//ntH+ysf6tOrVy/HvA2RnJzMq6++CsCrr77q8qD3pb3kEydOUFFRQa9evRzTrtzey/u+PjZv\n3syHH37I/v37GT58OL/97W/p378/W7duZcaMGVe1nz17NkVFRbz33ns1riwNDAx0PGvjknbt2jFx\n4kQGDRpEmzZtWLx4Mdu2baOsrIwzZ85wxx13sGzZMs6fP09ubi4ZGRn8+c9/rve2lJWVERgYWO/5\n3aWgl2brj3/8I/v27SM7O5uTJ0+SlZXl1rjtpVA9c+aM470rz+C5fEgjKCgIf39/bDab44OmtLSU\nU6dO1bnmK4dLevToQXFxMeXl5Y73jhw5UiPcLx8aqKqqIi8vjx49enD48GHuv/9+XnjhBYqLiykp\nKeGmm26qsf1XPo/58OHD9OjRo841X+k//uM/2LVrF7t37+b9999n+vTpbi0rKCgIPz+/GqduXrm9\ndRluu5aNGzdiNpsJCAjggQceYOfOnRQUFLBmzRq6du16VfvWrVuzePFiHn/88Rp9Fx0dzcGDB2sc\ny4mOjq61nu+++46ysjJmzJiBj48PoaGh/PznP2fz5s1ubc+17NmzhwEDBtR7fncp6KXZKi8vx9/f\nn06dOlFcXOz0wN3lunXrRmhoKOvXr+fixYu88sorTp/M4+Pjw3333cdDDz3E8ePHgeoA3bJlS51r\nvvJDqGfPngwbNozHHnuM8+fPs2vXLl555ZUae547d+7knXfeobKykv/5n/+hbdu2xMbGcvr0aUwm\nE0FBQVRVVbFmzRp2795dY/nHjh1j5cqVVFRU8Oabb7J3717uuuuuOtVsNpv597//XeM9f39/Jk2a\nxLRp0xgyZEiNoHa23ZeGkBYtWkR5eTmHDx9mxYoV19zTdqW2D/T6nJI6c+ZMzp07R0ZGhuO9sLAw\nIiIiagyFzZo1i3feeYdvv/2WiooK/vu//5sRI0YQEBBAREQEFy5c4PXXX6eqqorCwkJSU1MdQX3p\nJIAjR4446j937hwVFRXY7XbOnz/PhQsXHOvKz8+nuLiY2NjYOm9PXSnoW7DAQM+eR9/Qv0gfeugh\nzp49S1BQEMOGDWPcuHG1/pJfuSf90ksv8fTTTxMUFITNZmP48OG1tgV48skniYiIIDY2lk6dOhEf\nH8++ffvcWperaa+//jqHDh2iR48eTJw4kWXLljmGekwmE4mJiaSmptKlSxdee+013n77bVq1aoXF\nYuHhhx9m6NChhISEsHv3bm6//fYa64qNjWX//v1069aNxx9/nE2bNl1zKODKui7/evbs2dhsNsex\nkUuSk5PZvXu3W9cqXL68P/3pT7Rv357evXszYsQIpk+fzqxZs1z2naua6+ryeX18fFi2bBklJSU1\n3n/ggQdqHD8YOXIkTzzxBOPHj3d8AG7YsAGoHup58803efrppwkMDCQmJobo6GjH+H9ubi7h4eGO\n4bSsrCzatWvH+PHjyc3Nxd/fnzvvvNOxrg0bNnDvvffi5+dX7210l+5e2UKoz6WucnNziYqKoqio\niA4dOjR1OR5x4cIFYmJi+OSTT2pcNFUff/jDHwgODua+++5z2fb8+fMMHDiQzz77jKCgoKumN/bd\nKxX0LYT6XOqiqqqK+fPnU15eXuO8f/GOxg56PXhERGo4ffo0ZrOZG2+8scaYtly/tEffQqjPRa4f\nevCIiIjUiYJeRMTgFPQiIgang7EtRGBgoEfvfS4ijaexb4ugg7HSLBw6BHFx1f96w4IF1Q9GWbDA\nO+sTaQw6GCsiItekoRvxuNOn4aOPqp86VZtjx7xXj0hLo6AXj/vwQ/jVr+CyB/lc02W3/RaRRqSg\nF4+z22HYMNi0qakrEWmZNEYvImJwCnoREYNT0IuIGJyCXkTE4BT0IiIGp6AXETE4t4L+4sWLxMTE\ncPfddwNQXFxMfHw8ffv2ZcyYMZSWljraLl++nMjISKKiour1cGUREWlcbgX9c889h8VicdwUy2q1\nOh6ePGrUKKxWKwA2m43U1FRsNhsZGRnMmzePqqoqz1UvIiIuuQz6vLw8Nm/ezJw5cxw300lPTyc5\nORmoflL8u+++C0BaWhpTp07Fz8+P8PBwIiIiyM7O9mD5IiLiisug/+1vf8vTTz+Nj88PTYuKihxP\nTDebzRQVFQFQUFBAWFiYo11YWBj5+fmNXbOIiNSB01sgvPfeewQHBxMTE0NmZuY125hMJqf3Oa9t\n2pIlSxxfx8XFERcX57JYEZGWJDMzs9bsrQunQb9t2zbS09PZvHkz586d49SpU8ycOROz2UxhYSEh\nISEcPXqU4OBgAEJDQ8nNzXXMn5eXR2ho6DWXfXnQi4jI1a7cCV66dGm9luP2g0eysrJ45pln+Nvf\n/sYjjzxC165dWbBgAVarldLSUqxWKzabjWnTppGdnU1+fj6jR48mJyfnqr16PXikZdm0CTZsaF43\nNVuwAD79FG65xXm7+HhITPROTSKu1Dc763T3ykuB/eijj5KUlMTq1asJDw9n48aNAFgsFpKSkrBY\nLPj6+rJq1So9vk6apVmz4IYbnLfZuRNef11BL9c/PUpQPK457tG7IzUV3n67+l+R5kCPEhQRkWtS\n0IuIGJyCXkTE4BT0IiIGp6AXETE4Bb2IiMEp6EVEDE5BLyJicAp6ERGDU9CLiBicgl5ExODqdFMz\nkSsVFcFXXzlv8/e/e6cWEbk2Bb00yPPPw8aNEBnpvN1Pf+qdekTkagp6qdX778OXXzpvs3Ur3HMP\nLFrknZpEpO4U9FKr55+HoCDo27f2NmPHwoQJ3qtJROpOQS9OTZsG48Y1dRUi0hA660ZExOAU9CIi\nBqegFxExOAW9iIjB6WCsSC38/CAjAwYMcN5u4EBISfFOTSL1oaAXqUVCAnz2Gdjttbc5fBgefth7\nNYnUh4JepBa+vhAd7bxNu3beqUWkITRGLyJicAp6ERGDU9CLiBicgl5ExOAU9CIiBqegFxExOAW9\niIjBKehFRAxOQS8iYnAKehERg1PQi4gYnIJeRMTgnAb9uXPnGDJkCAMHDsRisfDYY48BUFxcTHx8\nPH379mXMmDGUlpY65lm+fDmRkZFERUWxZcsWz1YvIiIuOQ36tm3bsnXrVr755ht27drF1q1b+fzz\nz7FarcTHx7Nv3z5GjRqF1WoFwGazkZqais1mIyMjg3nz5lFVVeWVDRERkWtzeZvidv9/H9YLFy5w\n8eJFAgMDSU9PJysrC4Dk5GTi4uKwWq2kpaUxdepU/Pz8CA8PJyIiguzsbGJjYz27FSJNqLQU/vIX\n5206dIDp08Fk8k5NIpdzGfRVVVXccsstHDhwgLlz59K/f3+Kioowm80AmM1mioqKACgoKKgR6mFh\nYeTn53uodJGmFxoKU6fCN984b7d6dfWDTDp29E5dIpdzGfQ+Pj588803nDx5krFjx7J169Ya000m\nEyYnuynOpolc79q1g5UrXbfbsMHztYjUxu0nTHXq1Inx48ezc+dOzGYzhYWFhISEcPToUYKDgwEI\nDQ0lNzfXMU9eXh6hoaHXXN6SJUscX8fFxREXF1e/LRARMajMzEwyMzMbvByT3V77EzFPnDiBr68v\nnTt35uzZs4wdO5bFixfz4Ycf0rVrVxYsWIDVaqW0tBSr1YrNZmPatGlkZ2eTn5/P6NGjycnJuWqv\n3mQy4WS10kyMGwe/+U31v9IwHTtCXp6GbqRh6pudTvfojx49SnJyMlVVVVRVVTFz5kxGjRpFTEwM\nSUlJrF69mvDwcDZu3AiAxWIhKSkJi8WCr68vq1at0tCNiEgTc7pH77GVao/+uqA9+sajPXppDPXN\nTl0ZKyJicAp6ERGDU9CLiBicgl5ExOAU9CIiBqegFxExOAW9iIjBKehFRAxOQS8iYnAKehERg1PQ\ni4gYnIJeRMTgFPQiIganoBcRMTgFvYiIwSnoRUQMTkEvImJwCnoREYNT0IuIGJyCXkTE4HybugBp\nGnl5UFLivE1ZmXdqERHPUtC3UEOGQMeO4OvkJ8DHB3r29F5NIuIZCvoWqqICsrIgOLipKxERT9MY\nvYiIwSnoRUQMTkEvImJwCnoREYNT0IuIGJyCXkTE4BT0IiIGp6AXETE4XTAl4iWVldUvZ5xdqSxS\nX/qxEvGC4GDXVyHb7fDVV3DLLd6pSVoOBb2IF+TkuG4TFwenTnm8FGmBNEYvImJwCnoREYNzGfS5\nubmMHDmS/v37c9NNN7Fy5UoAiouLiY+Pp2/fvowZM4bS0lLHPMuXLycyMpKoqCi2bNniuepFRMQl\nl0Hv5+fHihUr+O677/jyyy954YUX2LNnD1arlfj4ePbt28eoUaOwWq0A2Gw2UlNTsdlsZGRkMG/e\nPKqqqjy+ISIicm0ugz4kJISBAwcC0KFDB/r160d+fj7p6ekkJycDkJyczLvvvgtAWloaU6dOxc/P\nj/DwcCIiIsjOzvbgJoiIiDN1GqM/dOgQ//jHPxgyZAhFRUWYzWYAzGYzRUVFABQUFBAWFuaYJyws\njPz8/EYsWURE6sLtoC8vL2fSpEk899xzBAQE1JhmMpkwmUy1zutsmoiIeJZb59FXVFQwadIkZs6c\nSWJiIlC9F19YWEhISAhHjx4l+P+vBgkNDSU3N9cxb15eHqGhoVctc8mSJY6v4+LiiIuLa8BmiIgY\nT2ZmJpmZmQ1ejslut9udNbDb7SQnJ9O1a1dWrFjheP+RRx6ha9euLFiwAKvVSmlpKVarFZvNxrRp\n08jOziY/P5/Ro0eTk5NTY6/eZDLhYrXiYcHBsHu3nhnbnMTFwZIl1f+KXEt9s9PlHv0XX3zBq6++\nSnR0NDExMUD16ZOPPvooSUlJrF69mvDwcDZu3AiAxWIhKSkJi8WCr68vq1at0tCNiEgTcrlH75GV\nao++yWmPvvnRHr24Ut/s1JWxIiIGp6AXETE4Bb2IiMEp6EVEDE5BLyJicAp6ERGD0xOmRJqRjAw4\ndMh5mzFjoEcPr5QjBqGgN5jSUujfH86dc97u1Clo08Y7NYl7pk2DbdugsLD2NtnZUFwM8+d7ry65\n/inoDaa8HKqqYN8+5+1at4Yr7k0nTez++6tfzijgpT4U9AbUqhV07drUVYhIc6GDsSIiBqegFxEx\nOAW9iIjBKehFRAxOQS8iYnAKehERg1PQi4gYnIJeRMTgFPQiIganoBcRMTgFvYiIwSnoRUQMTkEv\nImJwCnoREYNT0IuIGJyCXkTE4PTgEZHrzJ498MEHzttYLNCrl3fqkeZPQS9yHRk+HF5+GVaurL1N\nQQFER8P69d6rS5o3Bb3IdWTSpOqXM6+9Bps3e6ceuT5ojF5ExOAU9CIiBqegFxExOAW9iIjBKehF\nRAxOQS8iYnAKehERg3MZ9L/4xS8wm83cfPPNjveKi4uJj4+nb9++jBkzhtLSUse05cuXExkZSVRU\nFFu2bPFM1SIi4jaXQT9r1iwyMjJqvGe1WomPj2ffvn2MGjUKq9UKgM1mIzU1FZvNRkZGBvPmzaOq\nqsozlYuIiFtcXhk7YsQIDh06VOO99PR0srKyAEhOTiYuLg6r1UpaWhpTp07Fz8+P8PBwIiIiyM7O\nJjY21iPFtzT798Ojj4LdXnubs2fBZPJeTSLS/NXrFghFRUWYzWYAzGYzRUVFABQUFNQI9bCwMPLz\n8xuhTIHqm1kdOQKPPea8Xffu3qlHRK4PDb7XjclkwuRkF7K2aUuWLHF8HRcXR1xcXENLaRG6d4eJ\nE5u6ChHxhszMTDIzMxu8nHoFvdlsprCwkJCQEI4ePUpwcDAAoaGh5ObmOtrl5eURGhp6zWVcHvQi\nInK1K3eCly5dWq/l1Ov0yoSEBFJSUgBISUkhMTHR8f4bb7zBhQsXOHjwIPv37+e2226rV2EiItI4\nXO7RT506laysLE6cOEHPnj1ZtmwZjz76KElJSaxevZrw8HA2btwIgMViISkpCYvFgq+vL6tWrXI6\nrCMiIp5nstudncPhoZWaTDTBaq976enVD51IT2/qSqQ5u3Q/+tdea+pKpLHVNzt1ZayIiMHpCVMi\nBvT997Bzp/M2wcHQs6d36pGmpaAXMZg+feDECbj//trbnD0LrVrBP//pvbqk6SjoRQwmNhb+/nfn\nbfbs0fUYLYnG6EVEDE5BLyJicAp6ERGDU9CLiBicgl5ExOB01k0z8fnnsHat8zaHD4O/v1fKERED\nUdA3E//7v3DsGCQk1N4mNhYGDvReTSJiDAr6ZmTQIJgzp6mrkJaitNT1X5EdO+p8eyNQ0Iu0QKGh\ncPfd4OqZFuvXQ3m5hgyvdwp6kRaoY0d48UXX7VJTPV+LeJ7OuhERMTgFvYiIwSnoRUQMTmP0IuLU\n/Png6yQpWrWC3/0OgoK8V5PUjYJeRGq1dm319R3OrFgBkyfD7bd7pSSpBwW9iNTq5z933UZn5jR/\nGqMXETE4Bb2IiMEp6EVEDE5BLyJicAp6ERGD01k3XpCeXn2/eWc+/xzGjvVOPSLSsijoveCll6Br\nV7BYam+TmOj8XvQiIvWloPeSSZOqbwsrIuJtGqMXETE47dGLiMdt2wbvv++63ejRMHKk5+tpaRT0\nItJgOTnVDzOpzfPPw8mTMGxY7W2+/hry8xX0nqCgF5EGufVWePZZ1+2eegruvLP26WvXun60odSP\ngl5EGmTFiqauQFxR0DfQ+fNw7pzzNhUV3qlF5Hp34UL1EI8z/v7QurV36jEKBX0DRUdDQQH4ODl/\nyccHFi/2Xk0i1yOzufqA7Q031N7m4kW4+WbYvt17dRmBR06vzMjIICoqisjISJ588klPrKLZKC+H\nvXur90Jqe5WUwNChTV2pSPM2bpzz36OTJ6sDvry8qSu9/jT6Hv3Fixd58MEH+eijjwgNDeXWW28l\nISGBfv36NfaqDCEzM5O4uLimLqNZUF/8QH3xgyv7orgY1qxxPk/nzjBhgmfrup40etBnZ2cTERFB\neHg4AFOmTCEtLe26DPrz56v/VHSmqqph69Av9A/UFz9QX/zg8r644QYYPx4+/dT5PCkp1X9J+/nV\n3sbHB9q2bbw6m7NGD/r8/Hx69uzp+D4sLIwdO3Y09moaJC8PnnjCeUiXlcGGDdUHfpxp2xbat2/c\n+kTk2jp1ghdfdN0uMxO6d3fe5uzZ6hsJ/v8+aYMkJMBddzV8OZ7S6EFvMpkae5GN7tCh6gOozs7p\nBZgyRfenEbkeHTzouk1mZvXxtYbauxd27mzeQY+9kW3fvt0+duxYx/dPPPGE3Wq11mjTp08fO6CX\nXnrppVcdXn369KlXLpvsdrudRlRZWcmPfvQjPv74Y3r06MFtt93G66+/fl2O0YuIGEGjD934+vry\n/PPPM3bsWC5evMjs2bMV8iIiTajR9+hFRKR58ej96F1dOLV3716GDh1K27Zt+eMf/+jJUpqcq754\n7bXXGDBgANHR0QwfPpxdu3Y1QZXe4aov0tLSGDBgADExMQwaNIhPPvmkCar0DncvLvzqq6/w9fXl\n7bff9mJ13uWqLzIzM+nUqRMxMTHExMTw+9//vgmq9A53fi4yMzOJiYnhpptucn0qbr2PurpQWVlp\n79Onj/2KdWwMAAAD2UlEQVTgwYP2Cxcu2AcMGGC32Ww12hw7dsz+1Vdf2RctWmR/5plnPFVKk3On\nL7Zt22YvLS212+12+wcffGAfMmRIU5Tqce70RXl5uePrXbt21fsAVHPnTl9cajdy5Ej7+PHj7W+9\n9VYTVOp57vTF1q1b7XfffXcTVeg97vRFSUmJ3WKx2HNzc+12u91+/Phxp8v02B795RdO+fn5OS6c\nuly3bt0YPHgwfs6uajAAd/pi6NChdOrUCYAhQ4aQl5fXFKV6nDt90f6yCxPKy8sJCgrydple4U5f\nAPzpT39i8uTJdOvWrQmq9A53+8LeAkaa3emLDRs2MGnSJMLCwgBc/o54LOivdeFUfn6+p1bXrNW1\nL1avXs1dzfqk3Ppzty/effdd+vXrx7hx41i5cqU3S/Qad/oiPz+ftLQ05s6dC1wf16nUhzt9YTKZ\n2LZtGwMGDOCuu+7CZrN5u0yvcKcv9u/fT3FxMSNHjmTw4MGsX7/e6TI9dvdKo/5A1kdd+mLr1q28\n8sorfPHFFx6sqOm42xeJiYkkJiby2WefMXPmTP71r395uDLvc6cvHnroIaxWKyaTCbvdbtg9Wnf6\n4pZbbiE3N5d27drxwQcfkJiYyL59+7xQnXe50xcVFRV8/fXXfPzxx5w5c4ahQ4cSGxtLZGTkNdt7\nLOhDQ0PJzc11fJ+bm+v4M6Olcbcvdu3axX333UdGRgaBgYHeLNFr6vpzMWLECCorK/n+++/p2rWr\nN0r0Gnf6YufOnUyZMgWAEydO8MEHH+Dn50dCQoJXa/U0d/oiICDA8fW4ceOYN28excXFdOnSxWt1\neoM7fdGzZ0+CgoLw9/fH39+fH//4x3z77be1Br3HDsZWVFTYe/fubT948KD9/PnztR5ostvt9sWL\nFxv6YKw7fXH48GF7nz597Nu3b2+iKr3Dnb7IycmxV1VV2e12u33nzp323r17N0WpHleX3xG73W6/\n99577Zs2bfJihd7jTl8UFhY6fi527Nhh79WrVxNU6nnu9MWePXvso0aNsldWVtpPnz5tv+mmm+zf\nffddrcv02B59bRdO/fWvfwXggQceoLCwkFtvvZVTp07h4+PDc889h81mo0OHDp4qq0m40xfLli2j\npKTEMRbr5+dHdnZ2U5btEe70xaZNm1i3bh1+fn506NCBN954o4mr9gx3+qKlcKcv3nrrLf785z/j\n6+tLu3btWvTPRVRUFHfeeSfR0dH4+Phw3333YbFYal2mLpgSETE4j14wJSIiTU9BLyJicAp6ERGD\nU9CLiBicgl5ExOAU9CIiBqegFxExOAW9iIjB/R+nrtf7E6EeowAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a0034d0>"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}