Poisson Distribution and Poisson Regression

Poisson Distribution and Poisson Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson Distribution and Poisson Regression\n",
    "[ここ](https://gist.github.com/shotahorii/6e9275c01a5206075028e30018a12e56#file-discreteprobabilitydistribution-ipynb)でポアソン分布について簡単に説明し、[ここ](https://gist.github.com/shotahorii/c4f793a200f6b56683e383eae4352ea8#file-poisson-regression-ipynb)でポアソン回帰について書いた。  \n",
    "おさらいの意味も兼ねてまとめて簡単に整理しておく。\n",
    "\n",
    "**ポアソン分布**  \n",
    "ポアソン分布は、確率変数が0以上の整数であることを前提とした確率分布である。つまり、何らかの事象Xが0,1,2,3,...と非負の整数となり得るとき、ポアソン分布に従う可能性が考えられる。\n",
    "分布の形を決めるパラメータはλ一つだけで、ポアソン分布の平均と分散はλに等しくなる。\n",
    "\n",
    "この性質(平均=分散)を考えると直感的にポアソン分布の形を理解しやすい。つまり、λが小さい時=分布の平均が小さい時、分散も小さいので、確率変数Xがλ付近の値になる確率がとても高い。一方、λが大きい時、つまり平均の大きなポアソン分布では分散も大きいのでλ付近の値をとる確率は他と比べてそこまで高くはなくなる。\n",
    "\n",
    "また、λ>0である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "ポアソン分布の確率質量関数は以下である。<br><br>\n",
       "$$p(X=k) = \\frac{\\lambda^ke^{-\\lambda}}{k!}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "ポアソン分布の確率質量関数は以下である。<br><br>\n",
    "$$p(X=k) = \\frac{\\lambda^ke^{-\\lambda}}{k!}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "実際に複数のλでポアソン分布をグラフに書いてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import numpy.random as rd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define factorial function & poisson mass function\n",
    "factorial = lambda x: 1 if x==0 else x*factorial(x-1)\n",
    "poisson = lambda l,k: (l**k)*math.exp(-l)/factorial(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAGxCAYAAABY7ANPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcXHWV9/HPaUIQEnYCmAbakCjLDKAIMcpWIROmWcPj\nRqADE3CJzzNER2UEH226ezKKaIYBwwMaJ9o6xAmirCqSjFAOAQMtsggmJnQ6ZTZZDIQsLCF9nj9u\ndbpSqe6u6q66S9X3/XrVi7p1t1NF18mp3z33XnN3RERERKSy6qIOQERERKQWqOgSERERCYGKLhER\nEZEQqOgSERERCYGKLhEREZEQqOgSERERCYGKLqk4M+s2syOLXLYhu7z+NkUkFpTDpFz0R1HjzKzL\nzM6s8G5KvRhcWS8eZ2a7m9kd2ffabWanl3P7IhKdWshhuczs2mweq/R7lgpQ0SVhsKgDAB4GmoD1\nUQciIokThxxGdrTto8C6qGORwVHRJbsws/3M7D4ze9HM/pp9Xp8z/yEzm2Vmj5jZJjO7x8wOMLPb\nzGyjmT1mZkfkbfZcM+vMbvObOduqM7PZZvaSmT0PnJsXy3Qz+6OZvWZmz5vZp0t9P+6+zd2/7e6P\nAt2lri8iyVJtOSzH/wO+BGwbwjYkQiq6pJA64PvA4cARwFbg5rxlLiIYORoNjAMeBeYB+wPLgJa8\n5S8ETsw+ppjZFdnXPw2cA5wAnETwKy7XC8A57r4PcDnw72b2XgAzO9zMXjGzDdn/5j7fYGZTh/Qp\niEhSVV0OM7OPAW+4+68G95FIHAyLOgCJH3ffANyVnXzTzK4Dfp232A/cfRWAmd0PHOPuD2Wn7wD+\nJW/5b7j7RmCjmd0IXEyQFD8G3Oju67LrXgeckRPL/TnPHzazhcBpwFPuvpogQYqI7FBtOczMRgJf\nAyYV8fYlxjTSJbswsz3N7LtmtsrMXgV+A+xnZrl9DS/kPH+9wPTIvM2uyXmeIfh1Sfa/q/Pm5cZy\ntpn9NnuI4BXgbOCgkt+UiNSMKsxhrcCPskWaJJiKLinki8C7gZPdfT+g52y/oTSTHp7zvIHeRtD1\nBeYFOzMbDvwU+CYwyt33B+7viSM7NL8p2yuR++h57eIhxCsiyVVtOWwS8FkzW29mPfv7iZn98xDe\nj0RAhxcFYLiZ7ZF9bgTD3a8Dr5nZAQS/sobqn83scWBv4LPA7OzrPyFIJr8g6Lu4Ojeu7ONld+82\ns7OBs4A/AGR/9e1dzM6zya/nR8YeZraHu785xPckIvFQ7TnsTGD3nOnfAf8EqL8rYTTSJQA9yeJ1\nYAuwL7An8DJBc+kv85YfzDVr7gGeAH4P3EfQCwHwPeAB4GmCRPKzHSu5byZIbneY2QZganY7g/En\ngvc2miBRbS1wdpKIJFNV5zB3f8XdX+x5AG8Dr7r71lK3JdEy94H/9sysEbiRoEib5+7X97HcyQR/\n4Be5+52lrCsiUgnKXyISFwMWXRbcymA5wTHldUAHMNXdlxVYbhHBL43vu/udxa4rIlIJyl8iEifF\nHF4cD6xw94y7bwMWAFMKLDeToGHwxUGsKyJSCcpfIhIbxRRd9ex8Ouya7Gs7mNlo4EJ3v5Wdzw4Z\ncF0RkQpS/hKR2CjX2Ys3svMZGyUzs4rdIFRE4snd43BPO+UvERmUUnNYMSNdawluo9DjsOxruU4C\nFphZF8EtEG4xswuKXHcHd4/9o6WlJfIYFKdirYY4Q6L8lcC/jaTEmaRYFWf5H4NRzEhXBzDOzBoI\nLgI3leD2B7nJ5sie52b2A+A+d7/XzHYbaF0RkQpS/hKR2Biw6HL37WZ2JbCQ3tOml5rZjGC2z81f\nZaB1yxe+iEjflL9EJE6K6uny4K7mR+W99t0+lr0ib3qXdZMslUpFHUJRFGf5JSXWpMQZFuWvXkn5\n20hKnJCcWBVnPBR1cdQwmJnHJRYRqTwzw+PRSD9kyl8itWcwOUy3ARIREREJgYouERERkRCo6BIR\nEREJgYouERERkRCo6BIREREJQayKrmnT2ujqyuz0WldXhmnT2pg4saXgfBEREZEkiNUlI2AzY8e2\nsGjRTMaMaaCrK8PkyXPo7GwDRgBbdpovIsmlS0aISJJVwSUjRtDZ2caXv9zOli3w5S+35xRcvfOb\nm9sjjFFERESkdEVdkT5cI7jjjm7uuw/eeKOb3oKrd/66dd1RBCYiIiIyaDEb6QLYwsUX17FlC1x8\ncR2wZZf5o0fHMGwRERGRfsSsegl6tmbNmg7ArFnTGTu2hd7Ca+f5IiIiIkkRq0b6pqZWZs2avlOT\nfFdXhpkz23nwwW4+/OG6XeaLSDKpkV5EkmwwOayoosvMGoEbCUbG5rn79XnzLwBmAd3ANuDz7v5I\ndt4qYGPPPHcf38c++kxaGzfCYYfBpk1FvisRib2wiq6o85eIVKeKFF1mVgcsByYB64AOYKq7L8tZ\nZi9335p9fhzwE3c/Jju9Eni/u78ywH76TFrusPfesG4d7LNP0e9NRGIsjKIrDvlLRKpTpS4ZMR5Y\n4e4Zd98GLACm5C7Qk7CyRhL8KtwRV5H76ZMZjB4dFF0iIiWIPH+JiPQoJpnUA6tzptdkX9uJmV1o\nZkuB+4ArcmY5sMjMOszsU4MNtL4e1q4d7NoiUqNikb9ERKCM1+ly97uBu83sVOBfgcnZWae4+3oz\nG0WQvJa6++JC22htbd3xPJVKkUqldkxrpEsk2dLpNOl0OuowCqp0/hKR5CtHDiump2sC0Orujdnp\nawDPb0bNW6cTONndN+S93gJscvcbCqzTb0/El74EBxwA11zTb7gikhAh9XTFIn+JSPWpVE9XBzDO\nzBrMbDgwFbg3b8djc56fCAx39w1mtpeZjcy+PgI4C3i2lAB7aKRLRAYhFvlLRASKOLzo7tvN7Epg\nIb2nXC81sxnBbJ8LfMTMLgPeAl4HPp5d/RDgruBm1gwD5rv7wsEEWl8PDz88mDVFpFbFJX+JiEDM\nLo7aXyyPPAJXXQW//W2IQYlIxejiqCKSZJU6vBgLOntRREREkiwxI11vvhlcIPWNN6AuMaWiiPRF\nI10ikmRVPdK1xx6w777w0ktRRyIiIiJSusQUXaAzGEVERCS5ElV0qa9LREREkipRRZdGukRERCSp\nElV0aaRLREREkipRRZdGukRERCSpElV0aaRLREREkipRRZdGukRERCSpElV0aaRLREREkipRRdeo\nUbBxY3B1ehEREZEkSVTRVVcHhx4K69dHHYmIiIhIaRJVdIH6ukRERCSZiiq6zKzRzJaZ2XIzu7rA\n/AvM7Gkze9LMHjezU4pdt1Tq6xKRUsQpf4lIbTN3738BszpgOTAJWAd0AFPdfVnOMnu5+9bs8+OA\nn7j7McWsm7MNHygWgJkzYdw4+NzninyHIhJLZoa7W4X3Eav8JSLVYzA5rJiRrvHACnfPuPs2YAEw\nJXeBnoSVNRLoLnbdUmmkS0RKEKv8JSK1rZiiqx5YnTO9JvvaTszsQjNbCtwHXFHKuqVQ0SUiJYhV\n/hKR2jasXBty97uBu83sVOBfgcmlbqO1tXXH81QqRSqV2mUZNdKLJFM6nSadTkcdRkFh5S8RSa5y\n5LBieromAK3u3pidvgZwd7++n3U6gZOB9xS7brE9EcuWwQUXwPLlAy4qIjEWUk9XrPKXiFSPSvV0\ndQDjzKzBzIYDU4F783Y8Nuf5icBwd99QzLql6hnpUn4TkSLEKn+JSG0b8PCiu283syuBhQRF2jx3\nX2pmM4LZPhf4iJldBrwFvA58vL91hxLwPvuAGbz2Guy771C2JCLVLm75S0Rq24CHF8NSyvD8UUfB\n3XfDMcdUOCgRqZgwDi+GRYcXRWpPpQ4vxo7OYBQREZGkSWTRpTMYRUREJGkSWXRppEtERESSJpFF\nl0a6REREJGkSWXRppEtERESSJpFFl0a6REREJGkSWXRppEtERESSJpHX6XrrLRg5Et54A+oSWTaK\niK7TJSJJVjPX6Ro+HPbbD158MepIRERERIqTyKIL1NclIiIiyZLYokt9XSIiIpIkiS26NNIlIiIi\nSZLYoksjXSIiIpIkRRVdZtZoZsvMbLmZXV1g/iVm9nT2sdjMjs+Ztyr7+pNm9ni5AtdIl4gUI475\nS0Rq07CBFjCzOuBmYBKwDugws3vcfVnOYiuB0919o5k1AnOBCdl53UDK3V8pZ+Aa6RKRgcQ1f4lI\nbSpmpGs8sMLdM+6+DVgATMldwN2XuPvG7OQSoD5nthW5n5JopEtEihDL/CUitamYZFIPrM6ZXsPO\nSSnfJ4H7c6YdWGRmHWb2qdJD7CMojXSJyMBimb9EpDYNeHixFGY2EbgcODXn5VPcfb2ZjSJIXkvd\nfXGh9VtbW3c8T6VSpFKpPvd10EHw2mvw5puwxx7liF5EKimdTpNOp6MOo09h5i8RSZ5y5LABbwNk\nZhOAVndvzE5fA7i7X5+33PHAz4BGd+/sY1stwCZ3v6HAvJJvo9HQAOk0jBlT0moiEgNh3AYozvlL\nRJKtUrcB6gDGmVmDmQ0HpgL35u34CIKEdWluwjKzvcxsZPb5COAs4NlSAuyP+rpEZACxzV8iUnsG\nPLzo7tvN7EpgIUGRNs/dl5rZjGC2zwWagQOAW8zMgG3uPh44BLjLzDy7r/nuvrBcwauvS0T6E+f8\nJSK1Z8DDi2EZzPD8Zz8LRx4J//RPFQpKRComjMOLYdHhRZHaU6nDi7GlkS4RERFJikQXXerpEhER\nkaRIdNGlkS4RERFJikQXXRrpEhERkaRIdNHVM9Kl/lURERGJu0QXXXvvDXV1wZXpRUREROIs0UUX\nqK9LREREkiHxRZf6ukRERCQJEl90aaRLREREkiDxRZdGukRERCQJEl90aaRLREREkiDxRZdGukRE\nRCQJEl90aaRLREREkqCoosvMGs1smZktN7OrC8y/xMyezj4Wm9nxxa47VPX1GukSkb7FOX+JSG0x\nH+By7mZWBywHJgHrgA5gqrsvy1lmArDU3TeaWSPQ6u4Tilk3Zxs+UCyFvPUWjBwJr78Ou+1W8uoi\nEhEzw92twvuIdf4SkeQaTA4rZqRrPLDC3TPuvg1YAEzJXcDdl7j7xuzkEqC+2HWHavhw2H9/eOml\ncm5VRKpErPOXiNSWYoquemB1zvQaepNSIZ8E7h/kuoMyerT6ukSkoNjnr4F0dWWYNq2NiRNbmDat\nja6uzKCWEZHoDSvnxsxsInA5cOpg1m9tbd3xPJVKkUqlilqvp6/r/e8fzF5FJAzpdJp0Oh11GH2K\nKn/1p6srw+TJc+jsbANGAFtYsqSFRYtmMmZMQ9HLiMjQlSOHFdPTNYGgx6ExO30N4O5+fd5yxwM/\nAxrdvbOUdbPzBt0T8elPw4knwmc+M6jVRSQCIfV0xT5/9WfatDbmz7+KoJjqsYUTTpjNRRe1AHD7\n7W08/fSuyzQ1zea221rKHpOIBCrV09UBjDOzBjMbDkwF7s3b8REECevSnoRV7LrloDMYRaQPsc9f\n/Vm7tpudiymAEbzySjevvQavvQavvFJ4mXXrusMJUkSKNuDhRXffbmZXAgsJirR57r7UzGYEs30u\n0AwcANxiZgZsc/fxfa1b7jcxejQsWVLurYpI0iUhf/Wnvr4O2EL+KNZpp9Vx3XXB1OrVdcyfv+sy\no0cn/jKMIlVnwMOLYRns8HxXV4YrrmjnmWe6OfvsOmbNmq4+BpEECOPwYlgqdXixqyvD6afPYc2a\n3n6tsWMH7unKX0ZEym8wOSzRRZeSjUhyqegqTlNThmefbefAA7sZPbrwD8uurgzNze10dnbz+9/X\n8dhj03nve5UDRSqp5oquvppM1UAqEn8quga2ZQscfjg88wwcdlhx61x+OYwZA9deW/ZwRCRHpRrp\nY6uvJlM1kIpINfjpT+FDHyq+4AL48pdhzpygyV5E4iXRRVdvk2kuNZCKSHX4j/+AT36ytHXe8x6Y\nPBluuaUyMYnI4CX68KJ6ukSSS4cX+/enP0EqBX/+M+y+e2nrPvccnHkmrFwJI/IPBohIWdRcTxf0\nNpA+8kg3I0bUcd99OntRJAlUdPXvS18CM7h+l0uxFuejH4VTToHPf76sYYlIVk0WXT1uvx3uuCPo\ngRCR+FPR1be33goa6B9+ODhcOBhPPgnnnQednfCOd5QtNBHJqrlG+lxHHw3LlkUdhYjI0P3850FO\nG2zBBfC+9wW3R5s3r3xxicjQVM1I19atcOCBsHkz7LZbGQMTkYrQSFffzjkHLr4YLr10aNt5/PHg\nMOPzz8Pw4eWJTUQCNT3StddecMghsGpV1JGIiAze6tXw2GPwkY8MfVvjx8Mxx8CPfjT0bYnI0FVN\n0QU6xCgiydfeDhddFPyQLIevfhWuuw7efrs82xORwRvwhtdJ0lN0nXtu1JGIiJSuuzvowbrzzvJt\n87TT4KCDMpx+ejt77NFNfb3uUSsSlaoquo46Cn7/+6ijEBEZnF//Gg44IGiAL5eurgyrV89h/fre\n6xkuWaLrGYpEoajDi2bWaGbLzGy5mV1dYP5RZvaomb1hZl/Im7fKzJ42syfN7PFyBV7I0UcHFxQU\nEemRlPwFg7sC/UCam9tzCi6AEXR2ttHc3F7eHYnIgAYc6TKzOuBmYBKwDugws3vcPbd76q/ATODC\nApvoBlLu/koZ4u2XerpEJFeS8tfLL8MDD8B3v1ve7eoetSLxUcxI13hghbtn3H0bsACYkruAu7/s\n7k8AhVo1rcj9DNmhh8Kbb8Jf/xrG3kQkARKTv267DS64APbbr7zb1T1qReKjmG9dPbA6Z3pN9rVi\nObDIzDrM7FOlBFcqs6CvS4cYRSQr9vmrqytDU1MbX/lKC+vXt9HVlSnr9mfNms7YsS30Fl7BPWpn\nzZpe1v2IyMDCaKQ/xd3Xm9koguS11N0XF1qwtbV1x/NUKkUqlSp5Zz19XR/60CCjFZGKSKfTpNPp\nqMMoVUXzV1dXhsmT59DZGfRc/fd/b2Hy5PI2uY8Z08CiRTNpbp7N8uXdPPdcHQ88oCZ6kVKVI4cN\neEV6M5sAtLp7Y3b6GsDdfZfbsJpZC7DJ3W/oY1t9zi/XFZ2//nXYuHHwN4kVkXCEcUX6uOevadPa\nmD//KnbuudpCU9NsbrutpeTtFeOEE+Cmm2AQv2lFJEelrkjfAYwzswYzGw5MBe7tL46cgPYys5HZ\n5yOAs4BnSwmwVGqmF5Ecsc5fUTS5X3aZrlAvEpUBDy+6+3YzuxJYSFCkzXP3pWY2I5jtc83sEOB3\nwN5At5l9DjgWGAXcZWae3dd8d19YqTcDKrpEpFfc81dvk/vOI12VbHK/5BI49li4+ebyXfVeRIpT\nNTe87vHmm7DvvrBpE+y+exkCE5GK0A2vg56uD3xgDi+91Hvh0rFjK3/h0sZG+Id/CG6qLSKDM5gc\nVnVFF8C73w333ReMeolIPKnoCpx7boYNG9rZc89uRo8O5xY9P/4x/Od/wv33V3Q3IlVNRVfW+efD\nJz4BFxa61KGIxIKKruBei4ceCh0d0BDiyYRbt0J9PSxdGuxfREpXqUb6xFFfl4gkwbPPwj77hFtw\nQdDLdeGFwYiXiISnKosuXSBVRJLgwQdh0qRo9n3ppTqLUSRsVVl0aaRLRJLg17+OruhKpWDDBnjm\nmWj2L1KLqrroikm7mojILt5+Gx5+GCZOjGb/dXXQ1BQ01ItIOKqy6DroINhtN3jxxagjEREp7He/\nC3q5Ro2KLoZLL4X582H79uhiEKklVVl0gfq6RCTeHnwQzjwz2hiOPTY4i/HXv442DpFaUbVFl/q6\nRCTOouznyqWGepHwqOgSEQnZ66/DY4/B6adHHQlMnQo//3lwFw8RqSwVXSIiIfvtb+G444JrdEXt\n4IPhtNPgzjujjkSk+lVt0aWeLhGJq7gcWuxx2WU6xCgShqq8DRDAtm2w997w6qvwjneUbbMiUia1\nfBugD34Qvva16Bvpe7zxBhxySIYzz2zn1Ve7qa8P5x6QIklWsdsAmVmjmS0zs+VmdnWB+UeZ2aNm\n9oaZfaGUdStl993hyCNhxYqw9igicRS3/PXaa8Htfz70oXJsrTzWr88Ac7j77qtIp9uYP/8qJk+e\nQ1dXJurQRKrKgEWXmdUBNwN/D/wNcLGZHZ232F+BmcC3BrFuxaivS6S2xTF//c//wPjx8RqBb25u\n57XX2oAR2VdG0NnZRnNze4RRiVSfYka6xgMr3D3j7tuABcCU3AXc/WV3fwJ4u9R1K+noo9XXJVLj\nYpe/4tbPBbB2bTe9BVePEaxb1x1FOCJVq5iiqx5YnTO9JvtaMYay7pAddZRGukRqXOzyVxwuipqv\nvr4O2JL36hZGj67ac61EIjEs6gBytba27nieSqVIpVJD2t7RR8OcOUOLSUTKI51Ok06now6jYorJ\nXy++CJkMnHRSeHEVY9as6SxZ0kJnZ88hxi2MHdvCrFkzI45MJD7KkcMGPHvRzCYAre7emJ2+BnB3\nv77Asi3AJne/YRDrlvXsRQjOXDz88KBx1ariHCmR6hHG2Ytxy1+33x7c6/Dee4fyriqjqytDc3M7\nDz7YzUEH1XHPPTp7UaQ/g8lhxYx0dQDjzKwBWA9MBS7uL44hrFtW++0HI0fCunXB/cVEpObEKn/F\n8dBijzFjGrjtthb+8Af4+79XzhSphAEP2Lv7duBKYCHwHLDA3Zea2Qwz+zSAmR1iZquBzwNfMbM/\nm9nIvtat1JspRH1dIrUrbvkrjk30+Y47Dt7zHl2hXqQSqvbiqD0+85kgifzjP5Z90yIyBLV2cdRM\nBk4+Gf7yF6iLeX/6nXfCv/0bPPJI1JGIxFfFLo6aZLpWl4jEQc+hxbgXXAAXXABr1sATT0QdiUh1\nScDXf2h0rS4RiYM493PlGzYM/vf/1tnfIuVW9YcXu7rgjDPgz38u+6ZFZAhq6fCie9CY/vDDMHZs\niIENwcsvw7hxwa3URo2KOhqR+NHhxQKOOCJIHps3Rx2JiNSqP/0Jhg8P7gebFAcdBB/+MHzve1FH\nIlI9qr7o2m03ePe7YfnyqCMRkVrT1ZVh2rQ2Lrighd12a2PVqmTdQHrmTLj1Vng7/wZJIjIoVV90\ngZrpRSR8XV0ZJk+ew/z5V7FiRRsrV17F5Mlz6OpKTuH1vvfBu94Fd98ddSQi1aEmiq6jjlIzvYiE\nq7m5Pee2OgAj6Oxso7m5PcKoSjdzphrqRcqlJooujXSJSNjWru2mt+DqMYJ167qjCGfQ/tf/gs5O\nePrpqCMRST4VXSIiFVBfXwdsyXt1C6NHJyvt7r57cJFpjXaJDF3VXzICgjMXDz44+G8SLkwoUguq\n/ZIRXV0Z/u7v5rByZc8hxi2MHdvCokUzE3cj6RdfDNo0nn8eDjww6mhE4mEwOawmiq6urgx/+7ft\nHH98N2PH1jFr1vTEJT2RalPtRRfAD3+Y4YtfDHLP6NHJzj0f/nCGrq529tuvm/r6ZL8XkXJQ0VVA\nzxlEvQ2tyf21KVJNaqHouvxyOOEE+Kd/iiCoMurqynD66XNYs0Z5VKSHLo5aQLWcQSQiyfLGG3DP\nPfDxj0cdydA1N7fnFFygPCoyOFVfdFXLGUQikiy/+lUwyjV6dNSRDJ3yqEh5FFV0mVmjmS0zs+Vm\ndnUfy3zbzFaY2VNm9r6c11eZ2dNm9qSZPV6uwItVLWcQicjgRJW//uu/4OKLhxp9PCiPipTHgD1d\nZlYHLAcmAeuADmCquy/LWeZs4Ep3P9fMPgDc5O4TsvNWAu9391cG2I96ukRqSBg9XVHlr82bgxtc\nr1xZHWf7FcqjRxzRQjqtPCq1azA5bFgRy4wHVrh7JruTBcAUIPfKV1OAHwG4+2Nmtq+ZHeLuLwBG\nhIcxx4xpYNGimTQ3z2b16m4ee6yOefOUKERqRCT565574NRTq6Pggp3z6Lp13WzcWMfuu8/kXe9S\nHhUpRTFFVz2wOmd6DUEi62+ZtdnXXgAcWGRm24G57h76PevHjGngtttaALjqKvjlL+GMM8KOQkQi\nEEn+WrCgeg4t9sjNo9u3wwc/CN//PnziExEHJpIgxRRdQ3WKu683s1EEyWupuy8utGBra+uO56lU\nilQqVfZgPvMZ+NCHoK0N3vGOsm9eRPqQTqdJp9NRh1GqkvPX66/Dgw+m+PGPU+FFGbLddoPvfQ8m\nT4ZzzoF3vjPqiEQqrxw5rJierglAq7s3ZqevAdzdr89Z5jvAQ+5+e3Z6GXBGdng+d1stwCZ3v6HA\nfip2cdR8jY3Q1ASXXhrK7kSkgJB6ukLPX9/7HixcCHfcUal3FR9f+QosX14b71UkX6Wu09UBjDOz\nBjMbDkwF7s1b5l7gsmwQE4BX3f0FM9vLzEZmXx8BnAU8W0qAlfB//g/cckvUUYhICELPX9V01uJA\nmpvhmWeCHjYRGVhRV6Q3s0bgJoIibZ67f8PMZhD8YpybXeZmoJHgvOLL3f33ZjYGuIugL2IYMN/d\nv9HHPkIb6dq+HY48Eu66C048MZRdikiesK5IH2b+Wr8ejj0W1q+vnfaF3/wGpk2DZ5+FffeNOhqR\n8Og2QCX4+teD07n/4z9C26WI5KjG2wDddBM8+SS0t0cdUbg+/WkYNkxHEKS2qOgqwQsvwFFHQVcX\n7L9/aLsVkaxqLLomTIDW1qBvtJa8+ir8zd/A7bcHl8oQqQUqukp0ySUwfnzyb0YrkkTVVnStXOl8\n4AOwdi3svnvUEYXvzjvhn/85w/jx7fzlL93U19cxa9Z0XRNRqpaKrhItXgxXXAHLlkGd7mYhEqpq\nK7q+/nVn9eraPcTW1ZXhb/92Dlu36u4fUhsqdfZi1TrllKDZ9cEHo45ERJKuls5aLKS5uT2n4AIY\nQWdnG83N7RFGJRIvNV10menyESJSHq+8EvyQq1Vr13bTW3D1GMG6dd1RhCMSSzVddEFwkdR0Gtas\niToSEUmyiy6q7TaF+vo6gitu5NrCoYfW8Icikqfmvw177x001M+dG3UkIpJkzz3XRldXJuowIjNr\n1nTGjm2ht/Dawp57trB+/XQ2b44wMJEYqelG+h4PPJDhwgvb+cAHujnsMJ1xIxKGamukh8013zje\n1ZWhublFaXDgAAAgAElEQVSddeu6GT26jmuvnc63vtXAE0/Az38Oo0dHHaFI+ejsxUHo6sowefIc\nOjt1xo1ImKqv6HJgC01Ns7nttpaoQ4oNd7j+erj1VvjOdzLMn9/O2rW6pIQk32By2LBKBZMUzc3t\nOQUX9J5xo8QpIqVS43g+M7jmGthzzwznnTeH7u7eH7hLlugHrtSWmu/p0hk3IlI+Wxg9uubTakEd\nHe05BRfokhJSi2p+pKv3jJvcwmsLnZ11vPgiHHxwRIGJSMIErQmzZs2MOpBY6usH7tKl3bgHI2I9\nPWE6/CjVquZ/khU642bMmBYaG6dz3HHwox8FPQldXRmmTWtj4sQWpk2r7bOURGRXTU2zdaisH31d\nUmL58jpOPBGuuy7DpElzmD//KtLpNubPv4rJk+co10p1cfcBH0AjsAxYDlzdxzLfBlYATwHvLWXd\n7HIelZUrV3lTU6tPnHitNzW1+sqVq9zd/Ykn3N/7XvdTT13lRxzxRYfNDg85bPaxY7+4Y7k4euih\nh6IOoShJidM9ObEmJc7sd76oHDSUR7Xnr1JE+bexcuUqHzu2J4/6jjz6/POr/P773evrW3PmPbRj\nmaam1l2209TU6qnUzvm62PnllpTvm+Isv8HksGISVh3wPNAA7J5NSkfnLXM28Ivs8w8AS4pdN2cb\nIXxEpXvrLfcTTshNBi0Fk0ExX/RyJItit9HQcHrF4khKnOXeRl+xhhFHUuIsdhn3cIquWs9f+Vpa\nWiLdf18/cN3dU6lrs7k1N8+6jx59rX/nO+6PPOL+9NOFC7ee7fRV2NXy961a4gw7jmKK9koVXROA\n+3Omr8n/xQd8B7goZ3opcEgx6+bM6/ONRa2vZHDIIdd6S4v77NmrvL6+/y/6QMmgmGRR2jZaKhJH\nUuKszDZ2jTWMOJISZ7HL9Aip6Kr5/JUr6qKrP01NhX/cnnhiq19xhftJJ7nvtlvuMr5jmQkTWv2e\ne9z/7u8Kz8//gVwr37dqiTOaOPrPX+6VK7o+AszNmZ4GfDtvmfuAD+VMLwJOLGbdnHkF31Qc9JUM\nJkxo9a9+1b2hofAXfY89Wv2ww9yPPNJ9n30KL3Pgga3+wQ+6H3RQ4fmjRrX6aae5n3aa+6hRhZc5\n+OBWP/1094MPLhxnz/xdlxloG6XPL30blYmzMtvYNdYw4khKnP0tk3+IyH1wCavUh/LXzuJcdBXz\nD+8ZZ1yb97cVPA455Fo/7zz3ffYpPN/sWt97b/cDD3R/xzsK/43uu2+rn3CC+/7795+rd83Xvd+3\nUaNa/dRT+87Vg8/nlcm1UeeFpPybUCh/uQ8uh1Uqaf33YJKWHnroUVuPUhNWyQlO+UsPPfSo4KPU\nnFTMJSPWAkfkTB+WfS1/mcMLLDO8iHUhiLwqrkwtIrGi/CUisVHMJSM6gHFm1mBmw4GpwL15y9wL\nXAZgZhOAV939hSLXFRGpFOUvEYmNAUe63H27mV0JLCQo0ua5+1IzmxHM9rnu/kszO8fMnie4EMvl\n/a1bsXcjIpJD+UtE4iQ2N7wWERERqWaRX5HezBrNbJmZLTezq6OOpy9mtsrMnjazJ83s8ajjyWVm\n88zsBTN7Jue1/c1soZn9ycweMLN9o4wxG1OhOFvMbI2Z/T77aIwyxmxMh5nZg2b2nJn9wcw+m309\nVp9pgThnZl+P42e6h5k9lv3+/MHMWrKvx+ozLVVS8hfEN4clJX9BMnJYUvJXH7HGMoeVM39FOtJl\nZnUEV3qeBKwj6KGY6u7LIguqD2a2Eni/u78SdSz5zOxUYDPwI3c/Pvva9cBf3f2b2X8M9nf3a2IY\nZwuwyd1viDK2XGZ2KHCouz9lZiOBJ4ApBIedYvOZ9hPnRcTsMwUws73cfauZ7QY8AnyW4AzB2Hym\npUhS/oL45rCk5K9sXLHPYUnJXwPEGrscVq78FfVI13hghbtn3H0bsIDgA48jI/rPqyB3XwzkJ9Ip\nwA+zz38IXBhqUAX0EScEn21suPtf3P2p7PPNBBfLPIyYfaZ9xFmfnR2rzxTA3bdmn+5B0E/qxOwz\nLVGS8hfENIclJX9BMnJYUvIXJCuHlSt/Rf0FrAdW50yvofcDjxsHFplZh5l9KupginBw9gws3P0v\nwMERx9OfK83sKTP7jzgMeecys3cB7wWWAIfE9TPNifOx7Eux+0zNrM7MngT+Aixy9w5i/JkWIUn5\nC5KVw5KUvyCG3zdITv6C+OewcuWvqIuuJDnF3U8EzgH+MTvMnCRxPWPiFuBId38vwR9znIaTRwI/\nBT6X/RWW/xnG4jMtEGcsP1N373b39xH86h5vZn9DTD/TKpXkHBbnv4tYft+Skr8gGTmsXPkr6qKr\nmAsXxoK7r8/+9yXgLoJDC3H2gpkdAjuOm78YcTwFuftL3ttY+D3g5Cjj6WFmwwiSwH+6+z3Zl2P3\nmRaKM66faQ93fw1IA43E8DMtQWLyFyQuhyXm7yKO37ek5C9IXg4bav6KuuhKxMUHzWyvbCWOmY0A\nzgKejTaqXRg7HwO/F5ieff4PwD35K0Rkpzizf6g9Pkx8PtfvA39095tyXovjZ7pLnHH8TM3soJ5D\nBGa2JzCZoH8jjp9psRKRvyAROSwp+QuSkcOSkr8gATmsnPkr8ut0ZU8FvYneiw9+I9KACjCzMQS/\nDJ2ggW5+nOI0sx8DKeBA4AWgBbgbuIPg9iYZ4OPu/mpUMUKfcU4kOI7fDawCZvQcI4+KmZ0C/A/w\nB3rvsfV/gceBnxCTz7SfOC8hfp/pcQSNpnXZx+3u/jUzO4AYfaalSkL+gnjnsKTkL0hGDktK/oLk\n5LBy5q/Iiy4RERGRWhD14UURERGRmqCiS0RERCQEKrpEREREQqCiS0RERCQEKrpEREREQqCiS0RE\nRCQEKrpEREREQqCiSyrCzLrN7Mgil23ILq+/RxGJnPKXVIr+SGqQmXWZ2ZkV3k2pV90t61V6zWx3\nM7sj+167zez0Astcb2Yvm9lLZhaLq3OLSP9qJH99wMwWmtlfzewFM7s979Y4yl8JpaJLKsUGXqTi\nHgaagPX5M8xsBnABcBxwPHC+mX063PBEJKaizl/7A98FGrKPzcAPemYqfyWXii4BwMz2M7P7zOzF\n7K+r+8ysPmf+Q2Y2y8weMbNNZnaPmR1gZreZ2UYze8zMjsjb7Llm1pnd5jdztlVnZrOzv9CeB87N\ni2W6mf3RzF4zs+cHk0zcfZu7f9vdHyW4f1e+y4B/c/f17r4emE3vjUtFJEGqMH/9yt1/5u6b3f0N\n4GbgQzmLKH8llIou6VFHcLf3w4EjgK0EX/RcFxGMHI0GxgGPAvMIfpUtI7j5a64LgROzjylmdkX2\n9U8D5wAnACcBH81b7wXgHHffB7gc+Hczey+AmR1uZq+Y2Ybsf3OfbzCzqUW+378Bns6Zfjr7mogk\nT7XnrzOA53Kmlb8SaljUAUg8uPsG4K7s5Jtmdh3w67zFfuDuqwDM7H7gGHd/KDt9B/Avect/w903\nAhvN7EbgYoLE+DHgRndfl133OoKk0hPL/TnPHzazhcBpwFPuvpogSQ7VSGBjzvRr2ddEJGGqOX+Z\n2fFAM3B+zsvKXwmlkS4BwMz2NLPvmtkqM3sV+A2wn5nl9ja8kPP89QLT+V/6NTnPMwS/MMn+d3Xe\nvNxYzjaz32YPE7wCnA0cVPKb6t9mYJ+c6X2zr4lIwlRr/jKzccAvgZnZVokeyl8JpaJLenwReDdw\nsrvvB/Sc7TeUhtLDc543AOuyz9cXmBfszGw48FPgm8Aod98fuL8njuzw/KZsv0Tuo+e1i4uM7TmC\nwwM93svOw/cikhxVl7/MrAFYBLS5+4/zYlP+SigdXqxdw81sj+xzIxjyfh14zcwOAFrLsI9/NrPH\ngb2BzxI0ewL8BPismf2CoPfi6ty4so+X3b3bzM4GzgL+AJAdnt+7mJ1nE2DPD4s9zGwPd38zO/0j\n4AvZwwwGfAG4cXBvU0RCVtX5K3sSwK+BOe7+vQKLKH8llEa6aldPwngd2EIwPL0n8DJBg+kv85Yf\nzHVr7gGeAH4P3EfQDwHwPeABgubP3wE/27GS+2aCBHeHmW0Apma3Mxh/Inhvo4FfAVt7zlBy9+9m\nY/pDNo57+0huIhI/1Z6/PgGMAVpzR8Jy9qP8lVDmPvDfopk1ElTRdcA8d7++j+VOJviDv8jd7yxl\nXRGRSlD+EpG4GLDosuDWBsuBSQTHtDuAqe6+rMByiwh+eXzf3e8sdl0RkUpQ/hKROCnm8OJ4YIW7\nZ9x9G7AAmFJguZkEDYQvDmJdEZFKUP4SkdgopuiqZ+fTY9dkX9vBzEYDF7r7rex8tsiA64qIVJDy\nl4jERrnOXryRnc/gKJmZlfWGoSISf+4e9T3uQPlLRAap1BxWzEjXWoLbKvQ4LPtarpOABWbWRXBL\nhFvM7IIi193B3WP/aGlpiTwGxalYqyHOkCh/JfBvIylxJilWxVn+x2AUM9LVAYzLXqhtPcEpsDtd\ngNLdj+x5bmY/AO5z93vNbLeB1hURqSDlLxGJjQGLLnffbmZXAgvpPW16qZnNCGb73PxVBlq3fOGL\niPRN+UtE4qSoni53/xVwVN5r3+1j2SvypndZN8lSqVTUIRRFcZZfUmJNSpxhUf7qlZS/jaTECcmJ\nVXHGQ1EXRw2DmXlcYhGRyjMzPB6N9EOm/CVSewaTw3QbIBEREZEQqOgSERERCYGKLhEREZEQqOgS\nERERCYGKLhEREZEQqOgSERERCYGKLhEREZEQqOgSERERCYGKLhEREZEQqOgSERERCYGKLhEREZEQ\nqOgSERERCUFRRZeZNZrZMjNbbmZXF5h/gZk9bWZPmtnjZnZKzrxVufPKGbyIyECUv0QkLszd+1/A\nrA5YDkwC1gEdwFR3X5azzF7uvjX7/DjgJ+5+THZ6JfB+d39lgP34QLGISPUwM9zdKrwP5S8RqYjB\n5LBiRrrGAyvcPePu24AFwJTcBXoSVtZIoDs3riL3IyJSbspfIhIbxSSTemB1zvSa7Gs7MbMLzWwp\ncB9wRc4sBxaZWYeZfWoowYqIlEj5S0RiY1i5NuTudwN3m9mpwL8Ck7OzTnH39WY2iiB5LXX3xYW2\n0drauuN5KpUilUqVKzwRiVg6nSadTkcdRkHKXyIykHLksGJ6uiYAre7emJ2+BnB3v76fdTqBk919\nQ97rLcAmd7+hwDrqiRCpISH1dCl/iUhFVKqnqwMYZ2YNZjYcmArcm7fjsTnPTwSGu/sGM9vLzEZm\nXx8BnAU8W0qAIiJDoPwlIrEx4OFFd99uZlcCCwmKtHnuvtTMZgSzfS7wETO7DHgLeB34eHb1Q4C7\nzMyz+5rv7gsr8UZERPIpf4lInAx4eDEsGp4XqS1hHF4Mi/KXSO2p1OFFERERERkiFV0iIiIiIVDR\nJSIiIhICFV0iIiIiIVDRJSIiIhICFV0iIiIiIVDRJSIiIhICFV0iIiIiIVDRJSIiIhICFV0iIiIi\nIVDRJSIiIhICFV0iIiIiIVDRJSIiIhICFV0iIiIiISiq6DKzRjNbZmbLzezqAvMvMLOnzexJM3vc\nzE4pdl0RkUpS/hKRuDB3738BszpgOTAJWAd0AFPdfVnOMnu5+9bs8+OAn7j7McWsm7MNHygWEake\nZoa7W4X3ofwlIhUxmBxWzEjXeGCFu2fcfRuwAJiSu0BPwsoaCXQXu66ISAUpf4lIbBRTdNUDq3Om\n12Rf24mZXWhmS4H7gCtKWVdEpEKUv0QkNoaVa0Pufjdwt5mdCvwrMLnUbbS2tu54nkqlSKVS5QpP\nRCKWTqdJp9NRh1GQ8peIDKQcOayYnq4JQKu7N2anrwHc3a/vZ51O4GTgPcWuq54IkdoSUk+X8peI\nVESlero6gHFm1mBmw4GpwL15Ox6b8/xEYLi7byhmXRGRClL+EpHYGPDwortvN7MrgYUERdo8d19q\nZjOC2T4X+IiZXQa8BbwOfLy/dSv0XkREdqL8JSJxMuDhxbBoeF6ktoRxeDEsyl8itadShxdFRERE\nZIhUdImIiIiEQEWXiIiISAhUdImIiIiEQEWXiIiISAhUdImIiIiEQEWXiIiISAhUdImIiIiEQEWX\niIiISAhUdImIiIiEQEWXiIiISAhUdImIiIiEQEWXiIiISAiKKrrMrNHMlpnZcjO7usD8S8zs6exj\nsZkdnzNvVfb1J83s8XIGLyIyEOUvEYmLYQMtYGZ1wM3AJGAd0GFm97j7spzFVgKnu/tGM2sE5gIT\nsvO6gZS7v1Le0EVE+qf8JSJxUsxI13hghbtn3H0bsACYkruAuy9x943ZySVAfc5sK3I/IiLlpvwl\nIrFRTDKpB1bnTK9h56SU75PA/TnTDiwysw4z+1TpIYqIDJryl4jExoCHF0thZhOBy4FTc14+xd3X\nm9koguS11N0XF1q/tbV1x/NUKkUqlSpneCISoXQ6TTqdjjqMPil/iUh/ypHDzN37X8BsAtDq7o3Z\n6WsAd/fr85Y7HvgZ0OjunX1sqwXY5O43FJjnA8UiItXDzHB3q/A+lL9EpCIGk8OKObzYAYwzswYz\nGw5MBe7N2/ERBAnr0tyEZWZ7mdnI7PMRwFnAs6UEKCIyBMpfIhIbAx5edPftZnYlsJCgSJvn7kvN\nbEYw2+cCzcABwC1mZsA2dx8PHALcZWae3dd8d19YqTcjIpJL+UtE4mTAw4th0fC8SG0J4/BiWJS/\nRGpPpQ4vioiIiMgQqegSERERCYGKLhEREZEQqOgSERERCYGKLhEREZEQqOgSERERCYGKLhEREZEQ\nqOgSERERCYGKLhEREZEQqOgSERERCYGKLhEREZEQqOgSERERCYGKLhEREZEQFFV0mVmjmS0zs+Vm\ndnWB+ZeY2dPZx2IzO77YdUVEKkn5S0Tiwty9/wXM6oDlwCRgHdABTHX3ZTnLTACWuvtGM2sEWt19\nQjHr5mzDB4pFRKqHmeHuVuF9KH+JSEUMJocVM9I1Hljh7hl33wYsAKbkLuDuS9x9Y3ZyCVBf7Loi\nIhWk/CUisVFM0VUPrM6ZXkNvUirkk8D9g1xXRKSclL9EJDaGlXNjZjYRuBw4dTDrt7a27nieSqVI\npVJliUtEopdOp0mn01GH0SflLxHpTzlyWDE9XRMIehwas9PXAO7u1+ctdzzwM6DR3TtLWTc7Tz0R\nIjUkpJ4u5S8RqYhK9XR1AOPMrMHMhgNTgXvzdnwEQcK6tCdhFbuuiEgFKX+JSGwMeHjR3beb2ZXA\nQoIibZ67LzWzGcFsnws0AwcAt5iZAdvcfXxf61bs3YiI5FD+EpE4GfDwYlg0PC9SW8I4vBgW5S+R\n2lOpw4siIiIiMkQqukRERERCoKJLREREJAQqukRERERCoKJLREREJAQqukRERERCoKJLREREJAQq\nukRERERCoKJLREREJAQqukRERERCoKJLREREJAQqukRERERCoKJLREREJARFFV1m1mhmy8xsuZld\nXWD+UWb2qJm9YWZfyJu3ysyeNrMnzezxcgUuIlIM5S8RiYthAy1gZnXAzcAkYB3QYWb3uPuynMX+\nCswELiywiW4g5e6vlCFeEZGiKX+JSJwUM9I1Hljh7hl33wYsAKbkLuDuL7v7E8DbBda3IvcjIlJu\nyl8iEhvFJJN6YHXO9Jrsa8VyYJGZdZjZp0oJTkRkiJS/RCQ2Bjy8WAanuPt6MxtFkLyWuvviQgu2\ntrbueJ5KpUilUiGEJyJhSKfTpNPpqMMolfKXiADlyWHm7v0vYDYBaHX3xuz0NYC7+/UFlm0BNrn7\nDX1sq8/5ZuYDxSIi1cPMcHer8D6Uv0SkIgaTw4o5vNgBjDOzBjMbDkwF7u0vjpyA9jKzkdnnI4Cz\ngGdLCVBEZAiUv0QkNgY8vOju283sSmAhQZE2z92XmtmMYLbPNbNDgN8BewPdZvY54FhgFHCXmXl2\nX/PdfWGl3oyISC7lLxGJkwEPL4ZFw/MitSWMw4thUf4SqT2VOrwoIiIiIkOkoktEREQkBCq6RERE\nREIQxnW6pIy6ujI0N7ezdm039fV1zJo1nTFjGoqeLyIiItFQI32CdHVlmDx5Dp2dbcAIYAtjx7aw\naNFMxoxpGHC+SJyokV5EkkyN9FWuubk9p6ACGEFnZxszZrTzyCMwY0bh+c3N7ZHEKyIiIr10eDFB\n/vznbnoLqh4jeOyxbr70JXj22cLz163rDidAERER6ZNGuhLiwQfh97+vA7bkzdnC+efX8cgjcP75\nhed3d+t/s4iISNT0r3HMbdwIM2bAP/wD3HTTdMaObaG3sAp6tmbNmg7ArFm7zn/nO1vIZKZz/vnQ\n2Rm82tWVYdq0NiZObGHatDa6ujJhviUREZGapEb6GMk/83DSpOm0tDTQ2Ajf+hbsu2/vMuvWdTN6\ndN9nL+bOHz26gX//d5g9Gy6+OMMvfjGHri4120u01EgvIkk2mBymoismCp15OGxYC9///kwuvbQ8\nxdDq1XDqqW38+c9XsXPv1xaammZz220tZdmPSDFUdIlIkunsxQQrdGbi22+38cAD7WXbx+GHw5FH\nqtleREQkCiq6YmLt2nCKofr6ws32o0frT0FERKSSivqX1swazWyZmS03s6sLzD/KzB41szfM7Aul\nrCuBgw8Opxgq1Gx/wAG9zfgi1Ub5S0TiYsCeLjOrA5YDk4B1QAcw1d2X5SxzENAAXAi84u43FLtu\nzjZqtidi+3aYPDnDE0/M4bXXKt/gnttsv99+dTz11HSuuKKBr361rLsR6VcYPV3KXyJSKYPJYcVc\nHHU8sMLdM9mdLACmADsSj7u/DLxsZueVuq7AVVeBWQMdHTP5l3+ZnXPmYWXOKBwzpmGnpvn16+HM\nM6G7G669tuy7E4mS8peIxEYxRVc9sDpneg1BMirGUNatCd/5Dtx/P/z2t7D//g2RnEH4znfCQw8F\nhZc7tOgkRqkeyl8iEhuxug1Qa2vrjuepVIpUKhVZLGFYtAhaW2HxYth//2hjOfTQoPCaOBE2bMjw\n17/2Xi8s/1pgIoORTqdJp9NRh1ExtZa/RGpNOXJYMT1dE4BWd2/MTl8DuLtfX2DZFmBTTk9EKevW\nVE/E0qVwxhnw05/C6adHHU2vjo4Mp5wyh23bdPFUqayQerqUv0SkIip1na4OYJyZNZjZcGAqcG9/\ncQxh3arWc/udU09tYfz4Nq66KhOrggvgppvacwougBF0drbR3NweYVQig6b8JVVtoNu66bZv8TLg\n4UV3325mVwILCYq0ee6+1MxmBLN9rpkdAvwO2BvoNrPPAce6++ZC61bs3cRYoSvOz53bwsc+Fq8R\npLCuFyYSBuUvqWaF/l1ZsqT3yMRA8yUC7h6LRxBK9WpqanXY7EGres9jszc1tUYd2k6SEqckX/Y7\nH3nuKcej2vOXRGPlylXe1NTqqdS13tTU6itXrtppfl/5+oILWn3lSvcLLhg4nw+0D+nbYHJYrBrp\nq1lSRpBmzZrOkiUtO/0yqqtr4dJLZ0YcmYhI7ehrlOr++2eyaVMDixfDwoWF/1154IFuzjwT1q8v\nPP+JJ7r5059g990znHWWRsLCpHu/hGS//ZJx+50xYxpYtGgmTU2zmTixhaam2XzrWzOZMaOB9euj\njk5EpDYUuh9vZ2cbxx7bzqWXwnPPwbhxhf9d+ehH6+jqgo9+tPD811+v4+//Ho45pvA+1MNbORrp\nCkF3N/zlL9M54IAWNmzY+azAWbPiN4KUf/FUgK1b4bzz4De/gZEjIwpMRKRGrFpVeJRqwoRuHn44\nmOrqms7kyTsfmcj9d6XQkYves9HhlFO6efTRXffR1RWvIzDVREVXCG64AYYNa+Cxx2bS2lr5K85X\nwle+Al1dMHUq3H03DNNfjojIkPTcki33mohbtzZw443w+OM9o1S5RdEWGhp6j470HJlobi7878rA\n8+t49NFd99HRUcdHPwozZsCkSZDJ7BpnUv7tipsBr9MVlmq9zs1TT8HkydDRAe96V9TRDM22bXDu\nuXDwwRncg3s36gsogxXGdbrCUq35SyqnUM/Wnnu2MGLETD772QbOPjvD1Klz+hilKk++LRTD2LEt\n3HnnTB59tIHvfCe4WPaWLXN2OUqjvq/B5TAVXRX0+utw0klwzTVw6aVRR1MezzyT4aSTdPFUGToV\nXVLLpk1rY/78q8gfZZo6dTb/9V9Be0fPSFjvKFX5f+D2tw93OPvsNh54YNc4m5pmR3Lbujip1A2v\nZZCuvhqOOw6mTYs6kvL55jf7uniqvoAiIsVaurRwz9YLL/T2UxXqry23/vZhBm++WTjOoOdMShWv\nU+eqyK9+FfQ+3Xpr8IdbLZJy6QsRkTh6+WX41Kfgj39Mxhnt9fWF4+zoqKOlBTZu1FXvSxGv/7tV\n4qWX4BOfgB/+MPobWZdbX1/AuCUKEZEo5Rcizz+f4dZb4dhjgzPAH3tsOmPHttCbT3vOPJweXdAF\nzJpVOM5Fi6azahWMGRO0nMyffxXpdHDIdPLkOSq8+qCerjLJPQvl+efrOPvs6cydW309ToUaL3fb\nLejpmjix+t6vVI56uqRaFcqTw4e3cPzxM/n+9xs47rje5Srds1UO/cV53nlt/OIXtdnzpUb6iBT6\ngh15ZAv//d/V2Vye/wU8+ujp/OAHwRWS3/nOqKOTpFDRJdWqryb5Sy6Zzfz51VWITJzYQjrdVvD1\nBx/c9fVqokb6iBS6cvDKldXbXN5X4+XZZ0M6DfvtF35MIiJx0Vfva3BbnurS23Kyc4G5cmUda9bA\nYYdFFFhMqRGnDNRcHlw89fTTYcqU4FIZIiK16NVXIZOpnd7XQj1fY8a0cP750znhBPja1+CNN9Rs\n30MjXWXQV6VfjV+wvpjBjTcGl8eYOhV+9jNdtV5EasvChfDJT8Lpp08HWujqiv9t34aqv6vef/7z\ncNVV8O53Z3j77Tn85S+6sXZRPV1m1gjcSDAyNs/dry+wzLeBswmqj8vd/cns66uAjUA3sM3dx/ex\nj3C3XigAAA/2SURBVMT2RNxwQ4arr57D22/rgqFvvQUXXAB7751h+HBdtV76FlZPl/KXVELuyVMH\nH1zHbrtNZ/HiBubNC+5CkpQm+TBMmtTGgw9WX7N9RXq6zKwOuBmYBKwDOszsHndflrPM2cBYd3+3\nmX0AuBWYkJ3dDaTc/ZVSAkuKP/4RrruugXvumcmPf5zM+yqW0/Dh8G//luH975/Dm2/qV41ES/lL\nKqHQyVN7793C4sUzOf743vseJrmgKKfubrXg9Cjm+Nd4YIW7Z9x9G7AAmJK3zBTgRwDu/hiwr5kd\nkp1nRe4ncbZsgY99DL7xDTjnnOAL9uCDbdx2W0tNFxfXXdeeU3BB71Xr2yOMSmqU8peUXaGTpzZt\nauOb32yPMKr46uv6jlu31tFdY3XX/2/v/IOrKs88/nlC0AZSKAsJAQIXBFHrj7gwsswaMGlNjYjF\nWUZR4mjqTEtnLN0Z8oeMrbkX0Zl1xzrTsrNtRSt1B0bZ7Sywgi6sa0ZLJyUVqYpYiXhjVQztuswS\nLYjw7B/nJrkJ59zcG27ueU/yfGbOcM8577n34cn7fs/z/s5GTKYBf0w7/yB1LVOaD9PSKLBHRNpE\n5NuDNdRF7r0X5s+He+4J2xK3sIkFhkOYfhl5J2gLH9M4f/wG21dWxvn880YWLoTW1pEz0L4QQ52v\nVdWjIlKGJ16HVPXXfgkTiUTP55qaGmpqagpg3uB46inYtw/a2obXNj/5IGhiweTJ1mAwkmlpaaGl\npSVsM3JlWOqXMThOnYJ169K38Bm5k6dyIWiwfSwWY/Nm+OY3Ozh5cgMnTrg9JCUvGqaqGQ+8sQ0v\npJ2vBe7rl+ZnwIq087eByT7fFQfWBPyORoXXX1edNEn14MGwLXGTI0eSOnt2k0KXevvUd+mYMU1a\nXZ3UEyfCts5whVSZH1CDzucw/TLyxf79qldeqbpsmeq+fedq3OzZTXrkSDJsMyPJbbcl0nypPT5t\naEiEbVpGBqNh2bR0tQFzRCQGHAVuB+7ol2YHcC/wrIgsBI6raqeIjAGKVLVLRMYC3wAiuURt90yU\n998/y4EDRTQ3N/LVr7oTgbuEX60mHl/NI4/E+NrXYOdOKCsL20pjhGD6ZeRM+szEKVOKKC9vZMuW\nGD/6kbcsjkjwMglG7hw75t9d6w1VGWZkE5kB9cAfgMPA2tS1VcB30tL8E9AO/B6Yl7o2CzgAvAa8\n0f1swG8MbUh6Hvi13FitJnfOnlW9/37VuXNVX3klqQ0NCa2padaGhoT5cgRCAVq61PTLyBE/vS8p\nadK9e02jhoqGBv+WrvHjE/rcc967w0UGo2FDLnhZG+KwaAVlCNebPl3lgQeSOmqUBbEjnUIFXYU4\nXNYvIzdM7wtPUMPGz3+e1MsuU73uOtXWVi+dS5X1wWiYrRmeBTYbL78cObKJM2f8lpSI9kJ5hmFE\nnzffNL0vNJlWtb/nHnj6aVi2rIOurg18+qnbg+0HwoKuLPjLX2ymSj6xINYwDNf48ENYswba203v\nwyBoMdniYm9Zpt27N/Hss9GvrFsuGoAnn4T33mtk+vS+a4x4+2g1hmdYhAlaKK+01LKjYRhDT/qa\nUCtXruMHP+igqgrmzoV9+85dU8r0Pnw6O/0r696aadHBWroysHEjPPggvPJKjNGjbaZKvli/vpHW\n1nifLTQmTozz8sur+elP4bvftbXPDMMYGvy28CkpibNjx2quvz4G2MxEFwla//Gdd4qoroa1a+Gy\nyzqIx71Zp67u+ZvVhteFwLUNYx9/HB56CF58ES6+OGxrhh9+m8GeOhXjrrtgwgSvhfH06d5p264W\nIGPwFGrD60Lgmn4Zwdx55zo2bx5+my8Pd/yC5dmz4zz//GpefTXGgw920N6+gdOn+94fyjFfg9Ew\nC7rouybLtGlFXHppIxs3xnjxRZgzJxSTRixffOHtZfnYYx2MHr2BY8cKV4CMwmJBl1FoDh+GRYvi\ndHaeu9xaba23d67hLn6V9e73QUPDOrZsKWwwPRgNG/Hdi37Rc3Gx93KfM8de7oWmuBh++EPYu3cT\nL7wQ/UGThmEUlv6V6PXrGzl9OsZDD8GuXTB1ahGdnTZQPooEDbYHUhOxzh3z9dxzZ3nqKbj1Vigt\n9c8fhazIj/igy2+3+C++WMcTTzxKTY293MPi5Emb4WgYRm74VaK3b49TXLyaNWtivPsufPJJI3V1\n8XO6qdavXx2u8cZ5ETTmq6qqiG3bvJmp11/fwW9+s4GPPgpv2YkRH9rb8gVuEjTDsb29iLY272yk\n7EpvGEZ2+FWiu7rWUVe3iQcegPHje9eEamh4lNraOA0Nj9qwhWHA+vX+s043bWpk+3Y4dAjefXdT\nWsAFvT0om3q+Z6jfKyO6pevzz+H4cVuTxUX8ZjjOmhWnoWE1y5dDZWUHyeQGjh6N9kJ5hmHkRlD3\nUDIJe/f6V6L//Oe+lehM3VRGNMm0wCpARQWMH++fP/bsOcuWLd7sx1tv7dtSmvf3Sq5L2A/VwRBu\no+G3dcDu3aqXXKJaU5PUGTNsSxoX6f671db23fLh1CnVhQttq46og20DZOSI33YxZWVNumBBUidO\nVL34YtMFI5igLZ4WLEjozTerFhfnln8Go2Ghi1WPIUMkWn6FdOzYJq2sTOqOHd5GmkEvd8Ndamqa\n+xUM75gzp1kPH+5N59peXUYvFnQZ/RmovAa9NBctSujJk8F7+Fm5N1QHzh+LF/u/V2prm32/bzAa\nNuy7F/36+D/9dB1Llz7KzTd7zcvW1Bw9ggZNqhaxaJE3dmPx4g527sw8aDLsmSyGYXj4DYJvbY3z\nq1+t5v33Y7S0wPbt/t1DxcVnufDCgbuYjJHNQPlj+vQCDDfKJjID6oG3gXeA+wLS/AQ4DBwArs7l\nWT2PmmKmmtF776nOnZtb5DoQL7300qCeKzTD3c5MNZYzZ1Tb2lSvuMK/Vlxfn9DPPsu+Vtydx6qq\n7vKtfbvWmhaVvz0FaulyWb8KzVDnjYHKQqb7fVuxXuopk8XFCa2rU334YdW6Ove6D6NS3szOgcm1\npXQwGpaNYBUB7UAMGJ0SpUv7pbkR2Jn6/DdAa7bPpn1HzoXUz0EVFU16991JnTtXtbxcdebM/BbS\neDw+qOcKzUiwc6Bu4aAuyLFjm/VLX1L98pf988aKFYk+v9Gbx+LnFMJcA7fBvIyyuZ+eJhZbfN7f\nkQ87BgpCCxF0uaxfhfT5QHkjH3YMVBaCxmOtWpXUW25RHTMmvbzGez5fd11z1r8RBiNBawtJ2Hbm\nMtxoqIKuhcDzaedr+9f4gJ8BK9LODwGTs3k27V5WhXTGjCZ95pmkbt2qOn++/0uzqiqh+/ernjmT\n/0IadobIFrMzePxHQ0NCT51SveYa/6BMpFnHjfMmWpSXp39HvOc7qqsTumtXcM37jjuCArfsXka5\nBnb5CA7zb0fm8lagoMsp/QrD5wPljcHaMWtWk/7ud0nt6FB96y3V+nr/snDRRQldulR14kT/+1dd\nldCtW1WXLPEva/0ryK6NwTWtzS9RsVN16IKu5cDjaed3Aj/pl+Y/gL9NO98DzMvm2bR7PYVswoSE\nXn55cEvEpEkJXb5cddq07LoO81lIo5IhzM6BXyRBQdnKlQn95BPVgwdVq6r8a99lZc16ww2q48b5\n50Fo1lGjVEtKVEeP9v+dceMSesUVql/5iv/9srKELlrUP/DrvV9R4XW71NWpVlT4v7CmTEnojTeq\nTpni/x1Tpyb0pptUb7pJderUzGky3V+6VHXp0uA0fi3LBQq6Cq5f3f44H39l4/MlS1SXLAlOM2VK\nQm+4of/fvjdvlJcndPFi1epq1UmT/L+jtDShM2eqlpT437/ggoRWVnoVlNJS/7JwySXNum2b6rx5\nmfU6m+DQRUxr80tU7FQdnIYNlWj91+BEyw477BhJR66ClbPAmX7ZYYcdQ3jkqknZzF78EJiRdl6Z\nutY/zXSfNBdk8Sx4lg+LjW8Nw3AK0y/DMJwhm3mQbcAcEYmJyAXA7cCOfml2AHcBiMhC4Liqdmb5\nrGEYxlBh+mUYhjMM2NKlqmdE5HvAbrwg7UlVPSQiq7zb+riq7hKRJSLSjrfIxbcyPTtk/xvDMIw0\nTL8Mw3AJSY1HMAzDMAzDMIaQ0Hd1FpF6EXlbRN4RkfvCticIEUmKyO9F5DUR2Re2PemIyJMi0iki\nr6ddmyAiu0XkDyLynyIyPkwbUzb52RkXkQ9EZH/qqA/TxpRNlSLy3yJyUETeEJHvp6475VMfO1en\nrrvo0wtF5Lep8vOGiMRT153yaa5ERb/AXQ2Lin5BNDQsKvoVYKuTGpZP/Qq1pUtEivBWev468BHe\nGIrbVfXt0IwKQESOAPNV9X/DtqU/IlINdAFPq+pVqWuPAP+jqv+YehlMUNW1DtoZB06o6mNh2paO\niFQAFap6QERKgVeBZXjdTs74NIOdK3DMpwAiMkZVPxORUcBe4Pt4MwSd8WkuREm/wF0Ni4p+pexy\nXsOiol8D2OqchuVLv8Ju6VoAHFbVDlU9DTyD53AXEcL3ly+q+mugv5AuA36Z+vxL4JaCGuVDgJ3g\n+dYZVPVjVT2Q+tyFt1hmJY75NMDOaanbTvkUQFU/S328EG88qeKYT3MkSvoFjmpYVPQLoqFhUdEv\niJaG5Uu/wi6A04A/pp1/QK/DXUOBPSLSJiLfDtuYLChPzcBCVT8GykO2JxPfE5EDIvKEC03e6YjI\nTOBqoBWY7KpP0+z8beqScz4VkSIReQ34GNijqm047NMsiJJ+QbQ0LEr6BQ6WN4iOfoH7GpYv/Qo7\n6IoS16rqPGAJcG+qmTlKuDpj4p+Bi1T1arzM7FJzcinwb8Dfp2ph/X3ohE997HTSp6p6VlX/Gq/W\nvUBELsdRnw5ToqxhLucLJ8tbVPQLoqFh+dKvsIOubBYudAJVPZr690/Av+N1LbhMp4hMhp5+82Mh\n2+OLqv5JewcWbgSuCdOebkSkGE8E/kVVt6cuO+dTPztd9Wk3qvp/QAtQj4M+zYHI6BdETsMiky9c\nLG9R0S+Inoadr36FHXRFYvFBERmTisQRkbHAN4A3w7XqHIS+feA7gMbU57uB7f0fCIk+dqYyajd/\nhzt+/QXwlqr+OO2aiz49x04XfSoik7q7CESkBKjDG7/hok+zJRL6BZHQsKjoF0RDw6KiXxABDcun\nfoW+TldqKuiP6V188B9CNcgHEZmFVzNUvAF0m12yU0S2ADXARKATiAPbgH/F296kA7hNVY+HZSME\n2lmL149/FkgCq7r7yMNCRK4FXgbeoHePrfuBfcBWHPFpBjtX4p5Pr8QbaFqUOp5V1YdF5K9wyKe5\nEgX9Arc1LCr6BdHQsKjoF0RHw/KpX6EHXYZhGIZhGCOBsLsXDcMwDMMwRgQWdBmGYRiGYRQAC7oM\nwzAMwzAKgAVdhmEYhmEYBcCCLsMwDMMwjAJgQZdhGIZhGEYBsKDLMAzDMAyjAPw/EKxUofjfFaIA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108c6af98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xrange = range(30)\n",
    "fig, ax = plt.subplots(2,2,figsize=(10,7))\n",
    "axes = [ax[0,0],ax[0,1],ax[1,0],ax[1,1]]\n",
    "ls = [1,4,10,20]\n",
    "for i in range(4):\n",
    "    axes[i].set_ylim(top=0.4)\n",
    "    chart = axes[i].plot(xrange,[poisson(ls[i],x) for x in xrange],'o-')\n",
    "    title = axes[i].set_title('Lambda='+str(ls[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ポアソン回帰**  \n",
    "さて、ではポアソン回帰とはどんな回帰か。これは、モデル化したい目的変数がポアソン分布に従うだろう、と考えられるデータに対して用いられる回帰分析である。「モデル化したい目的変数がポアソン分布に従う」というのはどういうことだろうか。  \n",
    "\n",
    "統計初学者であれば、\"何らかの説明変数を横軸、目的変数を縦軸にとったグラフを思い浮かべて、その形が上のポアソン分布のどれかに近い形をしている\"というイメージをしがちだと思う。が、これは間違っている。\n",
    "  \n",
    "ポアソン回帰を用いる時目的変数に想定しているのは、例えばデータセットの中から特定の説明変数,目的変数のペア(y<sub>i</sub>,x<sub>i</sub>)を取り出したとき、y<sub>i</sub>の値が、対応する説明変数の値x<sub>i</sub>から求められるλ<sub>i</sub>をパラメータとしたポアソン分布に従って生成されている、ということである。これはすなわち、説明変数の値が決まったとき、目的変数は(一意に決まるのではなく)ポアソン分布を誤差構造とした確率的な値として決まる、ということ。\n",
    "\n",
    "例えば一説明変数の場合であれば、ものすごくざっくり言うと、説明変数xの各値でデータを切り分けていって、xの各値ごとに目的変数yの値のヒストグラムを作ってみる。そのどれもが、何となくλの異なるポアソン分布の形をしていそうな感じであれば、その説明変数xを使ってyをポアソン回帰できるかもしれない、と考えられる。\n",
    "\n",
    "**ポアソン回帰できそうなデータを人工的に作ってみる**  \n",
    "ここでは簡単に、説明変数はx一つだけで、λ=x+1とする"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+QldWZ5z+nf9y+TZOG7oiQ8DNRjBB1kcxYLkLooBJ0\nJ5KZnZKQ7I7MOnEzMxBnnNpNnJRAWVOZmK1xyzGLNTGukMFGXavUWJVENKE1aEwTDZFJAzoxoMSA\nJJhEke4GfPaPcy/vvfDevud97u33dt9+PlVv3e6377fPuadfzns43/d5HiciGIZhGPVDQ607YBiG\nYVQXm9gNwzDqDJvYDcMw6gyb2A3DMOoMm9gNwzDqDJvYDcMw6ozEE7tzrsE594Jz7lu57zucc1ud\nc3udc4875yZUv5uGYRhGKJoV+41AX8H3XwSeFJEPAd8Hbq5GxwzDMAwdiSZ259w04GrgGwWnlwOb\ncl9vAj5Zna4ZhmEYGpKu2P838D+AwnDVySJyCEBEDgJnV6lvhmEYhoLgid0595+AQyKyE3BDvNVy\nFBiGYdSQpgTvvQy4xjl3NdAKvMc596/AQefcZBE55JybArwRJ3bO2YRvGIahQESGWkyfQfCKXUT+\nXkRmiMgHgU8B3xeR/wo8BqzKve064NEhfocdIqxbt67mfRgph42FjYWNxdCHhmo8x/4V4Ern3F7g\n8tz3hmEYRo1IshVzChF5Cngq9/UR4IpqdsowDMPQY5GnNaCrq6vWXRgx2FhE2FhE2FhUhtPu4SRu\nyDlJqy3DMIx6wTmHDJd5ahiGYYwObGI3DMOoM2xiNwzDqDNsYjcMw6gzbGI3DMOoM2xiNwzDqDNs\nYjcMw6gzbGI3DMOoM2xiNwzDqDNsYjcMw6gzbGI3DMOoM2xiNwzDqDNsYjcMw6gzbGI3DMOoM5IU\ns25xzv3IOfcT59wu59y63Pl1zrkDzrkXcsey4euuYRiGUY5E+didc+NE5B3nXCPwDPB54CrgLRG5\nvYzW8rEbhmEkZNjzsYvIO7kvW/Bl9fIzdaJGDcMwjOEj0cTunGtwzv0EOAg8ISI7cj9a7Zzb6Zz7\nhnNuQtV7aRiGYQSTqJi1iLwLXOycawceds7NBTYAt4qIOOf+AbgduD5Ov379+lNfd3V1WV1DwzCM\n0+jp6aGnp6ei36GueeqcuwU4Wri37pybCTwmIhfFvN/22A3DMBIyrHvszrmz8tsszrlW4Epgj3Nu\nSsHb/gT4tyQdMAzDqHcOHz7Mjh07OHz4cCrtJdljfx+wzTm3E/gR8LiIfBv4qnPuxdz5xcDfDkM/\nDcMwRiVbtjzAzJnnc+WVn2PmzPPZsuWBYW9TvRWTuCHbijEMY4xx+PBhZs48n2PHtgEXAS/S2vox\n9u/fw6RJk4J+x7A/7mgYhmGEs2/fPjKZWfhJHeAimptnsm/fvmFt1yZ2wzCMYWLWrFkMDu4DXsyd\neZHjx/cza9asYW3XJnbDMIxhYtKkSdxzzwZaWz9Ge/t8Wls/xj33bAjehtFie+yGYRjDzOHDh9m3\nbx+zZs1KPKlr9thtYjcMwxjBmHlqGIZh2MRuGIZRb9jEbhjGmCXtiNC0sIndMIwxSS0iQtPCzFPD\nMMYc1YgITQszTw3DMAKoVURoWtjEbhjGmKNWEaFpYRO7YRhjjlpFhKaF7bEbhjFmqSQiNC0s8tQw\nDKPOMPPUMAzDSFQar8U59yPn3E+cc7ucc+ty5zucc1udc3udc4/ny+cZhmEYtSF4YheRAeBjInIx\nMA+4yjl3CfBF4EkR+RDwfeDmYempYRhGldFEno6GaNVEWzEi8k7uyxagCRBgObApd34T8Mmq9c4w\nDGOY0ESejpZo1UTmqXOuAXgeOAf4PyJys3PuTRHpKHjPERHpjNGaeWoYxohAE3laq2hVjXnalOTN\nIvIucLFzrh142Dn3YfyqvehtpfTr168/9XVXVxddXV1JmjcMw6gK+cjTY8fOjDwtNUlrNBp6enro\n6emp6HeoH3d0zt0CvAP8BdAlIoecc1OAbSIyJ+b9tmI3DGNEUO8r9iRPxZyVf+LFOdcKXAnsBr4F\nrMq97Trg0SQdMAzDSBtN5Glek80upq3tQ2Szi4OjVdM2XJNsxbwP2JTbZ28AHhCRbzvnngMedM79\nN2A/cO0w9NMwDKOqrFy5giuuWJI48tRPga251/Js2fIA11//V2QyPj/NPfdsYOXKFRX0PKCPFnlq\nGIZRnlpt31jkqWEYxjChSfVbq/TANrEbhmEEoEn1W6v0wDaxG4YxZkliamrM01qlB070HLthGEa9\noDU1k5qnWpO2Esw8NQxjzGHPsRuGYdQZo8kI1WATu2EYY47RZIRqsIndMIwRhSZKc/fu3WzatInd\nu3cHvb/Q1Gxr+w+JIk+TaCr5TJVgE7thGCMGTVrcNWv+hrlzP8KqVV9m7tyPsGbNjcHt+byGA7nX\n4dHUJNWviKRy+KYMwzDieeONN6S1tVPgpwIi8FNpbe2UN954o6Smr69PoLVIA63S19dX9bbS0pxO\nbu5MNN/ait0wjBGBxpzs7e0FphdpYFrufHXbGk2Gq03shmGMCDTm5CWXXAK8VqSBA7nz1W1rVBmu\nSZf42gPbijGMMcUbb7whvb29ibYdurvvl9bWTmlvv1haWzulu/v+sprVqz+f246ZLdAqq1d/ftja\n6u6+X1paJkpLy1xpaZkYrGluHi+Nje+X5ubxQZpCUGzFWICSYRhVp5JUtYcPH04cpbl79256e3u5\n5JJLmDPnjDo/VWtr6dKreeKJbcAk4DBLl36Mxx//9pCa6dM/yIEDB4FpwAGmT5/Mq6/+IriPmgAl\nW7EbhlFVqmEYjkS2b98ea9Ru3769pGbz5s2xms2bNwe3i5mnhmHUmtEUoZmErVu34lfdhUbt1Nz5\neB588MFYjT8/fCQpjTfNOfd959zPnHO7nHNrcufXOecOOOdeyB3Lhq+7hmGMdEZThGYSli5dChyg\n2Kj9Ze58PNdee22sxp8fPoL32HOFqqeIyE7n3HjgeWA5sAJ4S0RuL6OX0LYMwxjd5PfYm5tncvz4\n/mHfY9doAJ555hm2bt3K0qVLueyyy8q+/+Mfv5qtW78PTAB+x9KlS8rusc+Y8QFee+0Q8H7g9ZG9\nxw48AlwOrAP+LuD9wXtKhmGMfip5KmbChPmJnlRJqhERufLKq4qeplm69KqymmnTPlCkmT59VlnN\nBRdcLJAVmC6QlQsvnBfUvzwo9ti1k/osYB8wPjex/wLYCXwDmFBCk+jDGIYxtkgzsjMtI/Sxxx6L\n1Tz22GPlBySHZmJPXGgjtw3zEHCjiLztnNsA3Coi4pz7B+B24Po47fr160993dXVRVdXV9LmDcOo\nU/Km67FjZ5qupbZXNBoY2ggttSUzlBH6mc98JlbzyCOPxGoeeeQR/uiP/ihW09PTQ09PT8m+B5Hk\nLoCvuPRd/KQe9/OZwIslfhZ8hzIMY+xhK/Z4GO6tGOCbwO2nnZtS8PXfAt0ltMEfxDCM0Y92j72l\npV2y2RnS0tIevMeeyYyXTGayZDLhkZ1Ll14l0CJwtkBL0B779OmzchP1ucF77BdeOC+nOUegNZU9\n9iSPO14GfAZY4pz7ScGjjV91zr3onNsJLM5N7oZhjGG0qWqfffaHDAwcp78/y8DAcZ599tmymnvv\n3cTg4EkGB9sZHDzJxo2bgtry680G/BMuYVPh+efPAd4Ffg+8GxTl+vOf7wMEGASEV155Jaitikh6\nJ9Ae2IrdMMYE2u0RTQpezZaKVqfR3HbbbbGa2267bcj+FYJFnhqGUWu0kaeaFLyaaFCtTqPZsmVL\nrMafHz5sYjcMo6poI081KXg10aBanUazcuXKWI0/P4wkXeJrD2wrxjDGDN3d90s2O1Gy2XMlmw1L\nbyuiS8HrTdCswDSBbJAJGukiIzRE5zWNOdO1MUjT1tZe1E5bW3tQ//JgWzGGYYwEnn32h/T3D9Df\n7+jvHwgyQQEWLFhAS0sz2ewALS3NLFiwoKzGrxsdkM29hrFjxwt4U7MfEH784x+X1Tz55DYgA8wA\nMnzve98L6F8DhYZrKiS9E2gPbMVuGGOCNOuQas3TDRs2xOo2bNhQUnPTTTfFam666aaSGjNPDcOo\nC9KsQ6o1TzWm5kMPPRSr8eer1041sIndMIyqkmYdUq15qjE1//RP/zRW489Xr52qkHSJrz2wrRhj\nlKGJnEwTbf/6+vpk48aNZbdGKqGSOqTZ7ERpazsv2HSNTNBzgk1QEZHOzklFus7OSWU1jY2ZIiO0\nsTFTVuPN0yjC1cxTw6gR2sjJtND2b82av2Hu3I+watWXmTv3I6xZc+Ow9O/OO++gr+95Nm78En19\nz3PnnXcEa51rAFpzr+V5/fWDFEZ2/upXvwrS9fcfL9INDAyU1cyZ8+GcZgAQ5s6dW1Zz7NgghRGu\n/f39Qf2riKR3Au2BrdiNUcJIr9mZZmRnmmg+lzbJlsbU1LSlMVxPB1uxG0bljPSanWlGdqaJ5nMN\nlRZ3KDSmpqYtjeFaDWxiN4zTGOk1O9OM7EwTzef65Cc/SZw56c+XRmNqatrSGK5VIekSX3tgWzHG\nKCJfbq29/eJE5dbSQts/b2pGUZqhpqbGcNVourvvl+bmdmlsPE+am8PS9vq0uJE5GZoW15ua0ViE\nmJqatrzhmhHoEMgEGa6FkFZpPM1hE7sx2qjHp2JWr74x8dMqkea8YdWIiHR0TJbC+qAhT6qMGzex\n6DOFPnUCTUVPuEBDWU1zc1tRW5lMa4J2Zge3U6y3id0wjBJozNO0NCK6aFBtZOeKFStidStWrCip\nueWWW2I1t9xyS1XbOZ1hndjxDsD3gZ8Bu4DP5853AFuBvcDjWDFrwxiRbNy4MbeCloJjtmzcuLHm\nGhGRRYsW5Va1hbpzZdGiRSU18+bNi9XMmzf0Fkl7e3usrr299Gr/vPPOi9Wcd955VW3ndDQTexLz\n9ARwk4h8GPiPwF87584Hvgg8KSIfyk38N+t3/A3DGC405mlaGtAZmtrIzquuuipW58/Hs2LFiliN\nP1+9dqpC0jtB/gAeAa4A9gCTc+emAHtKvD/4DmUYYw3Nfrlujz0fEXpugv1yraZZ4D0CzcF77D4a\nNEqLG7LHrk2LCw2J99gzmdYiTdgee/J2ivXDu2I/hXNuFjAPeC43qR/KzdwHgbMrutMYxhhDE0Wq\njTzdu/dlCiMnX3rp5QCVK9KE0NPzA6AJv9Zr4qmnng7SvfXWOxSmxX377bfLajKZ1qL+tbS0BLW1\nevUafBrd3wHvsnr16rKaxYu7KEzB29XVVVYzblx7Uf/a2sYH9a8ikt4JgPHAj4Hlue+PnPbz35TQ\nJbpLGcZYQBNtqY081aS41Rih2mhQjTmpMVy1n6tua54655qAh4B/FZFHc6cPOecm534+BXijlH79\n+vWnjp6eniRNG0Zdoom21EaealLcaqJVtdGgDzzwQKzOn49HmxY3rfqqmv719PQUzZUqktwFgG8C\nt5927jbgC7mvvwB8pYQ2+A5lGGMFW7FH2Io9Hob5ccfLgJPATuAnwAvAMqATeBL/uONWYGIJffAH\nMYyxhCaKtLv7fslkxksmM1kymfHBkaeaOp+aaFUfoRlpQqNBvTkZ6ULMSW+4RpoQw1Uk/7miKNJk\n9VWnSmh91VrUPG1KsLJ/Bmgs8eMrEvwnwTCMAlauXMEVVyxh3759zJo1i0mTJpXV3HvvJgYHTwLt\nwO/ZuHETK1eWfuwuj19jnQR+k3stjzdcHdAKuCDDdc+evMbXId27d29QW/4zNZ7SHT9+vKzG1xRN\nXvP0kUceI0qn+xaPPvqtsumFt217OtfGOOBI0Jay/0xC3jwN+UwVk/ROoD2wFbthVAVtnU+NTqPR\nbKmIiHzuc5+L1X3uc58rqdFuxWzevDlWt3nz5qp+Lu1YFIKl7TWM+kdb51Oj02g0JijAww8/HKvz\n5+PRmqcPPvhgrM6fj0fzubRjUSk2sRvGKENb51Oj02g0EZoAf/zHfxyr8+fj0UaeXnvttbE6fz4e\nzefSjkXFJF3iaw9sK8YYA2gzQm7fvl3Wrl1bdjslT2TiTQs28SJdMvNUo9FEaIqIOFeccdG5prKa\nqHap14Sap9OnzyrSTZ8+q6xGY+56TWTSho5FHmwrxjBqhzYadOnSq1m48EpuvXULCxdeycc/fnVZ\njV8rJTcMfX3QKHIypD5oX98eCs2/3bt3l9WcOCE5TT8gnDwZZtS2to4v6t+4cePKao4cOUKhIXzk\nyK+D2jp48HBRW4cOHSqr0Ri1XhPVPE2FpHcC7YGt2I06Js1ny7Xmqeb5co3JqK3zqXnme/ny5bGa\n5cuXD9lWWkZorcxTm9gNowr09vbKhAnzpTA9a3v7xdLb2zukbu3atRKX1nXt2rVV1YiIXH/99bG6\n66+/vqTmmmuuidVcc801JTUzZsyI1cyYMWPI/mlS8La1tcVq2trahmxLk4I3Lc3paCZ224oxjCqg\nrUOalqEJupqdGpNRW+dTY4ReccUVsRp/vjRpGaFmnhrGCCOpEdrdfb9ksx3S2nqBZLMdCaNBkxmh\nGkNTJB8RGqXTDYkI9SZjpAkxGX2dz+gzhdb5jOqQ+sjOkChNcEVtgQtqy5uaUS3ScCO0VeAcCTWF\nNYZrIdiK3TCqg8YIffbZH9Lff4xjx35Pf/8xnn322aC2/LonmSG3Y8cLFBp/P/7xj4N0b775OwrT\n6f72t78tq3nttdeLNAcOHCir0UaDHj16lCiy03H06FtlNc3N4yiMjM1kskFt+T42AmdROqi+mMHB\nE3hTeBAQBgcHy2q8kRz1L9RIroikdwLtga3YjVGCxgjV1vnUGKFpRltqokG15qnGCNWakxrddddd\nF6u57rrrqj4WhWDmqWFUjsYI1db51BihmtqgIjojdPLkybGayZMnl9RozVONEao1JzW6zs7OWE1n\nZ2fVx6IQzcRuWzGGcRoaI1Rb51NjhKYZbamJBtWapxojVGtOanSf+MQnYjX+fDzasaiYpHcC7YGt\n2I1RhDdCJ0pb23mSzU4MMkJ9GtjIjAut86kxQn20ZYN4U7MhYbRlZOSFGKE+GjTShESDevM0GotQ\n81RjhGrNSY3O12Mt7F9jWY0fi+jvGzoWUZu2YjeMquFcA9Caey2PT28bmXFh9UTz5mlkhIbga4O2\nALOAlqDaoAATJnRQaORNnDgxsH+RESpSvo/eIIzG4uTJsFS1kRHaQqgRGpmTvn+h5uR5551fpPvQ\nhz5UVtPQ0FykaWwsb7p6kzZZquSKSXon0B7Yit0YJWjM0zRT6WoNQ03k6YoVK2I1K1asKKlZsmRJ\nrGbJkiVD9k/zubTmpGYsNG2NePMUuAc4BLxYcG4dfgPphdyxbAh98AcxjFqiMU+10aAandYw1ESe\ntre3x2ra20s/X57JZGI1mczQWxCaz6U1JzVjoWlrNJin9wIfjzl/u4jMzx3fVf7HwTBGDBrzNM1U\nulrDUBN5etVVV8Vq/Pl4Fi5cGKvx50uj+Vxac1IzFpq2RoV5CszkzBX73wVqg+9QhlFr8uZpNntu\nsHmqqYcZ6aLIzvC0uI3i08E2BhuGPvI0MvJCIk+9SVtoGDYEaFxRO8miQZNFdmrNyWgszgkeC01b\no9U8Xe2c2+mc+4ZzbkIVfp9h1BwfRTpAf7+jv38gKIrUR4NGUZOh0aBPPrmNwsjO733ve2U1PgIy\nA8wAMkERkAC7dv0bFKTg3bXrxTIKaGjIUGwYli+VPG7chKJ22treE9S/qD6oj+wMqQ/qDc3IfA4x\nNAF27fpZUVu7du0qq/FGaDKj+33vmwqcAA4CJ3j/+98fpKuIJHcBzlyxT8LfigH+AbhnCG2iu5Rh\n1ApNFKk2GlRjrmkiIEV0Rqimf5r0uyK6KFetkZxWFKkm2vd0GO7I09Mn9tCf5X4u69atO3Vs27Yt\n+IMZRppooki10aAac00TASmiM0I1/dOk3xXRRblqjeS0okg10b7btm0rmivTmNhnAbsKvp9S8PXf\nAt1DaIcYZsMYOdiKvbL+2Yo9YsSv2IFu4HX8ptmrwJ8D38TbvDuBR4DJQ+iDP4gxNujr65ONGzeW\nTZRViLamaFKdjyKNDMOQKFJt7U1NlKY3TgvNyfIRkF6X3AjVpOD16Xej/oWk3xXJR7lG9UFDoly1\nNUXTiiLVRPsW9zP5xB5snorIp0Xk/SLSIiIzROReEfkzEblIROaJyCdFpHzRQMMA1qz5G+bO/Qir\nVn2ZuXM/wpo1N5bVaGuKanQ9PT+gMELzqaeeLqvx0aCRsRYaDerNvyhKM8T8a27OQkFN0UwmE9hW\nciPUR3IWRnaWNzR9+t2ofyHpd4FczdOoPmhIzdOmppYiTXNzc1BbHR1nUfi5Ojs7y2q8ERpFkYYY\nocuXL6fwWlq+/Jqg/lVE0juB9sBW7EYOzVaHtqaoRqeJSkwzfay2Lc1WwrJly2I1y5YtK6lZuHBh\nrGbhwoVD9k+zhaPd9tFsnWm2VbTpnAthuPfYKzlsYjfyaMxJbU1RjU4TlZhm+lhtWxrzL5vNxmqy\n2WxJTVNTU6ymqWnobRWN6ao1ajVmt8YI1aZzLkQzsVsSMCN1NClutTVFNTpNVGKa6WO1bWmiILu6\numI1/nw8l156aazGny+NJh2xNoWxRqdJe6xN51wxSe8E2gNbsRsFeHOyNbcCag0yJ7u775fW1k5p\nb79YWls7g2uKamqRaiI0NTU0I13UVrLam+EaEZ0Rqoki1UaeakzXqE6q/0yhRq03uyPTNcTs1hih\n0bV+bvC1Xgi2YjdGC3feeQd9fc+zceOX6Ot7njvvvKOsZuXKFezfv4cnn/wX9u/fw8qVYZXeNbVI\nfW3QyAgNqQ169tlTKDRBJ0+eHNQ/n3Y2aisk7azXnMDn5TsRnKr25EkhaYrbKIrUG6EhUaTd3VuA\nY8C/A8dy35fHm67RWISYrpmMNyXzn6mlpSWorSNHfkeh6frmm2+W1URpj31bIWmPFyxYQDbbQjbr\nyGZbWLBgQVD/KiLpnUB7YCt2owZozCuNSaZ9XjnNVLCaZ7c15qTW6NbUPNXGD2iemdeY6tqxKAQz\nTw2jGI15pTHJNBqRdFPBaqItNeak1ujW1DzVRvxqolw1prp2LArRTOy2FWPUjMOHD7Njxw4OHz48\nbJrIvNqET0a6iXLmVWSS3Qn8Re51aJMs0jySa+ORshqoNBVscVvlUsFGNTuLdUPV7IxMxmLNUCZj\nZFh/Cfgo8KUgozuqeVrc1lA1T7XmqaaWq8ZU15r+FZP0TqA9sBW7UUDeCJ0wYX6wEarRiIg0Nmal\n0KhNFtk5W8INw4bEGq/TmJOVtjW8n0sz5pX3L5lRq9F58zTShJinWtM/6qdtxRijAM2+o3avUrNH\nrNmL1uzZihTul28T6M29Dk/+EU2umDQDqDR72NHft3j8ygUoRSUJHxbYmHsduiRhdA3+s8D1Av8c\nvF+uTYUhYhO7MUrQ7Dtq9yo1e8SavWjNnq1IehkDRXTZHdMMoNLsYWsDlDQlCauxX65BM7HbHruR\nOpp9R+1epWYPNtqLLtYMtRet2bMF3R67JlAGdGXu0gyg0uxha/fYNSUJa7ZfriHpnUB7YCv2ukXz\n38zu7vulpaVdMplp0tLSHrzH3tIyQVpaPigtLROC9yo1ASyazH8+M2G0/xqSmVBEl90x2uv1Zd1C\nMwZqsjumGUClCQzTBij5koRRWyElCbu775fm5jZpauqU5ua24Gtw+/btsnbt2iG3ekqBrdiNtNFm\nXLz33k0MDBxncLCVgYHjbNy4KVAzyMBAIwMDg0Ea0AWwXHDBRTlNC+C48MILy2p8lsEouCY0y2BL\nyzgKA5uy2WxZzS9/eRAKyrq9/vrrQW35dqKx8MfQaAKofAm5KAtiKIsXL6YwGGrx4o+W1XR0vJek\nWRp9HyFpmbsvf/k2jh9/lxMn3svx4+/yj//4lbKapUuvZuHCK7n11i0sXHglH//41UFtVUTSO4H2\nwFbsdYfW0IyMq2KTbKjVjEYjogtgSTO7o8bcTTNASdOWdizSCiYT0V1PmutCe90WgpmnRppozSSN\ncaXRiOgCWNLM7qgx/9IMUNK0pR2LtILJRHTXk+a60F63hWgm9uCtGOfcPc65Q865FwvOdTjntjrn\n9jrnHnfOTaj2/yiMkYvWTNIYVxoN6My1NLM7avqnMVxBZwpr2tKOhSYTotZI1lxPmutCe91WTOgd\nAFgIzKOgYDVwG/A/c19/AfjKEPrgO5QxevCZEydKNnuuZLMTg80kjXEVac4J1ojks/hF5lpIFj9v\n4mUFpgpkE2R3bBBoFmgINgy9+dci4LMNhph/mhJtIiKaknqatvxYRGMeOhaaTIiaoCER3TWoMXc1\n7RTCcK7YRWQ7cHr6s+X42F9yr6VvXUZd4jMnDtDf7+jvHwjKnAiwatV1tLQ0k8kM0NLSzKpV15XV\nbN36BIWG4datjwe1deTIEQrNtSNHfl1Ws2vXrpxmHODYteunZTWDg4N4s3UW0MLgYH9Q/3xGwwZg\nItAQlNFQa05GJfUGILCk3pIll1NoMl5++eVlNTfccAOF5vMNN3w2qH97975c1L+XXnq5rOb88+dQ\nOBZz5swJastfgxlaWt6lpSUTdA16czcaixBz9/HHv8327U+wdu2n2b79CR5//NtB/auIJHcBYCbF\nK/Yjp/38yBDaRHcpY+SjLfulMV01UZMiujJtmnJw8+fPj9XMnz9/yP5p2tKapxpTU2P+aa+LNE11\nzTVYjTJ3GhgBjzvKUD9cv379qaOnp6fKTRtp09vbC0wHLsqduQiYljtfmn379pHJzCrSNTfPZN++\nfSU13/nOd4Bpp7U1NXe+NM8991yszp+Px1+bZ2qGumZ37twZq/HnS6Np66GHHorV+POleeCBB2J1\n/nw8W7dujdX48/ForwtNWxoN6K5B7edKSk9PT9FcqSLJXYAzV+y7gcm5r6cAu4fQDutdzUgfW7FH\n2Io9wlbs1YXhftwRv3m4q+D724Av5L4283SEoYkITarxZldUXiy07Fd39/3S1DROGhomSFPTuCDT\n9cxsfOVfnjHYAAATSUlEQVSjJr0undJu2nJwGp0mWlVEFxHqzb/G3N+5Mcj805aDq8xUT2ZO5qOf\ns9kZwdHPmpKOlTKsEzvQDbyOdzVeBf4c6ACeBPYCW4GJQ+iHfQCMiLTS4l555VVFF3roP6qOjslF\nupAnVZqb20TzpIXXudzk5IJ00CRJQ+/90zDJ0hB4XX7SPDv3Wn5i958pGr9kYxG1FaJraEiegtdf\nF1mB6QLZxJNtJjN92FNNrF59Y+5znZdoku7r65ONGzcO+0o9z7Cv2Cs5bGJPj7TS4qYZDaqNZtTo\nNNs+mqhOEV05uDTHQrPtk+b2iDb6uVbbKhpsYjdEJL20uGlGg2qjGTU6TXpbTVSniK4cXJpjoYk8\n1V4XaaZz1kS51grNxG5JwOqQtNLiphkNqo1m1Og06W01UZ1QWA6uWDdUObg0x0ITeaq9LtJM56yJ\nch1VJL0TaA9sxZ4q+YjQtrbzgiNCu7vvl0xmvDQ1nS2ZzPggjTeuogjN4YwG9cZfsv3hSNcs8B6B\n5sA99uRGrSaq0+ucFO/Nl99j16bF1eg0kafRdTEt0XWhiWTWXOsieSM08gGG0whNu4KSrdjrGOca\ngNbca3m+/OXbGBw8yYkTExgcPBmUknT79h9SGKH5zDPPBLV15MhvKY4GPVJW4yM7G4AJQENwZKfX\nNeGfyG0K1DWQNKWrj+qMIiBDojpzSgrHwrnGsgqfSlfIp7cNSaULsHhxV5Guq6urrMZHng4A+4GB\noMhTv5ZzQD5lchjaSOak1zrAggULyGZbaW2dQDbbyoIFC4K1SdCmtq6IpHcC7YGt2FNDYyhVVm8y\nPOWsiM5o1JiMWl1atUFFdLVStc+xp/WcuNY81RiaWvNUq0tKNdrBzFNDRGcopVlvUmM0akxGrS6t\n2qAiulqp2rS9aaVL1pqnGkNTa56mVb+0Gu1oJnbbiqlDNIZSmvUmNUajxmTU6tKqDQq6WqnatL1p\npUvWmqcaQ1NrnqZVv7RmdVKT3gm0B7ZiTxVvKHVIa+sFks12BBlKUapab3gNZ73JyGg8R0KNxsoi\nOyODNyyKtEGKg4ZCa4MmN3c1tVK1aXv1kZ1RlOtwRoNqIjvzgXXt7RcHB9YV6traLkqkS4q2f3mw\nFbuRx5tQxzh27Pf09x8LMqF+/vN9FJp4r7zySlnN0aPHijQ+BW0oQj4Fb4hBOW7cBOAEcBA4QVvb\newLbaaDQ4A3DUWjUMnR+OyBf8zTShNY8nTp1OoVG7bRp08pq5sz5MH4sDgEnmDt3blBbq1ZdRybT\nSCbzFplMY1CqWr8ui2qyhraTNCUuwJ133kFf3/Ns3Pgl+vqe58477yirWblyBfv37+HJJ/+F/fv3\nsHLl0P9LKkTkXWAg9zo8VNI/NUnvBNoDW7GnhsaE0hih2mjLJUuWxOqWLFlS1f5p29IYrtr+aWp2\naoxukfSii9MyJithNPQxD2aeGiI6E0pjhGqjLTOZTKwukym9naA1ajVtaQxXbf80NTs1RrdIetHF\naRmTlTAa+phHM7HbVkwdojGhNEaoNtpy4cKFsTp/vnr907alMVy1/dPU7NQY3ZBedHHNDMMEjIY+\nVkTSO4H2wFbsarZv3y5r164t+xxwIZp0uhojtLJoy2QZDX3/8kZjUwKjVpuCN1k0qO9f1E5o/3zN\nzqitkJqdGqNbRBelqTFCKzUM02A09FFEt2JvqvF9xSjD0qVX88QTPcA0br31f7F0aVdQzcSenh8Q\nGXlv8dRTT5fVHD3ajzfHfLTgO++8E9hLIR/JGBql6fsW9S8Eb8xm8BVzDgTVBo3aEvJ1NENobh7H\n8ePvkjeFM5lsWc3g4Mmido4fPx7U1oQJHbz22qFTbU2cOLGsZvHixezatfeUJqT2Zp6kUZqrVl3H\nU089g8gAzoXVp125cgVXXLGEffv2MWvWLCZNmhTcv7QYDX1Uk/ROoD2wFXtitBF8GnNNEwGpqf6j\n1c2ZMydWM2fOnCHb0lRQ0kSRaiNPNX+rNCtXjSaTsV7BzNP6QhvBpzHXNBGQ2Ww2VpPNZofsn0YH\nxGrKXVdNTU2xuqam0s+Ka6JItZGnmr+VNuVsmmlxjeqhmdirYp465/Y5537qnPuJc666lV3HMNoI\nPo25pomA9AmkztSUSyyl0c2ZMydW48+X5tJLL43V+fPxaKJItZGnmr+VNuVsmmlxjRqT9E4QdwCv\nAB1l3jOsd7V6RRvB5821yJwMMdd8BGRkyIVEQHqTMSMwMfeaNBo03Jz0mij97nDWFNWkt9Wm0vV/\nq0gX8reKUs768UtSazZptOVoMRnrFWr4uGM+TM+oMrl7IkmMP4A9e16mMApy7969AW0JhalWRUJS\nwTbgDddJhEYlAjQ0tFAYsdrYGBKl2UBh+t1QmpvHAcfxJXuPBxmhPoo0igYNiSK94YYbKEzbe8MN\nnw3q3803f5FMppGmpt+TyTRy881fLKvZu/dlCv9WL730clBbAEmjLWsSOWlURtI7QdyBX7G/AOwA\nPlviPcN7W6tDtOapxsjTmIwXXnhhrObCCy8csn+atLOaCFLtWGiiSNM0NNOsKWrUHmplngLvy71O\nAnYCC2PeI+vWrTt1bNu2bXhHow7QmqcaI09jMjrnYjXODb3VoUk7q4kg1Y6FJoo0TUMzzZqiRvps\n27ataK6s2cRe9AthHXBTzPlhHYx6xFbsEbZij7AV+9iiJhM7PmXe+NzXbcAzwNKY9w335x8V9PX1\nycaNG8v+g8/jzdMWgbMEWoLNU2/kNQmME2gKrPOpjdBMnkrXp52NzN2QtLMaw1UkPxaRLmQsfBRp\n1L+QKFJtDU2NOVlJTVEzQkcXmom9GpGnk4GH/X/LaQLuE5GtVfi9dceaNX/D1772dWA68BqrV3+2\nbFrSrVufwNfE7ACOsnXr40FtDQ4ez+mmAgcYHBwIUOUjNPNRpCE0UGgYhuJrdGbIR56ePBlSh7SR\nQsM1NAVvW1s7g4O/J280jh8/vqymo+O9HD168FT/Ojs7A1rK9ynfvzA0EZA+GvSHQDvwdnBa3LqO\ntjQikt4JtAdjfMWu+a+6pvamVjd79uxYzezZs0tqtFsxGp12LDZs2BCr27BhQ0mNJpWuditGg22p\njC0YCXvsJRsa4xO7xlzT1N7U6iB5ZKfWPNXotGOxaNGiWN2iRYtKajSpdLXmqQYzQccWmondnj1P\nCU20oKb2plY3e/bsWI0/H88FF1wQq/HnS6PRacdCk05Xk0pXGw2qwaJBjbIkvRNoD8b4il2kMJXu\nJAlNpetrb0YmaEjtzUiXrGZnWult9W01nKYJG4vOzklFus7OSWU1PpVupAlJpRvV6/SaUPNUxG+v\n9Pb2Bm+nmAk6dsC2YkY2PuQ+KtQbFrLfWKQJnziT67wmeqojTFNciDl0su3omCyF6QFCJttx4yYW\naULznXd33y9NTeOkoWGCNDWNC54EN2/eLNdcc82Qe+unt5PNTpRsdnZwrvO8rrW1UyZMmJ9okk56\nMzBGJzaxj2A0aXE1hqaIyJQpU2J1U6ZMKamZP39+rGb+/PklNWkamtqaomkZjdp2zAg1ymET+whG\nkxYXhaGp1TU0NMRqGhpKr8DTNDS1NUXTMhq17ZgRapRDM7GbeZoSmrS4GkMTYMqUKbE6fz6eefPm\nxWr8+XjSNDS1NUXTMhq17ZgRagwLSe8E2oMxvmIXyafFzQh0CGQSpMVtFTgneK+8WJfUCI36F25o\nFpqnw2doamuKaup8atAammaEGkNBjSJPjUD836gROAvoRySkJmY+GnSQ8GjQvC5KOxuuifrnU90O\nzbhx7bzzzgD5aMu2tvJRnQD9/ceL+jcwUD4y9u67v86f/dlf8O67gzQ0NHL33V8PaguS1/nUoI3q\ntGhQo+okvRNoD8b4il0Tbamt8zl16tRY3dSpU0tqNEm2tIamRmfmpDFWwczTkYsm2pIUzVNNWlyt\noanRmTlpjFU0E7uZpymhibbU1vmcOnVqrM6fj2fhwoWxGn8+Hq2hqdGZOWkY4Th/Q0ihIb8yTaWt\nkYrf383iMy7+Er/PPvT+t0YT6ZrxmZSPAseHpa3x4ydw9Oggfl/+17S1ZXj77d+V7Z/XDeCzJ/6O\ntraWsrotWx7g+uv/iubmmRw/vp977tkQVKZNqzOMkYBzDhEJTxeKTeyp4lwTfrKdBBwmpO6kRhPp\nMsA0/Oq4/CTd2TmFN9/8NT778gk6Ozv5zW/eGFIzffoHOXDgIPmbwfTpk3n11V+U7Z9v6/en+tfZ\nOb5sWwCHDx9WmYxanWHUGs3EblsxKdHR0YGfaH8EvJp7zebOx+MnoDM15SamSy+9NKd7Dngp95rN\nnY/nrrvuyk20L+CfiHmBI0fe5q677iqpue+++3KT+nPAy8BzvPbaIe67774h+xe1FfWvXFt5Jk2a\nxB/+4R8mnpy1OsMYjdiKPSWcc8Bs/ESWZzbw75QaF40GoKmpiZMnP3iGrrHxF5w4cSJW89GPfpQf\n/ODgGZpFi97H008/HatZvnw53/rW7jM011wzl0cffbRk/zRtGcZYpWYrdufcMufcHufcS865L1Tj\nd9YbEydOJM4w9OfjOeuss2I1/nxp/uAP/iBW58/Hk1Z6W21bhmGEU/GK3XnH7SXgcuB1YAfwKRHZ\nc9r7xvSKHSoxNFvIm5Phe+zJjdD3vvdsjhx5+5QmZN97xowP8Nprh0i6x65pyzDGIrVasV8CvCwi\n+8WHUt4PLK/C761DCiM7GxNoHPl6naF0dJwNHAN+Dhyjs3PoVT7Apz/9aSioefrpT5dfQb/66i/Y\nvPlurrlmLps33x00qQP85jdvsGHDP7Fo0fvYsOGfbFI3jGqS9MH30w/gPwNfL/j+vwD/HPO+6j2x\nPwppbMznRy+OtmxsbCypWbZsWaxm2bJlQ7alSYubZs1OwzDCYaTnilm/fv2pr7u6uujq6kqz+Zpy\n8uRJ/KN9F+XOXARM5eTJfy+p6enpAaafofHnS7Nly5bYtrZs2cJf/uVfxmp6e3tj2ppGb29v2YAo\nwzCqR09PT9l/42VJeic4/QAuBb5b8P0XgS/EvG+Y72sjG1uxG4ahAcWKvRrmaSOwF2+e/groBVaK\nyO7T3ieVtjXaSTPyVGNOrllzI1/72t3kg4ZWr/4sd955R9m2DMMYPmoWeeqcWwbcgXf67hGRr8S8\nZ8xP7JB/xvwkjY2NJZ8pP52rrrqKnp4eurq6+M53vhPc1l133cWWLVtYuXJlyS2Y09m9eze9vb1c\ncskltgVjGCMASylgGIZRZ1hKAcMwDMMmdsMwjHrDJnbDMIw6wyZ2wzCMOsMmdsMwjDrDJnbDMIw6\nwyZ2wzCMOsMmdsMwjDrDJnbDMIw6wyZ2wzCMOsMmdsMwjDrDJnbDMIw6wyZ2wzCMOsMmdsMwjDrD\nJnbDMIw6o6KJ3Tm3zjl3wDn3Qu5YVq2OGYZhGDqqsWK/XUTm547vVuH31T0VF6qtI2wsImwsImws\nKqMaE3uiyh6GXbSF2FhE2FhE2FhURjUm9tXOuZ3OuW845yZU4fcZhmEYFVB2YnfOPeGce7Hg2JV7\n/QSwAfigiMwDDgK3D3eHDcMwjKGpWjFr59xM4DERuajEz62StWEYhoKkxaybKmnMOTdFRA7mvv0T\n4N+q1THDMAxDR0UTO/BV59w84F1gH/DfK+6RYRiGURFV24oxDMMwRgapRp6O9YAm59wy59we59xL\nzrkv1Lo/tcQ5t88591Pn3E+cc7217k/aOOfucc4dcs69WHCuwzm31Tm31zn3+Fh5yqzEWIy5ucI5\nN805933n3M9yD6l8Pnc+8XWR6ordObcOeEtExtzTM865BuAl4HLgdWAH8CkR2VPTjtUI59wrwEdE\n5M1a96UWOOcWAm8D38w/cOCcuw34jYh8NXfj7xCRL9ayn2lQYizG3FzhnJsCTBGRnc658cDzwHLg\nz0l4XdQiV8xYNVEvAV4Wkf0ichy4H/9HG6s4xnCuIhHZDpx+U1sObMp9vQn4ZKqdqhElxgLG2Fwh\nIgdFZGfu67eB3cA0FNdFLf5hjdWApqnAawXfH8idG6sI8IRzbodz7rO17swI4WwROQT+Hzlwdo37\nU2vG6lyBc24WMA94Dpic9Lqo+sRuAU1GIJeJyHzgauCvc/8dN4oZy082jNm5IrcN8xBwY27lfvp1\nUPa6qPRxxzMQkSsD33o38Fi12x/B/BKYUfD9tNy5MYmI/Cr3etg59zB+q2p7bXtVcw455yaLyKHc\nfusbte5QrRCRwwXfjpm5wjnXhJ/U/1VEHs2dTnxdpP1UzJSCb4cMaKpDdgDnOudmOucywKeAb9W4\nTzXBOTcutyrBOdcGLGVsXQt5HMX7yN8CVuW+vg549HRBHVM0FmN4rvi/QJ+I3FFwLvF1kfZTMd/E\n7xudCmjK7x2NBXKPbN2Bv6HeIyJfqXGXaoJz7gPAw/j/UjYB9421sXDOdQNdwHuBQ8A64BHg/wHT\ngf3AtSLy21r1MS1KjMXHGGNzhXPuMuBpYBf+34YAfw/0Ag+S4LqwACXDMIw6Y8w+bmYYhlGv2MRu\nGIZRZ9jEbhiGUWfYxG4YhlFn2MRuGIZRZ9jEbhiGUWfYxG4YhlFn2MRuGIZRZ/x/MvIU1I8cXwsA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bf62748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = list(range(20))*100\n",
    "y = []\n",
    "for i in x:\n",
    "    y.append(rd.poisson(i+1)) # random generator following poisson distribution. λ=x+1 here.\n",
    "df = pd.DataFrame({'x':x,'y':y})\n",
    "\n",
    "fig,ax = plt.subplots(1,1)\n",
    "chart = ax.scatter(df.x,df.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "このデータは各xでスライスしてy値のヒストグラムを作るとそれぞれλ=x+1のポアソン分布の形に近くなっているはずである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAGxCAYAAABGJTP8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X24XHV56P3vTcD4Au5ELNAKBjyeqvVosbXAuSzpBsX6\nwktSn4OplQM21r6KBk4rKJqg7VG8GoL12HPaGjkxtSdan8Nr7QP6JLuh+iTSAkoVUY+AWiXaQDZS\nlMZwP3+stclkZ2bPzF7zsmbP93Ndc+2ZNWvWuvfaM/e+5/f7rd+KzESSJEnzc8iwA5AkSRplFlOS\nJEkVWExJkiRVYDElSZJUgcWUJElSBRZTkiRJFVhMCYCIeEJEfCkiji4fXx0R7+7Tvh6LiGd1uO6y\ncv2BvVcj4o8j4rcGtT9J1ZnDDtinOWzALKYWqDKxfCQipiPiOxGxps1L3gT8XWbuGkB43U5u1tPJ\n0CLisIj464i4p0xyy2et8sfA2yPi0F7uV1LnymLo0Yh4KCJ+UP6MOV4yNjkMICKeFBF/GhHfj4gH\nI2Kq4Wlz2IBZTC1clwP/DjgOOB34g4h4+Rzr/xaweRCBAXMlxEG5Bfg14Luzn8jM+4G7gLMHHZSk\nA1yRmU/NzCPKn3MVJeOWw/4CWAI8B3ga8PgXZnPY4FlM1UxEPCsidkfEieXjn4qI7zVpPWnnPwPv\nzsyHMvMrwJ8DF7TY53HACcDOFs8viYgbyjh2l/ef0fD8toh4T0R8tvwGeV1EPC0i/rJsGdsZEc+c\ntdlXR8T/Kbf5/oZtHVI2UX8/Ir4OvHpWLBdExJfLb6lfj4g3dXlcyMy9mfknmfk54LEWq/3d7H1L\naq+HOaybfY5VDouI5wBnAm/KzAeycPus1cxhA2QxVTOZ+Q3gD4C/jIgnAVcDV2fmdoCI+FDZpPtA\nw8+Z+3eU6ywBfhL4YsOmvwA8v8VuXwB8IzNbFRaHAB+haOV6JvAI8N9mrfNaipaenwKeDXwO2Ags\nBb4CrJ21/grg58rbORHx6+XyNwGvAn4WeDHwf8163S7gVZn5VOANwIaGpH3crGMz+zitavH7NXNX\nGYOkLvQihzX4nYj4l4i4NSJ+ZY7djlsOOwm4D3h3WbR9ocnxMYcNUmZ6q+ENuJaiGLoDOKzL1x4L\n7AOe0LDsZRTJptn6rwM+N2vZ1RQtW83WPxHY3fB4G3Bpw+M/Bv6m4fGZwG0Njx8Dzmh4/NvAp8v7\n/y/Ft62Z584of5dDWsRyDfDmCsf5W8DyJstfBnx92O8Db95G9VYlh5WvP5GikDkEeCXwEPAfW6w7\nVjkMuLSM4Z3AocBy4AfAcxrWMYcN8GbLVH19mKIl6YOZubfL1z5c/nxqw7IJig9bMw8CR7TaWDnQ\n8c8i4t6I2EPRfLwk4oDBoI2DPn/Y5PHhszb77Yb791F8G6T8+a1ZzzXG8sqI+P/KpvoHKZLs01vF\nXsERwJ4+bFcaF1VyGJl5R2Y+mJmPZebfAh8DWrVOjVsO+yHwb8AfZuaPs2j12wY0jos1hw2QxVQN\nRcRTgKsompjXld12M8/999h/Zkvj7QcRcSdAZu6hGFjd2MT7s8CXWuzyi8AJ0frU3f8C/HvgFzJz\nCcW3IKg2CPO4hvvLgO+U97/b5LliZxFPAD4JvB/4icxcCvztTBxlE3mrY/NQRPxqF/E9j6JrVFKX\nquawFpLWOWfcctjMEI7G+GcPzjeHDZDFVD39CfD5zHwT8Cngz2aeyMzfzv1ntjTejsjMFzRsYzNw\nWRQDL58H/AZFs/dBMvOfga9T9MM3czjFN6GHIuJpwLqqvyDw+2VsxwEXAlvK5Z8ALoyIZ0TEUuBt\nDa95Qnn7l8x8LCJeScM3scz81hzH5qmZ+b9m1o1i6ognlg8XR8TiWfH9EkWSk9S9yjksIl4TEU+J\nwsspxjNd12xnY5jDtgPfBC6NiEUR8RJgEripYV/msAGymKqZiDib4sP1O+Wii4AXddmqAsVgyW9Q\nNDFvBd6XmZ+eY/0/ozgDsJmrgCcD/0IxKPNTs56fz5wr1wH/CNwG3EAxOBSK031vovhG9Q/A//34\nizIfpkhafx0RDwCraJFcO3A38K8UTfL/D/BIlGfrRMRPUnyru3ae25bGVg9z2FsoutIeBK4A3piZ\nt8yx/tjksMz8MXAOxdl6eyh+9/My86tgDhuGyGz/HoqIe4FpigFvezPzpLLi/jhFE+a9wLmZOd2/\nUNVPZfPzbcBLczCT3tVWRPwxxcDN/zHsWFSd+Ws8mMP2M4cNXqfF1DeAn8/MBxuWXUFxNsT7I+Jt\nwNLMvKR/oUpS98xfkvqt026+aLLuOcCm8v4mijk3JKluzF+S+qrTYiqBT0cxcdoby2VHzzSlZjF1\n/VH9CFCSKjJ/SeqrTi+C+JLM/G5E/ARwc0TczcED9pr2F0ZEzy/wKKn+MrMO1y8D85ekLnWbvzpq\nmcrM75Y/v09xdsBJwK6IOBogIo4BvjfH62t/W7t27dBjWEhxjlKsxtn7W53kGOSvUXp/GOd4xjlK\nsc5H22IqIp4cEYeX959CccrrncD17L9w7vnM/xR1SeoL85ekQeikm+9o4JqyuftQ4GOZeXNE/APw\niSgu7ngfcG4f45Sk+TB/Seq7tsVUZt5DcVHI2csfoLiQ4oIwOTk57BA6MipxwujEapwL17jkLxid\n94dx9taoxAmjFWu3OppnqtIOIrLf+5BULxFB1mcA+ryZv6TxM5/85eVkJEmSKrCYkiRJqsBiSpIk\nqQKLKUmSpAospiRJkiqwmJIkSarAYkqSJKkCiylJkqQKLKYkSZIqsJiSJEmqwGJKkiSpAospSZKk\nCiymJEmSKrCYkiRJqsBiSpIkqQKLKUmSpAospiRJkiqwmJIkSarAYkqSJKkCiylJkqQKLKYkSZIq\nsJiSJEmqwGJKkiSpAospSZKkCiymJEmSKrCYkiRJqqDjYioiDomI2yLi+vLx0oi4OSLujoibImKi\n1WtXrryM7dt39CJeSepalfwlSe100zL1FuDLDY8vAT6Tmc8BtgKXtnrhtdeu5eyzv8bFF181vygl\nqZp55y9JaqejYioijgVeBXy4YfE5wKby/iZgRestHMb09Hls3Hgkt9yyc36RStI8VM9fkjS3Tlum\nNgC/D2TDsqMzcxdAZt4PHNVuI9PTq7jyyhu6DlKSKuhJ/pKkVg5tt0JEvBrYlZl3RMTkHKtm66fW\nPX7v3nu/1XFwkkbD1NQUU1NTww7jIL3IX+vWrXv8/uTkJJOTc21G0qjpRf6KzDlqICAi/ivweuDH\nwJOAI4BrgBcDk5m5KyKOAbZl5vOavD7356m9rFhxOddc84eVgpZUbxFBZkYN4qicv9rlSEkLy3zy\nV9tuvsx8e2Y+MzOfBawCtmbmecANwAXlaucD17Xb1sTEFi666Kxu4pOkeetl/pKkVqrMM/U+4IyI\nuBt4afm4hb0sXryZ1at3c+qpJ1fYpST1RBf5S5Lm1rabr/IOIvK5z30Hp5xyFldfbSEljYO6dPNV\nZTefNH7mk7/aDkDvhde+1jFSkiRpYfJyMpIkSRVYTEmSJFVgMSVJklSBxZQkSVIFFlOSJEkVWExJ\nkiRVYDElSZJUgcWUJElSBRZTkiRJFVhMSZIkVWAxJUmSVIHFlCRJUgUWU5IkSRVYTEmSJFVgMSVJ\nklSBxZQkSVIFFlOSJEkVWExJkiRVYDElSZJUgcWUJElSBRZTkiRJFVhMSZIkVWAxJUmSVIHFlCRJ\nUgUWU5IkSRW0LaYiYnFE7IyI2yPizohYWy5fGhE3R8TdEXFTREz0P1xJ6pz5S9IgtC2mMvNR4LTM\nfBFwIvDKiDgJuAT4TGY+B9gKXNrXSCWpS+YvSYPQUTdfZj5S3l0MHAokcA6wqVy+CVjR8+gkqSLz\nl6R+66iYiohDIuJ24H7g05l5K3B0Zu4CyMz7gaP6F6YkzY/5S1K/HdrJSpn5GPCiiHgqcE1EPJ/i\n290Bq7V6/dTUOgDWrYPJyUkmJyfZvn0HGzbcyJ49i1iyZB9r1pzJ8uWnzOuXkDRcU1NTTE1NDTuM\npqrmr3Xr1j1+fyZ/SVo4epG/IrNlDmn+goh3Ao8AbwQmM3NXRBwDbMvM5zVZP9euLfYxk5Muvvgq\nNm48kunpVcBhwF4mJrawevVu1q9/a4VfR1IdRASZGcOOY7b55K9uc6Sk0Taf/NXJ2XxPnznTJSKe\nBJwB3AVcD1xQrnY+cF0nO9y+fUdZSJ1HUUgBHMb09Hls3Hgkt9yys5v4JamlXucvSWqmkzFTPwls\ni4g7gJ3ATZn5KeAK4IyIuBt4KfC+Tna4YcONZYvUwaanV3HllTd0FLgkdaCn+UuSmmk7Zioz7wR+\nrsnyB4CXdbvDPXsWsb9FarbDmJ5e1O0mJampXucvSWpm4DOgL1myD9jb4tm9TEzsG2Q4kiRJlQy8\nmFqz5kwmJrY0fW5iYgsXXXTWgCOSJEmav4EXU8uXn8Lq1btZvHgz+1uo9rJ48WZWr97NqaeePOiQ\nJEmS5q2jeaZ6bf36t/LAAzvZseNyfvSjRTzxifs45ZSzWL/eQkqSJI2WoRRTAMuWncyyZRZPkiRp\ntA28m0+SJGkhsZiSJEmqwGJKkiSpAospSZKkCiymJEmSKrCYkiRJqsBiSpIkqQKLKUmSpAospiRJ\nkiqwmJIkSarAYkqSJKkCiylJkqQKLKYkSZIqsJiSJEmqwGJKkiSpAospSZKkCiymJEmSKrCYkiRJ\nqsBiSpIkqQKLKUmSpAospiRJkiqwmJIkSaqgbTEVEcdGxNaI+FJE3BkRF5bLl0bEzRFxd0TcFBET\n/Q9Xkjpn/pI0CJ20TP0YuCgznw/8R+B3I+K5wCXAZzLzOcBW4NL+hSlJ82L+ktR3bYupzLw/M+8o\n7z8M3AUcC5wDbCpX2wSs6FeQkjQf5i9Jg3BoNytHxPHAicAO4OjM3AVFwoqIo3od3PbtO9iw4Ub2\n7FnEkiX7WLPmTJYvP6XXu5E0Bgadv3rBHCiNho6LqYg4HPgk8JbMfDgictYqsx8/bmpqHQDr1sHk\n5CSTk5Nt93fxxVexceORTE+vBQ4D9rJt2xZWr97B+vVv7TRsSQMwNTXF1NTUsMNoqUr+Wrdu3eP3\nO81fvWAOlAajF/krMlvmkP0rRRwK3Aj8bWZ+oFx2FzCZmbsi4hhgW2Y+r8lrc+3aYh8NOemA+7OX\nbd++g7PP/hrT0+cdtM7ExGZuuOGnOfXUk9vGLWk4IoLMjGHHAdXzVyc5stfMgdLwzCd/dTo1wkeA\nL88kotL1wAXl/fOB67rZ8Vw2bLiR6elVTZ+bnl7FlVfe0KtdSVr4Bpq/esEcKI2Wtt18EfES4NeA\nOyPidorm8LcDVwCfiIhfB+4Dzu1VUHv2LKJo1m7mMKanF/VqV5IWsGHkr14wB0qjpW0xlZmfBVp9\ncl/W23AKS5bsA/bSPJnsZWJiXz92K2mBGUb+6gVzoDRaajkD+po1ZzIxsaXpcxMTW7joorMGHJEk\nDY45UBottSymli8/hdWrd7N48WaKb2cAe1m8eDOrV+924KWkBc0cKI2WruaZGqT169/KAw/sZMeO\ny/nRjxbxxCfu45RTzmL9epOIpIXPHCiNjtoWUwDLlp3MsmUmDknjyRwojYZadvNJkiSNCospSZKk\nCiymJEmSKrCYkiRJqqDWA9AlSf21ffsONmy4kT17FrFkyT7WrDmT5ctPGXZY0kixmJKkMXXxxVex\nceORTE+vpZhtfS/btm1h9eodrF//1mGHJ40Mu/kkaQxt376jLKTOY/9law5jevo8Nm48kltu2TnM\n8KSRYjElSWNow4YbmZ5e1fS56elVXHnlDQOOSBpddvNJ0gLWakzUnj2LaH4hZShaqFpdH1rSbBZT\nkrRAzTUmasmSfRTX/WtWUO1lYmLfQGOVRpndfJK0ALUbE3X66ScwMbGl6WsnJrZw0UVnDSxWadRZ\nTEnSAtRuTNTWrfewevVuFi/eTNFCBbCXxYs3s3r1bk491WsCSp2ym0+SFqBOxkStX/9WHnhgJzt2\nXM6PfrSIJz5xH6ecchbr11tISd2wmJKkBajTMVHLlp3MsmUWT1IVdvNJ0gK0Zs2ZjomSBsRiSpIW\noOXLT3FMlDQgdvNJ0gLlmChpMCymJGkBc0yU1H9280mSJFVgy5QkaU6tLkkjqWAxJUlqaa5L0qxf\n/9ZhhyfVgt18kqSm2l2S5pZbdg4zPKk22hZTEbExInZFxBcbli2NiJsj4u6IuCkiJvobZnPbt+9g\n5crLOO20taxceRnbt+8YRhiSaqzOOazu2l2S5sorbxhwRFI9ddIydTXwy7OWXQJ8JjOfA2wFLu11\nYO1cfPFVnH3217j22rVMTV3Otdeu5eyzv8bFF1816FAk1Vstc9go6OSSNJI6KKYy8++BB2ctPgfY\nVN7fBKzocVxzsulZUqfqmMNGxf5L0jSz/5I00rib75ipozJzF0Bm3g8c1buQ2rPpWVJFQ81ho8JL\n0kid6dXZfDnXk1NT6wBYtw4mJyeZnJystDObnqV6mZqaYmpqathhVNEyh61bt+7x+73IX73Wz2kL\nikvS7OBDH9rMo4+uYuZsvsWLt3hJGi0Yvchf8y2mdkXE0Zm5KyKOAb4318qTk+uAopjqhU6vhi5p\nMGYXGZdffvnwgulMxzlsXa8SVx8MYtoCL0mjha4X+avTYirK24zrgQuAK4Dzgeu63nMFa9acybZt\nW8oxUwey6VlSE7XKYb1w4NjRGTNjRzezYsXOnrUceUkaaW6dTI3wV8DngJ+OiG9GxBuA9wFnRMTd\nwEvLxwPj1dAldaqOOawXHDsq1UfblqnMfF2Lp17W41i6YtOzpE7UNYdV5dhRqT5G+nIyNj1LGleO\nHZXqw8vJSNIIctoCqT4spiSpD/p9uSvHjkr1MdLdfJJUR4OYsgAcOyrVhcWUJPXQIKcsAMeOSnVg\nMSVJ89Bq5vFiyoK1TV9TTFlw+ePFVD9nLx+khfJ7SPNlMSVJXZqrG6/TKQsG1RXYbwvl95CqcAC6\nJHXhwG68maJpphvvSH784++zf0D4bMWUBe22ccstO/v9a/TEQvk9pKospiSpC+1mHodH205ZsFBm\nL18ov4dUlcWUJHWhXTfeYYcd23bKgoUye/lC+T2kqhwzJUld6GTm8XZTFiyU2csXyu8hVWUxJUld\nWLPmTLZt2zJr6oNC48zjc01Z0Ok26m6h/B5SVXbzSVIXejHz+EKZvXyh/B5SVbZMSVKXejHz+EKZ\nvXyh/B5SFRZTkjQPvZh5fKHMXr5Qfg9pvuzmkyRJqsCWKUnS0HlJGo0yiylJ0lB5SRqNOrv5JElD\n4yVptBBYTEmShsZL0mghWPDdfPbDS1J9DfqSNP5PUD8s6GLKfnhJqrdBXpLG/wnqlwXbzWc/vCTV\n35o1ZzIxsaXpc728JI3/E9RPC7aYsh9ekupj+/YdrFx5GaedtpaVKy9j+/YdQHeXpGm1jU74P0H9\ntGC7+QbdDy9Jaq5d91onl6Sp2kXn/wT104ItpgbZDy9Jau7A7rUZM91rm1mxYiennnrynJek6XQb\nc/F/gvppwXbzDaofXpLUWi+613qxDf8nqJ8qtUxFxCuAqyiKso2ZeUVPouqBoh9+Bx/60GYefXQV\nM83CixdvOaAffuY02Xvu+RYnnHBc09Nk251K28mptlW3MSpxdhLroOOY7zEdlTgHHcdCOqW8kxy2\ncuVlI/G73nvvFMcfPznsMA5ycPfaFDBZ3u+se63TLrq53qN1+p8wbnmhDrm27/krM+d1o0g+XweW\nUbwr7wCe22S9XLs2c+3aPMDMssbbbO3W6WQbF1ywI5/73Hfk8ce/K5/73HfkBRfsePy5iy7akBMT\nH034t4S1Cf+WExMfzYsu2tBinTxonXbP92IboxJnJ7EOJ47uj+moxFmnv32jIrXML7cM6tZJDgNa\n/q69yE+93MYv/dLaeW2j33GuWPGOhvdNlu/jfPx9tGLFO+axjTxoG52+R4f9P2Hc8kIdcu0g8leV\nlqmTgK9l5n0AEbEFOAf4SoVt9lyrfvhO+uAzc851jj9+Y9+3MYh9uI3RjLNO22g3XqWmOsxhC+J3\nHZo1a85k27Yts94/hU6719pt4/TTT+Cd73xCR+/RYf5PGKXP9ELZxqDyV5Vi6hnAtxoef5siOY2E\nog9+bdPnij74y8v7rdd573tfy/T0x/u6jUHsw22MZpx12saIFhhd5bAR/12H5uDuNWjWvdbdNg7s\notu69V/a5vN2+xnE/4RR+kwvlG0MKn9F0aI1jxdGvAb45cx8U/n49cBJmXnhrPXmtwNJIy0zY9gx\nzKWTHGb+ksZTt/mrSsvUPwPPbHh8bLmsUkCSNCBtc5j5S1InqkyNcCvw7IhYFhFPAFYB1/cmLEnq\nO3OYpJ6Yd8tUZu6LiN8Dbmb/acV39SwySeojc5ikXpn3mClJkiT1cQb0iHhFRHwlIr4aEW/r1356\nISLujYgvRMTtEfH5YcczIyI2RsSuiPhiw7KlEXFzRNwdETdFxMQwY5zRIta1EfHtiLitvL1iyDEe\nGxFbI+JLEXFnRFxYLq/dMW0S65vL5XU7posjYmf52bkzItaWy2t3TLs1KjmsrvkLRieHjUL+KmMa\niRw2jvmrLy1TEXEI8FXgpcB3KMYmrMrMWs1BNSMivgH8fGY+OOxYGkXELwIPAx/NzBeWy64Admfm\n+8sEvzQzLxlmnGVczWJdC/wgM68canCliDgGOCYz74iIw4F/pJhX6A3U7JjOEetrqdExBYiIJ2fm\nIxGxCPgscCHwGmp2TLsxSjmsrvkLRieHjUL+gtHJYeOYv/rVMvX4ZHiZuReYmQyvroIaXqcwM/8e\nmJ0gzwE2lfc3ASsGGlQLLWKF4tjWQmben5l3lPcfBu6iOIOrdse0RazPKJ+uzTEFyMxHyruLKcZh\nJjU8pl0apRxWy/wFo5PDRiF/wejksHHMX/36ADabDO8ZLdatgwQ+HRG3RsRvDDuYNo7KzF1QvGGB\no4YcTzu/FxF3RMSHh9303CgijgdOBHYAR9f5mDbEurNcVKtjGhGHRMTtwP3ApzPzVmp+TDswSjls\nlPIXjFYOq9VnrdGo5LBxyV+1/DYzBC/JzJ8DXgX8btnkOyrqfAbBnwLPyswTKd6otWjaLZudPwm8\npfzWNPsY1uaYNom1dsc0Mx/LzBdRfEM+KSKeT42P6QI0yvkL6vveqN1nbcao5LBxyl/9KqY6mtCz\nLjLzu+XP7wPXUO/L4uyKiKPh8X7p7w05npYy8/u5f1DeXwC/MMx4ACLiUIoP9+bMvK5cXMtj2izW\nOh7TGZn5EDAFvIKaHtMujEwOG7H8BSPy3qjrZ21Ucti45a9+FVMjMxleRDy5rJ6JiKcALwf+abhR\nHSA4sI/5euCC8v75wHWzXzBEB8Ravgln/Ar1OK4fAb6cmR9oWFbXY3pQrHU7phHx9Jmm+oh4EnAG\nxfiIuh7TTo1EDhuB/AWjk8NGIX/B6OSwscpffZtnqjzl8QPsnwzvfX3ZUUURcQLFt7mkGHz2sbrE\nGhF/BUwCRwK7gLXAtcBfA8cB9wHnZuaeYcU4o0Wsp1H0lT8G3Av85kw/9DBExEuA7cCdFH/vBN4O\nfB74BDU6pnPE+jrqdUxfQDFA85Dy9vHM/KOIeBo1O6bdGoUcVuf8BaOTw0Yhf8Ho5LBxzF9O2ilJ\nklSBA9AlSZIqsJiSJEmqwGJKkiSpAospSZKkCiymJEmSKrCYkiRJqsBiSpIkqQKLKT0uIn45Iv53\nw+PHIuJZfdjP+RFxSxfrr42Izb2OY74i4syI2DLsOCTtZ/7qjPmrPyymRlhE/KeI+GxE/GtEbG3y\n/IkR8Q/l87dGxM+22eQfAu9teNzPGV273XZPY4mIyYjYGhF7IuIbTZ5/d0R8MSL2RsS7Dggk80bg\nZyLiP/QyJmmcdJC/HouIH5S3hyLiz9ts0vx14DpviYhvRMTDEfGliHg2mL/6xWJqtO0GNnBgAgEg\nIg6juGzDR4El5c/ryotPHiQiXgw8NTNvbVzc84jr41+BjcB/afH814DfB25s8fwW4Df7EJc0Llrm\nr1ICL8zMIzLzqZn5plYbMn8dKCLeCLwBeGVmHg6cCfxLwyrmrx6zmBqCiHhWROyOiBPLxz8VEd+L\niOXdbCczt2bmJ4HvNnl6EliUmX+SmXsz84MUyeX0Fpt7JfB3c8T8qoi4LSKmI+K+iFjb8Nyy8lvk\nBRHxzfJ3+82IeHFEfCEiHoiID87a5CER8cHym9WXI+L0hu0dHxFT5b5uAp4+K5ZPRMR3I+LBcr2f\nmes4NZOZt2bmx4B7Wjy/OTNvAh5usYkp4NXd7lcadQPKX1Dkq07/R5m/9m8/gHcBazLz7nL9e2Zd\nW24K81dPWUwNQWZ+A/gD4C+juFL11cDVmbkdICI+VH7QHmj4OXP/jg5383zgi7OWfaFc3swLgLvn\n2N7DwHmZOUHxIfytiDh71jonAc8GXgtcRXFhy9OB/wCcGxGnNqx7MkXrz5HAOuB/R8SS8rm/Am6l\nSEJ/SHHV7kafAv4dcBRwG/CxmSci4m2zjt0Bx3GO369bdwHLIuLwHm5Tqr0B5a8ZfxcR34mIT0bE\nsjnWM3/td2x5e0FZHP6fiFg3ax3zV49ZTA1JZm4Evg7sBI4GLmt47nczc2lmPq3h58z9EzvcxeHA\n9KxlDwFHtFh/CfCDOeLdnplfKu//E0Uz8S81rgK8OzP/LTM/Q9EM/b8yc3dmfge4BXhRw/q7ylaz\nfZn5CYpE+OqIOA54MfCuskXtFuCGWbH8z8x8JDP3Au8GfjYijiifu2LWsTvgOLY8Wt37AcU35yXt\nVpQWmgHkL4DlwPHAcylar26MiFb/s8xf+x1b/jyD4svz6cCvRsTqhnXMXz1mMTVcH6Z4s3+w/GD1\n0sPAU2ctm6B1wnmQ1oUWEXFyFAMevxcReyj6258+a7XvNdz/IbBr1uPGb0H/POu19wE/Vd4ezMwf\nznpuJo5DIuJ9EfH1Mo57KBLh7Fj67Yhyv3varSgtUP3MX2Tm32fmjzPzIeAtFIXV81qsbv46MFaA\nKzLzB5lVMrLxAAAXtUlEQVR5H/BnwKsa1jF/9ZjF1JBExFMompI3AusamoiJiP8e+89gabz9ICLu\n7HAXXwJeOGvZC8vlzXwR+Ok5tvcxigHtz8jMJRQfzioDPJ8x6/Ezge9QfANdWnYfND4349eAs4DT\nyziOL+MIgIi4dI5j91CFeGd7HnBvZrYaUyUtWAPIXwftctbP2cxf+90N/NusZbPPJjR/9ZjF1PD8\nCfD58gyVT1F8uAHIzN9uOIOl8XZEZr5gZr3yW85i4DBgUUQsjv1n600B+yLizRHxhIi4EHgMOOgU\n5NKnKAatt3I4xTeuvRFxEvC6Wc93m5iOLmM7NCL+E0VT/t9k5jeBfwAuj4jDIuIXKZJPYxyPAg+W\nCf29NCSKzHzvHMfu8Za6KCwGnkAxmHRxFGdAzjx/aEQ8keIzclj5fOPn5ZeAv+3yd5YWir7mr4j4\nmYj42XKdw4ErgW9TjPVpxvxV5q+yVWwL8AcRcXhEHAu8iQO7G81fPWYxNQTlwMeXA79TLroIeFFE\n/GqXmzqPokn3Q8AvAo8Afw5QNruvoBj8+CDwn4FzMvPHzTaUmbcDeyLiFxoXN9z/HeA9ETFNMT7i\n47M30eXjHcC/pzhd9z3Aa3L/2SavA06hOHX6ncCmhtd9FPgmRTP7PwGfa/b7dGA5xbG7ETiO4tjd\n1PD8X5TLVlEMRH0EeH3D879Kwz8QaVwMIn9RjMP6OMW4z69TfEbPzMx9zTZk/joof72ZYtzXd4DP\nAn+Zmf+z4XnzV49FZvu5xCLiXoo39WPA3sw8KSKWUrwhlwH3Audm5uwBzxohEXEG8NuZ+SvDjqXO\nIuJM4PWZuWrYsag989d4MH91xvzVH50WU98Afj4zH2xYdgWwOzPfHxFvA5Zm5iX9C1WSumf+ktRv\nnXbzNZs87Rz2N19uouhSkqS6MX9J6qtOi6kEPh3F9d3eWC47OjN3AWTm/RQTkElS3Zi/JPVV0+u0\nNfGSzPxuRPwEcHNE3E37AXoAREQ/LzYpqaYysy7XRjN/SepKt/mro5apzPxu+fP7FHN1nATsioij\nASLiGA6c8Gz262t/W7t27dBjWEhxjlKsxtn7W53kGOSvUXp/GOd4xjlKsc5H22IqIp5czvMxM1Hb\ny4E7geuBC8rVzgeum1cEktQn5i9Jg9BJN9/RwDVlc/ehwMcy8+aI+AfgExHx6xTT5Z/bxzglaT7M\nX5L6rm0xlZn3AAddnDIzHwBe1o+ghmFycnLYIXRkVOKE0YnVOBeucclfMDrvD+PsrVGJE0Yr1m51\nNM9UpR1EZL/3IaleIoKszwD0eTN/SeNnPvnLy8lIkiRVYDElSZJUgcWUJElSBRZTkiRJFVhMSZIk\nVWAxJUmSVIHFlCRJUgUWU5IkSRVYTEmSJFVgMSVJklSBxZQkSVIFFlOSJEkVWExJkiRVYDElSZJU\ngcWUJElSBRZTkiRJFVhMSZIkVWAxJUmSVIHFlCRJUgUWU5IkSRVYTEmSJFVgMSVJklSBxZQkSVIF\nFlOSJEkVWExJkiRVYDElSZJUQcfFVEQcEhG3RcT15eOlEXFzRNwdETdFxET/wpSk+TN/Seqnblqm\n3gJ8ueHxJcBnMvM5wFbg0l4GJkk9ZP6S1DcdFVMRcSzwKuDDDYvPATaV9zcBK3obmiRVZ/6S1G+d\ntkxtAH4fyIZlR2fmLoDMvB84qsexSVIvmL8k9dWh7VaIiFcDuzLzjoiYnGPVbPXEunXrHr8/OTnJ\n5ORcm5E0aqamppiamhp2GAcxf0lqpxf5KzJb5pBihYj/Crwe+DHwJOAI4BrgxcBkZu6KiGOAbZn5\nvCavz3b7kLSwRASZGTWIw/wlqSvzyV9tu/ky8+2Z+czMfBawCtiamecBNwAXlKudD1zXZbyS1Ffm\nL0mDUGWeqfcBZ0TE3cBLy8eSNArMX5J6pm03X+Ud2EwujZ26dPNVZf6Sxk9fuvkkSZLUmsWUJElS\nBRZTkiRJFVhMSZIkVWAxJUmSVIHFlCRJUgUWU5IkSRVYTEmSJFVgMSVJklSBxZQkSVIFFlOSJEkV\nWExJkiRVYDElSZJUgcWUJElSBRZTkiRJFVhMSZIkVWAxJUmSVIHFlCRJUgUWU5IkSRVYTEmSJFVg\nMSVJklSBxZQkSVIFFlOSJEkVWExJkiRVYDElSZJUgcWUJElSBW2LqYhYHBE7I+L2iLgzItaWy5dG\nxM0RcXdE3BQRE/0PV5I6Z/6SNAiRme1XinhyZj4SEYuAzwIXAq8Bdmfm+yPibcDSzLykyWuzk31I\nWjgigsyMYccB5i9J3ZlP/uqomy8zHynvLgYOBRI4B9hULt8ErOhmx5I0COYvSf3WUTEVEYdExO3A\n/cCnM/NW4OjM3AWQmfcDR/UvTEmaH/OXpH47tJOVMvMx4EUR8VTgmoh4PsW3uwNWa/X6devWPX5/\ncnKSycnJrgOVVF9TU1NMTU0NO4ymzF+S5tKL/NXRmKkDXhDxTuAR4I3AZGbuiohjgG2Z+bwm6zvm\nQBozdRoz1cj8JamdvoyZioinz5zpEhFPAs4A7gKuBy4oVzsfuK6raCWpz8xfkgahbctURLyAYoDm\nIeXt45n5RxHxNOATwHHAfcC5mbmnyev9ZieNmbq0TJm/JHVrPvmr626+bpmMpPFTl2KqKvOXNH76\nNjWCJEmSmrOYkiRJqsBiSpIkqQKLKUmSpAospiRJkiqwmJIkSarAYkqSJKkCiylJkqQKLKYkSZIq\nsJiSJEmqwGJKkiSpAospSZKkCiymJEmSKrCYkiRJqsBiSpIkqQKLKUmSpAospiRJkiqwmJIkSarA\nYkqSJKkCiylJkqQKLKYkSZIqsJiSJEmqwGJKkiSpAospSZKkCiymJEmSKrCYkiRJqsBiSpIkqYK2\nxVREHBsRWyPiSxFxZ0RcWC5fGhE3R8TdEXFTREz0P1xJ6pz5S9IgRGbOvULEMcAxmXlHRBwO/CNw\nDvAGYHdmvj8i3gYszcxLmrw+2+1D0sISEWRm1CAO85ekrswnf7VtmcrM+zPzjvL+w8BdwLEUCWlT\nudomYEV34UpSf5m/JA1CV2OmIuJ44ERgB3B0Zu6CImEBR/U6OEnqFfOXpH45tNMVyybyTwJvycyH\nI2J223fLtvB169Y9fn9ycpLJycnuopRUa1NTU0xNTQ07jJbMX5Ja6UX+ajtmCiAiDgVuBP42Mz9Q\nLrsLmMzMXeW4hG2Z+bwmr3XMgTRm6jJmCsxfkrrTlzFTpY8AX55JRKXrgQvK++cD13WzY0kaEPOX\npL7q5Gy+lwDbgTspmsITeDvweeATwHHAfcC5mbmnyev9ZieNmbq0TJm/JHVrPvmro26+KkxG0vip\nSzFVlflLGj/97OaTJElSExZTkiRJFVhMSZIkVWAxJUmSVIHFlCRJUgUWU5IkSRVYTEmSJFVgMSVJ\nklSBxZQkSVIFFlOSJEkVWExJkiRVYDElSZJUwaHDDkCd2b59Bxs23MiePYtYsmQfa9acyfLlpww7\nLEmSxl70+4roXnW9uosvvoqNG49kenoVcBiwl4mJLaxevZv169867PCkg8znqut1ZP6Sxs988pfF\nVM1t376Ds8/+GtPT5x303MTEZm644ac59dSThxCZ1JrFlKRRNZ/85Zipmtuw4cayRepg09OruPLK\nGwYckSRJauSYqZpoNSZqz55FFF17zRzG9PSiQYYpSbVVl7GldYlDg2MxVQP7x0StZWZM1LZtW1i9\negdLluwD9tK8oNrLxMS+gcYqSXU0Vx4d5NjSusShwbKbb8i2b99RfvDOY3/BdBjT0+exceORnH76\nCUxMbGn62omJLVx00VkDi1WS6qhdHr3llp1jFYcGz2JqyNqNidq69R5Wr97N4sWbKVqoAPayePFm\nVq/e7eBzSWOvLmNL6xKHBs9uviHrZEzU+vVv5YEHdrJjx+X86EeLeOIT93HKKWexfr2FlCTVZWxp\nXeLQ4FlMDVmnY6KWLTuZZcssniRptrqMLa1LHBo8u/mGbM2aMx0TJUkV1CWP1iUODZ7F1JAtX36K\nY6IkqYK65NG6xKHBcwb0mnjDG3ayY8cNB4yJuvrq/R+8desOfk2zZVIdOAO6hqFdHh23ODQ/88lf\njpmqCcdESVI1dcmjdYlDg9O2mIqIjcCZwK7MfGG5bCnwcWAZcC9wbmZO9zHOkeeMuNJwmMPUjbrk\n6rrEoc500jJ1NfBB4KMNyy4BPpOZ74+ItwGXlsvUhDPiSkNlDlNH6pKr6xKHOtd2AHpm/j3w4KzF\n5wCbyvubgBU9jmvBcEZcabjMYepEXXJ1XeJQd+Z7Nt9RmbkLIDPvB47qXUgLizPiSrVkDtMB6pKr\n6xKHutOrAehznu6yruG0s8nJSSYnJ3u02/qr04y4nfTB96Kf3r7+8TM1NcXU1NSww6iiZQ4b5/w1\nTuqSq+sSxzjpRf6abzG1KyKOzsxdEXEM8L25Vl43xufw12VG3E764HvRT29f/3iaXWRcfvnlwwum\nMx3nsHHOX+OkLrm6LnGMk17kr067+aK8zbgeuKC8fz5wXdd7HhN1mBG3kz74XvTT29evGjOHaU51\nyNV1ikPdaVtMRcRfAZ8DfjoivhkRbwDeB5wREXcDLy0fq4k6zIjbSR98p/3027fvYOXKyzjttLWs\nXHkZ27fv6Go/0qCZw9SJbnL1XHlwUHH0MwZ1r203X2a+rsVTL+txLAvW+vVv5YEHdrJjx+UHzIi7\nfv1gJnXrpA++mOR57nXadeHZ1686MoepU53k6kEMZWgXh8Mp6scZ0AdkmDPidt4H33qdvXu/zcaN\nzy678GbMdOFtZsWKnfb1Sxp5c+XqA4cyzDgwD/aqt6FVHIOMQZ3zQsdjoJM++HbrwOK2XXj29Uta\nyOowlKEOMehgFlMLSKs+9E764Nutc+ihP0G7Lrw6jA+TpLlUGWvUq6EMdYhBvWU33wLRrg+9k7EA\nc62zcuVldNKFN+zxYZLUStWxRr0YylCHGNR7FlMLQKd96J2M22q1zpo1Z7Jt25ZZ+yjM7sLziumS\n6qYXY426yYN1jUH9YTffAjCIPnS78CSNsl7kyap5cNAxOH3C4NgytQAMqg/dLjxJo6pXebJKHhxk\nDE6fMFgWUwvAIPvQ7cKTNIp6mSfnmwcHFYPTJwye3XwLgFMSSNLc6pAnBxWD0ycMnsVUDwy7X9rx\nTJI0tzrkyUHF4PQJg2c3X0V16Zd2PJMkza0OeXIQMTh9wuBZTFVQt35pxzNJ0tzqkCf7HYPTJwye\n3XwV2C8tSaqbOnRpjhtbpjqwffsONmy4kT17FrFkyT7WrDmT5ctPsV9akgaoVS7WwerQpTlOLKba\nmGtMlP3SkjQYdRmfOkrq0KU5Luzmm8OBY6JmCqaZMVFHcvrpJwz9VFtJWuja5eJbbtk5zPAki6m5\ntBsTtXXrPfZLz0MnU0m0W2fY01FI6k6Vz6zjU1V3dvPNoZMxUfZLd6eTpvp269jcL42Wqp9Zx6eq\n7iym5tDpmCj7pTvTyVQSmTnnOscfv7FW01FImlsvppBxfKrqzm6+OdTh8gMLSSdN9e3Wee97/8bm\nfmmE9KKLzlysurOYmoNzdfRWJ0317db54Q8Pb7sNSb1VZbxTp110c+3DXFxfjl8t2M3XhmOieqfz\npvrW6zzpSQ+zZ4/N/dKgVB3v1MnnvpN9mIvrx/Gr+1lMdcAxUb3RySUOMnPOdS699NW8851eJkEa\nhF6Md2r3uT/99BN45zuf0NE+zMX1UbfLqQ3bSHfz9aJ50SbKwemkqb7dOm9+8+qeNff7t5fm1ul4\npypddFu33uM4yCGpMk1NN2PhxuF/9ci2TPWiedEmysHrpKm+3Tq9aO73by+118l4p6pddKedtrbt\nPtR7Vaep6XQs3Lj8rx7JYqoXzYs2UQ5PJ0317dap0tzv317qTLvxTnv3fpuNG59dqYvOaQ8GrxfT\n1LzgBV+l3d9tnP5Xj2Q3Xy9OtXVG3fHl317qTLspCWCx0x6MoF5MUwOPtv27jdP/6kotUxHxCuAq\niqJsY2Ze0Wy9lSsvO+jq3p1c/bvVOt3Mhtv5NqaAyabbqJN7753i+OMnhx1GR4Yda9W/fZX3aKfP\nd7ONe+75FieccNzQ45hrG6OmkxzWLH9BfY55u/dHlX0U45128KEPbebRR1cx08WyePEWVq/ezW23\n/QTddtHNzgvt9jGsVodh569OzSfOTv6HZjLnOocddiyrV++e8+/2rnd9iqq5tpupNYaZv+ZdTEXE\nIcB/A14KfAe4NSKuy8yvzF732mvXdn05kLnW6bRZuLttTLH/j1zfpuVR+YDDcGOt+rcf1GVvutvG\nH/GFL7yjBnHUc8xCtzrNYbPzF9TnmLd7f/RiH3ONd1q58jK67aJrlhfqOO3BqOTa+cTZi2lqJib2\ntf279SLX9mJqjUHkryotUycBX8vM+wAiYgtwDnBQMdXt5UDa9dW+5z1HtD3Fvl0/ayfb0Giq+rfv\n5FTtXlz2ZiFtow5jFuahwxxWz2NedRvd/F1bjXfqZLqTTjntweD0Ypqamb/tXH+3XkyLUXUbg8pf\nVcZMPQP4VsPjb5fLWur0ciDt+ki3br2n7enxvdiGRlPVv30np2r34rI3C2kbI6qrHFa3Y151G734\nuzoz+WjqxTQ1nfxtezEtRtVtDCp/RRYdo92/MOI1wC9n5pvKx68HTsrMC2etN78dSBppmRnDjmEu\nneQw85c0nrrNX1W6+f4ZeGbD42PLZZUCkqQBaZvDzF+SOlGlm+9W4NkRsSwingCsAq7vTViS1Hfm\nMEk9Me+WqczcFxG/B9zM/tOK7+pZZJLUR+YwSb0y7zFTkiRJ6uMM6BHxioj4SkR8NSLe1q/99EJE\n3BsRX4iI2yPi88OOZ0ZEbIyIXRHxxYZlSyPi5oi4OyJuioiJYcY4o0WsayPi2xFxW3l7xZBjPDYi\ntkbElyLizoi4sFxeu2PaJNY3l8vrdkwXR8TO8rNzZ0SsLZfX7ph2a1RyWF3zF4xODhuF/FXGNBI5\nbBzzV19apsrJ8L5Kw2R4wKpmE3rWQUR8A/j5zHxw2LE0iohfBB4GPpqZLyyXXQHszsz3lwl+aWZe\nMsw4y7iaxboW+EFmXjnU4EoRcQxwTGbeERGHA/9IMa/QG6jZMZ0j1tdSo2MKEBFPzsxHImIR8Fng\nQuA11OyYdmOUclhd8xeMTg4bhfwFo5PDxjF/9atl6vHJ8DJzLzAzGV5dBTW8TmFm/j0wO0GeA2wq\n728CVgw0qBZaxArFsa2FzLw/M+8o7z8M3EVxBlftjmmLWGfmQKrNMQXIzEfKu4spxmEmNTymXRql\nHFbL/AWjk8NGIX/B6OSwccxf/foAdj2h55Al8OmIuDUifmPYwbRxVGbuguINCxw15Hja+b2IuCMi\nPjzspudGEXE8cCKwAzi6zse0Idad5aJaHdOIOCQibgfuBz6dmbdS82PagVHKYaOUv2C0clitPmuN\nRiWHjUv+quW3mSF4SWb+HPAq4HfLJt9RUeczCP4UeFZmnkjxRq1F027Z7PxJ4C3lt6bZx7A2x7RJ\nrLU7ppn5WGa+iOIb8kkR8XxqfEwXoFHOX1Df90btPmszRiWHjVP+6lcx1dGEnnWRmd8tf34fuIai\nib+udkXE0fB4v/T3hhxPS5n5/dw/KO8vgF8YZjwAEXEoxYd7c2ZeVy6u5TFtFmsdj+mMzHyI4kqm\nr6Cmx7QLI5PDRix/wYi8N+r6WRuVHDZu+atfxdTITIYXEU8uq2ci4inAy4F/Gm5UBwgO7GO+Hrig\nvH8+cN3sFwzRAbGWb8IZv0I9jutHgC9n5gcaltX1mB4Ua92OaUQ8faapPiKeBJxBMT6irse0UyOR\nw0Ygf8Ho5LBRyF8wOjlsrPJX3+aZKk95/AD7J8N7X192VFFEnEDxbS4pBp99rC6xRsRfAZPAkcAu\nYC1wLfDXwHHAfcC5mblnWDHOaBHraRR95Y8B9wK/OdMPPQwR8RJgO3Anxd87gbcDnwc+QY2O6Ryx\nvo56HdMXUAzQPKS8fTwz/yginkbNjmm3RiGH1Tl/wejksFHIXzA6OWwc85eTdkqSJFXgAHRJkqQK\nLKYkSZIqsJiSJEmqwGJKkiSpAospSZKkCiymJEmSKrCYkiRJquD/B2JggN09PYevAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a2ef160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xrange = range(30)\n",
    "fig, ax = plt.subplots(2,2,figsize=(10,7))\n",
    "axes = [ax[0,0],ax[0,1],ax[1,0],ax[1,1]]\n",
    "vis_x = [0,5,10,15]\n",
    "for i in range(4):\n",
    "    hist = df.query('x=='+str(vis_x[i])).groupby(by='y').count().reset_index()\n",
    "    hist = pd.DataFrame({'y':xrange}).merge(hist,on='y',how='left').fillna(0)\n",
    "    axes[i].set_ylim(top=50)\n",
    "    chart = axes[i].plot(hist.y,hist.x, 'bo', ms=8)\n",
    "    chart = axes[i].vlines(hist.y, 0, hist.x, colors='b', lw=5, alpha=0.5)\n",
    "    title = axes[i].set_title('x='+str(vis_x[i])+' (lambda='+str(vis_x[i]+1)+')')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "と、いうことで、xの値からそのx点におけるポアソン分布のλがわかれば、だいたいデータの構造が記述できそうである。\n",
    "\n",
    "**回帰モデル**  \n",
    "さて、回帰である。手元にこのデータがあったとして、どうやら誤差がポアソン分布に従いそうだ、というのもわかった。どう回帰モデルを作ったら良いだろうか。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "上の方に書いた、ポアソン分布の確率質量関数を思い出してみる。<br><br>\n",
       "$$p(X=k) = \\frac{\\lambda^ke^{-\\lambda}}{k!}$$<br><br>\n",
       "少し書き方をかえる。<br><br>\n",
       "$$p(y|\\lambda) = \\frac{\\lambda^ye^{-\\lambda}}{y!}$$<br><br>\n",
       "これはすなわち、あるλが与えられた時、変数がyの値をとる確率を表している。\n",
       "<br><br>\n",
       "ここで、$$\\lambda$$を何らかの$$x$$の関数$$f(x)$$として表せたとしたら、こんな風に書き換えられる。\n",
       "<br><br>\n",
       "$$p(y|x)=\\frac{f(x)^ye^{-f(x)}}{y!}$$\n",
       "<br><br>\n",
       "これを全ての$$x,y$$のペアについて計算して、その確率の積の最も大きくなるような$$f(x)$$\n",
       "を見つけるのがポアソン回帰である。"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "上の方に書いた、ポアソン分布の確率質量関数を思い出してみる。<br><br>\n",
    "$$p(X=k) = \\frac{\\lambda^ke^{-\\lambda}}{k!}$$<br><br>\n",
    "少し書き方をかえる。<br><br>\n",
    "$$p(y|\\lambda) = \\frac{\\lambda^ye^{-\\lambda}}{y!}$$<br><br>\n",
    "これはすなわち、あるλが与えられた時、変数がyの値をとる確率を表している。\n",
    "<br><br>\n",
    "ここで、$$\\lambda$$を何らかの$$x$$の関数$$f(x)$$として表せたとしたら、こんな風に書き換えられる。\n",
    "<br><br>\n",
    "$$p(y|x)=\\frac{f(x)^ye^{-f(x)}}{y!}$$\n",
    "<br><br>\n",
    "これを全ての$$x,y$$のペアについて計算して、その確率の積の最も大きくなるような$$f(x)$$\n",
    "を見つけるのがポアソン回帰である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**一般化線形モデルとしてのポアソン回帰**  \n",
    "さて、データはある、それに誤差構造がポアソン分布であるということもわかっている。ここから分析を実行するために足りない部品は何だろうか。  \n",
    "\n",
    "それは、上の式でいうf(x)の形である。「λを何らかのxの関数f(x)で表せたとしたら」と上で書いたが、このf(x)の式の形をあらかじめ指定しなければならない。（でないと無限の可能性のあるf(x)を探索して最適な関数を探すことになってしまう...）\n",
    "\n",
    "ここで登場するのが一般化線形モデル(GLM)という概念である。[ここ](https://gist.github.com/shotahorii/cc2508accc0b2e6815663b0531b1edb9#file-glm-ipynb)で詳しく書いているが、「λを何らかのxの関数f(x)で表せたとしたら」について、どのように表すかを「線形予測子」「リンク関数」という概念で以下のように定義している。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "注:上の例では説明変数をxただ一つとしていたが、ここでは説明のためM個の説明変数があるとする。<br>\n",
       "まずこの説明変数$$x_1,x_2,...,x_M$$を線形結合して「線形予測子」なるものを作る。<br><br>\n",
       "$$\\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M$$<br><br>\n",
       "この部分が、全ての説明変数を線形に結びつけてパラメータ予測に使う部品である。\n",
       "<br><br>\n",
       "次に、この線形予測子がどのように予測したいパラメータ(ここでは$$\\lambda$$)を表すのか、という部分である。<br>\n",
       "以下の式の$$f(\\lambda)$$の部分である。これを「リンク関数」という。\n",
       "<br><br>\n",
       "$$f(\\lambda) = \\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M$$\n",
       "<br><br>\n",
       "一番単純なのは、$$f(\\lambda)=\\lambda$$とするもので、これを\"identity関数\"と呼んだりする。\n",
       "<br><br>\n",
       "ポアソン回帰の場合は、このリンク関数に対数$$f(\\lambda)=log\\lambda$$を用いることが多い。これは$$\\lambda>0$$\n",
       "であるためで、対数にすることによって$$\\lambda$$は0にならない。\n",
       "<br><br>\n",
       "$$log(\\lambda)=\\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M$$\n",
       "<br><br>\n",
       "$$\\lambda = exp(\\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "注:上の例では説明変数をxただ一つとしていたが、ここでは説明のためM個の説明変数があるとする。<br>\n",
    "まずこの説明変数$$x_1,x_2,...,x_M$$を線形結合して「線形予測子」なるものを作る。<br><br>\n",
    "$$\\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M$$<br><br>\n",
    "この部分が、全ての説明変数を線形に結びつけてパラメータ予測に使う部品である。\n",
    "<br><br>\n",
    "次に、この線形予測子がどのように予測したいパラメータ(ここでは$$\\lambda$$)を表すのか、という部分である。<br>\n",
    "以下の式の$$f(\\lambda)$$の部分である。これを「リンク関数」という。\n",
    "<br><br>\n",
    "$$f(\\lambda) = \\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M$$\n",
    "<br><br>\n",
    "一番単純なのは、$$f(\\lambda)=\\lambda$$とするもので、これを\"identity関数\"と呼んだりする。\n",
    "<br><br>\n",
    "ポアソン回帰の場合は、このリンク関数に対数$$f(\\lambda)=log\\lambda$$を用いることが多い。これは$$\\lambda>0$$\n",
    "であるためで、対数にすることによって$$\\lambda$$は0にならない。\n",
    "<br><br>\n",
    "$$log(\\lambda)=\\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M$$\n",
    "<br><br>\n",
    "$$\\lambda = exp(\\beta_0 + \\beta_1x_1 + ... + \\beta_Mx_M)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**実装**  \n",
    "さて、実際にポアソン回帰してみる。ここでは、真のモデルがλ=1+xであると知っているので、リンク関数をidentityとしてポアソン回帰をしてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "# add constant to the data\n",
    "x_c = sm.add_constant(df.x)\n",
    "# modeling with GLM (log link function is automatically selected for Poisson)\n",
    "model = sm.GLM(df.y, x_c, family=sm.families.Poisson(sm.families.links.identity))\n",
    "result = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>y</td>        <th>  No. Observations:  </th>  <td>  2000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>GLM</td>       <th>  Df Residuals:      </th>  <td>  1998</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>        <td>Poisson</td>     <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>      <td>identity</td>     <th>  Scale:             </th>    <td>1.0</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -4935.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Wed, 27 Mar 2019</td> <th>  Deviance:          </th> <td>  2114.9</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>20:21:50</td>     <th>  Pearson chi2:      </th> <td>2.03e+03</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>         <td>7</td>        <th>                     </th>     <td> </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>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    1.0200</td> <td>    0.070</td> <td>   14.501</td> <td> 0.000</td> <td>    0.882     1.158</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x</th>     <td>    1.0004</td> <td>    0.010</td> <td>   99.139</td> <td> 0.000</td> <td>    0.981     1.020</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   No. Observations:                 2000\n",
       "Model:                            GLM   Df Residuals:                     1998\n",
       "Model Family:                 Poisson   Df Model:                            1\n",
       "Link Function:               identity   Scale:                             1.0\n",
       "Method:                          IRLS   Log-Likelihood:                -4935.7\n",
       "Date:                Wed, 27 Mar 2019   Deviance:                       2114.9\n",
       "Time:                        20:21:50   Pearson chi2:                 2.03e+03\n",
       "No. Iterations:                     7                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [95.0% Conf. Int.]\n",
       "------------------------------------------------------------------------------\n",
       "const          1.0200      0.070     14.501      0.000         0.882     1.158\n",
       "x              1.0004      0.010     99.139      0.000         0.981     1.020\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "λ=1.02+1.0004xなので、ほぼ真のモデルを正しく捉えられている。"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}