Logistic regression prediction interval

var.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "pd.options.display.max_rows = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>survived</th>\n",
       "      <th>pclass</th>\n",
       "      <th>sex</th>\n",
       "      <th>age</th>\n",
       "      <th>sibsp</th>\n",
       "      <th>parch</th>\n",
       "      <th>fare</th>\n",
       "      <th>embarked</th>\n",
       "      <th>class</th>\n",
       "      <th>who</th>\n",
       "      <th>adult_male</th>\n",
       "      <th>deck</th>\n",
       "      <th>embark_town</th>\n",
       "      <th>alive</th>\n",
       "      <th>alone</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>male</td>\n",
       "      <td>22</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>7.2500</td>\n",
       "      <td>S</td>\n",
       "      <td>Third</td>\n",
       "      <td>man</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Southampton</td>\n",
       "      <td>no</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>female</td>\n",
       "      <td>38</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>71.2833</td>\n",
       "      <td>C</td>\n",
       "      <td>First</td>\n",
       "      <td>woman</td>\n",
       "      <td>False</td>\n",
       "      <td>C</td>\n",
       "      <td>Cherbourg</td>\n",
       "      <td>yes</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>female</td>\n",
       "      <td>26</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7.9250</td>\n",
       "      <td>S</td>\n",
       "      <td>Third</td>\n",
       "      <td>woman</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Southampton</td>\n",
       "      <td>yes</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>female</td>\n",
       "      <td>35</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>53.1000</td>\n",
       "      <td>S</td>\n",
       "      <td>First</td>\n",
       "      <td>woman</td>\n",
       "      <td>False</td>\n",
       "      <td>C</td>\n",
       "      <td>Southampton</td>\n",
       "      <td>yes</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>male</td>\n",
       "      <td>35</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8.0500</td>\n",
       "      <td>S</td>\n",
       "      <td>Third</td>\n",
       "      <td>man</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Southampton</td>\n",
       "      <td>no</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   survived  pclass     sex  age  sibsp  parch     fare embarked  class  \\\n",
       "0         0       3    male   22      1      0   7.2500        S  Third   \n",
       "1         1       1  female   38      1      0  71.2833        C  First   \n",
       "2         1       3  female   26      0      0   7.9250        S  Third   \n",
       "3         1       1  female   35      1      0  53.1000        S  First   \n",
       "4         0       3    male   35      0      0   8.0500        S  Third   \n",
       "\n",
       "     who adult_male deck  embark_town alive  alone  \n",
       "0    man       True  NaN  Southampton    no  False  \n",
       "1  woman      False    C    Cherbourg   yes  False  \n",
       "2  woman      False  NaN  Southampton   yes   True  \n",
       "3  woman      False    C  Southampton   yes  False  \n",
       "4    man       True  NaN  Southampton    no   True  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = (sns.load_dataset('titanic')\n",
    "         .dropna(subset=['survived', 'age']))\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.672429\n",
      "         Iterations 4\n"
     ]
    }
   ],
   "source": [
    "mod = sm.Logit.from_formula('survived ~ age', df)\n",
    "res = mod.fit()\n",
    "me = res.get_margeff()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>     <td>survived</td>     <th>  No. Observations:  </th>  <td>   714</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>               <td>Logit</td>      <th>  Df Residuals:      </th>  <td>   712</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>MLE</td>       <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Fri, 11 Sep 2015</td> <th>  Pseudo R-squ.:     </th> <td>0.004445</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>08:16:34</td>     <th>  Log-Likelihood:    </th> <td> -480.11</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>           <td>True</td>       <th>  LL-Null:           </th> <td> -482.26</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                      <td> </td>        <th>  LLR p-value:       </th>  <td>0.03839</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0567</td> <td>    0.174</td> <td>   -0.327</td> <td> 0.744</td> <td>   -0.397</td> <td>    0.283</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>       <td>   -0.0110</td> <td>    0.005</td> <td>   -2.057</td> <td> 0.040</td> <td>   -0.021</td> <td>   -0.001</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:               survived   No. Observations:                  714\n",
       "Model:                          Logit   Df Residuals:                      712\n",
       "Method:                           MLE   Df Model:                            1\n",
       "Date:                Fri, 11 Sep 2015   Pseudo R-squ.:                0.004445\n",
       "Time:                        08:16:34   Log-Likelihood:                -480.11\n",
       "converged:                       True   LL-Null:                       -482.26\n",
       "                                        LLR p-value:                   0.03839\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0567      0.174     -0.327      0.744      -0.397       0.283\n",
       "age           -0.0110      0.005     -2.057      0.040      -0.021      -0.001\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Marginal Effects</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th> <td>survived</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>          <td>dydx</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>At:</th>             <td>overall</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "   <th></th>      <th>dy/dx</th>    <th>std err</th>      <th>z</th>      <th>P>|z|</th> <th>[95.0% Conf. Int.]</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th> <td>   -0.0026</td> <td>    0.001</td> <td>   -2.081</td> <td> 0.037</td>      <td>   -0.005</td>    \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "        Logit Marginal Effects       \n",
       "=====================================\n",
       "Dep. Variable:               survived\n",
       "Method:                          dydx\n",
       "At:                           overall\n",
       "=========================================================================\n",
       "                dy/dx    std err          z      P>|z| [95.0% Conf. Int.]\n",
       "-------------------------------------------------------------------------\n",
       "age           -0.0026      0.001     -2.081      0.037             -0.005\n",
       "=========================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "me.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Greene's notation here (~page 736). We've got a few things\n",
    "\n",
    "- Exogenous: $X$, a single row is $x.T$\n",
    "- estimated coef $\\hat{\\beta}$\n",
    "- predicted prob $\\hat{F} = \\Lambda(\\bf(x)^T\\hat{\\bf{\\beta}})$\n",
    "where $\\Lambda$ is the logistic CDF\n",
    "- Variance of $\\hat{\\beta} = V$, from statsmodels\n",
    "\n",
    "We want a prediction interval for the predicted probability around the CEF. We have\n",
    "\n",
    "$Asy. Var[\\hat{F}] = [\\partial \\hat{F} / \\partial \\hat{\\beta}]\\bf{V} [\\partial \\hat{F} / \\partial \\hat{\\beta}]$\n",
    "\n",
    "And the vector of derivatives is\n",
    "\n",
    "$[\\partial \\hat{F} / \\partial \\hat{\\beta}] = [d\\hat{F}/dz][\\partial z / \\partial \\hat{\\beta}] = \\hat{f} x$\n",
    "\n",
    "where $\\hat{f} = \\hat{\\Lambda}(1 - \\hat{\\Lambda})$ (the logistic pdf I beleive).\n",
    "\n",
    "And so\n",
    "\n",
    "$Asy.Var[\\hat{F}] = \\hat{f}^2 x^T\\bf{V}x$\n",
    "\n",
    "Note that the $x$ here is a single *row* from the matrix. But vectors are always column vectors so $x$ is $Kx1$, where $K$ is the number of predictors (2 for us)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "Λ = lambda x: stats.logistic().cdf(x)\n",
    "λ = lambda x: stats.logistic().pdf(x)\n",
    "\n",
    "β_ = res.params.values.reshape(-1, 1)\n",
    "V_ = res.cov_params().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# this is the equation from above for Var(F)\n",
    "# only works for a single observation though\n",
    "def var_π(x, β, V_):\n",
    "    #      λ(z)**s * x.T @ V_ @ x\n",
    "    prob = λ(x.T.dot(β))**2 * x.T.dot(V_).dot(x)  # maybe this is it?\n",
    "    return prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# for multiple observations, but not quite vectorized fully\n",
    "\n",
    "def var_πs(xx, β, V_):\n",
    "    α = λ(xx.dot(β))**2\n",
    "    out = np.empty((500, 1))\n",
    "    for i, x in enumerate(xx):\n",
    "        out[i] = x.T.dot(V_).dot(x)\n",
    "    return α * out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Making some fake data.\n",
    "xx = sm.add_constant(np.linspace(df.age.min(), df.age.max(), 500).reshape(-1, 1))\n",
    "πs = Λ(xx.dot(β_))\n",
    "vv = np.sqrt(var_πs(xx, β_, V_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x108b62438>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAECCAYAAADq7fyyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwW2l63/nve87BlSDAG3i/iZT06tIttVqa6Z7uGcfe\nZJzsJlOZOKlKuVKx3ZvZcpxUypWkymtPZZ3KbipJrWu2yk52tpxxNmPHjv+YVCbZxJuJK44TZ3p6\n1PeWWmodtUjxBt6vIAkQBHBO/ngBElJ3i2I3SVz4fKpUIgCKeAQQPxy8532fV/m+jxBCiMZhVbsA\nIYQQR0uCXQghGowEuxBCNBgJdiGEaDAS7EII0WAk2IUQosE4T7pRa20B3wSuADnga67rjlXc/jng\nG4ACUsBPua67e3zlCiGEOMhBR+xfBYKu674E/CImxAHQWivgnwI/47rul4A/AM4cV6FCCCGezkHB\n/jLwPQDXdW8CNypuOw+sAH9La/1fgBbXdd3jKFIIIcTTOyjY40C64nKxNDwD0AG8BPxj4E8Af1xr\n/WNHX6IQQojDOCjY00Bz5fe7ruuVvl4BHrhGAXNkf+PxHyCEEOJkPfHkKfAq8BXgO1rrF4FbFbeN\nAzGt9WjphOqXgN940g/zfd9XSn2WeoUQ4jQ6VHCqJzUBK50gLc+KAXgFuA7EXNf9Vmno5R+V7vRV\n13X/5gH35y8tbR6mvqpIJpuROo+O1Hm06qHOeqgR6qrOQwX7E4/YXdf1gZ977Or7Fbf/IfDCYe5Q\nCCHE8ZIFSkII0WAk2IUQosFIsAshRIORYBdCiAYjwS6EEA1Ggl0IIRqMBLsQQjQYCXYhhGgwEuxC\nCNFgJNiFEKLBSLALIUSDkWAXQogGc+LB7vkeP5j9Pus7ayd910IIcSpU5Yh9p7DDe0vvcmf5fQpe\noRolCCFEw6raUIxjOazn1nht9lWm05PVKkMIIRpO1cfYLWUxsTHB63M32dhZr3Y5QghR96oe7ACW\nZZH3dnl36R0ZnhFCiM+oJoK9rDw888O5HzCTnqp2OUIIUZdqKtjLFIqHGw95Y16GZ4QQ4rBqMtjB\nDM/sFs3wzN3lOzI8I4QQT6lmg73MsRzWcqv8cO4HMntGCCGeQs0He5lC7c2ekcVNQgjxyeom2GF/\n9kx5cVO+mK92SUIIUXPqKtjLKmfPTKUn8X2/2iUJIUTNqMtgL7OUxWR6gjfmb7K2s1rtcoQQoibU\ndbCDCfe8l+fW0nu8v3SLXCFX7ZKEEKKq6j7YyxzLYWN3g5tzr/FwY1yGZ4QQp1bDBHuZbdnMbE5z\nc+41lrPL1S5HCCFOXMMFO5jhmaJf5O7y+7y3+A7ZQrbaJQkhxIlpyGAvsy2brfwWb8z/kA/X7uP5\nXrVLEkKIY+dUu4CTYCuHhe15ljKLnEmM0BPrrXZJQghxbBr6iL2SUgofn/trLm8tvEk6l652SUII\ncSxOTbCXOZbDTiHLO4tvS3MxIURDOnXBXuZYNmu5VV6bfZXJjQmZHimEaBinNtjLLGUxuTnB63M/\nlOmRQoiGcOqDHcBWNgW/wN3l93l38W0y+Uy1SxJCiE9Ngr2Cbdls57d5bfo17q18QNErVrskIYQ4\ntFMx3fGwbMtmObvEcnaJ4fgw/fHBapckhBBPTY7YP4FSCqXM3quvz/2Q1exKtUsSQoinIsF+ALO5\nR57by7d4b/EdGX8XQtQ8Cfan5FgOW/kt3py/yQcy/i6EqGES7IdkWw4r2SV+MPt9mf8uhKhJTzx5\nqrW2gG8CV4Ac8DXXdccqbv+bwF8BlkpX/azruvePqdaaoZRCoZjcnGBue5aRlrN0RjurXZYQQgAH\nz4r5KhB0XfclrfULwDdK15U9D/xl13XfOa4Ca5mtbIp+kXsrd5jenOJcyznioUS1yxJCnHIHDcW8\nDHwPwHXdm8CNx26/Dnxda/3ftNa/eAz11QW7ov+MbM8nhKi2g4I9DlS2QSyWhmfKfhf4WeB/AL6o\ntf7TR1xfXdnbnm/+B9xfc+UEqxCiKg4aikkDzRWXLdd1K3er+FXXddMAWuvfA64Bv/ekH5hMNtOS\njhJyQp+m3hPT0hL9TP9+19/kbuZtzrScYbhlGKXUEVX2qGSy+eBvqgFS59GqhzrroUaonzoP46Bg\nfxX4CvAdrfWLwK3yDVrrBHBLa30JyGCO2v/ZQXe4tLTJ+kaGkF27R7MtLVHW149mvvrb6+9zW7kM\nJ84c+QYfyWQzS0ubR/ozj4PUebTqoc56qBHqq87DOCjYvwt8WWv9aunyK1rrnwRirut+qzSu/oeY\nGTP/yXXd7x224EZnKxsPjw/X7jOzNc1I4iztkfZqlyWEaGBPDHbXdX3g5x67+n7F7b+LGWcXB7At\nm93iLneWbxMLxjjbcp54KF7tsoQQDUgWKJ0w27LJFrK8s/gWtxffkxYFQogjJ8FeJY7lkM6neWP+\nh9xdviNTJIUQR6YqwS6r8Pc5VoC13Co3517DXb0ne7AKIT6zE+/H/i/+432+f1sx0FdgeNBioE8R\nDB7PVMB6Uu4Bv5CZp7eplzOJUWzLrnZZQog6dOLBfmm4lbfHUkxM+UxMFbEs6OtRDA9aDA4owqHT\nHfK2spnfnmdue46+WD/DiTNYSkbMhBBP78SD/bpOko74ZDYDTE57TEx5TKd8plNF1A+hp0sxNKgY\nHrCIRk9nyJebjM1upZjdnqUv1sdQfFgCXgjxVKqyNZ5S0NaqaGu1uXbFJp32mZj2mJj0mZ03f157\n3aMrWQr5QYvm2OkL+fJq1ZnNaWa3UvTHBhiID0rACyGeqCb2PI3HFVcu21y5DFvbPpNTHhNTPvOL\nPgtLPq+/5dHeBsODFsODFi2J0xXy5SCf2pxkZmt6L+CFEOLj1ESwV4o1KS5ftLl8EbJZn8lpn4kp\nj9l5n5VVj7fe9WhJ7Id8WyvH1oel1lQG/PTWFFecC8T8DjmCF0I8ouaCvVIkorhwXnHhvEUu5zM1\nY0I+Nevz7m2Pd297NDfD8IDF8KAi2aFORciXg3xyY5KN9Xv0xHoZjp+RWTRCCKDGg71SKKQ4N6o4\nN2qRz/tMp/y9E6+373rcvgvRaCnkhxRdSYVlNXbIW8pCKcXc1iyzWym6mroZSYziWHXztAohjkFd\nJkAgoBgZVowMWxQKPqk5E/JTMz53XY+7LoRDMDRgTrz2dCtsu3FDvjyLZimzyNzWLJ1NXZyJjxB2\nwtUuTQhRBXUZ7JUcRzE0oBgasPA8M6NmYspnctrDfeDjPigSDMJgv/me/l6F4zRuyDuWw2p2hcXt\nedrCHQzHz9Acarx+00KIT1b3wV7JshT9vYr+Xnjp8xYLSz6TU+Zo/sG4z4PxIo4DA30m5Af6FcFA\nY4a8YwVI727w1uKbJIJxBuPD0i5YiFOioYK9kmUperoUPV3wwg2L5RW/tNrV4+Gkz8PJIrYFfb2K\noUGLoX5FqAFXvQYsh0whw53lW4SdKL2xXvpi/afiJLMQp1XDBnslpcyMmWQH3LhmsbYOE1Pe3rj8\n1EyR7yvo6VYMDyquXPYO/qF1xrYc8t4u4xtjTKQf0hXtZjh+hoAdqHZpQogjdiqCvZJSirZWaGu1\nef6qzUbaL4W8z+yc+fODm5t0dZqQHx60iDU1ztGtrcyUyMXMArNbKdrC7QzGB0mEWqpcmRDiqJy6\nYH9cIq64+ozN1Wdgc8ucdJ1JKVJzRRYWfW6+6dHRvh/yiXjjhLxjOaR3N3h38R2iThPdsW76Yv2y\n4EmIOnfqg71Sc0zxzEWbL34hyuzctmlSNukzt+CzvOLz5jserS37q15bWxpj1atjOex6OSY2HjKx\n8ZCOaJKB2CCxYKzapQkhPgUJ9k8QjSgunre5eB5yOZ/JGZ+JSY/UnM87tzzeueURby6HvKKjvf5X\nvZaP1FezKyxsL9AciNHV1E1vrE+O4oWoIxLsTyEUUpwfVZwftdjd3V/1OpPyuXXH49YdaIruh3xn\nA6x6DVgOO8UdHm6MM74xRnu4g96mXlojbdUuTQhxAAn2QwoGFaNnFKNnzKrXmdlSa4MZnzv3PO7c\ng0gYhkr9a3q66zvky0fq67k1VrJLBOwQyUiSvlg/kUCkytUJIT6OBPtn4Dj7J1WLRZ+5eRPyk9M+\n9z70uPchBIMw1G++p7dX4dRxawPbcvD8IguZeaY3p4kHY7RFOuiPDci0SSFqiAT7EbFtRX+for/P\n4qUXfBYWSwuipj0+HPf5cLxIwIH+PsWZIdPaIFDHq16DdoCdYo7U5gyT6QniwTgdkSS9sT5pQiZE\nlckr8BhYlhmC6emGFz9nsbT8Mate7Yq9Xut41atSioAKkC1kmUpPMrb+gEQoQWu4jd6mPkJOqNol\nCnHqSLAfM6XMydTOJHzueYvVtY+uelXK7PU6PGh62NTrXq9KKYJ2kGwhS2ZzhsmNCWLBGCOqn0Au\nRnMoXu0ShTgVJNhPkFKK9jZob7O5/pzN+sb+mHx5r9cfvO7RmTQdK4cHLOJ1uiBKKUXADpAr5pjb\nmmNpdYOgFaQl1EpbpI1kpFM2BhHimEiwV1FLQvHcszbPPbu/1+vktNnrdXHJ54239xdEDQ3U9zaA\nASuAj89abpXlnSXurX5APNBMPNRCMpIkHkrU7f9NiFojwV4jHtnrdcdnurQN4GzFgqjmmJlGOTSo\n6Oyo3xC0lY2tbLLFHbKZeVKbM1iWRTwYpzkYpyPSQTyUqHaZQtQtCfYaFAkrzp9VnD9rsZv3mUmZ\nHjbTKZ/3P/B4/wMIh+H8aJaeLq/ud4hybPNruJXfYiu/xfTmFJayiAViNAfjtIbaaI20yupXIZ6S\nBHuNC1ZsA1gslneI8pia9rl1Z5dbdyAYKG0eMlj/0yiBvemSmUKGTCFDamsG8IkGYsQCTcQCMdoj\nSaKBaHULFaJGSbDXEdtWDPQpBvosvBd8MtkQ79/NMjntMTbhMzaxP41yqDSNMlyn0ygrlYN+t5hj\ntZhjJbvC2PoDbMuhyYkSDcZoCjTRFmonGojKWL049STY65RlKfr7HGJNNi/csFhZxXSjfGwaZXeX\nYri0J2xTg/SVNzNuggBmnD67w3JmiTHvQyxlEQk0EXWiRJ0osWAzLaEWWRkrThUJ9gaglKKjHTra\nzTTK8uYhk9OmzcHcvM9rb3gk2xVDDdhXHh4N+91ijt1ijvXcGsXNIgW/SNAO0ORECdkRIk6EpmCM\nRChByJYFVKLxSLA3oMrNQ7Yz5sTr5JTpK79U6ivfkig3KrNob6vfaZQHsS0bGzNfPlPIkilkWctB\ncatIoVggYDmEnDBhJ0Iv7WTSRZqcGM2hZoKlNwoh6o0Ee4NriiouaZtLGnZy+9MoU3M+773v8d77\nHrGm/WmUXQ3Qcvhp2MrGdkzg5708+d08C9sF1jczFLwCRd8joGyCTpiQHSRshwnaIUJ2iHgoTsSJ\nSk8cUbPkN/MUCYcU50YV50Yt8nnTcnhy2ozJl1sOh0MwOFDqRlnn0yg/Lcdy9l4YBS9Pwcuznd8G\nwPd9Cl4BD38v+IN2gJAVImiHCNpBQk6YmBMj7IRlda2oCgn2UyoQUJwZMp0mi0UzTDM5ZYL+/gOf\n+w+KBALQ32tCvr9PEazzaZRHodwqoawc/Bkye9cV/SLFYgGUwrEcE/iWQ8AOEbACBK0gATtAU6CJ\niBMlaAdljr44UhLswrQc7lX098IXPl/qRjltWhzsdaO0oLeiG2U4LCH/SSqHeWA//ClkH/m+glfA\n8z3AfEoIWEGCtoNjBQlaAfO3HSBgBYg6UcKBCAFLZveIg0mwi0dYlqKrU9HVCZ8vdaOcnPb2Vr5O\np8w0yq5So7LBAYt4s4T8p/H4GH3RL5AtFICdR673fZ+CX8D3fQA6tuNspfMErQC25RCwTPibNweH\noBMm6kRKnxSCDXtiXHwyCXbxiSq7UT5/1Sa9+WijsvlFn5tvmUZlQwOmUVkjz7CplnLP+zJLWXh+\nkZ1iEYof/X7P9yh45o1AKXCsAI4VIGA55vyBFSCgHOyK68J2hLATJmgHcSxHnsM698Rg11pbwDeB\nK0AO+JrrumMf833/FFhxXfeXjqVKURPizYpnL9s8exmyWZ+pGTMmPzvn8+5tj3dvezRFyyGv6O46\nHTNsao2lrI9M1fT8IrlikVwx97H/pugVKfrmXUJhtkH8pDcDx7JxLIeQHSbiRPbeDOQ8Qe046Ij9\nq0DQdd2XtNYvAN8oXbdHa/2zwDPAfzmWCkVNikQU+pxCn3t0hs10yueu63HXNfu9DvQpLl3I05rw\n676HTSOrnO9fduCbgV+k6BX3PxlUBH8yl2AzvVt6M3D2hokCVpCQEyLsRPaGj8TRO+hRfRn4HoDr\nuje11jcqb9RavwR8Hvh14MKxVChqXuUMG8/zmV/wmZw2QT/20GfsYWbv5OvQgDn5GolIyNc7W9nY\n9qNvBj4eec8jU8iwuZv5yL/ZO1/geaAsFP7eUFH5k8D+p4PypwXz6SDsmCmlASsgnw4OcFCwx4F0\nxeWi1tpyXdfTWvcAvwz8OeAvHleBor5YlqK3R9HbY/Z7XVmFhSUb98Pc3slXePTka6O1NxCfbO98\nwWO5XPQLFIuFAz8dlNnK2jtHELBKnxSUXXqDkE8HB/1P00BzxWXLdV2v9PVfADqA/x/oBqJa6w9c\n1/2tJ/1AhSIZSeLh4/lFPM8zT5pv/vZLY30+4APK91HKwir9kZM69aPcw+bsaJjLFzxz8nXanHxd\nXPJZWPJ5/e399gZDA4qOdiXPsfiIj/t0AL5ZNezlgexH/k3504HneSZDANsO4Kjym4FDV7GF9Eau\n9GZgwj/iROr+zUCVp1B9HK31TwBfcV33Fa31i8D/5rrun/6Y7/tp4MJTnDz95Dt7TPlkTtErkvfy\n5Ao5dou7e3N/zbJvc3v5ex+5rnzZK+7NFVZKYStb3iBqQCbrMf6wwIPxPJNTBQqlg7FYk2J0JMDZ\nEYeBPudUrnwV1eP7/l6GKKVQlBeZBffCP2CXPhWU/g7awdI6gzABK3Bcq40P9UI4KNgV+7NiAF4B\nrgMx13W/VfF9Pw1o13W/fsD9+UtLm4ep70gUvSIFv0C+WHqD8HKlRSPlN48Ceb9QetMoEIuHWFpN\nU/TyFH1v7+fYltnSrVa0tERZX//oOGatOajOfN4nNeczVWpvkNs11wfKG4gMmA1EgsHjDflGeTxr\nQT3UCEdT5+PrDMrDREErgGPvnzNwKtYbhO0IkcD+jKKDJJOHWyzyxGA/BlUJ9sNKJpsp1+n5Hnkv\nz24xR6aQZbeYo1A0bwIFL196Q9gl7xXJF3fNWKHv7b3TH+dJnkZ88Xiez8Li/snXLdOiBcuC3u7S\nuHy/RTR69CHfiI9ntdRDjVC9Os1IQwEqhohMuwln/02h4o3g2r+4FPD/rl942p9fnwNIJ8hSFqFS\nV7/mYPzA7/d9M+6XK+6wnc+UeoPnKXi77Hp58sVdcsVdCt7uI28AMjRkWJaip1vR0w0v3Hh05evM\nrJlW+epNj2SHCfmhAYuWhDx2or7Ylv2RIZu91hMfc74Acx5z5ml/vgT7EVNKEbSDBO3ggW8E+WKe\nnWKWrd0tdgo77Hq75k8hx04xR97LgW8Wi5zG6V2Pr3zd3Cp1oyytfF1aNr3lE/H9k6/JDjn5KoQE\nexUFbHMi5pPeADzfY6eww2Z+k2w+Q66ww46XM38XsuQKOYp+sabG/Y9Tc0zxzEWbZy7u95YvH8nf\nuuNx6w5EwjDQb47ke7sVjiMhL04fCfYaZimLaCBKNBD92NvjrUHGUim281vsFLJki1kyeXMeQCnV\n0J0AK3vLFwrm5OvktMf0jL/XdnhvY+8Bi4E+WRQlTg8J9joWckIko0mSJB+5vugVSe9usJHbIFvI\nkC1kyOSz5L18Q/b0cJz98XbPM0M0ZkPv/Y29ATqTisHS0XwiLs3KROOSYG9AtmXTGm6jNdz2yPW5\nYo6VzDJb+S0y+W22C9vkirsErUDDhFxl2+HPPW829p6aMePyC0tmYdSb73jEm2Gw32Jw4PRsByhO\nDwn2UyRkh+ht7nvkunLYp/Nptnc32c5v4/vg2I3xq5GIK569ZPPspdK4fMrMl0/N+rz/gcf7H+w3\nK7uoTbOy454vL8Rxa4xXr/jUymHfiwl83/dJ5zZY2Vlha3eTdD5NoZgn8Fgb2HoUDinOjSjOjexv\nBzg1bY7oy83KLAt6uhSDpfnysSYJeVF/JNjFI5RSJMItJMIte9dt7W6xnF0indswQe8V6v7E7OPb\nAa6swuKyzf0HOVJz5mTsa697tLeVhmz6ZRMRUT8k2MWBYsEYsWBs7/JmLs1idpGN3Dqbu2mUsup6\nymVls7JL2mNre//E69y8z8qqxzu3zCYiA/0WQ/1mEZX0sRG1SoJdHFpzKE5zyMy993yPxcwiazsr\nrOXWyRdyOHZ9H83HmhSXtM0lDbu7FVMpUz737nvcuw8BB/p6TdvhgT5FOCQhL2qHBLv4TCxl0d3U\nTXdTN2CGbRYyc6zvrJHe3SJQp21Py4LBRzcRWVgy4/KTMx4TUz4TU/ube5fH5aW/vKi2+n7ViZpj\nhm3OAaZlwtz2LF4gy1Jxo+574liWoqdL0dMFn79usZE2fWymZvY39379LdNffqDf7BTV2SFTKcXJ\nk2AXxyZgBxiMD5FMNtPvrDG3lWIlu8J6br3uQ14pRUsCWhI2V58xm3tPp8yQTWrO5/Ydj9t3IBSC\n/l5zJH8SrYeFAAl2cUIcy2EgPsRAfIiCV2B2K1WaaZMmUOdj8mA29z5/VnH+rGlxMDtXmjOfKk+l\nNEM23Z2KgX4ZshHHS4JdnDjHchiMDzEYHyJXyDGzNc1KdoVsfrvuT7yCaXEwOKAYHICXfIuVNZgu\nz7JZMH9ef8t0pRzoM0M2XZ0yZCOOjgS7qKqQE2K05SyjLWfZzKWZ2Z5hJbuM7/lYVv33tFFK0dEG\nHW02165AJuszkzLTKVNzFatfA6VZNv0W/TLLRnxGEuyiZjSH4lwMXcL3fRYy88xvz7GeW6/7xVCV\nopVDNkWf+YX9hmUPJ30eTpohm86kYqDPBH1LQhZGicORYBc1RylFd1MP3U095It5pjYnWdxeYNfb\nrdtd4z+OU7n69XMWa+ulIZuU2R5wYdE0LGuOlWbZ9Cm6uyTgxcEa51UiGlLADuwN1axkV0htTbOa\nXW2IE66VlFK0tUJbq83VZyG7Y4ZsplNmI5G79zzu3jMLo4aHtunu8hnolR7z4uNJsIu60R5ppz3S\nXjqKn2Ahs0ixWGiIsfjHRcL7G4kUi+bofWrGzLL5cKzAh2Pm+5Idpsf8QJ9FW6sM2QhDgl3UHXMU\nf46RxFkWM4uktqZJ727W/SrXT2Lbit4eRW+P2eDbJ8ydu5m9IZulZZ+33t3vZTPYbxZSybaAp1dj\nvhLEqaCUoqupi66mLrbz20ylJ1nMLGAru2GPXJVStLbYPHvZ5tnLkMv5zMz6TM94TM/u97Ipbws4\n0Kfo75P2w6eNBLtoCE2BJi62X+J8q2Zqc5L5rTnyXh7bqt+uk08jFFKMnlGMnjG9bBaXzJDNdKq8\nLaAPeLS2mDnzA32KTtkxquFJsIuGYls2ZxIjnEmMML89R2pzhq3dzYZY+HQQyzKzZrq74PPXbdKb\nPjMp05Vybt7n1rrHrTtmxyhzNG/aHMgJ2MYjwS4aVnnKZDq3wWR6kpWd5YaaE3+QeLPi0gWbSxcw\nbQ7mTZuD6dT+nHmAZLuiv8+0OuhoUw07jHWaSLCLhhcPJXg2eYWdwg4T6YcsZhawsE5VgDmOmT0z\n2A++b7G+gZlKmTJdKZdWfN65BeEwDPQqBvotensUIWlaVpck2MWpEXbCXGi7yLmW80ymJ5jbnqXo\nF+t696dPw5yAhdYWmyuX9zcTKQf9h+M+H46X+sx3lk7A9lq0tsh0ynohwS5OHduyGWkZ5UxihNTW\nDDObM+wWc9Uuq2oqNxPxfZ+VVXM0P50yLQ/mF3zeeNujqWn/BGxPlyIQkJCvVRLs4tRSStHfPEB/\n8wDz2/NknFWWvXRDtS04rPL+rx3tpmlZdscnNbu/AnZvOqUFPd2lsfk+i3izhHwtOb2/wUJU6G7q\nJpk8R4s3wWR6gvTu6Q74skhYcXZEcXakNJ1y+dFWBzOzPj98w7Qg7i8dzXd3SshXm/zmClGhLdJO\nW6S9NJNmgpXsKgFbXiZQmk7ZqejuhBvXbLa3faZnzZTK1JzPnQ887nxg+tkMDW7TlfTo67VojknQ\nnzT5jRXiY5iZNFfZzm8zsfGQ5ewizimaKvk0mpoUF84pLpwz/WzmF32mS4ujHowXeDAOYPaA7eu1\nGOhVdHUpHFuC/rhJsAvxBE2BJi53PMNOYYeHG+MsZhZkiOZj2Lair0fR1wMvfs4GFeLuvQwzpcVR\n5aN5xy6NzfeamTZx2R7wWMhvqBBPIeyEudh+idHEWcbTYyxuLzR8u4LPoiVhc0nbXNLsHc3PzO6v\nhJ1OmVYH8Wbo7zW7RknjsqMjwS7EIQSdIBfaLnK25RzjG2PMb8+dusVOh1V5NP/CdZutbXMCdmbW\nY3be567rcdc1M226u0pH830WibjMm/+0JNiF+BQcy+F8qzZH8BtjzG/NopQE/NOINSkunFdcOG/G\n5sszbWZmzUnY1JzPzbc8Yk3maL6v17QtDsq8+acmwS7EZ2BbNudazzOSGGUi/ZDZrRQK6bfytGzb\nDMH0dMHnnrfZzph58+WQv/ehx70PQSno7jTz5mUV7MEk2IU4ArZlM9pyljOJESY2TMCDhM9hNUX3\nN/v2PLOJSHm+/NyC+fPG2x7RCHtDNtLT5qMk2IU4QpayGGkZZThxxgT89iy+72Gpxtu+77hZlqKr\nU9HVCdefg2zWDNPMpDxm5nzuj/ncHzM9bTo7FH29Zny+vU36zUuwC3EMKgN+Mj1BaislAf8ZRSKP\nroJdWd0/ml9c9llY8nn7PQiFoLfbDNn09SiaTuHuURLsQhwjS1mcSYwwFB+WgD9ClqVIdiiSHXDt\nCuzkfGbnfFKz5mi+st98eYFUf4/ZiOQ0TKl8YrBrrS3gm8AVIAd8zXXdsYrb/zzwvwI+8Duu6/7a\nMdYqRN0oYnpxAAAOdUlEQVT6aMDP4Pu+BPwRCYcUI8OKkWHToXJ9A1JzHqnS2Hx5gZRtmVbEfaUF\nUomEX+3Sj8VBR+xfBYKu676ktX4B+EbpOrTWNvAPgevANnBXa/3bruuuHmfBQtQzCfjjV9lv/pmL\nUCj6LC76zJSO6GfnzW5Sb7zt0RTdpKcb+nrMsE2jbBN4ULC/DHwPwHXdm1rrG+UbXNctaq0vuK7r\naa27ABvYPb5ShWgcMkRzchzbzIPv7QGet8mUTsKmZj3m5uHBuM+DcTNs095WCvleRVdSYddpX5uD\ngj0OpCsuF7XWluu6HkAp1H8C+CfAvwcyx1OmEI2pMuBlmuTJiEYU50YU50YsEokI4w8ze8M284s+\nK6tm02/Hhu5uRX+Poq+3vlbCHhTsaaC54vJeqJe5rvuvtdbfBb4N/FTp70+UTDY/6eaaIXUeLanz\nYF2dz+H5VxhfHWcqPfXEhU4tLdETru7w6qFGgNGRJkZHzNf5vM9MqsDElPkzU9ouEDyamxXDAw5D\ngwGGBhzC4doN+YOC/VXgK8B3tNYvArfKN2it48C/A77suu6u1nobKB50h0tLm5+h3JORTDZLnUdI\n6jycBF1cinYwsTHO7PbsRwK+pSXK+nptfziuhxrh4+s04/Nw7YrN1rZFatYnNecxO+dz+26e23fz\nKAUd7eUeOGaGTi0N2xwU7N8Fvqy1frV0+RWt9U8CMdd1v6W1/m3gj7TWeeA94LePsVYhTg3bshlt\nPcdwYoSJ9DizWx8NeHH8Yk0KfU6hzz06dz4157O4ZFbGvnsbHMc0MOvrUfR2V7/lwROD3XVdH/i5\nx66+X3H7t4BvHUNdQgjKrQrOMRwf4eHGWGkla2NO0at1j8+d3901Uyln58wR/UzK3xu2iUTMIqm+\nHtPyoCl6siEvC5SEqAO2ZXO29TxnEqOs2/OsrT+QGTRVFgwqhgYUQwMAZqvA1LzPbGnYZuyhz9jD\n/UVSvd1mSmV3lyJ4zL1tJNiFqCO2ZXOh4wKJYlepXfCcbPhRI5qaFOdHFedHzSKptXWYnfNIzfvM\nL+z3nVcKkh2lYZseRWfH0fe2kWAXog6V+8GPJEYZW3/A/PacbNlXQ5RStLVCW6vNM5fY6zs/W+o3\nv7RsxujfuWU2/+7uMiHf12PRkvjs4/PymyBEHXMsB912gTPxEcbTYyxsz0vA16DKvvPXn4PcrtkL\ndrY0dFO5XWA0QmlBlUVv96cbn5ffACEaQHnLvpH4KA82HrAkm27XtFBQMTyoGB4EMNsFlk/Czs75\nj6yGbUkc/ufLMy9EAwk6QS61XyKXGGV8Y4zFzDyOFah2WeIAsab9DUbK4/PlkJ9b8KHpcD9Pgl2I\nBhRyQlxsv8SZxAjjG2OlI3gJ+HpQOT7/bGl8/t/++8P9DJkvJUQDCzthLrVf5vM9XyARbKHg5atd\nkjikT7OiVY7YhTgFIk6EZ5LPkslnGF9/wHJ2hYAtL/9GJc+sEKdINBDlmeQVMvkMY2sfsrKzKgHf\ngOQZFeIUigaiPNt5la3dLR5ujLOys0JAZtE0DHkmhTjFYsEYzyavsJnbZCI9bo7gJeDrnpw8FULQ\nHGrm2eRVrnfeIBaIkfcK1S5JfAYS7EKIPc2hZq52XuP5zusS8HVMgl0I8RHxULwU8M/TFGiiIAFf\nVyTYhRCfKB5K8Fzn81yTgK8rEuxCiAOVA/65zmsm4Iuy0KmWSbALIZ5aItTCc53Pc7XzGhEnIkfw\nNUqCXQhxaC3hVp7vusHV5HNEnagEfI2RYBdCfGot4VaudV3navI5OYKvIbISQQjxmbWEW3k+fIP1\nnTXGN8YoFCXgq0mO2IUQR6Y8RHO997ocwVeRBLsQ4si1Rh4bg5dZNCdKgl0IcWz2xuA7rxENSMCf\nFBljF0Icu5ZwK9fC11nfWWMi/ZCN3IbsyXqM5JEVQpyYlnArz4VbSec2eLgxznpuXQL+GMhQjBDi\nxMVDCa52XuNa5/PSbOwYSLALIaqmHPDXO2/QHGyWgD8i8hlICFF1zaFmriSf29vwY3VnDceyq11W\n3ZJgF0LUjPKGH3tb9smm25+KPGJCiJpT3rIvk8/wcGOc5ewijhWodll1Q4JdCFGzooEolzueIVvI\n8nBjnKXMggT8U5CTp0KImhdxIlxqv8yLPS/TFmmXVgUHkGAXQtSNkBPiYtslvtDzMu2RDopesdol\n1SQJdiFE3Qk6QS60XeSlvi+SjHbi+R6+71e7rJohwS6EqFuO5XC+VfNS7xfpifVIwJdIsAsh6p5t\n2Yy2nOPlvi/R19wPgOd7Va6qeiTYhRANw1IWZxIjvNT7RYbiwyhlUfRP3zi8THcUQjQcpRSD8SEG\nmgdJbc0wsznDbjGHfUpWs0qwCyEallKK/uYB+psHmN+eYzo9Tbawjd3gHSUb+38nhBAl3U09dDf1\nsJRZYmpzgq3drWqXdGwk2IUQp0oymiQZTbKWXSVtL7HkpQk02BH8E/83WmsL+CZwBcgBX3Ndd6zi\n9p8Efh4oALeBv+a6rsw1EkLUvNZIG+eTQ7T5s0ykx1nJrjZMw7GDZsV8FQi6rvsS8IvAN8o3aK0j\nwP8B/Kjrul8EEsCfOa5ChRDiOJQ7Sn6+5wVaw20N0a7goGB/GfgegOu6N4EbFbftAF9wXXendNkB\nskdeoRBCnIBoIMql9st8oedlOiJJin6xbhc7HRTscSBdcblYGp7BdV3fdd0lAK313wCaXNf9T8dT\nphBCnIygE0S3XeDl3i/RE+vFp/4WOx00oJQGmisuW67r7v0PSyH/fwJngT//NHeYTDYf/E01QOo8\nWlLn0aqHOuuhRnhynd1d1/D8q0xtTDG1MUW+mK+LufAHBfurwFeA72itXwRuPXb7r2OGZP7c0540\nXVraPHSRJy2ZbJY6j5DUebTqoc56qBGevs4m2rkQadubC79TyNT0XPiDKvsu8GWt9auly6+UZsLE\ngDeB/xn4I+A/a60BftV13X9zXMUKIUS1KKXoifXSE+tlKbPE9OYkm7ubODUY8E+sqHQU/nOPXX2/\n4uva/0wihBBHrDwXPp3bYDI9UXNTJWunEiGEqDPxUIJnk1fJ5DNMpidYzCzUxBG8dHcUQojPKBqI\ncrH9Ei/3fYnOaBe+71d1Jo0EuxBCHBHHcjjXep6X+kzbYAsLzzv5gK/+ZwYhhGgwlrIYjA8xGB9i\nfnuOmc1ptvPbJzZMI8EuhBDHqNxVci27ytTmJKs7awTtwLHepwS7EEKcgNZIG62RNjL5DBPphyxl\nFo/tCF6CXQghTlC5J02hVTOxMc58Zh7f97HU0Z3ylJOnQghRBY7lcLb1PC/3fonh+DCOChxZZ0k5\nYhdCiCpSSjEQH2IgPrS3ojW9u/mZNv+QYBdCiBpRXtG6mdtkanOS5ezSpxqHl2AXQoga0xxq5nLo\nGfLFPBPphwDbh/n3MsYuhBA1KmAHONd6Hv/v+muH+XcS7EII0WAk2IUQosFIsAshRIORYBdCiAYj\nwS6EEA1Ggl0IIRqMBLsQQjQYCXYhhGgwEuxCCNFglO/71a5BCCHEEZIjdiGEaDAS7EII0WAk2IUQ\nosFIsAshRIORYBdCiAYjwS6EEA3mRHZQ0lpbwDeBK0AO+JrrumMncd9PQ2v9AvCPXNf9Ma31WeDb\ngAe8D/x113WrPidUax0A/l9gCAgBfx/4gBqrVWttA98CzgM+8Fcxz/m3qaE6AbTWncBbwB/H1PZt\naq/Gt4GN0sVx4B9Sm3X+EvAVIAD8E+BVaqxOrfVPAz9TuhgBrgJfBH6V2qrTAn4D8xrygP8FKHKI\nx/Okjti/CgRd130J+EXgGyd0vwfSWv8CJohCpav+L+Drruv+CKCAP1ut2h7zl4ClUl1/Cvi/MY9j\nrdX6ZwDPdd0vAn8H+AfUYJ2lN8pfx2w5pqjB511rHQZwXffHSn/+CrVZ548CXyi9vn8UGKEGn3PX\ndX+z/FgCbwJ/A/hlaqxO4MeBptJr6H/nU7yGTirYXwa+B+C67k3gxgnd79N4APwE5sECeN513T8q\nff0fgD9Rlao+6juYX0Iwz1ueGqzVdd1/C/xs6eIwsAZcr7U6gV8B/h9grnS55h5LzBFlVGv9H7XW\nf6C1fpHarPPHgdta638D/Dvg/6M2n3MAtNY3gEuu6/4GtVlnFkhorRWQAHY5ZJ0nFexxIF1xuVj6\nuFF1ruv+a6BQcZWq+HoL88BWneu6267rbmmtmzEh/3d49PmrpVqLWutvYz7i/g419phqrX8G8+nn\n90tXKWqsxpJt4Fdc1/2TmCGt33ns9lqpMwlcB/4Cps5/SW0+nmVfB/5e6etarPNVIAzcw3yq/DUO\nWedJhWsaaK68X9d1vRO678OqrKsZWK9WIY/TWg8A/xn4Ldd1f5cartV13Z8BNGasMFxxUy3U+Qrw\nZa31HwLPAb+JCaeyWqgR4D6lMHdd90NgBeiquL1W6lwGft913YLruveBHR4NnlqpE611C3Dedd3/\nWrqqFl9DvwC86rquxvx+/hbm3EXZgXWeVLC/CvxPAKWPk7dO6H4/jXe01n+s9PX/CPzRk775pGit\nu4DfB37Bdd1vl66uuVq11n+5dCINzEfKIvBmLdXpuu4fc133R0tjre8CPwV8r5ZqLHmF0vkorXUv\n5gX9+zVY5/cx533KdUaBP6jBOgF+BPiDiss19xoCmtgf4VjDTHI5VJ0nMisG+C7mCOnV0uVXTuh+\nD6N8hvlvA9/SWgeBu8C/ql5Jj/g65ijol7XW5bH2nwd+rcZq/VfAt7XW/xVzlPHzmI+UtfiYlvnU\n5vP+z4B/rrUuv4hfwRy111Sdruv+ntb6R7TWr2MOFv8aMEGN1VlyHqickVeLz/uvYJ73/4Z5Df0S\nZvbWU9cp3R2FEKLB1MQJTCGEEEdHgl0IIRqMBLsQQjQYCXYhhGgwEuxCCNFgJNiFEKLBSLALIUSD\nkWAXQogG898B5UTS4hVfe6cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108b62470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xx[:, 1], πs)\n",
    "plt.fill_between(xx[:, 1], (πs - 1.96*vv).ravel(), (πs + 1.96*vv).ravel(),\n",
    "                 alpha=.25, color='g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}