Utility Theory and Logistic Regression Blog Post

post.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:632fcaeac19ef24277da49c7dd3db4e32af99d70db1d530b08809f6b9fd7dfb7"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Logistic regression is perhaps one of the most commonly used tools in all of statistics.  While I have mathematically understood and used logistic regression for quite some time, it took me much longer to develop intution for it.  It was only a few months ago that I learned about the connection between logistic regression and utility theory.\n",
      "\n",
      "Suppose a person must decide between two options, $Y = 0$ or $Y = 1$.  Utility theory posits that each option has an associated utility, $U_0$ and $U_1$, respectively.  The person chooses the option with the largest utility.  My intuition for logistic regression arises from understanding that it arises from a specific utility model."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In logistic regression, we are interested in modeling the effect of a covariate, $X$, on the person's choice.  In this post, we will work with the logistic regression model\n",
      "\n",
      "$$\\begin{align*}\n",
      "P(Y = 1 | X = x)\n",
      "    & = \\frac{e^x}{1 + e^x}.\n",
      "\\end{align*}$$\n",
      "\n",
      "This model is shown graphically below."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "\n",
      "from matplotlib import pyplot as plt\n",
      "import numpy as np\n",
      "from scipy.stats import gumbel_r"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def logistic_probability(x):\n",
      "    exp = np.exp(x)\n",
      "    \n",
      "    return exp / (1 + exp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots(figsize=(8, 6))\n",
      "\n",
      "n = 100\n",
      "xs = np.linspace(-5, 5, n)\n",
      "\n",
      "ax.plot(xs, logistic_probability(xs));\n",
      "\n",
      "ax.set_xlim(-5, 5);\n",
      "ax.set_xlabel('$x$');\n",
      "\n",
      "ax.set_ylabel('$P(Y = 1 | X = x)$');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAF/CAYAAAC/oTuRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmclVXhx/HP+SG5hAuGK+KWioFrimKWjqm5pJmVuael\n5a8wrTQtzaTS0LLccl/7aQhppZImgjHihoKJuACCiAoigqyihjDn98cZFXEYhpk799z73M/79Xpe\nd+7M5fr1vpTvnPOc5zwhxogkSapu/5M7gCRJajsLXZKkArDQJUkqAAtdkqQCsNAlSSoAC12SpALI\nUughhBtDCNNDCM8085rLQggTQghPhxB2LGc+SZKqTa4R+k3A/sv6YQjhQGCLGOOWwPeAq8oVTJKk\napSl0GOMDwGzm3nJV4A/N772cWCtEMJ65cgmSVI1qtRz6F2BV5d4PgXYKFMWSZIqXqUWOkBY6rl7\n1EqStAwr5Q6wDFOBbks836jxex8RQrDkJUk1Jca49IAXqNwR+t3AtwBCCL2BOTHG6U29MMZYVce5\n556bPUPRDz9jP+MiHH7Gfs6LFkVmz4688krk2Wcjjz3W/Bg2ywg9hHAbsCfQJYTwKnAu0BEgxnhN\njPHeEMKBIYSJwALg2zlySpLUFu+9B2++CTNnwqxZHz1mz07HnDkfPs6dm4558+Dtt6FTJ1h99Q+P\n5mQp9BjjkS14zcnlyCJJUkvFmMp22jR4/fUPjzfe+Ogxc2Y6FiyAtdeGT30qHWuvnY7OndPjBhuk\nr9daKx1rrvnh0akT/M9S8+ihycn2pFLPoRdWXV1d7giF52fc/vyM25+fcXks+TnHmEr4lVfS8eqr\n6Zg69cPjtddgpZVg/fXTscEGsN566dh1V1h3XVhnnfTYpQusscbHS7m9hBird11ZCCFWc35JUnkt\nWAAvvggTJ6bHyZPhpZfS4+TJsOqqsPHG0K1betxoI+ja9aNHp0758ocQiMtYFGehS5IKZfHiVM7j\nxsH48elx3LhU4rNnw+abwxZbpMfNN4dNN4XNNoNNNln+eercLHRJUuHEmM5fjx4Nzz774TFuXJru\n3npr6N79w8ettkoj7HJNgbcHC12SVNViTKPukSPhySdTiY8enUbj228P220H22yTjh49Kn+k3VoW\nuiSpqsyZAyNGwKOPwhNPwKhRsPLK0KsX7LQT7Lgj7LBDGnE3t/K7aCx0SVJFmzIF6uth+PBU4i+/\nDDvvDLvtllaP9+oFG26YO2V+FrokqaJMnw5Dh8KwYanI586FPfeEPfaA3XdPU+gdO+ZOWXksdElS\nVv/9Lzz8MNx/PwwenEbge+0Fe+8NdXXwmc9U92K1crHQJUllN2MG3HMPDBoEDzyQSnu//eBLX4Jd\ndkkbtGjFWOiSpLJ45RX429/gjjvguedgn33g4IPhwAPTDmpqGwtdktRuXn0VBg5MJT5xIhxyCHzj\nG/DFL6aV6SodC12SVFJz56YCv/VWGDMGvvY1OOywdF7cxWztx0KXJLVZQ0M6F3799XDffWlB2zHH\nwJe/7Ei8XCx0SVKrTZsGN92UinzNNeG734Ujjki3/1R5NVforjGUJH1MjOkys8suS6Pyww6Dv/41\n7dJWSzuzVRMLXZL0gXffhQEDUpG/9RaccgrceGNx90YvEqfcJUnMnQtXXgmXXpr2SD/lFNh/fzd7\nqTROuUuSmjR9Olx8MVx3XbpWfMgQ2Hbb3KnUGv7uJUk16PXX4dRT0+5t8+enu5ndcotlXs0sdEmq\nITNnwplnQs+eaXHb88/DFVfAZpvlTqa2stAlqQbMnw/nngvdu8O8efD003DJJbD++rmTqVQsdEkq\nsMWL0/Xj3bvDiy+mqfWrroKNNsqdTKXmojhJKqghQ+C002CtteCuu6BXr9yJ1J4sdEkqmMmT04K3\n556D3/0ODj3UzWBqgVPuklQQCxdCv36w887pfuPPPZdummKZ1wZH6JJUAMOGwQ9+AJ/+NIwc6ar1\nWmShS1IVmzcPTj893f3s8svhK19xRF6rnHKXpCo1ePCHG8E88wwccohlXsscoUtSlZk7F37ykw/v\nTb7vvrkTqRI4QpekKvLww7D99tCxYxqVW+Z6nyN0SaoCixbBeefB1VenG6kcfHDuRKo0FrokVbjJ\nk+Hoo+GTn4SnnoINNsidSJXIKXdJqmB33ZWuKf/a19JKdstcy+IIXZIq0KJFcM450L8/DBoEu+6a\nO5EqnYUuSRVm+nQ48kjo0CHdTGWddXInUjVwyl2SKsiIEWnr1t13T1PslrlayhG6JFWI/v3hRz+C\nG25wFbtWnIUuSZk1NMC558Ktt8K//w3bbJM7kaqRhS5JGb39Nhx3HEybBo8/DuuumzuRqpXn0CUp\nk+nTYY89YLXV0jaulrnawkKXpAwmToTPfS7dHe3mm2HllXMnUrWz0CWpzEaNSiPzM8+EX/7SO6Sp\nNDyHLklldP/9cMwxcO218NWv5k6jIrHQJalMBgyAU0+Fv/8dPv/53GlUNBa6JJXBzTfDWWfB0KGw\n7ba506iILHRJamfXXJNuffrvf8PWW+dOo6Ky0CWpHV16KVx8MdTXw6c/nTuNisxCl6R28rvfpcVv\nDz4Im2ySO42KzkKXpHbwhz/A9denMu/aNXca1QILXZJK7Ior0mGZq5wsdEkqoRtvhAsvTGXerVvu\nNKolFroklUj//nDOOTBsGGy2We40qjUWuiSVwD/+Aaedlq4z32qr3GlUi0KMMXeGVgshxGrOL6kY\nHnwQDjsM7rsPPvvZ3GlUZCEEYoxN7v7vzVkkqQ2eeSaVef/+lrnystAlqZVefhkOOAAuuwz22Sd3\nGtU6C12SWmHmTNhvPzj9dDjiiNxppEyFHkLYP4QwLoQwIYRwZhM/7xJCuC+EMDqE8GwI4fgMMSWp\nSe+8AwcfDIccAj/6Ue40UlL2RXEhhA7AeGAfYCowEjgyxjh2idf0BVaOMf48hNCl8fXrxRgXLfVe\nLoqTVFYNDWlEvtJK8Je/QGhyeZLUPiptUdwuwMQY4+QY43vAAOCQpV4zDVij8es1gDeXLnNJyuGX\nv4SpU9MGMpa5KkmO69C7Aq8u8XwKsOtSr7kO+HcI4TVgdeCbZcomScv0f/+XVrOPGAGrrJI7jfRR\nOUboLZkjPwsYHWPcENgBuCKEsHr7xpKkZXvoobQA7p//hHXXzZ1G+rgcI/SpwJI7HHcjjdKX9Dng\nfIAY44shhJeA7sCopd+sb9++H3xdV1dHXV1dadNKqnkvvpiuNb/1VujRI3ca1ZL6+nrq6+tb9Noc\ni+JWIi1y2xt4DXiCjy+K+yMwN8b4qxDCesCTwHYxxllLvZeL4iS1q/nzYbfd4Pvfhz59cqdRrWtu\nUVyWrV9DCAcAlwAdgBtijP1CCCcBxBivaVzZfhOwMem0QL8YY/8m3sdCl9RuYkwj87XWguuucxGc\n8qu4Qi8VC11Se+rXD+68M+3V7iI4VYLmCt27rUlSE+67Dy6/HJ54wjJXdbDQJWkpEyfCccfBHXfA\nRhvlTiO1jHu5S9ISFiyAQw9NG8h84Qu500gt5zl0SVrCt78NixalTWRcBKdK4zl0SWqBm25Ku8CN\nHGmZq/o4Qpck4NlnYa+9oL4eevbMnUZqWqXdnEWSKspbb6XrzS+6yDJX9XKELqmmxQjHHguf+ES6\ng5pUyTyHLknLcOON8PTT8PjjuZNIbeMIXVLNGj8edt8dhg/3piuqDp5Dl6SlLFwIRx0F551nmasY\nHKFLqklnnAEvvAD/+IeXqKl6eA5dkpYwdCj07w+jR1vmKg6n3CXVlJkz4fjj0yYyXbrkTiOVjlPu\nkmpGjGmf9i22SNecS9XGKXdJAm6+GSZPhoEDcyeRSs8RuqSa8PLLsPPO8MADsN12udNIreNla5Jq\nWkMDfOc7cNpplrmKy0KXVHhXXglvvw2nn547idR+nHKXVGgvvACf+xw8+ihstVXuNFLbOOUuqSYt\nXgzHHQfnnmuZq/gsdEmF9Yc/wKqrQp8+uZNI7c8pd0mF9P6NV0aNgk03zZ1GKg2n3CXVlIYGOPHE\nNNVumatWWOiSCufKK1OpO9WuWuKUu6RCmTw5bSDz8MOw9da500il5ZS7pJoQI5x0UtpAxjJXrbHQ\nJRXGn/8MM2a4gYxqk1Pukgrh9dfTtq6DB8OOO+ZOI7WP5qbcLXRJhXDUUbDxxnDBBbmTSO3H26dK\nKrTBg2HECLj++txJpHw8hy6pqr39NvzgB3DFFbDaarnTSPlY6JKq2nnnQa9ecMABuZNIeXkOXVLV\nevZZ2GsvGDMGNtggdxqp/XkduqTCaWhI15z/5jeWuQQWuqQqdf31qdS/973cSaTK4JS7pKozcyb0\n6AFDhsD22+dOI5WP16FLKpQTT4ROneCSS3InkcrL69AlFcaIEXDvvTB2bO4kUmXxHLqkqrF4cbrm\n/Pe/hzXXzJ1GqiwWuqSqcfXVsPrqaZtXSR/lOXRJVeGNN6BnTxg2DLbZJncaKQ8XxUmqescfD126\nwEUX5U4i5eOiOElV7dFH0yVq48blTiJVLs+hS6poDQ1wyilw4YXp/LmkplnokiraTTfBJz4BRx+d\nO4lU2TyHLqlizZkDn/kMDBoEO++cO42Un4viJFWln/wE5s+H667LnUSqDBa6pKozdizssQc89xys\nu27uNFJl8PapkqpKjPCjH8HZZ1vmUktZ6JIqzqBB8Oqr0KdP7iRS9fA6dEkVZeFCOO00+NOfoGPH\n3Gmk6uEIXVJF+dOfoHt32G+/3Emk6uKiOEkVY8YM6NEDHnoItt46dxqp8rjKXVJV6NMHOnSAyy7L\nnUSqTO7lLqniPfcc3H57ulxN0orzHLqk7GJMm8icfTZ86lO500jVyUKXlN2//gUvvww/+EHuJFL1\nylLoIYT9QwjjQggTQghnLuM1dSGEp0IIz4YQ6sscUVKZvPdeukztoou8TE1qi7KfQw8hdAD+BOwD\nTAVGhhDujjGOXeI1awFXAPvFGKeEELqUO6ek8rjuOujaFb785dxJpOqWY1HcLsDEGONkgBDCAOAQ\nYMmlMEcBf4sxTgGIMc4sd0hJ7W/uXPj1r+G++yA0uW5XUkvlmHLvCry6xPMpjd9b0pbA2iGEYSGE\nUSGEY8uWTlLZXHABHHAA7LBD7iRS9csxQm/JheMdgc8CewOrAY+FEEbEGCe0azJJZfPKK3DttTBm\nTO4kUjHkKPSpQLclnncjjdKX9CowM8b4DvBOCGE4sD3wsULv27fvB1/X1dVRV1dX4riS2sPZZ6eN\nZLouPT8n6QP19fXU19e36LVl3ykuhLASMJ40+n4NeAI4cqlFcVuTFs7tB6wMPA4cHmN8fqn3cqc4\nqQo9+SQcdBC88AKsvnruNFL1qKid4mKMi0IIJwODgQ7ADTHGsSGEkxp/fk2McVwI4T5gDNAAXLd0\nmUuqTjGmy9R+9SvLXCol93KXVFZ33w0//zk8/TSs5ObT0gopyQg9hPBJ4GhgG9LIehXS6PktYARw\ne4yxoe1xJRXVokVw5plpExnLXCqtFv0vFULYF+gB/DPGeO1SPwukBWs/CSEMjTGOLn1MSUVw442w\n/vpw4IG5k0jFs9wp9xDCKsBGMcaJy32zELaNMT5TqnAt+Oc55S5ViQULYMst4a67oFev3Gmk6lTS\n+6GHEDYHpjVeUpaVhS5Vj9/8Bp5/Hm67LXcSqXqVutCvIJ0vrw8hfB6IMcZHSpBzhVnoUnWYPh16\n9ICRI2HzzXOnkapXc4Xemq1fnwA2CyFsFmN8GFi3TekkFd6vfw3HHmuZS+2pNetMuwGTSIvgtgEe\nAf5R0lSSCuOFF2DgQBg3LncSqdhaU+iTSHdC6994W9OvlTiTpAI56yw4/XTo4k2QpXbVmin3gUDP\nxq83A9YrXRxJRTJiBDz+OJx6au4kUvG1aqe4EMJaMcY5IYTOMcbZ7ZCrpTlcFCdVqBhhzz3h+OPh\nO9/JnUYqhlIvigM4rvHxW63885IK7t574c034Vv+LSGVRWsLXZKWafFi+NnPoF8/t3iVysVCl1Ry\nt94Ka64JBx+cO4lUO/zdWVJJvfsu/PKX0L8/hCbP9ElqD47QJZXUlVfCDjvA7rvnTiLVFkfokkpm\nzhy44AIYNix3Eqn2tHaEPmSpR0nid7+Dgw6Cnj2X/1pJpdWq69ArhdehS5Xjtddg221h9Gjo1i13\nGqmY2nwdeghhp9JGklQ0v/pV2kDGMpfyaNEIPYTw9xhjk3u2hxC6xhinljxZCzhClyrD+PFpEdwL\nL8Daa+dOIxVXKXaKmx9C+G4I4SOL6EIIawK/bWtASdXtnHPgtNMscymnFp9DDyGsRtry9V5gZ+Bo\nYCfg3Rhj93ZL2HwmR+hSZqNGwSGHwIQJsNpqudNIxdbcCL1Fl62FEI4BpgK9gH7As8B5wFBgqxLl\nlFSFfvaztJGMZS7l1dJz6AtJl6jdCtwNbA1sHGP8R/vGW24uR+hSRkOGQJ8+8Nxz0LFj7jRS8TU3\nQm9poZ8SY7xsqe+tCxzS+B7XliTpCrLQpXwaGqBXrzRCP+yw3Gmk2tDmQm/mjVcBHoox9mr1m7SB\nhS7lM3AgXHQRPPGEe7ZL5dLmc+jLEmN8N4RwXlveQ1L1ee89+MUv4OqrLXOpUrT55iwxxrtKEURS\n9bjhBthsM9h779xJJL1vRS5bOwpobtnLezHG/iVJ1UJOuUvlt2ABbLklDBoEO7mHpFRW7XYOPTcL\nXSq/fv3Sfu0DB+ZOItWedi30EMInY4wL2vQmrf9nW+hSGc2aBd27wyOPwFbuQCGVXSm2fm3Od0vw\nHpKqwAUXwNe/bplLlail16H/EdgTmNfEjz8TY1y/1MFawhG6VD5TpsD228Mzz8CGG+ZOI9WmUmws\n8z/Aj2KMf2ziZz+OMV7c9pgrzkKXyufEE6FLlzRKl5RHSc6hhxA6xxhnN/F9z6FLBTd2LOyxR7o9\naufOudNItctV7pLa5Gtfg9694YwzcieRalvJCz2EsFaMcc6yRu3lYqFL7W/EiLRX+wsvwKqr5k4j\n1bb2WOV+XOPjt1r55yVVgRjTzVfOPdcylypdKS5bk1RQgwfD9Olw/PG5k0haHgtdUpMaGtLo/Pzz\nYaU23cZJUjlY6JKaNGAArLwyHHpo7iSSWsLfuyV9zMKFcM45cP313h5VqhaO0CV9zLXXpjuq7bVX\n7iSSWqqlO8XtEWMcvsTzHjHG599/bNeEzefysjWpxObPT2U+eHDa6lVS5SjF1q+PAnUxxoWlDtcW\nFrpUen37wosvwi235E4iaWmlKPSvAWsCD8YYJ5U4X6tZ6FJpTZ8OPXrAk0/CppvmTiNpaSXbKS6E\nsDfQMcZ4X6nCtYWFLpXWySdDx45wcZbbLUlanlKM0DvEGBc3fr0n8E3gNuA/wF4xxntKmLfFLHSp\ndCZOTPu1jxuX7qomqfKUotCvBaYBXwVWBf4JfBL4LLBxjHG90sVtOQtdKp0jjoBtt4Wzz86dRNKy\nNFfoLb0O/YvA9cARMcaxS735qW3MJymzUaNg+HC44YbcSSS1VktH6PvEGIcu42crxxj/W/JkLeAI\nXWq7GGGffdId1f73f3OnkdScNt1tLYSwMjB6WT9fssxDCBu3KqGkbAYPhqlT4YQTcieR1BbLLfTG\nwu4dQjgqhNDkDRRDCJ1DCN8DNil1QEntZ/FiOOMMuOCCtLpdUvVq6ZT76sCPG592AxYDHRsf3wam\nANfFGOe2U85l5XLKXWqDm29O+7U/9JB7tkvVoBSr3K8G5pLKvCtwYIxxQUlTtoKFLrXeO+9A9+4w\ncCDstlvuNJJaohSr3J+JMV7R+GYbAIcDN5Yon6QMLrsMevWyzKWiaGmhf7DwLcY4LYQwr53ySCqD\nmTPh97+HRx/NnURSqbS00H8WQtiBtDPcU8AH89whhPVijNPbI5yk9nH++XD44bDVVrmTSCqVlp5D\nPwcYCfQGegE7Aq8AjwDrxBi/tUL/0BD2By4BOgDXxxgvXMbregGPAd+MMf69iZ97Dl1aQZMmpan2\n55+H9bLs8SiptUp2c5al3vTTwK7Ad2OMe63An+sAjAf2AaaSflE4sokd6DoAQ0ir6G+KMf6tifey\n0KUVdPjhaYvXX/widxJJK6oUi+I+Jsb4IvBiCGHKCv7RXYCJMcbJjeEGAIcAY5d63Q+BO0gzApJK\n4LHH4JFH4KabcieRVGrL3VhmeWKMw1fwj3QFXl3i+ZTG730ghNCVVPJXvf+PaXVASUDa4vX00+G8\n82C11XKnkVRqbS70VmhJOV8C/KxxPj00HpLa4O9/hwUL4NhjcyeR1B5aPeXeBlNJG9S8rxtplL6k\nnYABIW1d1QU4IITwXozx7qXfrG/fvh98XVdXR11dXYnjStVv4UI480y4+mro0CF3GkktVV9fT319\nfYte2+pFca0VQliJtChub+A14AmaWBS3xOtvAga5yl1qvUsugSFD4J57cieR1BbtsiiutWKMi0II\nJwODSZet3RBjHBtCOKnx59eUO5NUZLNnw29/C8OG5U4iqT2VfYReSo7QpeU77TR46y24xl+VparX\nLtehVwILXWrehAnQuzc89xysv37uNJLaqrlCz7HKXVKZnHEG/PSnlrlUC3KscpdUBsOGwejRcNtt\nuZNIKgdH6FIBLV4MP/4xXHghrLJK7jSSysFClwro5puhUyc47LDcSSSVi4vipIKZPx+6d4e77kp3\nVZNUHC6Kk2pIv36w776WuVRrHKFLBTJ5Muy0E4wZA127LvflkqqMI3SpRpx+eloMZ5lLtcfL1qSC\neOABePJJuOWW3Ekk5eAIXSqARYvg1FPhj3+EVVfNnUZSDha6VABXXZV2g/vqV3MnkZSLi+KkKjdj\nBvToAfX10LNn7jSS2pM3Z5EK7H//F1ZeGS69NHcSSe2tou6HLql0nnoK7rwTxo7NnURSbp5Dl6pU\nQwP88Ifw619D586500jKzUKXqtQtt8DChXDCCbmTSKoEnkOXqtCcOfCZz8CgQbDzzrnTSCoXF8VJ\nBXPyyekWqVddlTuJpHJyUZxUIP/5D9x+uwvhJH2U59ClKtLQAH36wG9/C2uvnTuNpEpioUtV5Kab\n0uO3v503h6TK4zl0qUrMmpUWwv3rX/DZz+ZOIykHF8VJBXDiienGK5dfnjuJpFxcFCdVueHD4b77\n4PnncyeRVKk8hy5VuP/+F046CS67DNZYI3caSZXKQpcq3O9/D1tuCYcemjuJpErmOXSpgk2YALvt\nlq4933jj3Gkk5dbcOXRH6FKFihG+/3046yzLXNLyWehShfrLX+DNN+GUU3InkVQNnHKXKtAbb8B2\n26Wbr/TqlTuNpErhdehSlTniCOjWLS2Ik6T3eR26VEXuuguefBJuvDF3EknVxBG6VEHmzIFttknn\nz/fcM3caSZXGKXepSpx4InTs6H3OJTXNKXepCgwdCvffD88+mzuJpGrkZWtSBXjrLfje9+Caa9ze\nVVLrOOUuVYA+fVKp//nPuZNIqmROuUsVbMiQdL35mDG5k0iqZk65SxnNmQMnnAA33ABrrZU7jaRq\n5pS7lNFxx0GnTnDFFbmTSKoGTrlLFejOO+GRR+Dpp3MnkVQEjtClDGbMSHu133EH7L577jSSqoUb\ny0gVJEb4+tdhiy3gd7/LnUZSNXHKXaogN9wAkybBbbflTiKpSByhS2U0fnyaYh8+HHr0yJ1GUrVp\nboTuZWtSmSxcCEcdBb/5jWUuqfQcoUtlcsYZaYR+550Qmvz9WpKa5zl0KbOhQ9MtUUePtswltQ+n\n3KV2NnMmHH883HwzrLNO7jSSisopd6kdNTTAl78M227rJWqS2s5FcVImF14I8+bB+efnTiKp6DyH\nLrWT4cPh0kth1Cjo2DF3GklF5whdagdvvJEuUbvpJthoo9xpJNUCz6FLJbZ4MRxwAPTq5VS7pNLy\nHLpURuefnzaR+dWvcieRVEs8hy6V0D33wDXXwMiRsJL/d0kqI//KkUpkwgT49rfhH/+ADTfMnUZS\nrck25R5C2D+EMC6EMCGEcGYTPz86hPB0CGFMCOGREMJ2OXJKLTF/Pnz1q/DrX3t/c0l5ZFkUF0Lo\nAIwH9gGmAiOBI2OMY5d4zW7A8zHGuSGE/YG+McbeS72Pi+KUXYxw2GHQuTNce61bu0pqP5W4l/su\nwMQY42SAEMIA4BDgg0KPMT62xOsfB7z4RxXpggtgypS0V7tlLimXXIXeFXh1iedTgF2bef0JwL3t\nmkhqhXvugcsvT4vgVl45dxpJtSxXobd4njyEsBfwHcAzk6ooY8akm67cfTd07Zo7jaRal6vQpwLd\nlnjejTRK/4jGhXDXAfvHGGc39UZ9+/b94Ou6ujrq6upKmVNq0rRpcPDBaXS+226500gqqvr6eurr\n61v02lyL4lYiLYrbG3gNeIKPL4rbGPg3cEyMccQy3sdFcSq7t9+Gujo46CD45S9zp5FUS5pbFJdt\n69cQwgHAJUAH4IYYY78QwkkAMcZrQgjXA4cCrzT+kfdijLss9R4WusqqoQEOPzydL7/lFhfBSSqv\niiz0UrDQVW5nnQUPPggPPACrrJI7jaRaU4mXrUlV509/gjvugEcescwlVR4LXWqB22+Hfv3goYdg\nnXVyp5Gkj7PQpeWor4c+fWDwYNh889xpJKlp3j5VasbTT8M3vwkDBsCOO+ZOI0nLZqFLyzBpEnz5\ny+nc+Re/mDuNJDXPQpea8OqrsPfecPbZaYQuSZXOQpeWMm1aGpH/8Ifw/e/nTiNJLWOhS0uYMQP2\n2Sft0f6Tn+ROI0ktZ6FLjWbNgn33hUMPTVPtklRN3ClOIpX5fvvBF74Af/iDW7pKqkzN7RTnCF01\nb8aMdM58jz0sc0nVy0JXTZs27cM7p110kWUuqXpZ6KpZU6bAnnvCkUfCeedZ5pKqm4WumjRpUppi\nP+kk+MUvcqeRpLaz0FVznnoqLX776U/htNNyp5Gk0vDmLKopDzyQptivugq+/vXcaSSpdByhq2YM\nGJDK/PbbLXNJxeMIXYUXI1x6aVrFPnQobLdd7kSSVHoWugrtvffg1FPhwQfh4Ydh001zJ5Kk9mGh\nq7Bmz4bHG0AbAAAJVUlEQVTDDoNPfAIeewzWWCN3IklqP55DVyFNmAC9e8O228KgQZa5pOKz0FU4\nQ4bA5z+fLkm7+GLo0CF3Iklqf065qzAaGuC3v02XpP31r2kXOEmqFRa6CmHWLDj2WJg3D0aOhA03\nzJ1IksrLKXdVvf/8B3beGbp3h3//2zKXVJssdFWthga45JJ0H/MLLoA//hE6dsydSpLycMpdVWn6\ndDj++HRp2uOPw+ab504kSXk5QlfVufde2GGHNM3+0EOWuSSBI3RVkXnz4Iwz4F//goED0+1PJUmJ\nI3RVhfvvT3uwL1oEY8ZY5pK0NEfoqmhz56YNYoYMgWuvTQvgJEkf5whdFSnGdJvTnj1hpZXgmWcs\nc0lqjiN0VZwJE+Dkk+G11+C22+ALX8idSJIqnyN0VYx33oFzz4XddoMvfSltGGOZS1LLWOjKrqEB\n/vKXtNPb2LEwenQ6b+4mMZLUck65K6uHHkrlDdC/f7pLmiRpxVnoyuL55+Gcc2DUKOjXD444Av7H\n+SJJajX/ClVZTZgAxxwDe+0FvXvDuHFw1FGWuSS1lX+NqiwmTYITTkgL3rbeGiZOhJ/+FFZdNXcy\nSSoGC13t6pln4OijYZddYIMN0gj9F7+A1VfPnUySisVCV8nFCA8/DAcfnC4/2247ePFFOO886Nw5\ndzpJKiYXxalk/vtf+Otf4dJLYc4cOP30tNvbKqvkTiZJxRdijLkztFoIIVZz/qJ47bW0z/o118C2\n28Ipp8CBB7rQTZJKLYRAjDE09TNH6GqVxYvhvvtSkQ8fDocfDg88AD165E4mSbXJQtcKmTABbr0V\nbrwRNtwQvvvdtMtbp065k0lSbbPQtVwzZ8LAganIJ01Km8D885+w/fa5k0mS3uc5dDVp1iy46y64\n4w545JF0TvyYY2Dffd1jXZJyae4cuoWuD7z+ehp533EHPPYY7LMPfOMbcNBBXjcuSZXAQleTYoQx\nY2DQoHS88EK6bvzrX08jcs+LS1JlsdD1gRkzYOhQuP/+dKyyStoA5uCD073HP/GJ3AklSctiodew\nOXPSLUrr62HYsLRj2157pZH4l74En/40hCb/05AkVRoLvYZMmQKPPpqOhx5K0+i9e0NdHey5J+y6\nq4vaJKlaWegF9dZb8J//wMiR8MQTaSHbO+/A5z6Xjt13TzdFcRpdkorBQi+AWbNg9Oh0PP00PPkk\nvPRS2mq1V6907LYbbLGFU+iSVFQWehV5+20YNw6effajx5w5aSOXHXZIx447wjbbOPqWpFpioVeY\nRYvg5ZfTArUJE2D8+FTi48enVehbbJFG3j17ptLu2RM228ybnUhSrbPQy6yhAd54A155BSZPTlPj\nL72Uvn7xxfT9DTZIxb3FFtC9O2y9dXrcZBPo0CH3v4EkqRJZ6CX07rswfTpMmwZTp6bjtdfS45Qp\nqaynToU11oBu3dLIetNNP3zcfPN0rLxyWWNLkgqg4go9hLA/cAnQAbg+xnhhE6+5DDgAeBs4Psb4\nVBOvaXOhL1oEb76Zjpkz0zFjRhphv//4xhtpW9TXX08ry9dbD9ZfH7p2/fixySaw0Uaw6qptiiVJ\n0sdUVKGHEDoA44F9gKnASODIGOPYJV5zIHByjPHAEMKuwKUxxt5NvFd8993I/Pkwfz7Mmwdz5354\nzJuXFpPNmQOzZ3/4OHt2WjU+axYsWACdO8M660CXLun41KdSaa+7bvr+uuum5xtskF7blnPZ9fX1\n1NXVtf4NtFx+xu3Pz7j9+RmXR7V9zs0Veo7bp+4CTIwxTgYIIQwADgHGLvGarwB/BogxPh5CWCuE\nsF6McfrSb9apU5reXmONdAORNddMxxprfPj1pz6VdkTr3BnWWgvWXjt9b+21058p52KzavuPpxr5\nGbc/P+P252dcHkX6nHMUelfg1SWeTwF2bcFrNgI+VugLF3rdtSRJOS6Eaukc/9I13eSfs8wlScpz\nDr030DfGuH/j858DDUsujAshXA3UxxgHND4fB+y59JR7CKF6l+hLktQKlXQOfRSwZQhhU+A14HDg\nyKVeczdwMjCg8ReAOU2dP1/Wv5QkSbWm7IUeY1wUQjgZGEy6bO2GGOPYEMJJjT+/JsZ4bwjhwBDC\nRGAB8O1y55QkqZpU9cYykiQpcXfwTEIIp4UQGkIIa+fOUkQhhN+HEMaGEJ4OIfw9hLBm7kxFEULY\nP4QwLoQwIYRwZu48RRNC6BZCGBZCeC6E8GwI4ZTcmYoqhNAhhPBUCGFQ7iylYKFnEELoBuwLvJw7\nS4HdD/SMMW4PvAD8PHOeQmjcGOpPwP5AD+DIEMJn8qYqnPeAH8cYewK9gT5+xu3mVOB5Wn71VUWz\n0PP4I3BG7hBFFmMcEmNsaHz6OGkfA7XdBxtDxRjfA97fGEolEmN8PcY4uvHrt0ibbm2YN1XxhBA2\nAg4Erufjl0lXJQu9zEIIhwBTYoxjcmepId8B7s0doiCa2vSpa6Yshdd4NdCOpF9KVVoXAz8FGpb3\nwmqR47K1wgshDAHWb+JHZ5Omfr+05MvLEqqAmvmcz4oxDmp8zdnAwhhj/7KGK65CTE1WgxBCJ+AO\n4NTGkbpKJIRwEPBGjPGpEEJd7jylYqG3gxjjvk19P4SwDbAZ8HRIW9xtBDwZQtglxvhGGSMWwrI+\n5/eFEI4nTantXZZAtWEq0G2J591Io3SVUAihI/A34NYY45258xTQ54CvNN4IbBVgjRDC/8UYv5U5\nV5t42VpGIYSXgJ1ijLNyZymaxlv0/oG0w+DM3HmKIoSwEuluiXuTNoZ6gqXulqi2Cem3/T8Db8YY\nf5w7T9GFEPYETo8xHpw7S1t5Dj0vf5tqP5cDnYAhjZelXJk7UBHEGBeRdnEcTFodPNAyL7ndgWOA\nvRr/232q8RdUtZ9C/F3sCF2SpAJwhC5JUgFY6JIkFYCFLklSAVjokiQVgIUuSVIBWOiSJBWAhS5J\nUgFY6JIkFYCFLklSAXhzFkktFkLoABwObE66jeouwB9ijJOyBpPkCF3SCtmedBewSaS/P24HpmVN\nJAmw0CWtgBjjf2KM/wV2A+pjjPUxxndy55JkoUtaASGEXiGELsA2McaXQgifz51JUuI5dEkrYn9g\nOvBICOFQ4I3MeSQ18vapkiQVgFPukiQVgIUuSVIBWOiSJBWAhS5JUgFY6JIkFYCFLklSAVjokiQV\ngIUuSVIB/D9maUVZKn4KDwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10b82e990>"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In terms of utility theory, we assume that each person's utilities are related to $X$ by\n",
      "\n",
      "$$\\begin{align*}\n",
      "U_0 | X\n",
      "    & = \\beta_0 X + \\varepsilon_0 \\\\\n",
      "U_1 | X\n",
      "    & = \\beta_1 X + \\varepsilon_1\n",
      "\\end{align*}.$$\n",
      "\n",
      "Here, $\\beta_i X$ represents the observed utility of option $Y = i$ ($i = 0, 1$), and $\\varepsilon_i$ represents btthe random fluctuation in the option's utility.  Logistic regression arises from the assumption that $\\varepsilon_0$ and $\\varepsilon_0$ are independent with [Gumbel distributions](http://en.wikipedia.org/wiki/Gumbel_distribution)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Gumbel distribution has the density function\n",
      "\n",
      "$$\\begin{align*}\n",
      "f(\\varepsilon)\n",
      "    & = e^{-\\varepsilon}\\ \\exp(-e^{-\\varepsilon}).\n",
      "\\end{align*}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To show that this utility model gives rise to logistic regression, we see that\n",
      "\n",
      "$$\\begin{align*}\n",
      "P(Y = 1 | X = x)\n",
      "    & = P(U_1 > U_0 | X = x) \\\\\n",
      "    & = P(\\beta_1 x + \\varepsilon_1 > \\beta_0 x + \\varepsilon_0) \\\\\n",
      "    & = P((\\beta_1 - \\beta_0) x > \\varepsilon_0 - \\varepsilon_1).\n",
      "\\end{align*}$$\n",
      "\n",
      "It is useful to note that the difference of independent identically distributed Gumbel random variables follows a [logistic distribution](http://en.wikipedia.org/wiki/Logistic_distribution#Related_distributions), so\n",
      "\n",
      "$$f(\\Delta) = \\frac{e^{-\\Delta}}{(1 + e^{-\\Delta})^2},$$\n",
      "\n",
      "where $\\Delta = \\varepsilon_1 - \\varepsilon_0$.\n",
      "\n",
      "Therefore\n",
      "\n",
      "$$\\begin{align*}\n",
      "P(Y = 1 | X = x)\n",
      "    & = P((\\beta_1 - \\beta_0) x > -\\Delta) \\\\\n",
      "    & = \\int_{-(\\beta_1 - \\beta_0) x}^\\infty \\frac{e^{-\\Delta}}{(1 + e^{-\\Delta})^2}\\ d\\Delta.\n",
      "\\end{align*}$$\n",
      "\n",
      "Making the substitution $t = 1 + e^{-\\Delta}$, with $dt = -e^{-\\Delta}\\ dt$, we get that\n",
      "\n",
      "$$\\begin{align*}\n",
      "P(Y = 1 | X = x)\n",
      "    & = \\int_1^{1 + \\exp((\\beta_1 - \\beta_0) x)} t^{-2}\\ dt \\\\\n",
      "    & = \\left.t^{-1}\\right|_{1 + \\exp((\\beta_1 - \\beta_0) x)}^1 \\\\\n",
      "    & = 1 - \\frac{1}{{1 + \\exp((\\beta_1 - \\beta_0) x)}} \\\\\n",
      "    & = \\frac{{\\exp((\\beta_1 - \\beta_0) x)}}{{1 + \\exp((\\beta_1 - \\beta_0) x)}}.\n",
      "\\end{align*}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Recalling that\n",
      "\n",
      "$$\\begin{align*}\n",
      "P(Y = 1 | X = x)\n",
      "    & = \\frac{e^{x}}{1 + e^{x}},\n",
      "\\end{align*}$$\n",
      "\n",
      "we must have that $\\beta_1 - \\beta_0 = 1$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The fact that there are infinitely many solutions of this equation is a subtle but important point of utility theory, that the difference in utility is both location and scale invariant.  We will prove this statement in two parts.\n",
      "\n",
      "1. For any $\\mu$, let $\\tilde{U}_i = U_i + \\mu$ for $i = 0, 1$.  Then\n",
      "$$\\begin{align*}\n",
      "\\tilde{U_1} - \\tilde{U_0}\n",
      "    & = U_1 + \\mu - (U_0 - \\mu)\n",
      "      = U_1 - U_0,\n",
      "\\end{align*}$$\n",
      "so $\\tilde{U_1} - \\tilde{U_0} > 0$ if and only if $U_1 - U_0 > 0$, and therefore the difference in utility is location invariant.\n",
      "2. For any $\\sigma > 0$, let $\\tilde{U}_i = \\frac{1}{\\sigma} U_i$ for $i = 0, 1$.  Then\n",
      "$$\\begin{align*}\n",
      "\\tilde{U_1} - \\tilde{U_0}\n",
      "    & = \\frac{1}{\\sigma}\\left(U_1 - U_0\\right),\n",
      "\\end{align*}$$\n",
      "so $\\tilde{U_1} - \\tilde{U_0} > 0$ if and only if $U_1 - U_0 > 0$, and therefore the difference in utility is scale invariant."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Together, these invariances show that the units of utility are irrelevant, and the only quantity that matters is the difference in utilities.  Due to the location invariance in utilities, we may as well set $\\beta_0 = 0$, so $\\beta_1 = 1$ for convenience.  Our utility model is therefore\n",
      "\n",
      "$$\\begin{align*}\n",
      "U_0\n",
      "    & = \\varepsilon_0 \\\\\n",
      "U_1\n",
      "    & = x + \\varepsilon_1\n",
      "\\end{align*}.$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To verify that this utility model is equivalent to logistic regression in a second way, we will simulate $P(Y = 1 | X = x)$ by generating  Gumbel random variables."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 500\n",
      "\n",
      "eps0 = gumbel_r.rvs(size=(N, n))\n",
      "eps1 = gumbel_r.rvs(size=(N, n))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U0 = eps0\n",
      "U1 = xs + eps1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "simulated_ps = (U1 > U0).mean(axis=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots(figsize=(8, 6))\n",
      "\n",
      "ax.plot(xs, logistic_probability(xs), c='k', lw=1.5, label='Exact');\n",
      "ax.plot(xs, simulated_ps, c='r', label='Simulated');\n",
      "\n",
      "ax.set_xlim(-5, 5);\n",
      "ax.set_xlabel('$x$');\n",
      "\n",
      "ax.set_ylabel('$P(Y = 1 | X = x)$');\n",
      "\n",
      "ax.legend(loc=2);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAF/CAYAAAC/oTuRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcjeX/x/HXZcYyzFhG9qRk38kSWUZSBtUXka1s+Skk\n2qRSSlKStEiSrSiJZMsSGiFr9t2MfaxjxCxmxsxcvz+OZBkazMx9zsz7+Xicx8w593XOvM+c9Jnr\nuq/7uoy1FhEREfFsmZwOICIiIrdPBV1ERCQdUEEXERFJB1TQRURE0gEVdBERkXRABV1ERCQdcKSg\nG2PGG2NOGGO23qDNZ8aYvcaYzcaYqmmZT0RExNM41UOfADS53kFjTFOghLW2JPB/wOi0CiYiIuKJ\nHCno1trlwJkbNHkMmHSx7RogtzGmQFpkExER8UTueg69CHD4svtHgDsdyiIiIuL23LWgA5ir7muN\nWhERkevwdjrAdYQCRS+7f+fFx65gjFGRFxGRDMVae3WHF3DfHvps4GkAY8z9wN/W2hNJNbTWetTt\n7bffdjxDer/pd6zfcXq46XecAX7PiYnYU6ewK1diJ0zgwiuvcL55c2JKlyY+a1Zi8uQhzseHQ1Wr\nsrRdOz4fMOCGhdORHrox5gegAXCHMeYw8DaQGcBaO8Za+6sxpqkxJhiIAro4kVNEROQKZ85ApkyQ\nK9d1m8TExHDq1CnCwsI4ffo0YWFhhIWFER4eTnh4OFkOHOB/q1dTISwMrCXYy4tdCQnsSkxkN7AH\nCAYiY2PxB5ps3EjzjRtpDzx/g2iOFHRrbbtktOmdFllERMSD7NwJwcHw6KNp+mPjo6KIGjyYHKNH\nQ1wccVmycCp3bg5nz84hICIqisioKKIiI4mNi2MzsBQ4f9lr5ADezZKFzvHxfFuoECMeeACTLx+5\n8+Qhd+7c5MqVi3q5c9M8Vy5y5sxJrotfc+bMiZ+fH9myZYPcua+b0V3PoadbAQEBTkdI9/Q7Tn36\nHac+/Y6vcu4cvPsuTJoEWbLAyZPQrdttv2xAQADWWk6ePMn+/fs5dOgQhw8fvvT18OHDlAoJYdCZ\nM+wA+gH7gcIxMZQ+d45qfn6UypYN3xw5KFCwINlz5MA3WzZ6nThBvsOHOVelCnGNG5Mjf378PvgA\nExAAw4bRt1Ah+t52+isZaz13XpkxxnpyfhER+Q/Wwvffw6uvwsMPwwcfuIp7QAB89BG0b5+sl4mL\ni+PAgQMEBwdfuoWEhLB//34OHDjA+fPnr2jv7+tLi7x56RYdTfHYWH5r1ozogAAKFy5M4cKFKViw\nIPny5SNz5szX/6F//w0LF8LcuXDkiOsPknr1buOXAcYY7HUmxaXLgm5Mku9V3IAn//cmImnsyBF4\n6ilXAf/iC6hd+99j27fDQw/BqFHQsuWlh8+dO8f27dvZvn07u3fvZteuXezatYv9+/eTkJAAuIam\nX8qaldbe3pzNm5eYokUxpUvjW7Eid0dGUnDdOjIvWwblykHr1tCzJ2TNmsZvPmkZsqB78vtKr/S5\niEiyzZvnGlJ//nl47TXw8rricHx8PIdmzaJQ16780Lgx06Oj2bZtG4cP/7smWdasWSlVqhRlypSh\ndOnSlCpVimpnz1J61Ci8ihbF9OsHx47Bnj2uW3AwlC7tOj8fGAj58qX1u/5PKujiFvS5iMh/iouD\n11+HadNgyhSoV4/4+Hh27NjBunXrWL9+PX/99Rdbt24lJiaGmsAcYMRdd3GkXj0qVKhAhQoVKFeu\nHMWKFcPrnz8EjhyBl1+GVavgk0+gRQvwwNFcFXRxC/pcROSGNmyAnj2J9fNj6dNPE7RtG3/++Scb\nNmwgOjoagJw5c3LfffdRtWpVqlSpQpUqVSgTE0PmDh2gQQP49FPInv3f14yNdRXwjz6CXr1cvf3L\nj3sYFXRxC/pcROQaiYmcmjiRhI8+IsuhQ3yWPTvvhIUBkCVLFqpVq0atWrWoUaMGNWrUoESJEmTK\nlMSaaBER8NxzsGmTq3dfrpxrQlqfPq5h9E8+gXvvTeM3l/JU0MUt6HMREYCwsDAWL17M2dGjefjP\nPzkdH8/HwNI8eahdvz5169alTp06VKtWzXXtdXJZCxMmQP/+ULkyHDjg6rE3a5ZabyXNqaC7ibvv\nvpuTJ0/+e04H6NKlC5999lmK/6ygoCCeeuqpKyaIOM1dPxcRSV3x8fGsWrWK+fPns2jRIjZs2MBD\n1jLRGEbXq0f+Vq0IaNiQ8uXLJ937vlnbt8OKFdCpE9zMHwQe4EYFXQvLpCFjDHPnzuXBBx90OoqI\nSKo6e/YsCxcuZM6cOfz666+Eh4fj5eVF7dq1GfHyy/QcPx6vH39kcKNGKf/Dy5d33TIYd92cJUN5\n7rnneOKJJy7d79+/Pw899BAAZ86coXnz5uTPnx9/f38effRRQkP/3XguPDycLl26UKRIEfz9/WnZ\nsiXR0dEEBgZy9OhR/Pz8yJkzJ8ePH0/z9yUiGUt4eDgTJ06kWbNm5MuXjyeffJL58+fTrFkzfvzx\nR06fPs3y33+n79q1ZHnhBbxSo5hnYOqhp7GkhpxHjBhBlSpVmDRpEsWLF2f8+PFs3rz5Uvtu3box\nffp04uPj6dq1K71792bmzJkAPPXUU+TMmZMdO3aQI0cOVq1aRfbs2VmwYAEdO3Z0qyF3EUl/IiIi\nmDFjBlOnTmXp4sUUSEigXv78TG7YkHu6daNaq1ZXnGZk4EDInNl1aZqkqAx3Dr1v375s2rTptn92\nlSpVGDly5E095+677+b06dN4e//7d9Tw4cPp1q0ba9eupUmTJuTMmZMPP/yQJ598MsnX2LRpEw8+\n+CDh4eEcO3aMO++8k/DwcHJdtfOPzqGLSGq5cOECv/32G9999x2zZs2i0vnzTMiShXuthVy5yFy+\nPKZwYVi0CBo2hJdegvvvh99+g86dXZenFSjg9NvwSDqH7iaMMcyaNSvJc+g1a9akePHihIWF0bp1\n60uPR0dH069fPxYuXMiZM2cAiIyMxFrL4cOH8ff3v6aYi4ikmBMnIDQUqlVj7969fPPNN0ycOJGT\nJ0/i7+9P586d+WDVKvy6dMF06QJ+fv8+NyICJk6EDh1cBfzAAddiMSrmqSLDFfSb7VWnlVGjRhEX\nF0fhwoUZNmwYr732GgAff/wxe/bsYe3ateTPn59NmzZRrVo1rLUULVqU8PBwzp49e01R13r2InLb\nLlwg8X//I27HDlpWqsT8FSvw8vLi0UcfpXPnzgQGBpJlzRpXT7xnT/C+qqT4+bmWbu3ZE2bNgvPn\nXT12SRUZrqA7Lakh5z179jBw4ECWLVuGj48PNWvWJDAwkMqVKxMZGYmPjw+5cuUiPDycd95559Lz\nChUqRGBgID179mTUqFGXzqHXr1+fAgUKcPr0ac6dO0fOnDnT8i2KiKc4eRIOHoQaNa45dPToUfY9\n8QTnN2zgWFwczbZupd7779O5c2cKFSr0b8MhQ1zXfV9dzC/n5XXFBiqSSqy1Hntzxb/W9R532t13\n3219fHysr6/vpVuLFi1szZo17Ycffnip3ejRo23FihVtXFycPXr0qA0ICLC+vr62dOnSdsyYMTZT\npkw2ISHBWmtteHi47dSpky1QoIDNkyePbdWq1aXX6dq1q82bN6/NkyePPXbsWJq/36u56+cikmG1\na2dtlizWTpp06aE1a9bY9u3b20czZbIHwbZ/+GEbNHmyTfT3t/bgwSufv26dtXfeaW1MTBoHz7gu\n/n80yZqY4SbFiXP0uYikoePH4dQpqFgx6eNbt7q2H507F9u6NcEPPkj3kBCW/fEHZX19WZ2YyLnx\n47nznwm6b70F+/bB5Mn/vkbLlq710194IfXfjwBaKU7chD4XkTQSFwd167qG07dvhzvuuLZNixYk\nPvAAM4oV45t332XYtm1szJGDc2+/Ta8ZM/D63/9cG5n8IzLStSb6zJlQs6brdRs1chV5D97sxNOo\noItb0Ocikkb69nXNKL/3Xte2oT/+eMXhxDVriGnalAfy5WPT7t2ULFmSN59/no4//0ymAwegTBnX\nfuRXL8M6frzrtnw5PPWUazW2AQPS7G2JCrq4CX0uImlg5kzo1w82bnStY16tGrzzDrRpg7WWuXPn\n4t+hA5MjIlhWtiyDBg2i1T+Lv8TFwccfwzPPQL581752QoLr9Tp0gGHDICQEdNlsmlJBF7egz0Uk\nle3fD7VqwZw5rq8Aa9bA44+zeswY+g4dSuY1a/je25sV33xDm44dr1zFLTkWL4bGjV0rvr37bsq/\nB7khFXRxC/pcRFLRP+fN27Vz9dAvCg4OZkuzZrBnDy8UKcK67NnJ178/Xt263frP+uILVy89T54U\nCC43QwVd3II+F5FU1K+fa4LaL7+AMZw5c4Z3332XUaNG4Zs5M7t8ffEPDMR79WrYtu3G142L29LS\nryIi6dmsWa5z5xs2YIFJEyfyyiuvEB4eTteuXRk8eDD5Q0Ndw/Dff69ink6phy5pRp+LSCo4cMBV\nqGfNYmuOHPTs2ZMVK1ZQu3ZtvvzyS6pUqfJv2z17oGRJ0NLQHutGPXTth+4mpkyZwiOPPJIqr925\nc2cGDhyYKq99tYkTJ1KvXr00+VkiGV5cHLRtS1zfvrwyYwZVq1Zl586djBs3jhUrVlxZzAFKlVIx\nT8dU0NPYihUrqFOnDrlz5yZv3rzUrVuX9evX06FDBxYuXJgqP9MYk+zNWgICAhg3blyq5BCRFDZg\nAGFeXpT75huGDx9Oly5d2L17N127diXT1deQS7qnEylp6Ny5czRv3pwxY8bQpk0bYmNjWb58OVmz\nZk31n53coW7t0ibihqy9pmd9/scfOf/115SKjCRP8eL8/vvvBAQEOJNP3IL+hEtDe/bswRjDk08+\niTGGbNmy0bhxYypWrHjNUHWmTJkYPXo0JUuWJGfOnLz11luEhIRQu3ZtcufOTdu2bblw4QKQ9DB3\npkyZ2Ldv3zUZzpw5Q/PmzcmfPz/+/v48+uijhIaGAvDGG2+wfPlyevfujZ+fH3369AFg165dNG7c\nmLx581KmTBl++umnS693+vRpHnvsMXLlykWtWrUICQlJ8d+bSIa2YYNr//AiRVxbj/boQUiPHkS1\nb8+jkZE8/cILbNmyRcVcVNDTUunSpfHy8qJz584sWLCAM2fO3LD9okWL2LhxI6tXr+bDDz+ke/fu\n/PDDDxw6dIitW7fyww8/3HQGay3dunXj0KFDHDp0CB8fH3r37g3AkCFDqFevHqNGjSIiIoLPPvuM\nqKgoGjduTMeOHTl16hRTp06lZ8+e7Ny5E4BevXqRPXt2jh8/zvjx45kwYYJ6+SIpZds2aNoURo+G\nVauI79+fWfv3M+/rr3k/Xz4+XL6ckSNHkiNHDqeTihtQQU9Dfn5+rFixAmMM3bt3J3/+/Dz++OOc\nPHkyyfavvvoqvr6+lCtXjooVKxIYGMjdd99Nzpw5CQwMZOPGjTedwd/fnxYtWpAtWzZ8fX15/fXX\nWbZs2RVtLh+enzt3Lvfccw+dOnUiU6ZMVKlShZYtW/LTTz+RkJDAzz//zLvvvouPjw/ly5enU6dO\nmskukhL27IFHHoFPPoFWrdgbG0vtgQP532+/sfWZZxgcEkLdunWdTiluJGMWdGNS5nYLypQpw4QJ\nEzh8+DDbtm3j6NGj9O3bN8lebYECBS597+Pjc8X9bNmyERkZedM/Pzo6mh49enD33XeTK1cuGjRo\nwNmzZ68owpdnOXjwIGvWrCFPnjyXbt9//z0nTpwgLCyM+Ph4ihYteqn9XXfdddOZROQq+/e7tjYd\nPBjbti0TJ06katWqhISEMH36dMaOHateuVwjYxZ0a1PmdptKly5Np06d2LZt200/9/KimyNHDqKj\noy/dP378+HXbf/zxx+zZs4e1a9dy9uxZli1bhrX2UkG/+g+Lu+66iwYNGnDmzJlLt4iICEaNGsUd\nd9yBt7c3hw4dutT+8u9F5BaEhrq2Je3fn+i2benSpQtdunShRo0abN68mVatWjmdUNxUxizoDtm9\nezcjRoy4NAnt8OHD/PDDD9SuXTtZz7+8F33595UrV2b79u1s3ryZmJgYBg0adM3z/mkfGRmJj48P\nuXLlIjw8nHfeeeeKtgUKFLhiYlvz5s3Zs2cPkydP5sKFC1y4cIF169axa9cuvLy8aNmyJYMGDeL8\n+fPs2LGDSZMm6Ry6yO0YMQIef5yQJk2oU6cOkyZN4q233mLx4sVXjIaJXE0FPQ35+fmxZs0aatWq\nha+vL7Vr16ZSpUp8/PHHwJW946SK4tXH/7lfqlQp3nrrLR566CFKly5NvXr1rtu2b9++nD9/njvu\nuIM6deoQGBh4RdsXXniB6dOn4+/vT9++ffH19WXRokVMnTqVIkWKUKhQIQYMGEBcXBwAX3zxBZGR\nkRQsWJCuXbvStWvXFPyNiWQw1sKsWSwrVoz77ruPQ4cOMW/ePN55552b3xVNMhwt/SppRp+LyI0l\nbtvGuQceIM+5c1SrVo3p06dzzz33OB1L3Ih2WxO3oM9F5PqioqKYef/9/L1tGxu6dOHLL78kW7Zs\nTscSN6O13EVE3FhoaCj169en+LZtFOrRg3HjxqmYy01TD13SjD4XkWutX7+exx57jGxnz7LbGDKH\nh0OWLE7HEjelHrqIiBuaPXs29evXJ0uWLPzx6qtkbtZMxVxumQq6iIgDvvnmG1q0aEGFChVYu3Yt\nd27YAI895nQs8WAacpc0o89FxLUuxJAhQxg4cCCBgYFMmzYN30yZoGBBOHgQ8uRxOqK4MQ25i4i4\ngYSEBHr16sXAgQPp1KkTs2bNwtfXFxYvhurVVczltqTb/dC1WpmIuJO4uDg6duzITz/9RP/+/Rk6\ndOi//5+aPVvD7XLb0uWQu4iIO4mNjaV169b8MWcOwwYO5P/effffgwkJULgwrFoFxYs7F1I8wo2G\n3NNtD11ExB2cP3+eFi1asHDhQrbVrUv5jz5y7dbYvz9kzw5r10L+/Crmctt0Dl1EJCVERsLff1/1\nUCTNmjVj0aJFTBo1ivLbt0NQEOzeDeXKwcyZMGuWhtslRaiHLiKSEgYPho0bYdEiACIiIggMDGT1\n6tV89913dIiJgbp1oVYtmDoVfv8deveGPXtg5UqHw0t6oHPoIiIp4f77XQV9xQqiy5cnMDCQlStX\nMnXqVJ544gnX8TffhObN/33OhQuwZAk88ohrGF7kP2S4zVlERNJUVJTrPPjbb5OwfDlNYmJYunQp\n33//PU8++SRs2QJNm8KBA+CtgVG5dZoUJyKSmlavhipViOvRg8i33uJ4bCzjJ050FXOAsWOhWzcV\nc0lV6qGLiNyut98mMSaGNiEh3DtjBs9Ur07Jdetcx6KjoWhR2LABihVzNqd4PK0UJyKSiuzy5YxY\nv54ZM2ZQbOhQSh44AMHBroPTp0PNmirmkurUQxcRuR1xccT6+pL/wgVeGTyYN998E95+G0JD4Ztv\noF49ePFFaNHC6aSSDugcuohIKvnhpZcofeECT/XqxRtvvOF6sE8fKFkS2rRx9dQvn9kukko05C4i\ncosmT57Mpi++4HiJEnz66af/rs2eNy888ww88QR06QKZMzsbVDIEDbmLiNyChQsX0rx5c5bnzEm1\nL74gS7t2VzY4fhwqVnTNgL/3XmdCSrrjdpPijDFNjDG7jDF7jTH9kzh+hzFmgTFmkzFmmzGmswMx\nRUSStHHjRlq1akXFcuWoFR9PlgcfvLZRwYJw9KiKuaSZNO+hG2O8gN3AQ0AosA5oZ63deVmbQUBW\na+0AY8wdF9sXsNbGX/Va6qGLSJoKDQ2lZs2aeHl58de4ceR7/nnYtcvpWJJBuFsPvSYQbK09YK29\nAEwFHr+qzTEg58XvcwKnry7mIiJpLSoqiscee4xz584xd+5c8u3c6ZrFLuIGnCjoRYDDl90/cvGx\ny40FyhtjjgKbgRfSKJuISJISExPp2LEjmzZtYurUqVSqVAmWL4f69Z2OJgI4U9CTM0b+OrDJWlsY\nqAKMMsb4pW4sEZHrGzBgAL/88gsjRoygWbNmYC388YcKurgNJ65DDwWKXna/KK5e+uXqAEMArLUh\nxpj9QGlg/dUvNmjQoEvfBwQEEBAQkLJpRSTD++6rr/hq2DCee+45+vTp43pw717ImlUrwEmqCgoK\nIigoKFltnZgU541rklsj4CiwlmsnxY0Azlpr3zHGFAD+AipZa8Ovei1NihORVLV69WpO1qlDgLc3\nOUaOxKtHD/Dycq0CFxQEkyc7HVEyELeaFHdxcltvYCGwA/jRWrvTGNPDGNPjYrP3gerGmM3AYuDV\nq4u5iEhqO3bsGKObNqWCtzf211/xmjrVtS77qlWu8+eaECduRAvLiIgkIS4ujqYNGjBuzRoSxoyh\nePfurvPmP/wAr74KYWGwcSOULet0VMlAbtRDV0EXEUlCz549uXP0aDo88ADFVqy48mBEBMyeDe3b\ng0ny/60iqUIFXUTkJowbN47hzzzDOh8ffENCoFAhpyOJANptTUQk2davX0/P555jbZ48ZH/7bRVz\n8RjqoYtIxrZpk2u2eqlSnC1QgPtataJ5ZCTDCxfGe+NG14x2ETehIXcRketp1Qqio7GJiZxcsYJc\n0dFkzpYNr2XLXDPaRdyIW122JiLiNmJiYPFi+O47PnnkEQpGRzN22DC8Dh5UMRePo3PoIpJxLV0K\nlSuzOjiY/v3706JFC3q//LJmrotH0pC7iGRczz5LVKFClB03Dm9vbzZs2EDu3LmdTiVyXZrlLiJy\ntcRE7Jw5vFS6NCdOnGDlypUq5uLRVNBFJGP66y/OJCQw5vffGTlyJNWrV3c6kchtUUEXkQzp1Lhx\nfBsWRmBg4L87qIl4MJ1DF5EMJyYmhoN58vBStmyM27WLAgUKOB1JJFl02ZqIyGU+fPZZ8sTE0PPb\nb1XMJd1QD11EMpT58+czv2lT2pcrx/3btzsdR+SmaKU4ERHg5MmTVKxYkVlRUVQbP54sbdo4HUnk\npmjIXUQyPGstPXr0gL//pgaQpVkzpyOJpCjNcheRDGHy5Mn88ssvzGnfHq9z5yBHDqcjiaQoDbmL\nSLp3+PBhKlasSMWKFVlWuDCZHnoIund3OpbITdM5dBHJsKy1PPLII6xcuZIt69dzb506sHMnFCzo\ndDSRm6alX0UkwxozZgy//fYbX375Jffu2wcVKqiYS7qkHrqIpFshISFUrlyZOnXqsHDhQkynTlCj\nBjz/vNPRRG6JhtxFJMNJTEwkICCALVu2sG3bNu684w4oVAi2b4fChZ2OJ3JLNOQuIhnO6NGjWb58\nORMmTODOO++E2bOhcmUVc0m31EMXkXTn4MGDVKhQgTp16rBgwQKMMdCxI9SpAz17Oh1P5JZpyF1E\nMgxrLYGBgaxYsYJt27Zx9913w/nzruH2Xbs0IU48mobcRSTD+Pbbb1m4cCGff/65q5gDLFgA1aqp\nmEu6ph66iKQbx48fp1y5cpQvX55ly5aRKdPF1a3btYMGDeDZZ50NKHKbNOQuIhnCE088wdy5c9m8\neTOlS5d2PRgd7ZoIt3cv5MvnbECR26QhdxFJ92bOnMmMGTMYOnTov8UcYP5817XnKuaSzqmgi4jH\nO3foEKefeoroTJnIEhYGERHg5+c6OG0aaJtUyQC0faqIeK7ERJgwAVu2LDYqiuDJk/E6dQrKloUf\nfoDISNeEuBYtnE4qkup0Dl1EPNPu3dC5M5GRkTTcto37e/fm888/dx1buRJ694Zz56BECVi40Nms\nIilEk+JEJP1p357EYsWoPn8+J06dYseOHeTKlevf4wkJMH48lCkD9eo5l1MkBWlSnIikL9bC0qWM\nf+YZNm7ezPTp068s5gBeXtrzXDIU9dBFxPNs386FwEBynz7Ngw8+yOzZs13Lu4qkczfqoWtSnIh4\nnqVL+f1iAf/iiy9UzEXQkLuIeKBj33/PhEOHGDRsGMWKFXM6johb0JC7iHiU85GRxOXKxWMlSrB4\n2zYyZ87sdCSRNKMhdxHxLHPnwvLlSR76tl8/jiQmMuirr1TMRS6jHrqIuJ+6dV2Xna1adcXDwcHB\njC9Thvr33kuT3bsdCifiHPXQRcRz/P03bNkCoaGwfv2lh6219OnTh0ZAzf79ncsn4qZU0EXEvSxZ\n4uqh9+oFo0ZdenjWrFksnj+fet7e+GspV5FraMhdRNxL9+5QoQJ06AAlS8LevURnz07ZsmUJ8PJi\ngr8/mS7ruYtkJFopTkQ8g7WuzVRefhnuuAMefxzGjeODqCgOHTrEO126kEnboIokSQVdRNzHjh3g\n7Q2lSrnu9+7Nhf/9j+GnTtG2bVvuDg6GJ590NqOIm9I5dBFxHwsWQJMm8M/Kb9Wrsz86mmbAR4MG\nwYYNrvPrInINFXQRcR8LF7oK+kVLlizh3TNn+LBoUe48eBCqVYMcORwMKOK+NClORNxDVBQULAhH\nj4KfH/Hx8VSpUoX4qCh2RkdjateGqlXh7bedTiriGE2KExH3t2wZ3Hcf+PkBMHr0aLZv387MmTMx\n69fDkCHw0ksOhxRxX+qhi4h76NMHCheG114jLCyMkiVLUr16dRYtWoQ5cgQefhg2b4YsWZxOKuIY\nrRQnIu7vnwlxwJtvvklERASffvqpa2vUokVdM+BVzEWuSwVdRJwXEgIREVC5Mlu2bGHs2LH06tWL\ncuXK/dtGe56L3JCG3EXEeV9+CWvXYidMoHHjxmzYsIHg4GD8/f2dTibiVjQpTkTc28KF0K4dc+fO\nZcmSJXz22Wcq5iI3ST10EXFWXBzky0fcrl1UaNAALy8vtmzZor3ORZKgHrqIuK81a6BkSb788Uf2\n7t3LvHnzVMxFboF66CLirHfe4fzp0xT+7jtq1arF/PnzXTPbReQa6qGLiPtasoTxuXMTERHBxx9/\nrGIucoscuWzNGNPEGLPLGLPXGNP/Om0CjDEbjTHbjDFBaRxRRNJCVBSJ69czYN48evToQfny5Z1O\nJOKx0nzI3RjjBewGHgJCgXVAO2vtzsva5AZWAo9Ya48YY+6w1oYl8VoachfxZIsWsa1dOx6Ijyc4\nOJh82uulffcOAAAgAElEQVRc5IbcbaW4mkCwtfaAtfYCMBV4/Ko27YEZ1tojAEkVcxHxfAcnTGB6\neDhvvPGGirnIbXKioBcBDl92/8jFxy5XEvA3xvxujFlvjHkqzdKJSJpITEwkctYstufPT58+fZyO\nI+LxnJgUl5wx8sxANaARkB1YZYxZba3dm6rJRCTNTBszhmbnz9Pqq6/Ili2b03FEPJ4TBT0UKHrZ\n/aK4eumXOwyEWWvPA+eNMX8AlYFrCvqgQYMufR8QEEBAQEAKxxWRlBYdHc1vAwdSLGdO2nTs6HQc\nEbcVFBREUFBQsto6MSnOG9ekuEbAUWAt106KKwN8ATwCZAXWAE9aa3dc9VqaFCfigYYMGYLvm2/y\naI8eFP/qK6fjiHiMG02Kc2RhGWNMIDAS8ALGWWuHGmN6AFhrx1xs8zLQBUgExlprP0vidVTQRTzM\niRMnKFGiBNuN4a6lS6F6dacjiXgMtyvoKUUFXcTzPPfcc8weO5ZDOXLgFR4OXl5ORxLxGCmyUpwx\nJgfQAaiAq2edDVfvORJYDfxkrU28/bgikl7t2rWLsWPHMrZRI7x8fFTMRVJQsgq6MaYxUA6Ya639\n+qpjBteEtReNMYuttZtSPqaIpAcDBgzAx8eHtvnzQ82aTscRSVf+8zp0Y0w2YL+19lNrbcjVx63L\nJmvtcCAhNUKKiAc7dAisZeXKlfzyyy/0798fnz//hAcfdDqZSLpy0+fQjTHFgWMXLylzlM6hi7i5\nLVugcmVsrVoMDA9nUkQEu5csIXvDhnD8OGgjFpGbkqKT4owxo3CdLw8yxtTF1UlfmQI5b5oKuoib\ne/NNOH+eVZkzc+HDD6nq749fgwaQNSv88IPT6UQ8Tkqv5b4WuMcYc4+1dgWQ/7bSiUj6ZC1Mm0Z8\n69Z0+vlnni1bFp+5c8HXF57Sas4iKe1WVoorCuzDNQmuAq5d0WamaCoR8XybN0N8PGM3bGDv3r3M\nnj0b79q1oXZtp5OJpEu3MuT+z05oscaYO4CWV898TysachdxYwMGEBsby11TplCmTBmCgoIwOmcu\ncltSesj9R6D8xe/vAQrcajARSacuDrdPjIri5MmTDBs2TMVcJJXd0kpxxpjc1tq/jTF5rLVnUiFX\ncnOohy7ijv76i/jWrcl94gSBTZvy008/OZ1IJF1I6R46QKeLX5++xeeLSHo2bRpL/P2JiY1lyJAh\nTqcRyRCc2D5VRNIza7nw/fe8fuwY3bt3p1SpUk4nEskQVNBFJGWtW8eps2fZlSULc996y+k0IhnG\nrQ65i4gk6fjnn/NNRAQvvvQShQoVcjqOSIahHrqIpBxryfTTTyzKlYtfX3nF6TQiGYoKuoikmLWf\nf45vbCytP/iAnDlzOh1HJEO51cvWyllrd/zzNRVyJTeHLlsTcROJ8fEszpeP7UDP48fJmjWr05FE\n0p0bXbZ2Sz30f4q4k8VcRNyItQQ3a0aOv/+m0NixKuYiDkjWpDhjzH2pHUREPJS1JLz8MrHLlvFq\n+fK06drV6UQiGVJyZ7m/cb0DxpgiKZRFRDzRu+/y99SpBMTG8sawYWTKpItnRJyQ3H95EcaY7saY\nK4bojTG5gPdTPpaIeIRhw0icMoUGcXFUqF+fwMBApxOJZFjJKujW2k7AFKC7MaaYMaaVMeZnYAtw\nf2oGFBE3NW8ejBrF5//7H9vDwvjwww+1AYuIg5I1y90Y0xEIBZ4CWgLbgPeAxUAppybHaZa7iEPC\nw6FiRc5++SVFn3qKRo0aMXPmTKdTiaR7N5rlntyCHgf8BkwGZgNlgLustY7+C1ZBF3FIx47g78+L\n3t58+umnbN26lXLlyjmdSiTdS4nL1l621n522f2/jDGHjTHdcf1R8PVtpxQRzzBzJqxZw6HZsxlV\npQqdOnVSMRdxA7e0sMylJxuTDVhura2RcpFu6uerhy6SlsLCoGJF+OknOn/zDVOnTmXv3r0ULVrU\n6WQiGUJq7IcOgLU2Bte5dBHJCHr1gvbt2ZY7N99++y29e/dWMRdxE7fVQ3eaeugiaWj6dHjzTdi4\nkcfbtiUoKIh9+/aRN29ep5OJZBgpsvSrMaY9kPkGTS5Ya7+/2XAi4iF++AEGDuTPjRuZPXs27733\nnoq5iBtRD11EkqdyZey4cTR48UX27NlDSEgIOXLkcDqVSIaS4puzXPXiOay1Ubf7OiLixqyFkBB+\n27eP5cuXM2rUKBVzETdz2z10Y0xfa+3IFMpzsz9bPXSRtHD8OLZiRaoULkxUVBQ7d+4kc+YbnYET\nkdRw2z10Y8wIoAFwLonDZQFHCrqIpJHgYMJy52bLli18//33KuYibijZC8sAfa21I64+YIzpl7KR\nRMTdxO/axYpjx6hSpQpPPvmk03FEJAnJ3ZwlEZhwncNaJU4kndv4009siopi6NCh2h5VxE1plruI\n3FBERASL8+VjR/HivL59u3ZUE3FQiq8UZ4zJffFrntsJJiLub8SIERSNjeXxl15SMRdxY7fUQzfG\nvGCt/fSfr6mQK7k51EMXSUUnT57k3uLFOR4XR45jx0ALyYg4KtXWcheR9O3999/HJzqabD4+4O/v\ndBwRuYHbXlhGRNKnAwcOMHr0aAY9+iheR4+ChttF3Jp66CKSpLfeeotMmTLR46GHoEQJp+OIyH9Q\nQReRa2zdupXJkyfTp08f/MPD4d57nY4kIv9BBV1ErjFgwABy5cpF//79IThYPXQRD5Csgm6MqX/V\nQ79d9VVE0olly5Yxb948BgwYgL+/vwq6iIdI1mVrxpg/gQBrbVzqR0o+XbYmkkKmTIHy5bGVK1O7\ndm2OHDnC3r178fHxgXz5YOtWKFjQ6ZQiGV5KbJ86HOhgjFlmrd2XctFExHGJifDaaxAYyMwmTViz\nZg3jxo1zFfO//4bz56FAAadTish/uKmFZYwxjYDM1toFqRcp+dRDF0kBf/4JLVpgM2emTPbseGfO\nzObNm/H29oa//oJu3WDTJqdTiggpsLCMMcYLwFq7BDhvjBlljKlrjMlujGmWgllFJK1Nmwa9enHu\nwgV89u7lgw8+cBVz0PlzEQ+S3HPoXwPHgP8BPsBcIAdQDbjLWuvIeJx66CK3KTERihYlevZsptSv\nj82Xj+779/+7ZvuQIRARAR984GxOEQFS5hz6g8A3QFtr7c6rXvyF28wnIk7580/Im5eRCxcSFB3N\ndD+/KzdgCQ6GBx5wLp+IJFtyr0N/1lr7wdXF/KKvUjKQiKShH38kqlkzPvzwQ3I2b07OQ4cgLOzf\n4yEhGnIX8RD/WdCNMVmB686IsdbGXtb2rhTKJSKpLSEBpk9nZGgoUVFRvPfRR9CwISy4bM6rzqGL\neIz/LOgXC/b9xpj2xhifpNoYY/IYY/4PKJbSAUUklaxYQWyePLwzdSrPPPMMZcqUgWbNYN481/Go\nKDhzBgoXdjaniCRLcifF+QH9Lt4tCiQAmS9+jQaOAGOttWdTKef1cmlSnMit6tWLqcuX0y0khODg\nYAoVKgShoVCxIpw8CTt2QLt2sH2700lF5KKUmBT3EXAWVzEvAjS11kalUD4RSWsJCcRNncqb4eG8\n/NZbrmIOUKQIFCsGq1bBqVMabhfxIMkt6FuttaMAjDGFgCeB8amWSkRSlQ0KYn9cHBH58/Pyyy9f\nefCfYXd/fxV0EQ+S3Fnulya+WWuPAedSJ46IpIVDw4czPjKSt99+Gz8/vysP/lPQNSFOxKMk9xx6\nMLAA2ABsBIpba2dcPFbAWnsiVVNeP5fOoYvcpPiYGM76+tKmaFEW7NlD5syZr2yQkODaiCV3bvjy\nS2jc2JmgInKN2176FZiEa3W4u4D3gM+NMauNMR/jOr9+s4GaGGN2GWP2GmP636BdDWNMvDGm5c3+\nDBFJQmIi+xo14s+EBHp9/PG1xRzAywuaNFEPXcTD3NTmLFc80Zh7gVpAd2ttw5t4nhewG3gICAXW\nAe2SWIHOC9d+69HAhH9GBK5qox66SHJZy4X/+z82TJrE69WqsXjVqitXhbvc1Knw9NMQHQ3eyZ1q\nIyKpLSVmuV/DWhsChBhjjtzkU2sCwdbaAxfDTQUeB65ehe55YDpQ41YzishF1sLLL3Py119pfOEC\ni0aOvH4xB2jaFPr3VzEX8SDJHXK/LmvtHzf5lCLA4cvuH7n42CXGmCK4ivzof37MLQcUEXj7bS4s\nWECtM2cIbNOG+++//8btc+aEwYPTJpuIpIjbLui3IDnFeSTw2sXxdHPxJiK34pNPYPp0Xq1ShZPx\n8QwdOtTpRCKSCpwYTwvFtUDNP4ri6qVf7j5g6sUhwTuAQGPMBWvt7KtfbNCgQZe+DwgIICAgIIXj\niniwTZtg6FB2Tp7Mp02a0K9fP4oXL+50KhFJpqCgIIKCgpLV9pYnxd0qY4w3rklxjYCjwFqSmBR3\nWfsJwBxr7c9JHNOkOJHriYuDGjXgxRdp8sMPrF27lpCQEPLkyeN0MhG5RakyKe5WWWvjjTG9gYWA\nFzDOWrvTGNPj4vExaZ1JJF0aPBiKFWNhgQIsXLiQTz75RMVcJB1L8x56SlIPXeQ61q+HZs1I+Osv\nqjZtSnR0NDt27CBLlixOJxOR2+BWPXQRSWUxMdCpE3zyCeN+/ZWtW7cyffp0FXORdE49dJH05rXX\nYO9ezo4bR8lSpShbtixBQUE3vu5cRDyCeugiGcW+fTBuHGzfzvtDhxIWFsaIESNUzEUyACeuQxeR\n1LJ0KTzyCPsiIxk5ciSdOnXivvvuczqViKQBDbmLpCedOkGdOjzx228sWLCAPXv2ULhwYadTiUgK\nSYnd1kTEE/zxB+t8fJgxYwavvfaairlIBqIeukh6cfgwtlo17rvzTsJOn2b37t34+Pg4nUpEUpAm\nxYlkBMuXc/Cuu9i4YQPff/+9irlIBqMhd5F0InbxYsbv2UPt2rVp27at03FEJI1pyF0knTh+xx00\nP32aMevXa2a7SDqlSXEi6dzeVavwOX2a+7p2VTEXyaDUQxfxcNZa3q1WjQe2baNSaCj58+d3OpKI\npBL10EXSsblz5+K3aRM+Dz+sYi6SgamHLuLBYmNjKV++PL+EhlLm11/xbtjQ6UgikorUQxdJpz75\n5BNOhIRQBvCuXdvpOCLiIPXQRTzUkSNHKFOmDC9XqsQgb2/44w+nI4lIKtPCMiLp0IsvvkhCQgJ9\nq1WD3LmdjiMiDtOQu4gHWrRoET/99BNvvPEGubdsgfr1nY4kIg7TkLuIp9i3D7p3J+6116jYuzeJ\niYlsXbeObEWKwPHj4OfndEIRSWUachdJDyZNAuB8mza8/fffFPruO7Jt3QrlyqmYi4iG3EU8grUw\nZQrHXniBe2NiyFy6NA379oUBA6BePafTiYgbUEEX8QRr10KmTPSaMIFoY6i1aBGsXg1FikDLlk6n\nExE3oCF3EU8wZQp7a9Zk5pQpvP/++9x1112ux3/80dlcIuI2NClOxN3Fx2MLF+bBLFk45uvL5s2b\nyZo1q9OpRMQBmhQn4smWLOGItzdBoaEsW7ZMxVxEkqRz6CJu7swXX/DJiRN07dqV+rreXESuQ0Pu\nIm4sMTKSqNy5uT9XLv7Ys4e8efM6HUlEHKQhdxEPtaRvXxITEnht5EgVcxG5IfXQRdzUsWPH2HTX\nXWwuWZL+27djTJJ/lItIBqLtU0U80MCePakTH0/rKVNUzEXkP6mgi7ihOXPm4P3LLxypUIF7q1Z1\nOo6IeACdQxdxM3+HhbG2Y0fe9/LCb/hwp+OIiIdQD13EnSxbRlSpUjQ4d46jP/xA5kcecTqRiHgI\nFXQRd/D339C+PTFt2vDCmTP89sorVGjd2ulUIuJBNMtdxGkREfDww1woU4Yqv//OhSxZ2Lx5Mz4+\nPk4nExE3o+vQRdxVdDQ89hhUqMBrvr7sOHiQP/74Q8VcRG6aCrqIU2JjXVufFinC6s6d+aRePXr2\n7Ek97W8uIrdAQ+4iTrhwAVq3Bm9vosePp2qNGsTExLBt2zb8/PycTicibkpD7iLuplcviI+HadMY\n8Mor7Nmzh6VLl6qYi8gtU0EXSWtr18LcubB7N7+vXMlnn31Gnz59aNiwodPJRMSDachdJC1ZC/Xq\nQZcunGvdmkqVKpE1a1Y2btxI9uzZnU4nIm5OQ+4i7uLnn12XqXXuzIs9enD48GFWrlypYi4it009\ndJG0EhsL5crBmDHMi42lefPmvP766wwZMsTpZCLiIW7UQ1dBF0krI0bA0qWcmjCBSpUqkT9/ftau\nXUvWrFmdTiYiHkJD7iJOCwuDoUOxy5bRtWtXzpw5w8KFC1XMRSTFqKCLpIV334U2bRgdFMTcuXP5\n9NNPqVSpktOpRCQd0ZC7SGrbvx9q1GDXzJlUffhhGjZsyLx58zAmyVEzEZHr0jl0ESd9+CHx+/ZR\nbdUqTpw4wZYtWyhQoIDTqUTEA+kcuoiTfv6Zr4sUYevWrcybN0/FXERShXroIqnp8GHiypcnR0QE\nzz3/PJ999pnTiUTEg6mHLuKQMxMmsCAujnKVKjFs2DCn44hIOqYeukgquXDhAlvvuIMP4uJ4b/Nm\nSpUq5XQkEfFwN+qhZ0rrMCIZxfv9+lH83DmeGDNGxVxEUp166CKpYO7cucx89FGeu/deqgcHOx1H\nRNIJnUMXSUMHDx7k6aefZnbOnFR6+22n44hIBqEeukgKio2NpUGDBoTu2MHBxEQyHTsGfn5OxxKR\ndELn0EVSg7UwaxYMHOjaSQ3o06cPa9asYUbXrmRq2FDFXETSjIbcRW7Fnj3wwgtw4ADccw8EBjLh\n8cf5+uuvGTBgADX37IGWLZ1OKSIZiGM9dGNME2PMLmPMXmNM/ySOdzDGbDbGbDHGrDTGaCcLcV5k\nJAwYAHXqwEMPwebNMGcOx/LmpXrfvnQICGDwgAHw22/w2GNOpxWRDMSRHroxxgv4AngICAXWGWNm\nW2t3XtZsH1DfWnvWGNME+Bq4P+3TilymeXMoWBC2boVChQA4duwY961YwQv+/kwKDsbr88+hRg3I\nm9fhsCKSkTg15F4TCLbWHgAwxkwFHgcuFXRr7arL2q8B7kzLgCLX2LYN9u6FxYvB2/VPJy4ujiee\neIJzERE0XbUKr02boGtXGDnS4bAiktE4VdCLAIcvu38EqHWD9t2AX1M1kch/GTvWVawvFnNrLT17\n9uTPP/9k2rRpVKxYESpWhGrV4N57HQ4rIhmNUwU92deaGWMaAl2BB1Ivjsh/OH8epkyB9esvPTR8\n+HDGjRvHm2++SevWrf9tW768AwFFJKNzqqCHAkUvu18UVy/9Chcnwo0FmlhrzyT1QoMGDbr0fUBA\nAAEBASmZU8RlxgyoXh3uvhuAX375hf79+9OmTRveeecdZ7OJSLoVFBREUFBQsto6srCMMcYb2A00\nAo4Ca4F2l0+KM8bcBSwFOlprV1/ndbSwjKSN+vWhb19o2ZINGzZQr149KlasyO+//46Pj4/T6UQk\ng7jRwjKOrRRnjAkERgJewDhr7VBjTA8Aa+0YY8w3QAvg0MWnXLDW1rzqNVTQJfXt2gUNG8KhQ4Se\nPEnNmjXx9vZm7dq1FChQwOl0IpKBuGVBTwkq6JImXn4ZMmcm4vXXadCgAXv37uXPP/90TYITEUlD\n2pxF5FbFxsK33xK3bBktW7Zky5YtzJkzR8VcRNyOCrrIjfzyC7ZSJToPHszixYuZMGECgYGBTqcS\nEbmGhtxFbsA2asS3WbPSef58hg4dymuvveZ0JBHJwHQOXeRWBAcTVaUK/lFRPPfCC3zyyScYk+S/\nIxGRNKGCLnILtjRtyqL58/mrbVumTJlCpkzabVhEnKX90EVu0rTJk8k/fz4769Rh4sSJKuYi4vbU\nQxe5ysyZM/n+iSd4w8+PUkePkj17dqcjiYgAumxNJNnmz5/Pk08+yR9+fpT+6CN8VMxFxENoHFEy\nplmzIDAQ4uMvPbR06VJatmzJw6VKUdPLC5+OHR0MKCJyc1TQJeM5dQqefRaOHoVhwwBYsmQJjz76\nKCVKlGBakyZk6tgRtEa7iHgQDblLxmItPPccPPUU9O4N993H8ty5afbii5QqVYrf5s8ne82asHCh\n00lFRG6KCrpkLD/+CDt2wOTJkC0bG598Et9evahctSq//vYbeVeuhGLFoEIFp5OKiNwUDblLxnH8\nOLzwAkyaBNmyMW3aNGp+9RXnc+UiqEkT8ubNC2PHQvfuTicVEblpKuiSMVjrOm/+zDNQowYTJ06k\nXbt21K5Th0pr1uDzzTcwZw6sXAlt2jidVkTkpqmgS8YweTLs24cdOJChQ4fSpUsXGjVqxPz58/Et\nXRqGD4cWLaBtW8iRw+m0IiI3TQvLSPoXGgpVq5L466/0/fZbPv/8c9q1a8fEiRPJkiWLq4218Npr\n0LUrlC7tbF4RkevQWu6ScVkLzZoRX60aHfbuZdq0afTr14/hw4drOVcR8ThaKU4yrgkTiD9yhKbR\n0fy2bBkfffQRL7/8stOpRERSnHrokn4dOkR8lSq0zJmTBUePMn78eDpq9TcR8WDqoUvGYy3hLVvy\n5fnzrPLyYunSpdStW9fpVCIiqUYFXdKl5U89hc9ffzGtbFnWzJ1L8eLFnY4kIpKqNCtI0pW4uDgG\nPv00ZaZMYWzduqxYvVrFXEQyBJ1Dl3Tj8OHDtG7dmn5r1pC7Th0aLVuGt7cGoUQk/dBla5LuLV68\nmHbt2lE5OprZPj5kP3QItJe5iKQzNyroGnIXj5aYmMiQIUN4+OGHyZ8vH3NKliT7iBEq5iKS4aig\ni8cKDQ2lcePGvPnmm7Rt25a/Xn0VHy8v0KVpIpIB6QSjeKSZM2fyzDPPEBMTw7hx4+jSrh2mbFnX\nTmpaAU5EMiD9n088SnR0NM8++ywtW7bknnvuYePGjXTt2hXz2WdQtSo0aOB0RBERR2hSnHiMP4KC\nWPnEEzxw+jQxdevS8OOPyVyjBpw6BeXKwapVULKk0zFFRFKNVooTjxYVFcV7fftS95tveDhrVrK8\n9x4VT5+GDh3g/HkoUACeekrFXEQyNPXQxa39/vvvfNWhA8OPHWP3ffdRe/FicuTO/W+DPXtg8WJo\n3x4uf1xEJB3SdejicU6ePMmrr77KHZMm8aqXFyeHDqXCK684HUtExFEq6OIxEhISGDt2LAMGDKB1\nRATv58xJjrVr8SlRwuloIiKO08Iy4hHWrVtHnTp1eO655+h89918mTs3d6xerWIuIpIMKujiuEOH\nDtGxY0dq1qzJwYMHmf3xx4w4ehTvH3+EUqWcjici4hFU0MUxERERvPHGG5QuXZoZM2bw+uuvs3ft\nWh4dMwbz7rvQqJHTEUVEPIYKuqQ8a6FrV9iyJcnDMTExjBw5khIlSvD+++/TqlUrdu/ezZBBg/Dr\n2hWaNoUePdI4tIiIZ1NBl5T3++8wcya0aQMREZcejouL46uvvqJEiRL069ePChUqsGbNGiZPnsxd\nRYpAly6QLRsMH+5geBERz6RZ7pLymjWDFi1g9WqIiSH2m2/49rvveP/99zlw4AAPPPAAgwcPpmHD\nhq72iYmuHnlwMMybp53SRESuQ5etSdrZuRMaNoQDB4g8d47zlSvzwfnzjDh7lurVqzN48GAeeeQR\njLn436O10KcP/PUXLFoEvr7O5hcRcWMq6HLrEhPhwgXImjV57bt3JyJPHob7+PD5559T4MwZVmXO\nzM4vvuD+7t3/LeTgKub9+8PSpbBkCeTKlTrvQUQknVBBl1s3YgR89RUsWwaFCt2w6aaFCyn52GOU\nTEzkWHw8jz/+OAMGDKDW7t3w/vuwfr3rHPmBA64lW+fNgz/+gKAgyJs3Td6OiIgnU0GXW2MtlC4N\n993nmrEeFAT58l3RJCYmhp9//plRo0bx8J9/cqe3N5ufe47nn3+ekpdvltK1K8ydC+fOQeHCruvL\ny5SBAQNcm6uIiMh/UkGXWxMUBM8/7yrmb74Jv/7qGh7Pk4ft27czduxYvvvuO8LDwyl3zz2sO3WK\nhN9/x6969WtfKz4eQkKgWDFXL11ERG6atk+VW/P119C9OxgD771HTHg4Z6tVo0O+fCxZt47MmTPT\nokULunfvzoN795Lp118hqWIO4O3t6u2LiEiqUA9dknb6NNx7L7E7dzJv1Sq+++475s2dy6fx8dT0\n8eFI69bU69YN/5o1IUsWKFsWxoyBgACnk4uIpFsacpekzZsHmzfD669f8XBMTAx7e/Xi7NKlNAsP\n59y5cxQoUID27dvTsX17qi5ejFm1yjWxbf9+8Pd3TZhbv97VmxcRkVShgi7XOnECKleGTJngiy84\n07AhCxYsYPbs2cydM4c1UVG86udHwTZtaN26NY0aNcLbO4kzNPHxcPAg5MgBBQum/fsQEclAVNAz\nosREV285qR6ztdgWLThdoAALs2Wj6ejR1EhIICQxkXz58vHS/ffTe+NGsoSEkDlLlrTPLiIiSdKk\nuIwkIuL/27u34KqqO47j378JlDiQhPs1lpMBjQkzhNQoAVTkNoGhSDs6GKc6Vh/ajrbI1LZeHnTG\nB+sVWhynrYpSa43VVtEZDTLYvIhSqRhFQkw0oYQEkQSSA2gkyb8P+8AkGDRIkk12fp+ZMyc7Z599\n/nvNSf57rbXXWvDUU7BmDaSkwOuvw3nnAdDQ0MDmzZtpfuwx5mzZwvRjx/gKuH/8eEqBvc8/z0Wz\nZpF0441w661B37iIiPQLqqFHxZ49sHYtrFsH8+bBqlUcefNNzlm9mrWFhTz30UeUlZUxzp0yMx64\n4gqmrljBokWLmPz978Py5ZCZCXffDbEYVFbCqFFhn5WIiHSgJvco27YNHnkELymhaflyNufksKmy\nki1btrBjxw6udmctcH9uLmnLl/PLjRtJmz+fc+69t/NxGhshLy/oV09JgeLiUE5HREROTQm9P2ts\nDCZ3GTIkmF0t8Tj07ru0PfQQSbW1vDhhAr///HM+OXAAgNTUVAoKCpg9ezZz587lkqNHGXzddcE6\n44jZCkgAAAeTSURBVGVlsHVr183p77wDl14KJSUwf34fn6iIiHwbJfT+qqkJX7CA5vPP55Nhwzha\nVsag6mpGNjSwv7WVNcBLwPnZ2eTn51NQUMCsWbPIzs4mKSmp87E+/DBYb/zJJ4Na+Kns2QOTJmn4\nmYjIWUgJPSxHj0JLS+ffpabCyckWaGtrY/fu3VRUVLBz504+3r6dX2zYwNYvv+Tnra0AJCUlkZWV\nRW5uLjNmzCA/P5+8vDyGaslREZEBQQm9r7nDo48GE7Z0GLvt7e0cy8xk2z33UH7gAFVVVVRWVlJR\nUUFlZSUtieQ/BHhj8GAOjxrFpquvJnvaNHJzc8nJySElJSWkkxIRkbApoZ+u1tZOibj7b2tlX3k5\nKTffDLW1vHjVVbwfj1NdXU1NTQ27a2q4s6WFa4BC4H/JycRiMbKysrjgggvIysoiKxYj/777GDx6\nNDzzTJe1eRERGZiU0LvLPRi/ffvtkJZ24ga09ilTOJSXx57x49m3bx/19fXU19dTV1fH3r17Tzyf\nV1fHs+3tvATcDnwFjBgxglgsxuTJk4nFYmRmZjJ31y6mFhfDhg0kz5wZfPaRI7B+ffD5ubnw7LMw\naFDPnZuIiPR7Z93EMmZWCKwBkoAn3P3+Lvb5I7AYOArc4O7be+Kz3Z14PE5jY+OJR0NDA/GaGi5d\nt47vNTWxetEi6g4eZGhdHSM++IAJ8Tg/Bg4DfwFeBtqB9PR0YuPGsXjIEBaeey4XpaTw1vXXc+Gy\nZbyXkUFGRgapqaldBzJvHixbFiTwHTvg8cdhzpxgHPns2bopTURETkuf19DNLAmoABYAe4F3gSJ3\nL++wzxLgFndfYmaXAH9w95ldHMtLSkqIx+M0NzcTj8dpamr62uPgwYMcOnToxHNbW1un48wCngP+\nATyQlkba6NGMGTOGMWPGMHbs2OB55EjyamvJKSkhpakJLypicHl5sGb4jBmwdClcey1MnPiN519a\nWsrc4yuSvf023HRTMERs5UqYMuW7Fqt00KmMpVeojHufyrhv9LdyPttq6BcDVe5eA2BmxcCVQHmH\nfZYB6wHcfauZpZvZWHf/7OSDFRYWMhVYChQADhwcPJjDw4bRkJ5O+/DhTBw+nJxYjLS0NEYOHUrM\nnYwvvmBsczOj9u8nvayMI2vWcOs113Dbt/WdP/hgMF77hRegqAiefjpYbaybOn15Cgpg585uv1e6\np7/9gfZHKuPepzLuG1Eq5zAS+kRgT4ftWuCSbuwzCfhaQj+akUFySwtfLliAXXYZKc3NJFVVBUt7\nVlRAXV3nNyQnB1OcHp+kZfFiWLiQ9PHju38GM2cGDxERkbNEGAm9u238JzcpdPm+lFdegenTGaQ+\nZxERGcDC6EOfCdzj7oWJ7TuA9o43xpnZn4BSdy9ObO8CLj+5yd3M+u8t+iIiIt/B2dSHvg2YamaT\ngTpgBVB00j6vALcAxYkLgENd9Z+f6qREREQGmj5P6O7eama3ABsJhq096e7lZvazxOt/dvfXzGyJ\nmVUBR4Cf9nWcIiIi/Um/nlhGREREAueEHcBAZWa/NrN2M+v+mDfpNjN70MzKzazMzP5lZmlhxxQV\nZlZoZrvMrNLMfhd2PFFjZhlm9m8z+8jMdpjZr8KOKarMLMnMtpvZq2HH0hOU0ENgZhnAQmB32LFE\n2BtAjrtPBz4G7gg5nkhITAz1KMFyBNlAkZldGG5UkXMMWOXuOcBM4GaVca9ZCeyk+6OvzmpK6OF4\nBPht2EFEmbtvcvf2xOZWgnkM5MydmBjK3Y8BxyeGkh7i7vvc/f3Ez4cJJt2aEG5U0WNmk4AlwBN8\nfZh0v6SE3sfM7Eqg1t0/CDuWAeRG4LWwg4iIriZ9+ub5juU7S4wGmkFwUSo9azXwG4KlOSIhlMVZ\nos7MNgHjunjpLoKm30Udd++ToCLoG8r5Tnd/NbHPXcBX7v73Pg0uuiLRNNkfmNlQ4EVgZaKmLj3E\nzJYC+919u5nNDTuenqKE3gvcfWFXvzezaUAMKLNgZrtJwH/N7GJ339+HIUbCqcr5ODO7gaBJbX6f\nBDQw7AUyOmxnENTSpQeZ2SDgn8Df3P3lsOOJoFnAssRCYEOAVDP7q7tfH3JcZ0TD1kJkZtXAD9y9\nMexYoiaxRO/DBDMMHgg7nqgws2SC1RLnE0wM9R9OWi1RzowFV/vrgQZ3XxV2PFFnZpcDt7n7D8OO\n5UypDz1cuprqPWuBocCmxLCUx8IOKArcvZVgFseNBHcHP69k3uNmAz8Brkh8d7cnLlCl90Tif7Fq\n6CIiIhGgGrqIiEgEKKGLiIhEgBK6iIhIBCihi4iIRIASuoiISAQooYuIiESAErqIiEgEKKGLiIhE\ngBK6iIhIBGhxFhHpNjNLAlYAmQTLqF4MPOzun4YamIiohi4ip2U6wSpgnxL8/3gBqA81IhEBlNBF\n5DS4+3vu3gIUAKXuXuruX4Qdl4gooYvIaTCzfDMbBUxz92ozmxN2TCISUB+6iJyOQuAz4C0z+xGw\nP+R4RCRBy6eKiIhEgJrcRUREIkAJXUREJAKU0EVERCJACV1ERCQClNBFREQiQAldREQkApTQRURE\nIkAJXUREJAL+DyTkemPbLilUAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10cc3df90>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The red simulated curve is reasonably close to the black actual curve.  (Increasing `N` would cause the red curve to converge to the black one.)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Aside from being an important tool in econometrics, utility theory helps shed light on logistic regression.  The perspective it provides on logistic regression opens the door to generalization and related theories.  If the random portions of utility, $\\varepsilon_1$ and $\\varepsilon_0$ are normally distributed instead of Gumbel distributed, the utility model gives rise to [probit regression](http://en.wikipedia.org/wiki/Probit_model).  For a thorough introduction to utility/choice theory, consult the excellent book [_Discrete Choice Models with Simulation_](http://eml.berkeley.edu/books/choice2.html) by Kenneth Train, which is freely available online."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}