GLM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import uniform, normal, poisson, binomial\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12c51e4e0>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcnFWV//HP6eo9aUIgCcSQEMUYCTu2iD8UXIYZQQwzccXREUUDDCgIqOAoSEYdXEBBNIjAICIiAwEiChoWg0RZOhgiBMOOCYR0NjqddHqpqvP741aHTqe6u5an1v6+feWVXp6u51ZHTt0699xzzd0REZHRpabUAxARkeJT8BcRGYUU/EVERiEFfxGRUUjBX0RkFFLwFxEZhRT8RURGIQV/EZFRSMFfRGQUqs33AcysEbgfaEg93s3ufsGgaxqA64C3ABuAj7r7C8M97oQJE3z69On5Dk9EZFRZunTpenefONJ1eQd/oAd4j7tvMbM64AEzu9PdHxxwzUnAJnd/o5l9DPgO8NHhHnT69Om0tbVFMDwRkdHDzF7M5Lq80z4ebEl9Wpf6M7hh0PHAz1Mf3wy818ws33uLiEhuIsn5m1nMzJYB7cAid39o0CVTgFUA7h4HOoDdo7i3iIhkL5Lg7+4Jdz8Y2As4zMz2z+VxzGyumbWZWdu6deuiGJqIiKQRabWPu78K3Ae8b9C3XgKmAphZLTCOsPA7+OevdPdWd2+dOHHE9QoREclR3sHfzCaa2a6pj5uAo4G/D7psIfCp1McfAu51HSQgIlIyUVT7TAZ+bmYxwovJTe5+h5nNA9rcfSFwNfALM3sG2Ah8LIL7iohIjvIO/u6+HDgkzdfPH/BxN/DhfO8lIpKJjq4+2ju7mdTSyLjmulIPpyxFMfMXESkbi1e2M++OFds/P/+4WRw1c1IJR5Slewl1kwXOj6i9g4hUjY6uPubdsYJYjdFcX0usxph3xwo6uvpKPbSRvQp8Dngv8D0gWdjbKfiLSNVo7+wGoKE2tv1v99e+XrYWAvsB1wBfAv5EwaOzgr+IVI1JLY0A9MQT2/82e+3rZac/vXM8YdvrQ8B3gebC31rBX0TKQkdXH0+v7cwrRTOuuY7zj5tFIuls7YmTSDrnHzer/BZ9HfglMAtYAMwD2oDW4g1BC74iUnJRLtIeNXMSC6aOL99qn1XAKcDvgMOBqwgpnyLTzF9ESqoQi7TjmuuYsUdLeQX+JDCfEOj/CPwQeICSBH5Q8BeREiv1Im0U6aYRPQW8G/hP4DDgb8AZQKxwtxyJ0j4iUlIDF2kbamNFXaQt+J6AOHAJcAHhuKurgU8DZdDQXjN/ESmpUi3SFnxPwGPA24CvEFpdrgA+Q1kEftDMX0TKQCEXaYdq9ZAu3bS1J057Z3d+9+8BvglcBOwG3ERoZ1kmQb+fgr+IlIVxzXWRz/aHS+sUJN30Z+CzwJPAfxBSPmV6bJXSPiJSlUZK6/Snm3rjzqtdvfTG80g3bSEs4L4D2ArcSTi4tkwDPyj4i0iVyryKyAf9naVFwAHAZYRqnsfZ+TirMqTgLyJVaaRWD/3vDOpra9i1uZ762prsFnw3ERZw/xmoB+4HLgdaIn4iBaLgLyJVaaQqorz2FywgtGa4DjiXUNnzzkI8i8LRgq+IVK3BVUQAT6/tZFJLY24Lvq8ApwO3AAcDvwUOLfCTKBAFfxEpmWKcuNVfRZSu8uf842Yx744VbO2JY8bQC75OmOV/EegCvg2cA5RR94hsWbmeo97a2uptbW2lHoaIFEgxT9zq6OpjzvwlxGps+yw/kXQWnHoEwPAvQC8AJwN/AP4fYZfumwsyzEiY2VJ3H7E/qHL+IlI0/X10Vm3oKuqJW8Pl94dsApcEfgTsDyxJffwnyjrwZ0NpHxEpioEz/d54kr5Ekomp/Hpku2uHkHV+/++EzVpLgH8BfgrsHfmwSirvmb+ZTTWz+8xshZk9YWZnpLnmXWbWYWbLUn/Oz/e+IlI5Bm+4qq+tYcPWXrp640DhT9warvJnh66efYR8/kGEXjzXEjZsVVngh2hm/nHgbHd/1MxagKVmtsjdVwy67k/uflwE9xORCjM47dJcX8vuYxrojSdxH2GxNSLp+gcNfDcyY/VYvvf7g2h5si704vkRsGfBhlNyeQd/d18DrEl93GlmTwJTCK+bIiJp0y5jGmJce+JhdMcTRTtxa2D/oP53I03xGJ+6ezpz7t2LV8f0svWGOGNOqP6MeKQLvmY2HTiEcAzxYG83s8fM7E4zK9HZNSJSCkOlXabu3lyyE7faO7s54LlxXHHJW/jw3VO5963tfPqcR3j5PduKPpZSiOzlzczGErY+nOnumwd9+1Fgb3ffYmbHArcBM9I8xlxgLsC0adOiGpqIlIGyOlu3E6ad38wPrzyEV3bbxtf/83EeeuMGEkkvyiEy5SCSOn8zqwPuAH7v7pdkcP0LQKu7rx/qGtX5i0hB3Emo218Nqz/Zxckzl7KtPrF93aFQew2KJdM6/7xn/mZmhG0PTw4V+M1sT2Ctu7uZHUZIN23I994iIhnbQNih+wtgX2AJ7PX2Zm7oOnyndyPF2HlcalGkfY4APgn8zcyWpb72VWAagLtfQVg7P9XM4sA24GNerluLRaS6OHAzoSfPRuBrqT8N4duDD5Ep5s7jUoqi2ucBRjigzN0vJzQ7FREpnpeB0wirjG8htGg4aOjLB+5H6K9KmnfHChZMHV917wDU3kFECmqHTVTF+lknJKNnAXcB3wUeZNjAD3m2ea4w1V/MKiIlk08KJeeffY5QM3gPcCTwM+BNmY23IOf6linN/EUqSD6z6GIb6QzdyH82AfyQcKTiw8B84D4yDvww8gEw1UQzf5EKUWkLkelSKJk2b2vv7CaRdGrMSCSd2poatvX28ey6LRy69/idf+AJ4CTC9tJjgSuAqbmNu6z2IxSQZv4iFSCfWXSpjHSG7nCeXtvJS69u44UNW3lqbSdPrd3MK5u7OeumZSxe2f7ahb3AfxP6CjwDXE/YcZRj4O83ZJvnKqLgL1IBKnEhMtcUSkdXHxcveoqJY+sxjHjSiSdhz10adzxk/RGgFTgf+CChm9i/M0LtofRT2kekAlTqQmQuKZT+F7TdxzbSWFfLqo1dYNBYF6OhNka8M0nyHA899vcEbgdmF/RpVCXN/EUqQCUvRGabQhn4QtdYF9s+k6+L1TDz72P52Q9aGT+/PuT4V6DAnyPN/EUqxGhZiOx/oZt3xwriiQS7j6lnbE8tp9y4D7MffB3bpiZCGed7Sj3SyqbgL1JBBrciqFYDX+he90ATTV+IYa9Az+cTNF0Ug+ZSj7DyKfiLSE4K3fxs3NY6xn2xDn4F7AcsgIa3xSK/z2il4C8iWSvongMHbgS+AHQA3wDOA+qjeXgJtOArIlkp6J6D1YQF3I8DbyAcA3UBCvwFoOAvIlkpyJ6DJKF0cxZhMfdi4M/A/jtfWkktLsqZ0j4ikpXI9xw8A3wO+CPwbkIjtn3SXxpFumk0HNSSCc38RSQrke05iAPfJzRie5QQ9O9hyMAfRbpp8cp25sxfwinXL2XO/CU7tooYZTTzF5Gs5b3nYDlhk1Yb8AFCB84pw//I4HRTbU0NXT19PLuuk0P33m3EW2Z6UMtoeWeg4C8iOclpz0EP8O3Un/GEqp6PkFE/noHppp6+JC93bAPgrJuWc+HskdM/mXQZrbTOqflQ2kdEiuNBSBziMA96P5QMrRk+SsaN2PrTTb3x1wL/lF2bqK/NLP0zUpfRSuycmg8Ff5EqVTZVMVuBs8D/n7Ph5R6++pnlHNN6P4s3ZJ9vP2rmJC75yMHsuUsjMya10NJYN2y10cDfwUhrFZXYOTUfSvuIlJkocs5lk764h1DJ8zwsfPvL/O8HnicxFuiFr932ODd89nB2aarL6vnuM3Es9bU1xJNJYjVDVxsN9TsYaq2iUjun5srcvdRjSKu1tdXb2tpKPQyRosokaI/04tDR1cec+Ut2WNhMJJ0Fpx6R1YtJXi9CrwLnEA5RnwGrvtvFic8+THN9La929fLK5m7cnV2b6qirjdFcHxvy+abT/3tyB7Odfy7X38FIj1sJzGypu7eOdF3eM38zmwpcB+xB2Jh9pbtfOugaAy4lHLDWBZzo7o/me2+RapJJNUomLw75HJ/YL693DrcB/wmsBb4MfAN28TqYD+s6u1m7uYf+KefGrX3EYnFmTGohnkymrb5JZ6Rqo2x/B/0vdAdPHc+CU48YFdU+UeT848DZ7j4LOBw4zcxmDbrmGGBG6s9cQmGXSNXJJ88+Us450wXJfI5PzOY+O1lLqNz5N2Ai4Tzd7wBNYbH27KPfxPotvThhjbe2JmzsBehLJLPOsQ93TsCklkYSSdi8rY9E0of9HQyu/V+2alPVH+EIEQR/d1/TP4t3907gSXau2D0euM6DB4FdzWxyvvcWKSf5biAaKWhnuiCZ7yasrBc+HfgFoTXD7cA3CfX7gxIPM/Zo4XW7NlEXM+pra6iNhfDj7tTFaiLNsS9btYnuvjgvd2zj6fZOOrvjaX8Ho63CZ6BIF3zNbDrhKOWHBn1rCrBqwOerU19bM+jn5xLeGTBt2rQohyZSUJluIBrOwENMtvbEt+ec+38+mwXJfDZhZbXw+Q/gFOBO4O2EHP++Qz9uXcyY1NLA+i29JJJOrMbYfWw93X2JnZ5vrvr/LXZpqmO3MQ1s7Y1jwMFTx+90bRQpskoVWfA3s7HALcCZ7r45l8dw9yuBKyEs+EY1NpFCiyqIDBe0R3pxGCzXg18yuk8SuAL4SurjS4HTgAHt9gcvGA983MnjGkkknfOO2Zcj3jgh0hz74H+LXRrrhvy3GG0VPgNFEvzNrI4Q+H/p7gvSXPISMHXA53ulviZSFaIMIsMF7WId5TjsfVYCnwUeAI4mdON8/Y4/P9SC8VCPG+XzyObfItsX1GqSd6lnqpLn58BGdz9ziGveD5xOqPZ5G3CZux823OOq1FMqTTWUCQ6rvxHbN4Am4AfAp9hph242ZZaF6qOT7b9FNfXzKVqpJ3AE8Engb2a2LPW1rwLTANz9CuB3hMD/DKHU89MR3FekrFT1AevLCI3YHiVU8/wYGKJkI9MUWCE3omX7bzFazkYeKO/g7+4PMEJ3Dg9vL07L914i5a7qgkg38N+Eks0JwM3AB4f/kUzSLlEskI+k6v4tIqbePiKS3hLgYEIHzk8QGrGNEPihv6Z/Jtt6E3R296UtNR1tfXTKkYK/iOxoC+Hw9HcC24C7gGuBkVvmAyGdc/GilcRqjETSOfvoN+2Uzsl3I5rkT8FfRF7ze8K5uZcTErWPA/+S+Y8PTOe0NNbRVB/j4kVP7bRpKrLTwCRn6uopUsUyrmLZCJxFqNubCfyJUMqRpWz2O1T1AnkFUPAXqVIZV9PcQpjlryfU6X0dyDH7ku1+By3Klo7SPiJVKKOeNWsIC7gfAl5H6MfzLXIO/KB0TiXRzF+kCg2bfmmqCwu4ZxEWdC9KfRxRfFY6pzIo+ItUocHpl67eOH2JJGNerg1bMhcB7wCuIuT4I6Z0TvlT2kekQmRzVsDA9Mu6zh5Wr+/i3/40hfGH1xNfkgw7dBdTkMAvlUEzf5EKkEsrhKNmTuLaCWM579vLOe+WVvZ7cRyPvHkjl8xZyRUntjKuRjPz0UzBX6TMDC7PzLkVQh80fT/GNZceRndDgks+sZL73rqOrb2jo1+9DE/BX6SMpJvhv27XJiDLswKWAp+BCcsbuO+gdubPeYZt45PaSSvbKecvUgLp8vdDlWc2poJ+Rq0QthEOWHkb0A7cCjW/hg1jelV6KTvQzF+kyIbK3w9VntkdT2R24Mj9hENWnia0X/4eMB6OQqWXsjMFf5ECyiZ/P9zu2Bl7tAwdwDcD5wLzCSdq3Q28d8dxqPRSBlPwFymQbPP3M/ZoGXaGnzaA/w6SJzv2Emw5Oc7aL3czcc9GxkW1Y0uqloK/SAEMNcO/9sRweulQvW8y3h27HvgicD38Y48uvnnyChbvto7db6xnTEPtDqWg1XREIVTf8ykVBX+RAsgnfz9sisaBm4DPg29yfnH0i9zw3hd5vrMLHDZ29TGmoXZ7KmnZqk0FOyqxFAp59ONoo+AvkoORZp855++H8zJwKrAQaIV//LqLnz/6AjUWTlGtrakhkXRqLByi8uy6zoIflVhMxTj6cTRRqadUpGxaHURt8cp25sxfwinXL2XO/CUsXtm+0zUjdbcc11zHjD1aMgtaTujBMwv4A6GK5y+w69vqAUi6AxBPJsHC5+H1ILwoVMtRiTr6MVqa+UvFKeVb/2xmn5F0t3wW+BxwH3AU4UXgjeFb42rrtqeQxjfXs2FrD+Obwj3OP24W+0wcC2TeW7/cZXtWgAwvkuBvZtcAxwHt7r5/mu+/C7gdeD71pQXuPi+Ke8voUuq3/tmcVAV5lFgmgEuBrxH+K/0poYZ/0Hv1gS8wjbUxuuOJHV5oMtofUCH6301Vy/Mptahm/tcSTv28bphr/uTux0V0Pxmlsg2+USvK7PNxwiath4H3A1cAew19+XAvMNXWW7/ank8pRZLzd/f7CaeAihTUwODb/3cx3/oX9KSqXuBC4FDgOeAG4DcMG/gzkdX6QgWotudTKsXM+b/dzB4j1Cyc4+5PFPHeUiXK4a1/QWafDxNm+48DHwd+CEzM/2FBdfGSXrGC/6PA3u6+xcyOBW4DZgy+yMzmAnMBpk2bVqShSaUph7f+kbVL6ALOB34ATCaUcX5g6MuzDeSqi5ehmKfKxPJ+ILPpwB3pFnzTXPsC0Oru64e6prW11dva2iIZm0gh5D2jvo+wiPsccDLwHWDc0JdnG8g7uvqYM3/JDovjiaSz4NQj9A6gipnZUndvHem6osz8zWxPYK27u5kdRlhr2FCMe4sUQl4z6g7gS8DPgH0ILwLvGuFHcqhyKvXiuJS3SBZ8zexXwF+AmWa22sxOMrNTzOyU1CUfAh5P5fwvAz7mUb3lECmyofruZ7Th7DeEzVpXA+cAyxkx8ENuG5xKvTgu5S2Smb+7nzDC9y8nlIKKVLycZtTtwBnAjcABhFWvt2Z+z1xKTMthcVzKl3b4imQpq0DshJLNMwh99+cRTtqqz27NINdAXg6L41KeFPxFspRxIF5FaMT2W8KxilcD+4Vv5bJmkGsg10Euko6Cv0gOhg3ESeBK4MuENg0/AD4PhCxR1ou3g98hKJBLFBT8RXKUNhA/TSjfvJ9wlOKVwBt2vCSbNQPV6UuhqKWzjCoFawUdB74LHAg8Rui+uYidAj9kXoWTV1WRyAg085eiK1W7gYLNoh8jtGZYChwP/AR43dCXj2uu4+yj38T/3Pl3evqS1MYs7ZqB6vSlkBT8pahKlcYoSCvoHuCbwEXAboTjFT9E/xkqQ1q8sp2LFz1FrCacuHX20W9O+ztQ/3opJKV9pGhKmcaI/BSovwCHEIL/CcAK4MOMGPgH/g5aGutoqo9x8aKVaX8HBe0gKqOeZv5SNKVMY0Q2i95COGDlMkKr5d8Bx2T+49n+DlSnL4Wimb8UTSnbDUQyi15E2J17KfCfwBNkFfght9+B+tdLIWjmL0VT6nYDOc+iNxH68FxDaER+P/DO3MZQ6t+BSL/IWjpHTS2dq1exq33yut+thFn+OkInzvOBpvzvqQNWpFDKqqWzyEDF3KWac3XRK4RduTcDBxNaNBwa3T21U1dKTTl/qVq5VBd1bO3jlcu2kZzl4VStbxGOWMww8GtjllQKBX+pWtmWdz5473pWvqWTPc9oYsW4zTx82wb4KpDFBD3yklKRAlHaR4DqzEFnWt7ZsaWPzZf0cdC3dgWMKz70LLcevpr4CmfBUdkdeaiNWVIpNPMXFq9sZ878JZxy/VLmzF/C4pXtpR5SJDIp73z4zg3844Aupl7QTNuUTZz4xYf57ZFrqK/PbcaujVlSKTTzH+UK0vagjAxZ3tkH3d9OcPA3x9Ndl+Dij6/kir2fBYMZyRbiyWTOM3ZtzJJKoOA/ypVj87CoU1A7Vdb8FfgMNC6LsfjAdq76yPO8uksfU7qbeOnVbWze1kdDXU1eM3ZV80i5U/Af5colR90f8B9b9SqX3vM0sRojVmPRNn7rBi4EvgdMgK03xJn3cupdDzHqa2uYsmszl3zkQPaZqB21Ut2U8x/lyiFH3b/m8ImrHuKcm5fz0qvbWNPRTXdfIrIyyS13x+ndPxk6cP4H8CSMOaF2p+d+4exZHLr3bgV9/gU7U0AkC5HM/M3sGuA4oN3d90/zfSN0RDkW6AJOdPdHo7i35K+UOer+NQeAjV29ALiDu7N+Sy+TxzXml4LqhJdO6WLKDc2sGb+NH3zuKd5/9mSOGh/eTRT7uetkLikXUaV9rgUuB64b4vvHELqizCAcZT0/9beUiVLlqPvXHGrMCP8L7UbMjKQ7iaTnnoK6C5Jzncmrm7j1nau5YfY/6Ij1sfSOTTssaBf6ufentBprY1W9uC6VJZLg7+73m9n0YS45HrjOQyOhB81sVzOb7O5rori/VK7+wJ50xwxiNRBPQiKZxMw475h9sw+MG4CzgOsgPiPJ2ac9xgtv6gKggeIuaA+c6ffGk8STzoSxDWEsZbC4LqNXsXL+U4BVAz5fnfqajHL9aw4A45vrMTMmjq1n8rgmLv3owRx74OTMH8yB/wNmATcAX4NtDyZZMX1zTm2kO7r6ePTFTTz64sac8vODWz3U19awfksPXb3xrMciErWyqvYxs7nAXIBp06aVeDRSLAPz7o21Mbrjiezz72sI3TdvA94C/AE4CMaRWwvlxSvbOXfB31i/pQeA3cfU850PHphVfn5wGW1zfS27NdfTG0/irnbOUlrFCv4vAVMHfL5X6ms7cPcrgSshtHQuztCkHOScd3fgfwlpnh7gO6mPB/w/O9tF3Y6uPi5YuIL1W3qosXAu44atvVyw8AluPy3z/HxjbYzeeBKI01xfS088wdjGWq498bDcXuBEIlSs4L8QON3MbiQs9HYo3y95ew44GbibcLjKVcCb0l+azYtLe2c3iWSSGrPtwd8d4gnPOD/fn+vvSyR5ZXM3u49pYExDjPOPm8XU3ZszGodIIUVV6vkr4F3ABDNbDVxAqheiu19BOOn0WOAZQqnnp6O4r5RWyZrBJYAfAf8FxICfEF4EsljBGm7sk1oaidXUkBxw0JHj1MYs47WC/lz/xJZGxjTU0htPcu2JhynwS9mIqtrnhBG+78BpUdxLykM29eqRvkisAE4CHiQUEP+UHROKGRhp7OOa67hw9qydcv4Xzt4vo/Gny/W7x+lOLTqLlIOyWvCVypBNM7jINjX1EvL53wRagOuBjwNWmLEfNXMSd51xJM+u2wJ4Vu0eyqVlhshw1N5BspbpgSWRnWrVBryVcH7uHMLs/98ZMvAP1z4hm8NWxjXXceje47Nu91AOLTNERqKZv2Qt05lt3h1Du4BvABcDewK3A7OH/5GR3mkUa1auts5S7jTzl6xlOrMdGGj7/8440C4GDiJ04DwJeIIRA38m7zSKOSsf11zHjD3UHVTKk2b+kpNMZrb9gTarDVYdwFcIC7lvAO4B3pPZmDJ9p6FZuYiCv+Qhk9r5rALtbwklm2sIG7XmAWMyrxbKJqWjw1ZktFPwl4IbMdCuA84k9OPZD7iF7T1fs6kWyumdhsgopeAvBZHRbN2BXwOfJ6R7LgC+CtS/9hjZtkBWSkckMwr+krfBgT6j2fpqQiO23xDKOK8GDtjxklyrhZTSERmZgn+VKXbLhcGB/uyjZ3LxopVDz9aThB48XwL6CGWcZxDaNAyizVIihaPgX0WKfURgurTM/9z55PbPu/sS1MVqiCcSYbb+ch18Dvgj8G7gZ8A+Qz/+4Bx+0p25R76hYM9HZDRRnX+VSFfjfsHCJ3j0xU0FOyg8XVomVmNs6e7j6fZOXtzQxdPtnfRsS7DXNc1wIPAooWn3PQwb+PsdNXMSC049gs++8/WAcfUDzzNn/hJ+t3yNDkEXyYNm/lVicCDu6Uvycsc2zrjxr9TX1hTkXUC6tAxALJaaUxi8aW0L37vrQJpWx+ADhNObczjD7eoHnqe+Nryj2Li1hzN+/Vem7NpErMZ0CLpIDjTzrxIDA3Ei6bzcsQ2Alsa63HvqjCDdbtm5R76BsQ217LvbLnz94X25/ZojmLypiTVXbAvtGXII/ANf2BJJZ/2WXgDqYjUFe24i1U4z/yoxMD/e1RMCYf/MOFYzfJVMPovEg0srAZb8cj1f+r83M/2VMdxz6FouP/4Zrv3kYVl34Ow38IXNvf+wd9se/HUIukj2FPyrSH8gfnZdJ2fdtJz62hBt01XJ9Af8p9du4eJFK7d/PZcUyvbSyq3A1+Hyyw9l3S49nPeZ5Tw8a2POG60Gvij1v7AlkuGAlQljwjsaVQCJ5Mbcy/Oo3NbWVm9rayv1MCpWf+WPO9t3uvYH9f7vJZLOS69uY1JLA7uNadieMlpw6hHZB+t7CJU8zwOnQsf5fbRb7iWn6SqXDk69w3h6bScXL3oq7XMTGe3MbKm7t454nYJ/9UqXzuno6mPO/CXEagx3eGHDVmrM2Gfi2O0plO9/+EBaGusyC9yvEmr2rwLemPr7qPzH3T/G/oXkwS9K+aSqSnb8pEgRZBr8lfapYul2ug5ePDUzku70JZLEk9DVG+eMGx8DnFiNceHs/YaeVd8OnAqsBb5M6L3flP+4M9nZm+su3mLvhRApV6r2GWUGLp7GaowJY0IA7Usk6Y0n6YkneWXzNtZu7uGlV7fxlVuW71xJsxb4KPCvwETgIcIRixEE/sFj7P87irx+ZCeLiVQBBf9RZnB5ZlN9LZd+9GB+9h+tzDt+f7b0xKkxI1Zj1JixYWsvz67rDD/shLNzZwG3Af9NOGJxxDeY+Y0xqgNXsjnCUaTaKe0zCg3V+bKzOz7ETxj8AzgFuBM4nNCIbVbxx5gP9QoSeU0kM38ze5+ZrTSzZ8zs3DTfP9HM1pnZstSfz0ZxX8lduiMG95k4lgljG0i6k3An6c7E5gb2vXWX0Gd/MXAp8AAFDfzDjTHfx9PB6iJB3jN/M4sBPwaOJjTqfcTMFrr7ikGX/trdT8/+3AmvAAAMcklEQVT3flI445rruGjOAVywcAWJZJK9N4zh8tsOoaktBv9E6Mnz+lKPMj/q9y8SRJH2OQx4xt2fAzCzG4HjgcHBXyrAUTMncfvJ4+n7bpLdL6vHmgyuAU4k5x265Ub9/kWiSftMAVYN+Hw16Tu4fNDMlpvZzWY2Nd0DmdlcM2szs7Z169ZFMDTJ2jIY9+46JnyrATvWwkv4p8k68Hd09anrpkgZK9aC72+AX7l7j5mdDPwceM/gi9z9SkJygdbW1vLcfVatugnVO98BJgA3Ax/M7aFUSy9S/qKY+b8EDJzJ75X62nbuvsHde1KfXgW8JYL7SlT+DBwCfBv4BGG2n2PgVy29SGWIIvg/Aswws9ebWT3wMWDhwAvMbPKAT2cDT0ZwX8nXFuALwDuALuAu4Fpgt6F/ZKR0jmrpRSpD3mkfd4+b2enA7wknsV7j7k+Y2Tygzd0XAl8ws9lAHNhIWD6UPOXVo+YPwFxC/f5phFl/y/A/kkk6R7X0IpVBjd0qVM559Y3A2YQZ/kxCEu4dI/9YJs3WBo9NXTdFik+N3apYuoPT592xggVTxw//DuAWwix/PXAecD6QZkKe7h1FJs3W+qmWXqT8KfhXmI6uPh55YQOJpNNcH/75hgvEALwCnE4I/ocQWjQckv7xh3pHkW06Z3Atvdooi5QXBf8KsvMhLMnth7CkDcROKKr9IrAN+B9CymeI2DvSO4r+07S29sS3p3MyCeQq/RQpPwr+FWJwCeXEsQnaO3u2n2O7UyB+gbCgu4iQ07+KkOMfxkipnVzSOTmnqESkoBT8K8TgwLz72EbqYjG+9v59eev03V8LpAngJ4ScvgGXEw5cyaCoN5PUTratEbJZKxCR4lE//wqR7oCT2pjtGPifBI4k1O6/E3icsMCb4b9yIbpeFupgFhHJj2b+FWLYnHsf8F1gHjAWuI6wUzeHRmxRV+rks1YgIoWjOv8Ks1PVzKPAZ4DHgA8DPwL2KOkQ01K1j0hxqM6/Sm3PuW8DzgW+TzhHdwHwbyUd2rDURlmkvCj4V6L7gc8CTwMnAd8Dxpd0RCJSYbTgW0k2ExZwjyLk+RcRSjgV+EUkSwr+leJ3wP7AfOBMQiXPP+X3kDpwRWT0Utqn3K0n7NC9nnBo+p+Bw/N/WO26FRndNPMvVw7cRAj4NwJfJ1T2RBD4deCKiCj4l6OXCZU7HwWmAUsJNfwN0Tx8OR+4olSUSHEo7VNOHLgaOAfoIVTxnEnk/0rleuCKUlEixaOZf7l4jrCA+zngYOBvhBeBArw8F6KNQ76UihIpLs38Sy0BXAb8F+Ff4wrCC0CBX5bL7cAVNYATKS4F/1J6nLBZ6yHg/YTAv1fxbl9Ou27LNRUlUq2U9imFXuBC4FDgWeAG4DcUNfCXm3JMRYlUM838i+0RQiO2x4ETgEsJvXmk7FJRItUskpm/mb3PzFaa2TNmdm6a7zeY2a9T33/IzKZHcd+K0kVYwD0c2AQsJMz4Ffh3MK65jhl7tCjwixRY3sHfzGLAj4FjCFuSTjCzWYMuOwnY5O5vBH4AfCff+1aU+4ADgYsJi7lPAB+I7uFVGy8i2Yoi7XMY8Iy7PwdgZjcCxwMrBlxzPPCN1Mc3A5ebmXm5HiYQlQ7gy8CVwD7AvcC7o72FauNFJBdRpH2mAKsGfL469bW017h7nBAWd4/g3uXrN4T3QVcR0j3LiTzwqzZeRHJVVtU+ZjbXzNrMrG3dunWlHk5u1gEfB2YTXt4eJOzUbY7+VuXcpkFEylsUwf8lYOqAz/dKfS3tNWZWC4wDNgx+IHe/0t1b3b114sQKWwl1wgLuvoTE1oVAG/DWwt1Sh6OLSK6iCP6PADPM7PVmVg98jFDLMtBC4FOpjz8E3FtV+f5VhAXcfwfeCPwVOB+oL+xt+2vje+NJNm3tpTeeVG28iGQk7wVfd4+b2enA74EYcI27P2Fm84A2d19IaFf2CzN7BthIeIGofEngZ8CXgDhwCfAFwm+hqAwzB6zYNxaRCmXlOgFvbW31tra2Ug9jaE8TyjYXA+8hvAi8obhD6OjqY878JcRqbHtLhETSWXDqEZr9i4xSZrbU3VtHuq6sFnwrQpywgHsgsIxQzXM3RQ/8oAVfEcmd2jtkYzlhu1obYefCT4DXlW44aoYmIrnSzD8TPYQF3LcALwK/Bm6lJIF/4G5eNUMTkVxp5j+SvxBm+08CnyQ0pyjR9rShdvOqGZqIZEsz/6FsJRyheASwBfgdcB0lC/zD7eZVMzQRyZaCfzp3A/sT2i2fSmi/fExJR6TFXRGJlIL/QJsIKZ6jgTpCGeePgV1KOahAu3lFJEoK/v1uJTRi+zlwLvAYcGRJR7QDLe6KSJS04LsW+Dzwf8BBwB2Eqp4yNNTibkdXH8+u6wSMfSaO1QuCiIxo9AZ/B35BWNTdCnyL0KahzOPm4EPXF69s5yu3LGfD1l4AJoxt4KI5B6inv4gMa3SmfV4kLOB+itCF8zHgq5R94B+so6uPCxY+wYatvdSYUWPG+i09XLBQPf1FZHijK/gnCQu4+wMPAJcBfwLeXMpB5a69s5tE0jFse/A3MxLJpKqARGRYoyftsxL4LCHo/zPwU2B6KQeUv0ktjcRqDMdJpvrzuTuxmhpVAYnIsKp/5t8HXERYzH0CuBa4i4oP/BDy/xfO3o/dx9STdCfpzoSxDVw4W1VAIjK86p75/5VQt/9X4IPA5cCeJR1R5I6aOYnfn3mUqn1EJCvVGfy7gXnAd4EJhGMVP1jSERXUuOY6Dt17t1IPQ0QqSPUF/+cJlTwrgU8D3wcUF0VEdlB9wX8K4RzdywgLuyIispPqC/71hF26IiIypOqv9hERkZ0o+IuIjEIK/iIio1Bewd/MdjOzRWb2dOrv8UNclzCzZak/C/O5p4iI5C/fmf+5wD3uPgO4J/V5Otvc/eDUn9l53lNERPKUb/A/nnD8Cam//zXPxxMRkSLIN/jv4e5rUh+/AuwxxHWNZtZmZg+a2ZAvEGY2N3Vd27p16/IcmoiIDGXEOn8zu5v0HXH+a+An7u5m5kM8zN7u/pKZvQG418z+5u7PDr7I3a8ErgRobW0d6rFERCRP5p57jDWzlcC73H2NmU0G/ujuM0f4mWuBO9z95hGuW0c4diVXE4D1efx8OdNzq1zV/Pz03MrD3u4+caSL8t3hu5BwHtZFqb9vH3xBqgKoy917zGwCcASh5dqwMhn8cMyszd1b83mMcqXnVrmq+fnpuVWWfHP+FwFHm9nTwD+lPsfMWs3sqtQ1+wJtZvYYcB9wkbuvyPO+IiKSh7xm/u6+AXhvmq+3Ec7Nwt3/DByQz31ERCRa1bzD98pSD6CA9NwqVzU/Pz23CpLXgq+IiFSmap75i4jIEKou+JvZ+8xspZk9Y2ZDtZuoSGZ2jZm1m9njpR5L1MxsqpndZ2YrzOwJMzuj1GOKipk1mtnDZvZY6rldWOoxRc3MYmb2VzOrutM0zOwFM/tbqjdZW6nHE5WqSvuYWQx4CjgaWA08ApxQLdVFZnYksAW4zt33L/V4opTaJzLZ3R81sxZgKfCv1fBvZ2YGjHH3LWZWBzwAnOHuD5Z4aJExs7OAVmAXdz+u1OOJkpm9ALS6e6XU+Wek2mb+hwHPuPtz7t4L3EjoP1QV3P1+YGOpx1EI7r7G3R9NfdwJPEk4lLPiebAl9Wld6k/VzLrMbC/g/cBVI10r5aPagv8UYNWAz1dTJQFkNDGz6cAhwEOlHUl0UmmRZUA7sMjdq+a5AT8EvgwkSz2QAnHgD2a21MzmlnowUam24C8VzszGArcAZ7r75lKPJyrunnD3g4G9gMPMrCrSdmZ2HNDu7ktLPZYCeoe7HwocA5yWSr9WvGoL/i8BUwd8vlfqa1IBUvnwW4BfuvuCUo+nENz9VcJO9/eVeiwROQKYncqL3wi8x8yuL+2QouXuL6X+bgduJaSXK161Bf9HgBlm9nozqwc+Rug/JGUutSh6NfCku19S6vFEycwmmtmuqY+bCAUJfy/tqKLh7ue5+17uPp3w39u97v6JEg8rMmY2JlWAgJmNAf4ZqIpqu6oK/u4eB04Hfk9YMLzJ3Z8o7aiiY2a/Av4CzDSz1WZ2UqnHFKEjgE8SZo79R34eW+pBRWQycJ+ZLSdMUBa5e9WVRFapPYAHUr3JHgZ+6+53lXhMkaiqUk8REclMVc38RUQkMwr+IiKjkIK/iMgopOAvIjIKKfiLiIxCCv4iIqOQgr+IyCik4C8iMgr9f3fQib3/JiW5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate simulation data\n",
    "np.random.seed(5)\n",
    "n_sample = 100\n",
    "a = 0.6\n",
    "b = -0.4\n",
    "sd = 0.3\n",
    "\n",
    "x = uniform(1, 5, size=n_sample)\n",
    "mu = a * x + b\n",
    "y = normal(mu, sd)\n",
    "\n",
    "# linear regression\n",
    "slope, intercept, r_value, p_value, std_err  = stats.linregress(x, y)\n",
    "xvals = np.array([0, 5.5])\n",
    "yvals = slope * xvals + intercept\n",
    "\n",
    "plt.scatter(x, y, s=20, alpha=0.8)\n",
    "plt.plot(xvals, yvals, color='magenta')\n",
    "# plt.savefig('linear_regression.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Poission regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzwAAADSCAYAAACPUMWBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFy9JREFUeJzt3XGQXeV53/HvryLgjIlTHHY8rYQtmciZyHUC7Ua0TUrdKdgidBCdQhGuM2SGDCVFtVu3M5YbCh6ltMSZOu5M1RgSq7hpsUqMk25jMZQ6dhrHxdFiEzsSIyMUbKRSs0bU2DEGC57+sYf2slpp7+7e3Xvvq+9n5s6e8573nH1e3d1H99n33PemqpAkSZKkFv2ZYQcgSZIkSSvFgkeSJElSsyx4JEmSJDXLgkeSJElSsyx4JEmSJDXLgkeSJElSsyx4JEmSJDXLgkd9S/JIkiNJ3jzAa25PMp3k+SR3Deq6kkbbCuWTzyT5bpJvd4+Dg7q2pNGw2q9Fkrw2yW8l+dMkX03yjkF9X60eCx4txl8AvgJcNcBr/i/gXwC7B3hNSaNvJfIJwPaqOrt7/MiAry1p+Fb7tcgu4AXgdcDfA351kMWWVocFj/pWVS8CnwV+bIDX/ERV/Tbw9KCuKWn0rUQ+kdS+1XwtkuTVwN8B/nlVfbuqPgtMAT8zqO+t1WHBo74l+X7gWuDHT3L8d5L8n5M8fmd1o5U0ylYwn/yrJN9I8gdJ3roSsUsanlV+LfIm4HhVfaWn7Y8AZ3jGzBnDDkBj5TbgCPDXk5xdVd/uPVhVf2s4YUkaQyuRT94LHGD29pNtwH9NckFVPbbsaCWNitV8LXI28Oyctm8CPzDA76FV4AyP+pLkrwBXMzu1+03gLcONSNK4Wql8UlWfr6pvVdXzVfVR4A+Anx7EtSUN3xBei3wbeM2cttcA31rh76sBs+DRgpK8Cvj3wI1VdYzZ6dwT7p1Ncl/P6khzH/etdtySRs8q55MCMrjoJQ3LkF6LfAU4I8nGnrYfB/YvdRwaDm9pUz92Ap+rqk92+w8zz72zVXXZYi+c5Axmfw7XAGu6hHa8qo4vI15Jo2tF8kmSPwtcBPwecBy4BrgYePeyopU0KobxWuRPk3wC2Jnk54ALgK3AX13iGDQkzvDolJJsZnb6+B/3ND/M4FZHuRl4DtgBvLPbvnlA15Y0QlY4n3wfs8vKzgDfAP4hcOWcNxtLGkNDfi3yD4DvB54CPgb8fFU5wzNmUlXDjkGSJEmSVoQzPJIkSZKaZcEjSZIkqVkWPJIkSZKaZcEjSZIkqVkWPJIkSZKaNXKfw3PuuefW+vXrhx2GpB4PPfTQN6pqYthxLJb5RBo95hNJg9JvPhm5gmf9+vVMT08POwxJPZJ8ddgxLIX5RBo95hNJg9JvPvGWNkmSJEnNsuCRJEmS1CwLHkmSJEnNsuCRJEmS1CwLHkmSJEnNGrlV2hZj/Y5PnvL447df3le/3r6SJEmS2uEMjyRJkqRmWfBIkiRJapYFjyRJkqRmWfBIkiRJalZfBU+SLUkOJjmUZMc8x29M8uUkDyf5bJJNXfv6JM917Q8n+fCgByBJkiRJJ7PgKm1J1gC7gEuBI8C+JFNVdaCn291V9eGu/xXAB4Et3bHHquqCwYYtSZIkSQvrZ4ZnM3Coqg5X1QvAHmBrb4eqerZn99VADS5ESZIkSVqafgqetcATPftHurZXSHJTkseADwDv6jm0IckXk/xekr823zdIckOS6STTMzMziwhfkl7JfCJpUMwnUhsGtmhBVe2qqvOB9wI3d81PAq+vqguB9wB3J3nNPOfeWVWTVTU5MTExqJAknYbMJ5IGxXwitaGfgucocF7P/rqu7WT2AFcCVNXzVfV0t/0Q8BjwpqWFKkmSJEmL00/Bsw/YmGRDkjOBbcBUb4ckG3t2Lwce7donukUPSPJGYCNweBCBS5IkSdJCFlylraqOJ9kO3A+sAXZX1f4kO4HpqpoCtie5BPge8AxwXXf6xcDOJN8DXgJurKpjKzEQSZIkSZprwYIHoKr2AnvntN3Ss/3uk5x3L3DvcgKUJEmSpKUa2KIFkiRJkjRqLHgkSZIkNcuCR5IkSVKzLHgkSZIkNcuCR5IkSVKzLHgkSZIkNcuCR5IkSVKzLHgkSZIkNcuCR5IkSVKzLHgkSZIkNcuCR5IkSVKzLHgkSZIkNauvgifJliQHkxxKsmOe4zcm+XKSh5N8NsmmnmPv6847mOTtgwxekiRJkk5lwYInyRpgF3AZsAm4treg6dxdVW+pqguADwAf7M7dBGwD3gxsAf5ddz1JkiRJWnH9zPBsBg5V1eGqegHYA2zt7VBVz/bsvhqobnsrsKeqnq+qPwEOddeTJEmSpBXXT8GzFniiZ/9I1/YKSW5K8hizMzzvWuS5NySZTjI9MzPTb+ySdALziaRBMZ9IbRjYogVVtauqzgfeC9y8yHPvrKrJqpqcmJgYVEiSTkPmE0mDYj6R2tBPwXMUOK9nf13XdjJ7gCuXeK4kSZIkDUw/Bc8+YGOSDUnOZHYRgqneDkk29uxeDjzabU8B25KclWQDsBH4w+WHLUmSJEkLO2OhDlV1PMl24H5gDbC7qvYn2QlMV9UUsD3JJcD3gGeA67pz9ye5BzgAHAduqqoXV2gskiRJkvQKCxY8AFW1F9g7p+2Wnu13n+Lc24DblhqgJEmSJC3VwBYtkCRJkqRRY8EjSZIkqVkWPJIkSZKaZcEjSZIkqVkWPJIkSZKaZcEjSZIkqVkWPJIkSZKaZcEjSZIkqVkWPJIkSZKaZcEjSZIkqVkWPJIkSZKaZcEjSZIkqVl9FTxJtiQ5mORQkh3zHH9PkgNJvpTkU0ne0HPsxSQPd4+pQQYvSZIkSadyxkIdkqwBdgGXAkeAfUmmqupAT7cvApNV9Z0kPw98ALimO/ZcVV0w4LglSZIkaUH9zPBsBg5V1eGqegHYA2zt7VBVn66q73S7DwLrBhumJEmSJC1ePwXPWuCJnv0jXdvJXA/c17P/qiTTSR5McuV8JyS5oeszPTMz00dIkjQ/84mkQTGfSG0Y6KIFSd4JTAK/3NP8hqqaBN4BfCjJ+XPPq6o7q2qyqiYnJiYGGZKk04z5RNKgmE+kNvRT8BwFzuvZX9e1vUKSS4BfAK6oqudfbq+qo93Xw8BngAuXEa8kSZIk9a2fgmcfsDHJhiRnAtuAV6y2luRC4A5mi52netrPSXJWt30u8JNA72IHkiRJkrRiFlylraqOJ9kO3A+sAXZX1f4kO4Hpqppi9ha2s4HfTALwtaq6AvhR4I4kLzFbXN0+Z3U3SZIkSVoxCxY8AFW1F9g7p+2Wnu1LTnLe54C3LCdASZIkSVqqgS5aIEmSJEmjxIJHkiRJUrMseCRJkiQ1y4JHkiRJUrMseCRJkiQ1y4JHkiRJUrMseCRJkiQ1y4JHkiRJUrMseCRJkiQ1y4JHkiRJUrMseCRJkiQ1y4JHkiRJUrP6KniSbElyMMmhJDvmOf6eJAeSfCnJp5K8oefYdUke7R7XDTJ4SZIkSTqVBQueJGuAXcBlwCbg2iSb5nT7IjBZVT8GfBz4QHfua4FbgYuAzcCtSc4ZXPiSJEmSdHL9zPBsBg5V1eGqegHYA2zt7VBVn66q73S7DwLruu23Aw9U1bGqegZ4ANgymNAlSZIk6dT6KXjWAk/07B/p2k7meuC+xZyb5IYk00mmZ2Zm+ghJkuZnPpE0KOYTqQ0DXbQgyTuBSeCXF3NeVd1ZVZNVNTkxMTHIkCSdZswnkgbFfCK1oZ+C5yhwXs/+uq7tFZJcAvwCcEVVPb+YcyVJkiRpJfRT8OwDNibZkORMYBsw1dshyYXAHcwWO0/1HLofeFuSc7rFCt7WtUmSJEnSijtjoQ5VdTzJdmYLlTXA7qran2QnMF1VU8zewnY28JtJAL5WVVdU1bEkv8hs0QSws6qOrchIJEmSJGmOBQsegKraC+yd03ZLz/Ylpzh3N7B7qQFKkiRJ0lINdNECSZIkSRolfc3wnE7W7/jkKY8/fvvlqxSJJEmSpOVyhkeSJElSsyx4JEmSJDXLgkeSJElSs3wPjySNkX7fZ7hQv96+kiS1zBkeSZIkSc2y4JEkSZLULAseSZIkSc2y4JEkSZLULAseSZIkSc2y4JEkSZLUrL4KniRbkhxMcijJjnmOX5zkC0mOJ7lqzrEXkzzcPaYGFbgkSZIkLWTBz+FJsgbYBVwKHAH2JZmqqgM93b4G/CzwT+e5xHNVdcEAYpUkSZKkRenng0c3A4eq6jBAkj3AVuD/FTxV9Xh37KUViFGSJEmSlqSfW9rWAk/07B/p2vr1qiTTSR5McuV8HZLc0PWZnpmZWcSlJemVzCeSBsV8IrVhNRYteENVTQLvAD6U5Py5HarqzqqarKrJiYmJVQhJUqvMJ5IGxXwitaGfW9qOAuf17K/r2vpSVUe7r4eTfAa4EHhsETFKklbQ+h2fPOXxx2+/fJUikSRp8PqZ4dkHbEyyIcmZwDagr9XWkpyT5Kxu+1zgJ+l5748kSZIkraQFC56qOg5sB+4HHgHuqar9SXYmuQIgyU8kOQJcDdyRZH93+o8C00n+CPg0cPuc1d0kSZIkacX0c0sbVbUX2Dun7Zae7X3M3uo297zPAW9ZZoyS1LSFbikDbyuTJGmpVmPRAkmSJEkaCgseSZIkSc2y4JEkSZLULAseSZIkSc2y4JEkSZLUrL5WaZMkSZJG2WJWvPQDl08vzvBIkiRJapYFjyRJkqRmWfBIkiRJapYFjyRJkqRmuWiBJEmSRtYwFxhYzEIIGl3O8EiSJElqVl8FT5ItSQ4mOZRkxzzHL07yhSTHk1w159h1SR7tHtcNKnBJkiRJWsiCBU+SNcAu4DJgE3Btkk1zun0N+Fng7jnnvha4FbgI2AzcmuSc5YctSZIkSQvr5z08m4FDVXUYIMkeYCtw4OUOVfV4d+ylOee+HXigqo51xx8AtgAfW3bkkqRV5b3skqRx1M8tbWuBJ3r2j3Rt/VjOuZIkSZK0LCOxaEGSG5JMJ5memZkZdjiSxpj5RNKgmE+kNvRT8BwFzuvZX9e19aOvc6vqzqqarKrJiYmJPi8tSScyn0gaFPOJ1IZ+3sOzD9iYZAOzxco24B19Xv9+4F/2LFTwNuB9i45SkiRJGmHD/LwgndqCMzxVdRzYzmzx8ghwT1XtT7IzyRUASX4iyRHgauCOJPu7c48Bv8hs0bQP2PnyAgaSJEmStNL6meGhqvYCe+e03dKzvY/Z29XmO3c3sHsZMUqSJKkhrvqo1TQSixZIkiRJ0kroa4ZHkrR43s8tSdLwOcMjSZIkqVkWPJIkSZKaZcEjSZIkqVkWPJIkSZKa5aIFkiRJ0ipxSe7V5wyPJEmSpGZZ8EiSJElqlre0SZIGzs8gkiSNCgseSZIkDYR/7NAosuBZBn+pJUmSpNHW13t4kmxJcjDJoSQ75jl+VpL/3B3/fJL1Xfv6JM8lebh7fHiw4UuSJEnSyS04w5NkDbALuBQ4AuxLMlVVB3q6XQ88U1U/nGQb8EvANd2xx6rqggHHLUmSJEkL6meGZzNwqKoOV9ULwB5g65w+W4GPdtsfB/5mkgwuTEmSJElavH7ew7MWeKJn/whw0cn6VNXxJN8Efqg7tiHJF4FngZur6veXF7IkSZLUPt8vPhgr/Tk8TwKvr6oLgfcAdyd5zdxOSW5IMp1kemZmZoVDktQy84mkQTGfSG3op+A5CpzXs7+ua5u3T5IzgB8Enq6q56vqaYCqegh4DHjT3G9QVXdW1WRVTU5MTCx+FJLUMZ9IGhTzidSGfgqefcDGJBuSnAlsA6bm9JkCruu2rwJ+t6oqyUS36AFJ3ghsBA4PJnRJkiRJOrUF38PTvSdnO3A/sAbYXVX7k+wEpqtqCvgI8BtJDgHHmC2KAC4Gdib5HvAScGNVHVuJgUiSJEnSXH198GhV7QX2zmm7pWf7u8DV85x3L3DvMmOUJEmSpCXpq+CRJGmluAqRNNoW+h0Ff0812ix4JEmSpDHnH49OzoJHkhbBv3RKkjReVvpzeCRJkiRpaCx4JEmSJDXLgkeSJElSsyx4JEmSJDXLRQskSZKk08TpuPiOBY8kaSycjv9JSyvJZYx1urDgWQX+Jy1JkiQNhwWPJEmSpBO0MgvoogWSJEmSmtXXDE+SLcC/AdYAv15Vt885fhbwH4C/BDwNXFNVj3fH3gdcD7wIvKuq7h9Y9JIkzaOVv0pKkpZvwYInyRpgF3ApcATYl2Sqqg70dLseeKaqfjjJNuCXgGuSbAK2AW8G/jzw35O8qapeHPRAJGk5fIEsqQXmMg3DqL9fvZ8Zns3Aoao6DJBkD7AV6C14tgLv77Y/DvzbJOna91TV88CfJDnUXe9/DiZ8SZKWbtT/k5ak1gyjKO+n4FkLPNGzfwS46GR9qup4km8CP9S1Pzjn3LVLjvY04F9mJGk0mZ81LBbm0vKkqk7dIbkK2FJVP9ft/wxwUVVt7+nzx12fI93+Y8wWRe8HHqyq/9i1fwS4r6o+Pud73ADc0O3+CHBwieM5F/jGEs8dVY5p9LU2HjhxTG+oqolhBbMY5pNTam1MrY0HTo8xmU/GX2vjAcc0LpaUT/qZ4TkKnNezv65rm6/PkSRnAD/I7OIF/ZxLVd0J3NlHLKeUZLqqJpd7nVHimEZfa+OB8R6T+eTkWhtTa+MBxzRqzCfza2084JjGxVLH1M+y1PuAjUk2JDmT2UUIpub0mQKu67avAn63ZqeOpoBtSc5KsgHYCPzhYoOUJEmSpKVYcIane0/OduB+Zpel3l1V+5PsBKaragr4CPAb3aIEx5gtiuj63cPsAgfHgZtcoU2SJEnSaunrc3iqai+wd07bLT3b3wWuPsm5twG3LSPGxVj2tPMIckyjr7XxQJtjWqwW/w1aG1Nr4wHH1KrW/g1aGw84pnGxpDEtuGiBJEmSJI2rft7DI0mSJEljqZmCJ8mWJAeTHEqyY9jxLFeSx5N8OcnDSaaHHc9SJNmd5Klu2fKX216b5IEkj3ZfzxlmjIt1kjG9P8nR7rl6OMlPDzPGxUpyXpJPJzmQZH+Sd3ftY/1cLVVruQTMJ6OqtXxiLjmR+WQ0mU9G36DzSRMFT5I1wC7gMmATcG2STcONaiD+RlVdMMZLCt4FbJnTtgP4VFVtBD7V7Y+TuzhxTAC/0j1XF3TveRsnx4F/UlWbgL8M3NT9/oz7c7VoDecSMJ+MortoK5+YS3qYT0baXZhPRt1A80kTBQ+wGThUVYer6gVgD7B1yDGd9qrqfzC7al+vrcBHu+2PAleualDLdJIxjbWqerKqvtBtfwt4BFjLmD9XS2QuGVHmk9FnLjmB+WREmU9G36DzSSsFz1rgiZ79I13bOCvgvyV5KLOf9NyK11XVk932/wZeN8xgBmh7ki91U8pjNQ3eK8l64ELg87T7XJ1Ki7kEzCfjZuzzibkEMJ+Mm1Z/Ts0ntFPwtOinquovMjsVflOSi4cd0KB1H07bwjKBvwqcD1wAPAn86+GGszRJzgbuBf5RVT3be6yh5+p0ZT4ZH2OfT8wlzTOfjA/zSaeVgucocF7P/rqubWxV1dHu61PAbzE7Nd6Cryf5cwDd16eGHM+yVdXXq+rFqnoJ+DXG8LlK8n3MJpT/VFWf6Jqbe6760FwuAfPJOBn3fGIueQXzyXhp7ufUfPL/tVLw7AM2JtmQ5ExgGzA15JiWLMmrk/zAy9vA24A/PvVZY2MKuK7bvg74L0OMZSBe/sXr/G3G7LlKEuAjwCNV9cGeQ809V31oKpeA+WTcjHM+MZecwHwyXpr7OTWf9FyvlQ8e7Zba+xCwBthdVbcNOaQlS/JGZv9qAnAGcPc4jifJx4C3AucCXwduBX4buAd4PfBV4O9W1di8ye4kY3ors9PFBTwO/P2e+0tHXpKfAn4f+DLwUtf8z5i9V3Zsn6ulaimXgPlklLWWT8wlJzKfjCbzyegbdD5ppuCRJEmSpLlauaVNkiRJkk5gwSNJkiSpWRY8kiRJkpplwSNJkiSpWRY8kiRJkpplwSNJkiSpWRY8kiRJkpplwSNJkiSpWf8X4WKWpn6I8oQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Poissson distribution\n",
    "fig, axes = plt.subplots(1, 3, figsize=(14, 3), sharey=True)\n",
    "\n",
    "xx = np.arange(20)\n",
    "lambdas = [1, 5, 10]\n",
    "\n",
    "for i, lam in enumerate(lambdas):\n",
    "    yy = stats.poisson.pmf(xx, lam)\n",
    "    axes[i].bar(xx, yy)\n",
    "    axes[i].set_title(r'$\\lambda={}$'.format(lam))\n",
    "\n",
    "# plt.savefig('Poisson_distribution.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "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>   100</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    98</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>log</td>       <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -213.95</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Tue, 24 Sep 2019</td> <th>  Deviance:          </th> <td>  120.86</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>20:35:37</td>     <th>  Pearson chi2:      </th>  <td>  115.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>5</td>        <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</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>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -0.2789</td> <td>    0.162</td> <td>   -1.720</td> <td> 0.085</td> <td>   -0.597</td> <td>    0.039</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.5646</td> <td>    0.042</td> <td>   13.508</td> <td> 0.000</td> <td>    0.483</td> <td>    0.647</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:                  100\n",
       "Model:                            GLM   Df Residuals:                       98\n",
       "Model Family:                 Poisson   Df Model:                            1\n",
       "Link Function:                    log   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -213.95\n",
       "Date:                Tue, 24 Sep 2019   Deviance:                       120.86\n",
       "Time:                        20:35:37   Pearson chi2:                     115.\n",
       "No. Iterations:                     5                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -0.2789      0.162     -1.720      0.085      -0.597       0.039\n",
       "x1             0.5646      0.042     13.508      0.000       0.483       0.647\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Y')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8XPV19/HPGa3W4l3yghfZ2DHIrK4DNNBAkoYAIaYlpCVNaUiTulDSJE+akqR9gNpZSpqnaROSQGlKyUKWJoHgUEIhIYFCFrAJGOMF23iVsSQvWkfrzHn+mJERYiSNRjNz70jf9+ul18zcucu5P1/p+N7f/Z1r7o6IiMhoIkEHICIihUEJQ0RE0qKEISIiaVHCEBGRtChhiIhIWpQwREQkLUoYIiKSFiUMERFJixKGiIikpTjoALJp9uzZXldXF3QYIiIFY9OmTUfcvSadeSdUwqirq2Pjxo1BhyEiUjDMbF+68+qSlIiIpEUJQ0RE0qKEISIiaVHCEBGRtChhiIhIWpQwREQkLUoYIiKSFiUMEZE8ao32sbOxndZoX9ChjNmEGrgnIhJmj+1oYv0DW098vvnyei5cURtgRGOjMwwRkTxojfax/oGtFEWMitJiiiLG+ge2FtSZhhKGiEgeNLV3A1BWXHTi1f2V6YUgZ5ekzOwu4HKgyd1PS077HrAiOct0oMXdz0qx7F6gHYgB/e6+OldxiojkQ211OQA9/THKiovo6Y9h9sr0QpDLM4y7gUsGT3D3P3b3s5JJ4ofAvSMs/6bkvEoWIlLwplWUcPPl9cTiTmdPP7G4c/Pl9UyrKAk6tLTl7AzD3R83s7pU35mZAX8EvDlX2xcRCZsLV9Ry78IZNLV3U1tdXlDJAoK7S+r3gEZ33znM9w48bGYO/Ju735m/0EREcmdaRUnBJYoBQSWMdwPfGeH7C9y9wcxqgUfMbLu7P55qRjNbC6wFWLRoUfYjFRERIIC7pMysGLgS+N5w87h7Q/K1CbgPOGeEee9099XuvrqmJq2HRomISAaCuK3294Ht7n4w1ZdmVmlm1QPvgYuBLXmMT0REUshZwjCz7wC/AlaY2UEze3/yq6sZcjnKzOab2YPJj3OAJ8zsOeAp4L/d/aFcxSkiIunJ5V1S7x5m+rUpph0CLku+fwk4M1dxiYhIZjTSW0RE0qKEISIiaVHCEBGRtChhiIhIWpQwREQkLUoYIiKSFiUMERFJixKGiIikRQlDRCaF1mgfOxvbC+qRqGETVLVaEZG8eWxHE+sf2Hri882X13PhitoAIypMOsMQkQmtNdrH+ge2UhQxKkqLKYoY6x/YqjONDChhiMiE1tTeDUBZcdGJV/dXpkv6lDBEZEKrrS4HoKc/duLV7JXpkj4lDBGZ0KZVlHDz5fXE4k5nTz+xuHPz5fUF+5jUIKnTW0QmvAtX1HLvwhk0tXdTW12uZJEhJQwRmRSmVZQoUYyTLkmJiEhalDBERCQtuXym911m1mRmWwZN+wczazCzZ5M/lw2z7CVmtsPMdpnZJ3IVo4iIpC+XZxh3A5ekmP4v7n5W8ufBoV+aWRHwFeBSoB54t5nV5zBOEZHQCWMpk5x1erv742ZWl8Gi5wC73P0lADP7LnAFsHXEpUREJoiwljIJog/jg2a2OXnJakaK708CDgz6fDA5TURkwgtzKZN8J4zbgZOBs4CXgX8e7wrNbK2ZbTSzjc3NzeNdnYhIoMJcyiSvCcPdG9095u5x4N9JXH4aqgFYOOjzguS04dZ5p7uvdvfVNTU12Q1YRCTPwlzKJK8Jw8zmDfr4h8CWFLM9DSw3syVmVgpcDWzIR3wiIkELcymTnHV6m9l3gIuA2WZ2ELgFuMjMzgIc2Av8ZXLe+cDX3P0yd+83sw8C/wMUAXe5+wu5ilNEJGzCWsrE3D3oGLJm9erVvnHjxqDDEBEpGGa2yd1XpzOvRnqLiEhalDBERCQtShgiIpIWJQwRkaQwluMIEz0PQ0SE8JbjCBOdYYjIpBfmchxhooQhIpNemMtxhIkShohMemEuxxEmShgiMumFuRxHmKjTW0SE8JbjCBMlDBGRpGkVJUoUI9AlKRERSYsShoiIpEUJQ0RE0qKEISKhodIcr5WqTYJqJ3V6i0goqDTHa6VqEyCwdtIDlEQkcK3RPq68/UmKIkZZcRE9/TFicefe68+ftHctpWqT3n4HnNLiSNbaSQ9QEpGCotIcr5WqTWLxOLG4B9ZOShgiEjiV5nitVG1SFIlQFLHA2ilnCcPM7jKzJjPbMmja581su5ltNrP7zGz6MMvuNbPnzexZM9M1JpEJTqU5XitVm6xbU8+6NSsDa6ec9WGY2RuBDuAb7n5actrFwKPu3m9mnwNw94+nWHYvsNrdj4xlm+rDEClsrdE+leYYIlWbZLOdxtKHkbO7pNz9cTOrGzLt4UEffw1clavti0jhUWmO10rVJkG1U5B9GH8O/GSY7xx42Mw2mdnakVZiZmvNbKOZbWxubs56kCIikhBIwjCzvwf6gXuGmeUCd18FXArckLy8lZK73+nuq919dU1NTQ6iFREJF3dn9427OfDPB/K63bwnDDO7FrgceI8P04Hi7g3J1ybgPuCcvAUoIhJi7s7uv9nNgc8foGtPF/kcS5fXhGFmlwA3AmvcPTrMPJVmVj3wHrgY2JJqXhGZfAq5fMh4Y/e4s/OGnRz8l4Oc9KGTWH7bcswsy1EOL2ed3mb2HeAiYLaZHQRuAT4JlAGPJHfy1+5+nZnNB77m7pcBc4D7kt8XA99294dyFaeIFI5CLh8y3tg95uz4wA4O332YhTcuZOmtS/OaLEClQUSkQBRy+ZDxxh7vi7Ptmm00f6+ZunV1LL5pcdaShUqDiMiEU8jlQ8YTe7wnzgvveoHm7zWz9J+WUndzXd7PLAYoYYhIQSjk8iGZxh6Lxnj+iuc5ev9Rlt22jEV/uyjnsY5ECUNECkIhlw/JJPb+1n42X7KZ4w8fZ8XXVrDggwvyGHFq6sMQkYJSyOVD0o29t7GXzZdupvP5Tk791qnU/nHuOvZDURpERCQXCrl8SDqxd+3pYvPFm+k51MNpPz6NWZfMylN0o1PCEBEJiY7nO9j8ts3Eu+Oc+dMzmfa704IO6VWUMEREQqDliRa2vGMLkcoIZ//v2VSurAw6pNdQp7eITEgjjarOdMR1JsuNtkxrtI/N3zzAc2/dTEltCaueXJVWsghixLvOMERkwhlpVHWmI64zWW60ZR7b0cSGm7by9h/C4fkw/+6FlC8e/TbhoEa86wxDRCaU1mgf6x/YSlHEqCgtpihirH9gK63RvhG/y3SdmS7T0tnLYx/cyju+D/teZ9xzA3zqly/mJJZsUcIQkQllpFHVmY64zmS5kZaJ98Z58c938MafwpZzjR/9ZYRIZXHOYskWXZISkQll8KjqgbpNQ0dVj/RdputMd5lZXszzb3+e6E+P89jF8OTbnLIiy2ks2aIzDBGZUEYaVZ3paPFMlku1zE2/s5zdb9tCyy9aWPGfK7joS/XEnJzHki0a6S0iE9JIo6ozHS2eyXIDy1Tt6mfPu7YR64yx8gcrmfnWmXmPJRWN9BaRSW+kUdWZjhbPZLlpFSX0PtzCtvdso2R2CWc/eTZVp1UFEst46ZKUiEiOuDv7P7+fF658gcqVlaz6zapXJYtCozMMEZEciHXFeHHtizR+q5Gad9VwytdPoWhKUdBhjYsShohIlnUf6GbLH26hY1MHdZ+qY/HfLcYiwTz0KJtyeknKzO4ysyYz2zJo2kwze8TMdiZfZwyz7HuT8+w0s/fmMk4RSQii3EQhSad9Wv63hU2rN9H1YhenbTiNuv9bNyGSBeT+DONu4MvANwZN+wTwM3e/1cw+kfz88cELmdlM4BZgNeDAJjPb4O7HcxyvyKQVVLmJQpFO+zTc0cCuv95F+ZJyTvvFaVSeGr4CguOR0zMMd38cODZk8hXA15Pvvw78QYpF3wY84u7HkkniEeCSnAUqMskFWW6iEIzWPvHeODv+cgc7r9/JjLfOYNVTqyZcsoBg7pKa4+4vJ98fBuakmOck4MCgzweT017DzNaa2UYz29jc3JzdSEUmiSDLTRSCkdqn53APz77pWV6+82UWfWIRp//4dEqmF+YDnkYT6G21nhg1OK6Rg+5+p7uvdvfVNTU1WYpMZHIZXG5i4DVf5SYKwXDtU/liH5tWb6Lj2Q7qv1fP0n9cihVNjP6KVIJIGI1mNg8g+dqUYp4GYOGgzwuS00QkB4IsN1EIUpb56J/Pzrc+jxUbq365ito/mvj9PUHcVrsBeC9wa/L1/hTz/A/w2UF3UF0MfDI/4YlMTheuqOXehTOyUm5iIhpon8YjUXo//TJH/v0Q0940nfr/qqd0dmnQ4eVFThOGmX0HuAiYbWYHSdz5dCvwX2b2fmAf8EfJeVcD17n7B9z9mJl9Cng6uar17j6081xEsiyIchOFpOxIjPardtH+dDsLP7aQJZ9dQqRk8hTMUPFBEZE0HH3wKNuu2Yb3O6fcfQo1fzgx+kxVfFBEJEvifXH23LSHA587QNVZVdT/Vz0VyyuCDisQShgiIsPoPtDN1ndvpe3JNuZfN5+T/+VkisoLux7UeAx78c3MHjSzuvyFIiKpqFxHdoylHVujfWz+xgGePnMjHc91MOP2Jcz55yWTOlnAyGcY/wk8bGZfB/7J3XW0iuSZynVkx1ja8RfPH+Y3f7Wdc5+AhvnOt6+N0d25B27fM+nbf9gzDHf/PrAKmApsNLOPmdlHB37yFqHIJKVyHdkxlnZsfKaVw5ckksXGNxqf/tMuXpzSS1lxkdqf0Qfu9QKdQBlQPeRHRHJI5TqyI512dHca7mhg+wXPUt0K962N8D9XQLzEMIy+WFztzwiXpMzsEuALJAbarXL3aN6iEpFXlaMoKy5SuY4MjdaOvY29bH//do799zGq3zydL5zfQtcMpyQSIZ4cdlBSFFH7M/IZxt8D73L3TyhZiOSfynVkx0jteGTDEZ4+/WmO//Q4y764jFWPnMlH35OYt7svxuyqMmZVltLdF1P7o4F7IqHXGu1TuY4sGNyOVR5h10d38fKdL1N5ZiX199RTubIy5bzAhG5/DdwTmUBUriM7Btqx9ZetbHzvdrp2d7HwxoUsWb+ESFkk5byDP4sShohMEvGeOHtu2cOBzx+gbGEZZz56JjMuSvmEaBmGEoaITHjtm9rZ9t5tRF+IMu8D8zj5n0+meKr+/I2VWkxEJqx4T5y96/ay/5/2UzqnlNP/+3RmXTYr6LAKlhKGiBSMsdwA0PZUG9vft53o1ihz3zeXk79wckE9OrU12sfu5g7AObmmOhT9KEoYIlIQ0i3vEeuMsefmPRz814OUzS/j9J+czqxLCuus4rEdTXzi3uc50tEDwKzKUj73zjMCL0syeZ78ISIFK93yHscePsbTpz/NwS8cZP7a+bz+hdcXXLJojfZxy4atHOnoIWJGxIyjnb3csuGFwMuS6AxDREIvVXmPzp5+mtq7mVZRQm9TL7s+uoume5qYsmIKZ/3iLKZfOD3IkDPW1N5NLB4/kSwA3KE/5if2NyhKGCISesOV96ipKuPlu19m99/sJtYeY/Eti1n8ycWvGVdRSGqryykaVJYEwHGKiyzwsiR5b1UzW2Fmzw76aTOzjwyZ5yIzax00z835jlNEwiNVeY+b6k9m7+Vb2fG+HVTWV7L6udUs+YfXDsIrNNMqSli3pp7ZVWXE3Ym7M6uylHVrVgbe8R1oaRAzKwIagHPdfd+g6RcBH3P3y8eyPpUGEZnYWqN9HG6MEr+9maZ/baCoqoilty5l3gfmYRELOrysytddUoVUGuQtwO7ByUJEJBV3p++RVo59eCc9+3qY82dzOPnzJ1NaWxp0aDkxraKEVYvDNRI96IRxNfCdYb77XTN7DjhE4mzjhVQzmdlaYC3AokWLchKkiASrc1snuz6yi+MPH6eivqKgO7ULWWCXpMyslEQyWOnujUO+mwrE3b3DzC4Dvujuy0dbpy5JiUws/a397F23l4bbGohURliybgnz/2o+kZLC7qcIk0K5JHUp8MzQZAHg7m2D3j9oZl81s9nufiSvEYpIIDzuHP7Pw7z0yZfoO9LHvA/MY8lnllBaMzEvPxWKIBPGuxnmcpSZzQUa3d3N7BwSd3MdzWdwk81kfeZCGMsvZCqxL+2AcXJNVcHuS+svW9n5oZ10bOpg6humcsZDZ1C9KvVToSfrcRuUQBKGmVUCbwX+ctC06wDc/Q7gKuB6M+sHuoCrfSI96Slk0i25MNGEtfxCJh7b0cTHf7iZo529AMyuKuPWK08vqH3p3t/Nnr/fQ+O3GimdX8qp95xK7btrMUt999NkPW6DpCfuTXKt0T6uvP1JiiJ2YkBULO7ce/35E/p/bK3RPq74ypM0tERPjKaNu3PS9Cncf8MFBbXviX15goaWriH7UsH9N4T/37HveB/7P7ufg7cdBGDh3yxk0ScXUVw1/P9nJ+txmwuF0ochITBayYWJKszlF8YqsS+OMWhfgFg8Hup9iXXHaLitgf2f3U9/az9z/mwOS9YvoXzR6KOZJ+txGzQljEluuJILQZcgyLUwl18Yq8S+GI4TT+6Ou1MUiYRyXzzmNH6rkT037aHnQA8zL53J0luXUnVGVdrrmKzHbdB0b9okl6rkws2X10/4/6WFufzCWCX2ZSWzKktP7MvsqjLWrQnXv6O7c/TBo2w8eyPbr91O6ZxSznz0TM548IwxJQuYvMdt0NSHIcDkvdtEd0nlR9tTbbz08Zdo+UUL5SeXs/SzS6l5V82wHdrpmqzHbTapD0PGbFpFyaT8hQtj+YVMJfZlZtBhvEp0Z5Q9f7+H5u83U1JTwrLbljF/7Xwipdm5uDFZj9ugKGGISNZ1butk/z/up/HbjUTKIiy+aTELP7aQ4qn6k1PI9K8nIlnT8VwH+z6zj+YfNBOZEmHBhxew8G8XUja3LOjQJAuUMERk3NqeamPfp/dx9MdHKaouYtEnF7HgIwtUymOCUcIQCVghd9y2/G8L+z61j+OPHKd4ZjF16+s46YMnUTIjXPuRizbO5joL5RhQwhAJUCGWt3B3jv/0OPs+vY/Wx1spqS1h6eeWMv/6+RRXh+9PSi7aOJvrLKRjQOMwRALSGu1j/QNbKYoYFaXFFEWM9Q9spTXaF3RoKbk7R358hGfOe4bNF2+ma3cXy764jPP2nMeiGxeFMlnkoo2zuc5COwbC9y8sMkkUSnkLjzvN9zaz79P76Hyuk/K6cl73b69j7nvnhv752blo42yus1COgQFKGCIBCXt5i3h/nKbvNrH/s/uJbosy5XVTOOXuU6j9k9qCeYBRLto4m+sM+zEwVGH8q4tMQGEtbxHviXPoa4d4asVTbL9mO1Zk1H+3nnO2npM4qyiQZAG5aeNsrjOsx8BwVBpEJGBhuUOm72gfh/7tEA23NdB7uJfq1dUs/r+LmfWOWVhkfCU8gqa7pIan0iAiBSTo8hadL3Ry8LaDNH6jkXhXnBlvm8Ep3ziFGb8/Y9y1nsIiF22czXUGfQykSwlDZBLyeKJybMMXGzj+0+NYmTHnT+ew4CMLqDptbJVjZfJQwhCZRPpb+zn89cM03NZA164uSueXsuQzS5i3dh6lszUqW0YWWMIws71AOxAD+odeQ7PEufAXgcuAKHCtuz+T7zhl/MJyjX4y69jSwaGvHOLwNw8T74wz9byp1H2qjpp31rymE7vQr82H9XgLa1xjEfQZxpvc/cgw310KLE/+nAvcnnyVAlJIo1gnmnhPnOYfNnPojkO0/m9r4rLTn8zhpBtOovp3qlMuU+gjmMN6vIU1rrEK8/1xVwDf8IRfA9PNbF7QQUn6Cm0U60QR3RVl9427+dWCX7HtPdvoOdTD0n9ayhsa3sApd50ybLIo9BHMYT3ewhpXJoI8w3DgYTNz4N/c/c4h358EHBj0+WBy2suDZzKztcBagEWLFuUuWhmzQhvFWsjiPXGO3H+EQ3ceouVnLVAEs6+Yzfzr5jPjLTPSui220Ecwh/V4C2tcmQgyYVzg7g1mVgs8Ymbb3f3xsa4kmWjuhMQ4jGwHKZkrtFGshajzhU5e/o+XafxmI31H+ihbXEbd+jrmvX8eZfPH9gyKQh/BHNbjLaxxZSKwS1Lu3pB8bQLuA84ZMksDsHDQ5wXJaVIgCm0Ua6HoO95Hwx0NbDp3E0+f9jQNX25g2oXTOP0np3Pe7vOou6luzMkCCn8Ec1iPt7DGlYlARnqbWSUQcff25PtHgPXu/tCged4OfJDEXVLnAl9y96FJ5VU00jucJsLdIUHzmHPskWMcvvswR350BO9xKk+rZO775jLnmjlZfVCR7pLKjbDGVQgjvecA9yVHkRYD33b3h8zsOgB3vwN4kESy2EXittr3BRSrjFOhjGING3enc3Mnjfc00nhPI72HeimeWcz8tfOZe+1cqs6uyslI7EIfwRzW4y2scY1FIAnD3V8Czkwx/Y5B7x24IZ9xiYRB194umr7dROM9jUS3RrFiY+alM5l721xmvX1W6EuKy8QV9DgMEQF6j/TS/P1mGu9ppO3JNgCmnj+V5V9dTs27ajQKW0JBCUMkIP3t/Rz98VGavtPEsYeO4f1ORX0FSz6zhNo/qWVK3ZSgQxR5FSWMSSYMHW/pxDDcPGOZnu19HU/cA/o7+jn6wFGav9/MsQePEe+OU3pSKQs+soDa99RSdeZr+yVG2+fy4iK6+2Np72c6bQgEfpxkKgzH+ESlhDGJhKE8QToxDDfPWKYDWd3X8cTd3z4oSfwkmSTmljL3/XOp/eNapp0/bdiBdaPtc2dPP0c7e5ldVUZFadGo+5lOG0Z7Y7g7lWXFw+5rWIXhGJ/I9AClSaI12seVtz9JUcRODB6KxZ17rz8/r7c7jhbDcPPcfe05XHv3U2lN7+2PA0ZpcXb2NZO4rbWfk7fANe0zaP9ZC97jlM4rpeadNdS8qyaRJIpGvsNptLYAOHC8C3cnYsaCGYlLWMPtZzptWxyJsLOpHYDltdX0x+N5P04yFYZjvBAVwm21kmdhKE+QTgzDzbOjsS3t6dGeGGZOWXHpsNvJRdxTW+C07cayzTEW7oJIHDrndzD/uvnUXFXDtDcMfyYxlu0O7HPEDByKIxFiyaQRi/uw+5lO23b3xYiY4Q59sTjlJYVTxiIMx/hEp4QxSYShPEE6MQw3z4o5U9OeXlxkgGVtX4eLqaaqjPZN7Rz58RGO3X+EDz8LEOfoHPjVm2DbSvjyZ1/P9MrM7nAarS3i7mDQH48TMSPuPuJ+ptO2JUWRxHqBkqJIQZWxCMMxPtHphu5JIgzlCdKJYbh5Fs6qSHv6ujUrWbcme/s6OKaejn6WbnFuemo6W1+3kU2rN7Fv/T5KphTBx2r56o3wlY/Bo5fCX/xVfcbJIp22AJgxJbFPMyoS2xlpP9Np2+6+GLOryphVWUp3X6ygyliE4Rif6NSHMcmE4Q6SQrtLqnt/N8d+cozDG5pp/3kL3uUUVRUx420zmP2O2cy8dCaltaVZ3eZguktqbMJwjBeSsfRhKGGIDBHrjtH6eCvHHjrGsYeOEd0WBaC8rpxZl89i1jtmMf3C6RpxLROCOr1FxsDd6drZdSJBtPyihXhXHCszpr9xOvM+MI+Zl8yk4tSKnNRuEikUShgyKfUe6aXl5y0cf+Q4xx85TvfexB02U143hXl/kUgQ0y+cTlFFUcCRioSHEoZMCv0d/bQ+0UrLoy0c/9lxOn7bAQ5FU4uY/qbpLLxxITPfNpMpS1WOQ2Q4ShhpyEZJiLFsa3dzO2CcXFOVl0678cSeSbzZaKvR1hHrjtH267bEWcTPjtP+m3a837ESY+p5U6lbV8eMt86genU1keKR+yIyiTfRLh2Ac3JN9WuWy7TDOhP5PH7DZCLuU9CUMEYxnpIQmWzr4z/czNHOXgBmV5Vx65Wn57S0wXhizyTebLRVqnVcsGgWbb9qo+WxFloea6Ht1214j0MEqldXs/BjC5n+5ulMO3/amC4zZRLvYzua+MS9z3OkoweAWZWlfO6dZ5xYbmCd0d4YRzp6mFVZSmVZcU7KWOTz+A2TibhPYaC7pEYwnlIWYy1H0Brt44qvPEFDS1diBC+JgVknTa/g/htyU9pgPLFnEm822mpgHRU9sGR/hLk74yx4CRY3GN6XSBBVZ1cx/aLpiZ/fm07xtMz+X5RJvIl2eZKGluiQdpnC/TdcAMCVtz8JwMHjXcnBdsbCUcp65Cr+iVhOYyLuUy7pLqksGU8pi7GWI2hq7yYWdww78YfGgVg8nrPSBuOJPZN4M92eu9O9t5u2X7ax/5EjrH0Yag+DeZxYBF5eAJVra1lyaS3TLpiWcYLIRryJdkmMvD7RLg79MT+xPuDEd8WRCLH46GU9chX/RCynMRH3KSyUMEYwnlIWYy1HUFtdTlHEcJx48qTP3SmKRHJW2mA8sWcSb7rbi3XH6PhtB22/bqP1yVbaftlG78uJy16R6iI65sGOM+Dw8gh7F8TpKYZ7r1+W9T8GmbRPol1eKa8B4DjFRfaq5Qa+74/HsTTKeuQq/olYTmMi7lNY5H3kkZktNLOfm9lWM3vBzD6cYp6LzKzVzJ5N/tyc7zhhfKUsxvrHa1pFCevWrGRWZSlxd+LuzK4qY92a3JU2GE/smcSbans3vf1USg710XhPIzv/eiebztnEE1Of4Ldv+C27P7qbjk0dTH/zdJZ/ZTmrn13N7x2/gGUb6nnsrbB9cSJZ5Kr8Qybtk2iXemZXlZ1ol1mVpaxbs/LEM51PlPVIlvMYKO+R7f3I5/EbJhNxn8Ii730YZjYPmOfuz5hZNbAJ+AN33zponouAj7n75WNZd65GeusuqZGXHUu8fUf7OPzEcZqePI5t7ia6sZ3+o/0ARCojTH39VKrPrWbquVOZeu5UyuaXZT3msdJdUoVpIu5TLhRUaRAzux/4srs/MmjaRYQoYUhm+lr66Himg/ZN7bRvbKf96Xa69ySv4xtU1FckEsN5ieRQubJy1GdEiEh2FUynt5nVAWcDv0nx9e+a2XPAIRLJ44Vh1rEWWAuwaNGi3AQqo+o90kvHbzsldUfvAAANr0lEQVTo+G0iQXQ800HXrq4T35fXlVP9+mrmXzef6tXVVK+upniqutBECklgv7FmVgX8EPiIu7cN+foZYLG7d5jZZcCPgOWp1uPudwJ3QuIMI4chC+Bxp+ulLjqf66TjuQ46nk0kiZ6DPSfmKa8rp2pVFXPfN5fq1dVUraqidHbmZb5FJBwCSRhmVkIiWdzj7vcO/X5wAnH3B83sq2Y2292P5DPOya6/rZ+OzR10bu5MvD7XScfzHcQ744kZIlCxooJpb5xG1dlVVJ9dTdXZVZTM1PVikYko7wnDEuU+/wPY5u5fGGaeuUCju7uZnUPibq6juY4tX51kQ7eTzQ7zTNYT74kT3R6lc0vniZ/2zR307n/lrKF4ejGVZ1Yy78/nUXlGJVVnVlF5WiVFU149aro12kdTY3teOxrD/hyHsT7bQySsgjjDOB+4BnjezJ5NTvs7YBGAu98BXAVcb2b9QBdwtee4dz5fpQSGbucPzprPj549NO7tphN/rCtG14tddG7tpPOFTqJbo3S+0Jnoa0ieNFix4UtK2TKrh6aV0DgP3n3Ncs6/cP6opb2DKMcweJudPf2YGRWlRXnb/miGaxOVrpBCFPhdUtmU6V1S+SolMHQ70d5+9h+LsmhmBRWlxeMqKzJ4vdbWz4zDcEv9MvylHjq3dRLdFk3coTTwz10EU5ZNoXJlJZX1lVSsrKDytEr6FpRw1X/8asxtEUQ5hsHbLI5E2NnUDsDy2mr64/HAy0EM1yZ3X3sO1979lEpXSCgUzF1SYZGvUgJDtxMxA3+lTMRYthvvidP1UhfRHVEObGrh8oeh9ogzo7Gfis7EPPvZhZUZFSsqqH59NXOumUPlqZVU1FdQ8bqKlE+M29nY/qoY040piHIMg7fZ3RfDMDDoi8UpLwm+HMRwbbKjsS3ldJWukLBTwiB/pQSGbifuDvZKmYih2411x+je003Xzi66dnWdeI3ujNKzv+eVswVgeTUcq4FdZxhNNc6R2bD+/6yi5tTqMY1tyLQtgijHMHibJUURHAeHkqJIKMpBDNcmK+ZMTTldpSsk7HRJKmngmrI7mOW+D2NgO1cum8svHzvM9CMw4xj8fuV0ph+Frl1dr0kKxTOLmbJsClOWT0m8LptCxYoKpiyfwi8bj2ct/kzbIl9tONw2o72JPowpJUV52/5Y4hscUxBtJZJKQY30zqbxjvTOxV0rsc4Y3fu66d7bnXjd0037riidu7uIH+wldqz/VfMXzypmysmJpFCxvOJVCWK021WzGX+m6wrizh/dJSWSOSWMPHF3+lv76dnfk0gG+7rp2deTSA7JBNHX3PeqZazUKK8rp3xJOVOWTnnldWniNVuluUVE0qFO7yyJRWP0HOw58dO9v5ueAz30HEi+399DrD32qmWszChfXE55XTmzV81OJIfk5/K6ckrnlmIR1UsSkcKjhDHI3k/tpe1XbScSRP/x/tfMUzK7hLJFZUxZNoUZb5lB+aJyyhaVnXgtnaOEICITkxLGINHtUXobeylfUs6035tG2UlllC0se+V1YRlF5ek/D1pEZCJRwhjkpH9fPmzn5PbmDjr2d1BVVpLxcyqCKj2Sq/WP9CyHsHbo5iquMOxvGGKQiU0JI2mkEg6fuPd5mtp7ks9ehjlTy7n1ytPHdBtkUKVHsr2dgfVHe2Mc6ehhVmUplWXFr9pOWMte5CquMOxvGGKQiS/vj2gNo9ZoH+sf2EpRxKgoLaYoYqx/YCsHjka5ZcNWjrT3EE8+uNodmtu7uWXDVlqjfaOseeT1p7v8ePcjW9sZWD/Asc7EM7aPJdc9sJ187etY5SquMOxvGGKQyUEJg9QlHNxhR2MbsXgckn3YduLViMXjJ5bLdP3pLp+uXG9nYD0DpUyKI5ETpU0GtpOvfR2rXMUVhv0NQwwyOShh8OoSDgOvAyUcipJ/FCFxdpF4dYoikbRLOQy3/lyWHsnFdgbWM1DKpD+ZTOPuJ7aTr30dq1zFFYb9DUMMMjkoYQDTKkq4+fJ6YnGns6efWNy5+fJ6Fs6qYN2aemZXlxFJ3iprBjXV5axbU592x+Jw6892x2SutzOwfoAZFYkn6M2Yklj3wHbyta9jlau4wrC/YYhBJgeN9B5kpBIOu5s76Ojp011S6C6pfK630GKQwqPSICIikpaxJAxdkhIRkbQoYYiISFoCSRhmdomZ7TCzXWb2iRTfl5nZ95Lf/8bM6vIfpYiIDJb3hGFmRcBXgEuBeuDdZlY/ZLb3A8fdfRnwL8Dn8htlogNxZ2N7ysFPI3033m2ku+5MYzhwNMpPtx7mwNFoRjGLyOQVRGmQc4Bd7v4SgJl9F7gC2DponiuAf0i+/wHwZTMzz1MP/UhlFrJVgiHVeoC01p1pDLf97EW+9OiuE58/9OZl/PVbXjfm2EVkcgriktRJwIFBnw8mp6Wcx937gVZgVj6CG6nMQrZKMKRazy0bXuCWDaOvO9MYDhyNnkgWxZHEP/uXHt2lMw0RSVvBd3qb2Voz22hmG5ubm8e9vpHKLGSrBEOq9fTHnFg8Puq6M41hR2Mb8EqyKI5ETpQ/ERFJRxAJowFYOOjzguS0lPOYWTEwDTiaamXufqe7r3b31TU1NeMObqQyC9kqwZBqPcVFRlEkMuq6M41hxZypQLKcR/J1oPyJiEg6gkgYTwPLzWyJmZUCVwMbhsyzAXhv8v1VwKP56r8YqcxCtkowpFrPujUrWbdm9HVnGsPCWRV86M3LAOiLJZLGh968jIWzKsYUu4hMXoGM9Dazy4B/BYqAu9z9M2a2Htjo7hvMrBz4JnA2cAy4eqCTfCTZHOk9UpmFbJVgSLWedNedaQwHjkbZ0djGijlTlSxERKVBREQkPSoNIiIiWaeEISIiaVHCEBGRtChhiIhIWpQwREQkLUoYIiKSlgl1W62ZNQP7xrGK2cCRLIWTbWGNTXGNjeIaG8U1NpnEtdjd0yqTMaESxniZ2cZ070fOt7DGprjGRnGNjeIam1zHpUtSIiKSFiUMERFJixLGq90ZdAAjCGtsimtsFNfYKK6xyWlc6sMQEZG06AxDRETSMikThpndZWZNZrZlmO/NzL5kZrvMbLOZrQpJXBeZWauZPZv8uTlPcS00s5+b2VYze8HMPpxinry3WZpx5b3NzKzczJ4ys+eSca1LMU+ZmX0v2V6/MbO6kMR1rZk1D2qvD+Q6rkHbLjKz35rZAym+y3t7pRlXIO1lZnvN7PnkNl9Tojtnv4/uPul+gDcCq4Atw3x/GfATwIDzgN+EJK6LgAcCaK95wKrk+2rgRaA+6DZLM668t1myDaqS70uA3wDnDZnnr4A7ku+vBr4XkriuBb6c72Msue2PAt9O9e8VRHulGVcg7QXsBWaP8H1Ofh8n5RmGuz9O4sFMw7kC+IYn/BqYbmbzQhBXINz9ZXd/Jvm+HdgGnDRktry3WZpx5V2yDTqSH0uSP0M7C68Avp58/wPgLWZmIYgrEGa2AHg78LVhZsl7e6UZV1jl5PdxUiaMNJwEHBj0+SAh+EOU9LvJSwo/MbOV+d548lLA2ST+dzpYoG02QlwQQJslL2M8CzQBj7j7sO3l7v1AKzArBHEBvDN5GeMHZrYw1zEl/StwIxAf5vtA2iuNuCCY9nLgYTPbZGZrU3yfk99HJYzC8gyJYfxnArcBP8rnxs2sCvgh8BF3b8vntkcySlyBtJm7x9z9LGABcI6ZnZaP7Y4mjbh+DNS5+xnAI7zyv/qcMbPLgSZ335TrbY1FmnHlvb2SLnD3VcClwA1m9sZ8bFQJI7UGYPD/FBYkpwXK3dsGLim4+4NAiZnNzse2zayExB/le9z93hSzBNJmo8UVZJslt9kC/By4ZMhXJ9rLzIqBacDRoONy96Pu3pP8+DXgd/IQzvnAGjPbC3wXeLOZfWvIPEG016hxBdReuHtD8rUJuA84Z8gsOfl9VMJIbQPwZ8k7Dc4DWt395aCDMrO5A9dtzewcEv9+Of8jk9zmfwDb3P0Lw8yW9zZLJ64g2szMasxsevL9FOCtwPYhs20A3pt8fxXwqCd7K4OMa8h17jUk+oVyyt0/6e4L3L2ORIf2o+7+p0Nmy3t7pRNXEO1lZpVmVj3wHrgYGHpnZU5+H4vHu4JCZGbfIXH3zGwzOwjcQqIDEHe/A3iQxF0Gu4Ao8L6QxHUVcL2Z9QNdwNW5/qVJOh+4Bng+ef0b4O+ARYNiC6LN0okriDabB3zdzIpIJKj/cvcHzGw9sNHdN5BIdN80s10kbnS4OscxpRvXh8xsDdCfjOvaPMSVUgjaK524gmivOcB9yf8HFQPfdveHzOw6yO3vo0Z6i4hIWnRJSkRE0qKEISIiaVHCEBGRtChhiIhIWpQwREQkLUoYIjliiWq6e8xsZvLzjOTnumAjE8mMEoZIjrj7AeB24NbkpFuBO919b2BBiYyDxmGI5FCydMkm4C7gL4Cz3L0v2KhEMjMpR3qL5Iu795nZ3wIPARcrWUgh0yUpkdy7FHgZCEXFWpFMKWGI5JCZnUWiyN95wP/Jx4O4RHJFCUMkR5JVcm8n8ZyO/cDngf8XbFQimVPCEMmdvwD2u/sjyc9fBU41swsDjEkkY7pLSkRE0qIzDBERSYsShoiIpEUJQ0RE0qKEISIiaVHCEBGRtChhiIhIWpQwREQkLUoYIiKSlv8P439Hlzkdx3sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Poisson regression\n",
    "\n",
    "# generate simulation data\n",
    "np.random.seed(5)\n",
    "n_sample = 100\n",
    "a = 0.6\n",
    "b = -0.4\n",
    "\n",
    "x = uniform(1, 5, size=n_sample)\n",
    "mu = np.exp(a * x + b)\n",
    "y = poisson(mu)\n",
    "\n",
    "exog, endog = sm.add_constant(x), y\n",
    "\n",
    "# Poisson regression\n",
    "mod = sm.GLM(endog, exog, family=sm.families.Poisson(link=sm.families.links.log()))\n",
    "res = mod.fit()\n",
    "display(res.summary())\n",
    "\n",
    "y_pred = res.predict(exog)\n",
    "\n",
    "idx = x.argsort()\n",
    "x_ord, y_pred_ord = x[idx], y_pred[idx]\n",
    "plt.plot(x_ord, y_pred_ord, color='m')\n",
    "plt.scatter(x, y,  s=20, alpha=0.8)\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12e77b908>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH9BJREFUeJzt3Xt81PWd7/HXZ2ZyISGES7hoEgg3i1wsaEBR26LYFm1Xa29Ct1tLPWXbam9rt2vbc1zXnt112z27p1W7XXuR1rogVupSi9aeqq1bF7kpd5HIJRduAUIIhCRz+Zw/ZrCRAgkwyW9m8n4+HnnMzG9+Yd4PmLz55vf7zu9r7o6IiOSWUNABREQk/VTuIiI5SOUuIpKDVO4iIjlI5S4ikoNU7iIiOUjlLiKSg1TuIiI5SOUuIpKDIkG9cFlZmVdVVQX18iIiWWnNmjUH3H1oV/sFVu5VVVWsXr06qJcXEclKZrarO/vpsIyISA5SuYuI5CCVu4hIDlK5i4jkIJW7iEgO6rLczezHZrbfzDae5nkzs++aWY2ZrTezS9MfU0REzkZ3Ru4LgTlneP56YHzqawHwb+cfS0REzkeX89zd/fdmVnWGXW4CfurJ9fpWmNlAM7vA3fekKaOISLe4O9G40xFP0BFLfkXjJ778zfuxhBOLO/GEE0skUrdO4sStJ7/iCUikHsfdSXjyNZLbIOGOn7jlj48h+X0OuIPjqdvkhtkXD+ftlQN79O8iHR9iKgfqOj2uT237k3I3swUkR/eMHDkyDS8tItkunnAOt3bQ1NpBU2uUQ8c6aGmLcbQtmrxtj9HSHuNoW4xj7TGOR+PJr44/3rZF47THEnTEE2TDstDDBhRmRbl3m7s/BDwEUF1dnQX/BCJyrtydQ8c62HWolX3Nbew9kvza19zGnuY29re0c+hYB0faomcs5IJIiJLCCP0LIhTlRyjKD9O/IEJZ/wKK8sP0ywtTmBemIBIiPxJ68zY/HCI/EiYvbORHQkRCIfLCRl44RCRsRELJ23DICFvqNmREQkYoZIQsuT0UgnDqcfKLN+9bCIxOjy2Z+cR9A8wsdZu831vSUe4NQGWnxxWpbSLSB0TjCbbtO8rr+1rYceAYOw8eY8eB5FdLW+wt++aHQwwbUMCIAYVMvGAAQ/rnM6gon0FFeQwqzmdwcfLxgMI8+qcKPT+iSX3nIh3lvgy4w8wWA5cDzTreLpKbovEEr+9rYWNDMxsamtnQcIQte47QEUsAydFp+cB+jC4r5uZp5VQNKWbUkCJGlBYyYkAhg4vze3X02pd1We5mtgiYBZSZWT3wt0AegLt/H1gO3ADUAK3A/J4KKyK9y93ZfuAYL77eyIvbDrBi+0GOdcQBKCmIMKl8ALfOHMXk8lImXjCAkUOKKIiEA04t0L3ZMvO6eN6B29OWSEQC1RaN89xr+/nd1kZe3NbI7uY2AEYNKeLmS8uZXjWYSyoGMmpwEaGQRuGZKrBL/opI5nB31uxq4om19Ty1fg8tbTFKCiNcNbaM268t4x3jhjJySFHQMeUsqNxF+rC6Q60sXdvA0lfq2XWwlX55Ya6fPIIPXlrBFWMGEwnrZGa2UrmL9EFb97bw3d9uY/nG5NyHmWOG8Plrx3P95BEUF6gWcoH+FUX6kNf2HkmW+oa99C+I8LlZY/nY5aMoH9gv6GiSZip3kT5gy55kqT+9cS8lBRG+cO04PnX1aAYW5QcdTXqIyl0khx1pi/L3T23hsdV1yVKfPZ7brhpNaVFe0NGkh6ncRXLU81v38/WlG9h3pI2/fOcYPjdrnEq9D1G5i+SY5uNR/vdTm3l8TT0XDe/P9z9+VY9fpEoyj8pdJIc899o+vrZ0AweOdnD7NWP5wuzx+sRoH6VyF8kB0XiCe5Zt4tGXa3nb8BJ+8IlqLqnQaL0vU7mLZLmj7TE+9+hafv96IwveOYY733ORRuuichfJZvuOtDH/4VVs3dfCfR+cwtwZWgRHklTuIllq694W5j+8kubjUX50azWz3jYs6EiSQVTuIlnoDzUH+Mwja+iXH2bJZ2Yy6cLSoCNJhlG5i2SZpWvr+erP1zNmaDEPz5+hSwfIKancRbLIMxv3cOfj65g5Zgj/9vHLKO2nDyXJqancRbLE6p2H+OLiV5laOZAf3TqdfvmaESOnp4s1i2SBmv1Hue0nq7lwYD8Vu3SLyl0kw+0/0satP15JXtj4yfwZDC7WlRylazosI5LBjrbHmL9wFU2tHTy2YKaWupNuU7mLZKhoPMFnf7aG1/a28MNbq5lSoemO0n06LCOSgdydu57YwIvbDvCPH5zCNfqAkpwllbtIBlq0so4n1tbzpevG89HqyqDjSBZSuYtkmO2NR/nmU5u5elwZX7h2fNBxJEup3EUySDSe4MuPvUp+JMQ/f+TthEIWdCTJUjqhKpJB7v/tNtbVN/O9P7+UEaWFQceRLKaRu0iGWLPrEA88X8OHLq3ghikXBB1HspzKXSQDHG2P8eXH1nHhwH7cc+PEoONIDtBhGZEMcO8vN1Hf1MpjfzmTkkJdDEzOn0buIgF7ZuMelqyu53OzxjG9anDQcSRHdKvczWyOmW01sxozu+sUz480s+fN7BUzW29mN6Q/qkju2d/SxteWbmBKeSlfvE7THiV9uix3MwsDDwLXAxOBeWZ28kHB/wkscfdpwFzge+kOKpKL7nv6NY61x/nXW6aSF9Yv0pI+3Xk3zQBq3H27u3cAi4GbTtrHgQGp+6XA7vRFFMlN6+oOs3RtA5+6ejTjhvUPOo7kmO6cUC0H6jo9rgcuP2mfe4BnzezzQDFwXVrSieQod+fepzZT1r+AO64dF3QcyUHp+j1wHrDQ3SuAG4BHzOxP/mwzW2Bmq81sdWNjY5peWiT7LFu3mzW7mvjqe99G/wJNWpP06065NwCdr1xUkdrW2W3AEgB3/2+gECg7+Q9y94fcvdrdq4cOHXpuiUWy3PGOOPc9/RqTywfw4csqgo4jOao75b4KGG9mo80sn+QJ02Un7VMLzAYws4tJlruG5iKn8O+/f4M9zW3c/f5JunaM9Jguy93dY8AdwK+BLSRnxWwys3vN7MbUbncCnzazdcAi4JPu7j0VWiRb7T58nO//7g3eN+UCZozWnHbpOd062Ofuy4HlJ227u9P9zcBV6Y0mknu+9cxrJBzuun5C0FEkx2lirUgvWVvbxJOv7mbBO8ZQOVhroUrPUrmL9IJEwvm7X25mWEkBn501Nug40geo3EV6wbJ1u1lXd5i/mTOBYk19lF6gchfpYfGE853fbmPiBQO4eVp50HGkj1C5i/SwX23Yw44Dx/jC7HGa+ii9RuUu0oMSCeeB57Yxflh/3jNxRNBxpA9RuYv0oGc37+P1fUe541qN2qV3qdxFeoi788Dz26gaUsT7tCaq9DKVu0gPeeH1RjY2HOFzs8YR0bXapZfpHSfSA9yd+3+7jfKB/fiAZshIAFTuIj3gv7cfZG3tYT7zrjHkR/RjJr1P7zqRHvDAczUMLSngI9WVXe8s0gNU7iJptmbXIV564yB/+c4xFOaFg44jfZTKXSTN7n+uhsHF+Xzs8pFBR5E+TOUukkYb6pt5YWsjt109mqJ8XUNGgqNyF0mjB5+vYUBhhE/MHBV0FOnjVO4iaVLf1Mqzm/fy51eMoqQwL+g40sep3EXS5JEVuzAzPn6FRu0SPJW7SBoc74izeGUd7500nPKB/YKOI6JyF0mH/3y1gebjUW6dWRV0FBFA5S5y3tydhS/t5OILBjBj9OCg44gAKneR87Zi+yFe29vCJ68chZku6yuZQeUucp5+8tJOBhblcdNUXSBMMofKXeQ8nJj+OHf6SF1qQDKKyl3kPPxsRS0Af6EPLUmGUbmLnKO2aJzFq2p576QRmv4oGUflLnKO/vPVBg63Rrn1yqqgo4j8CZW7yDlwdx7+w04mjCjhck1/lAykchc5Byt3nJj+WKXpj5KRVO4i52Chpj9KhlO5i5ylfUfaeHbzPm6prqRfvqY/SmbqVrmb2Rwz22pmNWZ212n2+aiZbTazTWb2H+mNKZI5Hl9dRzzhzJuhlZYkc3W5VIyZhYEHgXcD9cAqM1vm7ps77TMe+Bpwlbs3mdmwngosEqREwlm8qo4rxw6hqqw46Dgip9WdkfsMoMbdt7t7B7AYuOmkfT4NPOjuTQDuvj+9MUUyw3/VHKC+6ThzNWqXDNedci8H6jo9rk9t6+wi4CIz+4OZrTCzOaf6g8xsgZmtNrPVjY2N55ZYJECLV9UyqCiP904aHnQUkTNK1wnVCDAemAXMA35gZgNP3sndH3L3anevHjp0aJpeWqR3NLa08+ymfXzo0goKIjqRKpmtO+XeAFR2elyR2tZZPbDM3aPuvgN4nWTZi+SMJ9bWE0s4c2dUdr2zSMC6U+6rgPFmNtrM8oG5wLKT9nmS5KgdMysjeZhmexpzigTK3XlsVR3TqwYxblhJ0HFEutRlubt7DLgD+DWwBVji7pvM7F4zuzG126+Bg2a2GXge+Gt3P9hToUV624rth9hx4Bhzp+tEqmSHLqdCArj7cmD5Sdvu7nTfgb9KfYnknEUraxlQGOF9l1wQdBSRbtEnVEW60HSsg2c27uXmaeVakEOyhspdpAtLX2mgI57Q3HbJKip3kTNwdxatrGVq5UAuvmBA0HFEuk3lLnIGa3Y1UbP/KPM0/VGyjMpd5AwWrayjOD/M+y+5MOgoImdF5S5yGs3Ho/xqw25unFpOcUG3JpaJZAyVu8hpLFu3m7ZoQodkJCup3EVOY8mqOiaMKGFKeWnQUUTOmspd5BQ27z7ChoZmbpleqTVSJSup3EVOYcnqOvLDIT6gNVIlS6ncRU7SHovz5KsNvHvScAYV5wcdR+ScqNxFTvKbzfs43BrllmqdSJXspXIXOcljq+ooH9iPq8aVBR1F5Jyp3EU6qW9q5b9qDvChyyoIh3QiVbKXyl2kkyfWJBcZ+8hlFQEnETk/KneRlETCeXxNHVeNLaNycFHQcUTOi8pdJOWlNw5S33Scj07XiVTJfip3kZQlq+so7ZfHeyYODzqKyHlTuYsAza1Rntm0lw9MvVCrLUlOULmLAE++2kBHLKFDMpIzVO4iJA/JTC4fwKQLdZEwyQ0qd+nzNjY0s2n3ET6qT6RKDlG5S5/32Ko6CiIhbnq7LhImuUPlLn1aa0eMJ19p4H1TLqC0KC/oOCJpo3KXPu1X6/fQ0h5j7oyRQUcRSSuVu/Rpi1fVMXZoMdOrBgUdRSStVO7SZ72+r4U1u5qYN2OkVluSnKNylz5r0cpa8sMhPnipLhImuUflLn1SWzTO0rUNvGfScAZrtSXJQSp36ZOe2biX5uNR5ulEquSobpW7mc0xs61mVmNmd51hvw+ZmZtZdfoiiqTfopW1jBpSxMwxQ4KOItIjuix3MwsDDwLXAxOBeWY28RT7lQBfBF5Od0iRdNreeJSXdxzilumVhLTakuSo7ozcZwA17r7d3TuAxcBNp9jvm8A/AW1pzCeSdotX1REJGR/WakuSw7pT7uVAXafH9altbzKzS4FKd/9VGrOJpF1HLMETa+q57uLhDCspDDqOSI857xOqZhYC/gW4sxv7LjCz1Wa2urGx8XxfWuSs/WbzPg4e62DuDF0kTHJbd8q9Aej8k1CR2nZCCTAZeMHMdgJXAMtOdVLV3R9y92p3rx46dOi5pxY5R4tW1lI+sB/vGK/3n+S27pT7KmC8mY02s3xgLrDsxJPu3uzuZe5e5e5VwArgRndf3SOJRc5R7cFW/qvmALdMrySsE6mS47osd3ePAXcAvwa2AEvcfZOZ3WtmN/Z0QJF0WbyqlpDBR6p1IlVyX6Q7O7n7cmD5SdvuPs2+s84/lkh6tUXjLF5Vx+yLh3NBab+g44j0OH1CVfqEX67bzaFjHcy/siroKCK9QuUuOc/dWfjSTi4a3p+ZY/WJVOkbVO6S89bsamLT7iPcemWVLu0rfYbKXXLewy/tZEBhhJunaY1U6TtU7pLT9jQf55mNe7lleiVF+d2aPyCSE1TuktMeXVFLwp1PzKwKOopIr1K5S85qi8ZZtLKW2ROGUzm4KOg4Ir1K5S4566n1ezh4rIP5V1UFHUWk16ncJSe5Oz95aSfjh/XnSk1/lD5I5S45aW1tExsamjX9UfoslbvkpIUv7aKkMMIHL9X0R+mbVO6Sc/YdaePpDXu4pVrTH6XvUrlLznl0xS7imv4ofZzKXXJKa0eMR1bsYvaE4YwcoumP0nep3CWnPLqilqbWKJ+dNTboKCKBUrlLzmiLxnnoxe1cNW4Il40aFHQckUCp3CVnLFldR2NLO3dcMz7oKCKBU7lLTuiIJfj+C28wvWoQV4wZHHQckcCp3CUn/OKVenY3t3HHteP1oSURVO6SA2LxBA8+/waXVJTyzvFlQccRyQgqd8l6v1y/m9pDrdxxzTiN2kVSVO6S1RIJ54HnapgwooTrLh4edByRjKFyl6z2zKa9vNF4jNuvGUcopFG7yAkqd8la7s79z9UwZmgxN0y5IOg4IhlF5S5Z67db9rNlzxFunzWOsEbtIm+hcpes5O7c/3wNlYP7cePUC4OOI5JxVO6SlZZv2Mu6usPcPmsceWG9jUVOpp8KyTpt0Tj/sHwLE0aU8JHqyqDjiGQklbtknR++uJ2Gw8e5+88m6li7yGmo3CWr7DvSxvdeeIP3ThrOlWP1aVSR01G5S1b51jNbicWdr99wcdBRRDJat8rdzOaY2VYzqzGzu07x/F+Z2WYzW29mvzWzUemPKn3durrDPLG2nvlXVzFqSHHQcUQyWpflbmZh4EHgemAiMM/MJp602ytAtbtfAvwc+Fa6g0rf5u7c+9RmyvoXcMc144KOI5LxujNynwHUuPt2d+8AFgM3dd7B3Z9399bUwxVARXpjSl/3y/V7WLOrib9+70WUFOYFHUck43Wn3MuBuk6P61PbTuc24OlTPWFmC8xstZmtbmxs7H5K6dOOd8S5b/kWJl04gA9fpqmPIt2R1hOqZvZxoBr49qmed/eH3L3a3auHDh2azpeWHPbQ77ezu7mNu9+vqY8i3RXpxj4NQOfhUkVq21uY2XXAN4B3uXt7euJJX1ff1Mr3f/cGN0wZweVjhgQdRyRrdGfkvgoYb2ajzSwfmAss67yDmU0D/h240d33pz+m9EXxhHPnknWEDL52vaY+ipyNLsvd3WPAHcCvgS3AEnffZGb3mtmNqd2+DfQHHjezV81s2Wn+OJFu+8GL23l5xyHuuXESlYOLgo4jklW6c1gGd18OLD9p292d7l+X5lzSx23a3cz/eXYr108ewYcv0+QrkbOlT6hKxmmLxvnS4lcZVJTPP9w8ReuiipyDbo3cRXrTfU+/xrb9R/npp2YwqDg/6DgiWUkjd8kov3u9kYUv7WT+VVW88yJNlxU5Vyp3yRiHjnXwlcfXMX5Yf/5mzoSg44hkNR2WkYzg7nx96QYOt3awcP50CvPCQUcSyWoauUtGeGTFLp7ZtJevvOdtTLqwNOg4IllP5S6Be3bTXu5ZtonZE4bxP94xJug4IjlB5S6BWlvbxBcWv8KUioHc/7FpunaMSJqo3CUw2xuPctvCVQwfUMiPbq2mKF+ngETSReUugWhsaefWh1cSMuMn82dQ1r8g6EgiOUVDJel1x9pjfGrhKg60dLBowRVUlWnJPJF008hdelU0nuD2/1jLpt3NPPCxaUytHBh0JJGcpJG79Jr2WJyvPL6eF7Y28o8fnMLsi4cHHUkkZ6ncpVc0t0ZZ8MhqXt5xiLuun8C8GSODjiSS01Tu0uPqDrUyf+Eqdh08xnfmTuWmqWdagldE0kHlLj1qQ30z8xeuoiMW56efupyZY7VUnkhvULlLj3nutX3c/ugrDC7OZ9GnL2f88JKgI4n0GSp3STt35ycv7eTepzYz8cIB/PiT0xlWUhh0LJE+ReUuadVw+Dh3PbGeF7cdYPaEYXx33jSKC/Q2E+lt+qmTtHB3/mNlLf/wqy048M0PTObPZ4wkpGvFiARC5S7nre5QK3ctXc8fag5y5dgh/NOHLqFycFHQsUT6NJW7nLNYPMGilbXc9/RrAPz9zZP52IyRWtBaJAOo3OWsReMJfvFKAw88V0PtoVauHlfGfR+aQsUgjdZFMoXKXbotGk/wi7UN3P/8NuoOHWdy+QB+8Ilqrrt4mEbrIhlG5S5daovGefKVBh58oYa6Q8eZUl7K335iErNV6iIZS+Uup+TurN7VxBNr6vnV+j20tMe4pKKUe/5sEtdOUKmLZDqVu7xF3aFWnlhbz9K1DdQeaqUoP8ycySP48KUVzBw7RKUukiVU7n1cWzTOml1NvLjtAC9ua2TT7iOYwcwxQ/ji7PHMmTxCH0ISyUL6qe1j2mNxXt97lJd3HOTFbQd4ecdB2qIJIiHjslGD+Oqct3HT1HLKB/YLOqqInAeVew5rj8XZureFDQ3NbGxoZkNDM1v3thCNOwBjhxYzd/pI3jG+jMvHDKG/RugiOaNbP81mNgf4DhAGfuju9530fAHwU+Ay4CBwi7vvTG9UOZW2aJw9zW3sPHiMnQeOsSP1tfPgMRqajpNI9jil/fKYUl7KbVePYUp5KdNGDuRCjc5FclaX5W5mYeBB4N1APbDKzJa5++ZOu90GNLn7ODObC/wTcEtPBM517s7xaJym1ihNxzo4dKyDptaON+/vPdLG3iPt7GtuY++RNpqPR9/y/SUFEarKiplWOYibp1UwYUQJU8pLqRjUTydDRfqQ7ozcZwA17r4dwMwWAzcBncv9JuCe1P2fAw+Ymbm7pzFrr0kknLg78YQTSzjxePJxLJEgFndicacjniCWSBCNOdFEgmgsQUc8QUcsQXvsxG2cjliCtmiC1o44x6Nx2qJxWjtiHI8maG2P0dIWo6U9xtH2KEfbYhxtj7152ORkZjC0fwEjSgsZOaSIGaMHM6K0kOEDCqkaUkRVWTFDivNV4iLSrXIvB+o6Pa4HLj/dPu4eM7NmYAhwIB0hO1uyqo5///0bOICDJ18zdQuOk0gk9024/3GbJ/dLeHJ7PJF87sT9P96mO/Ef5YdDFOaFKMqP0C8/TFF+mJLCCOUD+1FSWEL/ggglhRFKCvMYXJzHwKJ8BhfnM6gon0FFeZT2yyMSDvVcQBHJGb16Bs3MFgALAEaOPLcFkgcV5zNhxAAwsOSfmbpNPg6l7oTest0IhYxQanvIkt8XMiMcSm0LGWH7435hMyLhEJFQclskZIRTt5FwiLywkRcOkRcOEQkbeaEQBXkh8sOdb8Pkh0P0yw9TGAmpmEWk13Sn3BuAyk6PK1LbTrVPvZlFgFKSJ1bfwt0fAh4CqK6uPqcx8rsnDufdE4efy7eKiPQZ3RlKrgLGm9loM8sH5gLLTtpnGXBr6v6Hgeey9Xi7iEgu6HLknjqGfgfwa5JTIX/s7pvM7F5gtbsvA34EPGJmNcAhkv8BiIhIQLp1zN3dlwPLT9p2d6f7bcBH0htNRETOlc7wiYjkIJW7iEgOUrmLiOQglbuISA5SuYuI5CALajq6mTUCu87x28vogUsbpIFynR3lOnuZmk25zs755Brl7kO72imwcj8fZrba3auDznEy5To7ynX2MjWbcp2d3silwzIiIjlI5S4ikoOytdwfCjrAaSjX2VGus5ep2ZTr7PR4rqw85i4iImeWrSN3ERE5g6wtdzObamYrzOxVM1ttZjOCznSCmX3ezF4zs01m9q2g83RmZneamZtZWdBZAMzs26m/q/Vm9gszGxhwnjlmttXMaszsriCznGBmlWb2vJltTr2nvhh0ps7MLGxmr5jZU0FnOcHMBprZz1PvrS1mNjPoTABm9uXUv+FGM1tkZoU99VpZW+7At4C/c/epwN2px4Ezs2tIrin7dnefBPxzwJHeZGaVwHuA2qCzdPIbYLK7XwK8DnwtqCCdFoO/HpgIzDOziUHl6SQG3OnuE4ErgNszJNcJXwS2BB3iJN8BnnH3CcDbyYB8ZlYOfAGodvfJJC+h3mOXR8/mcndgQOp+KbA7wCydfRa4z93bAdx9f8B5OvtX4Ksk/+4ygrs/6+6x1MMVJFf6Csqbi8G7ewdwYjH4QLn7Hndfm7rfQrKoyoNNlWRmFcD7gB8GneUEMysF3klynQncvcPdDweb6k0RoF9qxboierC3srncvwR828zqSI6OAxvxneQi4B1m9rKZ/c7MpgcdCMDMbgIa3H1d0FnO4FPA0wG+/qkWg8+IEj3BzKqAacDLwSZ50/8lOWBIBB2kk9FAI/Bw6nDRD82sOOhQ7t5AsqtqgT1As7s/21Ov16sLZJ8tM/t/wIhTPPUNYDbwZXd/wsw+SvJ/6esyIFcEGEzy1+fpwBIzG9Mbyw52kevrJA/J9Loz5XL3/0zt8w2Shx8e7c1s2cTM+gNPAF9y9yMZkOf9wH53X2Nms4LO00kEuBT4vLu/bGbfAe4C/leQocxsEMnfBEcDh4HHzezj7v6znni9jC53dz9tWZvZT0ke6wN4nF78tbCLXJ8FlqbKfKWZJUheR6IxqFxmNoXkG2qdmUHy0MdaM5vh7nuDytUp3yeB9wOzA157tzuLwQfCzPJIFvuj7r406DwpVwE3mtkNQCEwwMx+5u4fDzhXPVDv7id+u/k5yXIP2nXADndvBDCzpcCVQI+UezYfltkNvCt1/1pgW4BZOnsSuAbAzC4C8gn4wkXuvsHdh7l7lbtXkXzzX9obxd4VM5tD8tf6G929NeA43VkMvtdZ8n/kHwFb3P1fgs5zgrt/zd0rUu+pucBzGVDspN7XdWb2ttSm2cDmACOdUAtcYWZFqX/T2fTgid6MHrl34dPAd1InJtqABQHnOeHHwI/NbCPQAdwa8Gg00z0AFAC/Sf1WscLdPxNEkNMtBh9ElpNcBfwFsMHMXk1t+3pqbWM5tc8Dj6b+k94OzA84D6lDRD8H1pI8BPkKPfhJVX1CVUQkB2XzYRkRETkNlbuISA5SuYuI5CCVu4hIDlK5i4jkIJW7iEgOUrmLiOQglbuISA76/wAtd5A+kATJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## logistic function\n",
    "def logistic(x):\n",
    "    return 1 / (1 + np.exp(-x))\n",
    "\n",
    "\n",
    "xx = np.linspace(-8, 8)\n",
    "plt.plot(xx, logistic(xx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "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>   100</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    98</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>        <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>         <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -31.508</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Tue, 24 Sep 2019</td> <th>  Deviance:          </th> <td>  63.016</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>20:35:37</td>     <th>  Pearson chi2:      </th>  <td>  109.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>6</td>        <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</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>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -7.0444</td> <td>    1.385</td> <td>   -5.087</td> <td> 0.000</td> <td>   -9.758</td> <td>   -4.330</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    2.4211</td> <td>    0.465</td> <td>    5.207</td> <td> 0.000</td> <td>    1.510</td> <td>    3.332</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:                  100\n",
       "Model:                            GLM   Df Residuals:                       98\n",
       "Model Family:                Binomial   Df Model:                            1\n",
       "Link Function:                  logit   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -31.508\n",
       "Date:                Tue, 24 Sep 2019   Deviance:                       63.016\n",
       "Time:                        20:35:37   Pearson chi2:                     109.\n",
       "No. Iterations:                     6                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -7.0444      1.385     -5.087      0.000      -9.758      -4.330\n",
       "x1             2.4211      0.465      5.207      0.000       1.510       3.332\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Y')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcXGWd7/HPr6p639JbErI2ZN9JjKgX7+CAzgRHg6Kj4ICCo3gVxRnnOi44COg4jI644sKoMKAMLsMgclFkgNFxAQmKgaSzkaSzd3c6Se9LLb/7R1VXKr0knaSrT3X39/165ZWqU0+d832erqpfnfOcqjJ3R0REBCAUdAAREckdKgoiIpKmoiAiImkqCiIikqaiICIiaSoKIiKSpqIgIiJpKgoiIpKmoiAiImmRoAOcrpqaGq+rqws6hojIuPLss88edvfaU7Ubd0Whrq6ODRs2BB1DRGRcMbOGkbTT4SMREUlTURARkTQVBRERSVNREBGRNBUFERFJU1EQEZE0FQUREUnL2ucUzOw7wOuAJndfPsTtBnwJeC3QBVzj7r/PVh6B1q4oTe09TC0rpKI4L+g4o2aofg1cNtK+n+4YtXZFebG5nY7eOKUFEebVlp5wv/7bwQbdNtDeli62NraxaFo5s6uLB2UqjITpicUpjIRp7ugZtM7M7MCw/ejP1NjWS18swZo5lSds73SdyZid7t9iYH8yxx2gtCBCbWkBzR096b9FbWkBPbH4oO0MHM+B6x/u8kiy9v9tBj4ehtpm//r2tnTxqx1NQIhXzq9J/y1O57EzmrL54bW7ga8C9wxz+6XAgtS/lwFfT/0vWfCLrU3c+vDm9PWbXreUixZNDTDR6BiqX8AJy95w/gwefO7ACW2G6vvpjtEvtjbxkf/YSFN7LwmHcMiYWlbAbZev4KJFU9O3t3T2AVBTevy2gb7y+Da+/MSO9PUbLp7PBy5ZmM7U2RujpbOP0oIIrd1RzIyQHV9nZp87e2OYGcX54UH96M90qK03va1wyPibS5LbO11nMmYjaZ/ZbmB/3nD+DL739J70uANYxn09dT2U+nsU54fT2xk4njWlBbh7ev2Z2zrZOA6VtbM3RlN7L+7JDP2Ph7ddMJsHnztwwjb7M23cd4zPP7Y9va6Qwd++egErZ00Z8WNntJm7Z2/lZnXAw8PsKXwT+G93//fU9a3Aq9z94MnWuXbtWtcnmk9Pa1eUy7/+a8IhoyASpjcWJ55wHnjvheN6j2GofvXFEoCRH0ku6+qLsedIF3OqiinOjwzb99Mdo9auKJfd8Sv2He0injj+QhQJGTMri7n3nRdw9XeeZv+xbkKWfMlKuDNzSjE/vv7Ede5t6eLi2/8bgEgoRCyRAODH77uQD9z/h2Sbo90kEgliieMZCiIhEu5MLy8iHEr2ORIKsb2pHYAFU8uIJRLpfgBcdsev2Huki/iAp30kZDz5d686rT2GMxmzkbTPbDewP72xOA0tnQAnjMVw8sMhZlcVAXD3NRdwzd2/A5Lj6e4Yyb8dwHk1pew83DHo8sBxHCorwJ4jXUQzBtaAsAFmzKos4kBrD+5OKHW9L5Zg/9FuBnYjbHBORRGN7T2nfOycDjN71t3XnqpdkHMKM4G9Gdf3pZYNYmbXmdkGM9vQ3Nw8JuEmkqb2HgAKIuH0/+7Hl49XQ/UrFnfiiUR6WcgMnPSTa7i+n+4YNbX3EE84lnqPav1vVQ3iiQRbG9vSt4cs+c/MiCcSg9a5tbENSBaE/v/dYUPDkRP6EE7dbhx/Z2xmROPxdJ+j8QSGpZYnTujH8cwn6r/en2OkzmTMRtI+s93A/vSPxUgkX/CTL8Lux/vXv45IKJRelZnRHY2nt5V5eeA4DpV1YC4bcCHhnLDNkFnqDcxgDvREYyN67GTDuPjuI3e/E7gTknsKAccZd/qPjfbG4ul3aGbHl49XQ/Urknxrll6WcAdLvtPqbztU3093jKaWFRIOGZ56JUhtJv3ivWhaefr2/kMc7k44FBq0zkXTygGIJRLpPQUzWDu3in/7bUO6D/HUHkTmE8DdyQuHCYeSfc4Lh5KZHPLCoUH9SGY6Uf/1/hwjdSZjNpL2me2G6s+gqjaM5N6bkXDH7Hj/+sczlkik9xTcnaK8cHpbmZeHGseBWRPpB0DmtoEEFMSg8hjkNRnF3VDSG2JG2Il0hogezaOo1yjsI/1/YdQoi4eIdEcoiBoFfXDbO3pIVAx+7GRDkEVhPzA74/qs1DIZZRXFedz0uqXp455myeOj4/nQEQzdr1vWLwM4YdkNF89PH9Mdru+nO0YVxXncsn5Zek7BPXkMu6asgFvWL2V2dXH69szjwresH7zO2dXF3HDxfL78xI7kO+NU5qUzK9KZKovyONLVx5SiCG09yTmFhDs1pQV8+g3H+xyLx6kuycfM6InGB/WjP9PAOYUPXjL/tCebz2TMRtI+s91Q/fngJQtGPKdQVZIPJLczu7p40HhWlRQAnt4j6N9W5uWB4+juxFpj9O3vI36gl3/om8ET/3OAi48WkncsQVmXUdpllHUbpT0QifenG/iCHgLy6cp3evKd7gLoKXCmTi+irKqA3x9upTUUozcf8qvz+dgQj51sCHJO4S+A95M8++hlwJfd/YJTrVNzCmdOZx/p7KPMTDr7aOizj+IdcQ5taqX0sNPb0Evri13ktcRJNEbpO9BH7/5eEt2DD/1YeZhIbR6JyhDxKWHyavKonFFE6bQCYmUhOgoTFFblEy01amcUU3lOIZGyCPuOdY/J2UcjnVPIWlEws38HXgXUAI3AJ4E8AHf/RuqU1K8C60ieknqtu5/y1V5FQUTOVqIvQff2bjrrO+ne1k3X1i66t3fT/WI30aboCW1DhSHyZ+ZTMLOAghkF5M9IXs6fkU/BjILk5XPyCReHA+rNyIy0KGTt8JG7X3mK2x24PlvbFxFxd/oO9NGxsYOOP3bQ+cdOOv7YQde2Logfb5c/M5/iBcXUrK+hcF4hRecVUXhuIYV1heTV5GE2womMCWBcTDSLiIxEX2MfbU+30f5MO23PtNHxbAfRw8ff+RfMLaB0ZSk1b6yhZGkJxUuKKVpYRKRUL4X9NBIiMi65O907ujn2i2O0/qqV1l+10vNi6pTNMJQsL6F6fTWl55dSuqqUkpUl5E2ZOHNp2aKiICLjRrQlypHHjnD00aMceewIffuTZ3bl1eRRfmE5M94zg/JXlFO2piznj/HnKhUFEclZnnDan22n5SctHPnZEdo3tCc/BFYZofLVlVReUknFRRUULyqeVMf9s0lFQURyisedY784RvOPmjn848P0HeiDEJS/vJy6m+uo+vMqytaWYWEVgWxQURCRnNC1rYtDdx/i0D2H6NvfR6g4RNWlVdRcVkP1a6vJq9Z8wFhQURCRwMRaYzT9oIlDdx+i7TdtEIKqdVVMv3061a+r1rxAAFQURGRMecI5+sRRDt19iMMPHCbRnaB4STHnffY8pl01jYJzCoKOOKmpKIjImEjEEjTd30TDpxvo3tpNZEqE6ddOZ/o105NzBJoozgkqCiKSVYlogsbvNbLnH/fQvaObkpUlLLlvCTVvrCFcqMNDuUZFQUSyItGXoPHeRho+00DPzh5KV5ey/MHlVL++GgtpryBXqSiIyKhK9CY4dPchGv6pgd6GXsrWljH/S/Op/otqHSIaB1QURGRUxHviHPrOIfbctofevb2UvayMhV9fSNW6KhWDcURFQUTOSqI3wYE7D7Dntj30Heij/MJyFn17EZWvrlQxGIdUFETkjHXWd1L/tno6nuug4qIKlty7hCl/OkXFYBxTURCR0+buHLzzIDv+dgfhkjDLf7ycmvU1QceSUaCiICKnJdoSZeu7tnL4wcNUvqaSxf+2WB84m0BUFERkxI4+cZT6q+uJNkeZ9/l5zPqbWTq9dIJRURCRU0r0Jdh10y72fnYvRQuLWPHwCspWlwUdS7JARUFETqprWxeb37aZjmc7OOe6c5h/+3zCJfok8kSloiAiQ3J3Dt11iO03bCdUEGLZA8uofWNt0LEky1QURGSQ6NEo296zjeYfNjPlT6ew+J7FFM4qDDqWjAEVBRE5wbFfHqP+qnr6DvZx3m3nMfv/ztavnE0iKgoiAiS/zbTh1gYaPtNA0XlFrP7NaspfWh50LBljKgoiQvfObur/qp62p9qYfu105n95PpFSvTxMRvqri0xyh757iO3v2w4hWHr/Uqa+dWrQkSRAKgoik1SsNca267fR9L0mKl5ZwZLvLqFwriaTJzsVBZFJqHt3NxvXbaR7Rzd1t9Yx52NzCEVCQceSHKCiIDLJtD/XzvOXPk+iN8H5j5/PlIumBB1JckhW3xqY2Toz22pmO8zso0PcPsfMnjSzP5jZRjN7bTbziEx2x355jOcueg6LGKt/tVoFQQbJWlEwszBwB3ApsBS40syWDmj2CeAH7r4auAL4WrbyiEx2R35+hI1/vpGCGQWs/s1qSpaWBB1JclA29xQuAHa4+0537wPuBy4b0MaB/hOhK4ADWcwjMmm1/LSF59c/T9GiIs7/5fkUztaEsgwtm3MKM4G9Gdf3AS8b0OZm4Odm9gGgBHh1FvOITEqHHz7MpjdtomR5CaseW0VeVV7QkSSHBX26wZXA3e4+C3gtcK+ZDcpkZteZ2QYz29Dc3DzmIUXGq8M/PsymyzdRurKUVf+lgiCnls2isB+YnXF9VmpZpr8GfgDg7r8FCoFBv+nn7ne6+1p3X1tbq29pFBmJ5gea2fTmTZSuKWXlYyvJq1RBkFPLZlF4BlhgZueaWT7JieSHBrTZA1wCYGZLSBYF7QqInKXmB5vZ9JZNlF1QxqqfryJvigqCjEzWioK7x4D3A48C9STPMtpkZrea2fpUs78D3m1mfwT+HbjG3T1bmUQmg6OPH2XzWzdTtraMlT9bSaRcH0eSkcvqo8XdHwEeGbDspozLm4ELs5lBZDJpfaqV5y97nuJFxax8ZCWRMhUEOT1BTzSLyCjp2NjB85c+T/70fFb+fKUmleWMqCiITAA9DT1s/PONhEpCrPqvVRRMLwg6koxT2rcUGeeix6JsfO1G4t1x1vx6DUV1RUFHknFMRUFkHEv0Jdh0+Sa6t3ez8tGVlCzTV1fI2VFREBmn3J2t79rKsSePsfiexVT+aWXQkWQC0JyCyDi1+5O7aby3kbpb65h+9fSg48gEoaIgMg4dvOsgDZ9qYPo7pzP3E3ODjiMTiIqCyDhz7BfH2HbdNipfU8nCbyzEzIKOJBOIioLIONK9s5sX3vQCRfOLWPbDZYTy9BSW0aVHlMg4EeuI8cJlL0AClv9kOZEKnScio0+PKpFxwBPOlndsoXNzJyt/tpLi+cVBR5IJSkVBZBzYc9seDj9wmHn/Mo+q11QFHUcmMB0+EslxRx49wq5P7GLqFVOZ9aFZQceRCU5FQSSH9TT0sPltmylZVsKiby3SmUaSdSoKIjkq0Zdg01s24VFn2QPLCJeEg44kk4DmFERy1M6P7qT9d+0s/eFSihdoYlnGhvYURHLQ4R8fZt8X9jHz/TOZ+uapQceRSURFQSTHdO/qZss1Wyh9SSnz/mVe0HFkklFREMkhib4Em9+6GU84y36wjFCBnqIytjSnIJJDdn5kJ+3PtLPsR8soOk8/liNjT29DRHJE8382s++L+5h5w0xq31QbdByZpFQURHJAT0MPW67dQtlLy5j3Wc0jSHBUFEQC5nGn/u31kICl9y/VPIIESnMKIgHb+/m9tP6ylUV3LdI8ggROb0lEAtT+XDu7PrGLmstrmP4O/aSmBE9FQSQg8Z449VfVk1edx8Jv6hfUJDfo8JFIQHZ9fBddm7pY8dMV5NfkBx1HBNCegkggjj5+lH1f2MeM62dQva466DgiaSoKImMsejRK/TvqKV5crNNPJefo8JHIGNv+vu1EG6Os+PEKwsX6OmzJLVndUzCzdWa21cx2mNlHh2nzFjPbbGabzOy+bOYRCVrjfY003d9E3c11lL2kLOg4IoNkbU/BzMLAHcBrgH3AM2b2kLtvzmizAPgYcKG7HzUzfUewTFg9e3rY9r5tlP+vcmZ/ZHbQcUSGlM09hQuAHe6+0937gPuBywa0eTdwh7sfBXD3pizmEQmMJ5wt12yBOCy5dwmhiKbzJDdl85E5E9ibcX1falmmhcBCM/u1mT1lZuuymEckMPu+sI9jTx5j/pfm61PLktOCnmiOAAuAVwGzgF+a2Qp3P5bZyMyuA64DmDNnzlhnFDkrnZs72XnjTqovq2b6tfrUsuS2bO4p7AcyD5zOSi3LtA94yN2j7r4L2EaySJzA3e9097Xuvra2Vl8pLONHIpZgyzVbCJeGWfTNRfrUsuS8bBaFZ4AFZnaumeUDVwAPDWjzIMm9BMyshuThpJ1ZzCQypvZ+bi/tz7Sz8GsLyZ+mTy1L7staUXD3GPB+4FGgHviBu28ys1vNbH2q2aNAi5ltBp4EPuzuLdnKJDKWOjd1svvm3dT+ZS1T36IT62R8MHcPOsNpWbt2rW/YsCHoGCInlYgl+MMr/kDP7h5euvml5NdqL0GCZWbPuvvaU7ULeqJZZELa9/l9tG9oZ+n3l6ogyLiik6VFRlnnlk52fTL5Gwm1f6kTI2R8UVEQGUUed7Zeu5VwSZiFX9NvJMj4o8NHIqNo35f20fZUG0u+u0RnG8m4pD0FkVHStb2LXTfuovr11Ux9m842kvFJRUFkFHjC2frOrYQKQyz8hg4byfilw0cio2D/Hftp/VUri+9eTMGMgqDjiJwx7SmInKXuXd3s/OhOqtZVMe3t04KOI3JWVBREzoK7s/XdW7GwsfCbOmwk458OH4mchUPfOcSxx4+x4OsLKJxTGHQckbOmPQWRM9R7oJcdf7eDiosqmHHdjKDjiIyKYYuCmT1iZnVjF0Vk/HB3tr13G97nLPrWIiykw0YyMZxsT+Eu4OdmdqOZ5Y1VIJHxoOn7TbQ81MK5nzqX4vnFQccRGTXDzim4+w/N7KfAPwAbzOxeIJFx++1jkE8k5/Q197HjAzsoe2kZMz848BdmRca3U0009wGdQAFQRkZREJmsdvzNDmKtMRZ9ZxGhiKblZGIZtiiY2TrgdpK/lrbG3bvGLJVIjjr88GGa7mui7uY6SpeXBh1HZNSdbE/hRuAv3X3TWIURyWWx1hjb/s82SpaXMOdjc4KOI5IVJ5tT+N9jGUQk17344RfpO9jH8v9cTihfh41kYtIjW2QEjj5+lIP/epDZH5pN+UvLg44jkjUqCiKnEOuIsfVdWylaUETdrXVBxxHJKn3Nhcgp7Pr4Lnoaejj/l+cTLgoHHUckq7SnIHISx/7nGPu/sp+ZH5jJlFdOCTqOSNapKIgMI94VZ+s7t1J4biHnfea8oOOIjAkdPhIZxu5P7qZ7RzerHl9FuESHjWRy0J6CyBDanm5j7+17Oec951B5cWXQcUTGjIqCyACJ3gRbrt1CwYwC5n12XtBxRMaUDh+JDLD7U7vpqu9ixU9XECnXU0QmF+0piGRo/0M7e27bw/RrplO9rjroOCJjTkVBJCXRl2DLNVvIr81n3u06bCSTk/aNRVIaPt1A58ZOlj+0nLxK/a6UTE5Z3VMws3VmttXMdpjZR0/S7k1m5ma2Npt5RIbT/mw7DZ9pYNrV06h5fU3QcUQCk7WiYGZh4A7gUmApcKWZLR2iXRnwQeDpbGUROZl4d5z6q+rJn57P/C/NDzqOSKCyuadwAbDD3Xe6ex9wP3DZEO0+Bfwz0JPFLCLD2vXxXXRt6WLxXYt12EgmvWwWhZnA3ozr+1LL0sxsDTDb3f/fyVZkZteZ2QYz29Dc3Dz6SWXSOvrEUfZ9cR8zrp9B1Wuqgo4jErjAzj4ysxDJn/v8u1O1dfc73X2tu6+tra3NfjiZFKLHomy5ZgtFC4v0ITWRlGwWhf3A7Izrs1LL+pUBy4H/NrPdwMuBhzTZLGNlxwd30HuglyX3LiFcrO82EoHsFoVngAVmdq6Z5QNXAA/13+jure5e4+517l4HPAWsd/cNWcwkAkDzA8003tPI3BvnUn6BfklNpF/WioK7x4D3A48C9cAP3H2Tmd1qZuuztV2RU+k91MvW67ZS+pJS5n5ibtBxRHJKVj+85u6PAI8MWHbTMG1flc0sIgDuztZ3bSXRmWDJvUsI5elD/SKZ9IlmmVQOfvsgR/7fEeZ/cT4lS0qCjiOSc/Q2SSaN7p3dvPi3LzLl4inM/MDMU99BZBJSUZBJweNO/dvrIQSL71qMhSzoSCI5SYePZFLY87k9tP26jcX3LKZwTmHQcURylvYUZMJre6aN3f+wm9q/rGXaVdOCjiOS01QUZEKLtcXYfMVm8mfks/CbCzHTYSORk9HhI5mw3J1t79lGT0MPq3+xWl92JzICKgoyYR266xBN9zdx7qfPpeLCiqDjiIwLOnwkE1JnfSfb37+dKRdPYc5H5wQdR2TcUFGQCSfeHWfzWzcTLg2z5N4lWFjzCCIjpcNHMuFsv347nc93suKRFRTMKAg6jsi4oj0FmVAOfvsgh+46xNxPzKX60uqg44iMOyoKMmG0/76dbddvo/LVldTdXBd0HJFxSUVBJoTo0Sib3rSJ/Np8ltyneQSRM6U5BRn3POFsefsWevf3cv4vzye/Nj/oSCLjloqCjHt7bttDy8MtzP/KfCpers8jiJwNHT6Sca3lkRZ2fWIXU6+cyszr9XXYImdLRUHGrc4tnWy+cjOl55ey6FuL9L1GIqNARUHGpeixKC9c9gKhghDLH1xOuDgcdCSRCUFzCjLueNypv7Kenl09rHp8lX4fQWQUqSjIuPPi37/IkZ8dYeE3FzLlf08JOo7IhKLDRzKuHPjmAfbdvo+ZH5jJjOtmBB1HZMJRUZBx48jPj7Dt+m1UXVrFvNvnBR1HZEJSUZBxoWNjB5vevImSpSUsvX8poYgeuiLZoGeW5Lyehh42XrqRcFmYFY+sIFKuqTCRbNGzS3JatCXKxnUbiXfGWf0/qymcpTONRLJJRUFyVrwrzvOvf57uXd2senQVpStKg44kMuGpKEhOSvQl2PTmTbQ91cayHy5jykU69VRkLKgoSM7xuFN/dT1HfnqEhXcupPZNtUFHEpk0sjrRbGbrzGyrme0ws48OcfuHzGyzmW00s8fNbG4280ju87iz5dotNP+gmfM+ex4z3q3PIoiMpawVBTMLA3cAlwJLgSvNbOmAZn8A1rr7SuBHwGezlUdyX39BaLy3kbpP1THnw3OCjiQy6WRzT+ECYIe773T3PuB+4LLMBu7+pLt3pa4+BczKYh7JYR53tlyTLAjnfvpc6j5RF3QkkUkpm0VhJrA34/q+1LLh/DXw0yzmkRzlcaf+HfU0freRc//xXObeqKOIIkHJiYlmM7sKWAtcNMzt1wHXAcyZo0MKE0l/QWj6XlOyIHxcBUEkSNncU9gPzM64Piu17ARm9mrgRmC9u/cOtSJ3v9Pd17r72tpanYkyUSR6E9RflSoIn1FBEMkF2dxTeAZYYGbnkiwGVwBvy2xgZquBbwLr3L0pi1kkx0SPRHnh8hdo/UUr5/3zecz5e+0BiuSCrBUFd4+Z2fuBR4Ew8B1332RmtwIb3P0h4HNAKfDD1E8p7nH39dnKJLmh+8VuNr52Iz27e1jyvSVMe9u0oCOJSEpW5xTc/RHgkQHLbsq4/Opsbl9yT+uvW3nhDS/g7qx6fBVTXqlPKovkEn1LqoyZxvsbee6S54hURljz2zUqCCI5SEVBss7dafjHBuqvrKf8gnLW/HYNxQuKg44lIkPIiVNSZeJK9CXY9p5tHLr7EFP/aiqLv72YUIHei4jkKhUFyZro0Sib3rSJY08eY+4n51L3yTpSJxSISI5SUZCs6NrexQvrX6D7xW4W37OY6VdPDzqSiIyAioKMusb7Gtn2nm1YvrHqsVX6LQSRcURFQUZNvCvO9hu2c+jbh6h4ZQVL7ltC4Wz9fKbIeKKiIKOi/fft1F9dT1d9F3NunEPdzXWEIppQFhlvVBTkrCRiCfb80x4abm0gb2oeKx9dSdVrqoKOJSJnSEVBzljnpk62vHML7b9rZ+qVU1nw1QXkVeUFHUtEzoKKgpy2RG+Chs80sOef9hAuD7P0+0uZ+papQccSkVGgoiCn5ejjR9l2/Ta6t3Yz7appzLt9Hvm1+UHHEpFRoqIgI9Kzp4cXP/wizT9opnBeISt+uoLqddVBxxKRUaaiICcVa4+x93N72fu5vWBQd3Mdsz8ym3BhOOhoIpIFKgoypERvggP/eoCGWxuINkepfWst8z47j8I5+tyByESmoiAniHfFOfivB9nzuT307e+j4qIK5n12HuUXlAcdTUTGgIqCABBri7H/a/vZd/s+os1RKi6qYPFdi6l8daW+xE5kElFRmOSiR6Ls+9I+9n95P7FjMarWVTHnxjn6ARyRSUpFYZLq+GMH+7++n8bvNpLoTFDzxhrmfHwO5Wt1mEhkMlNRmEQSvQmaf9TM/q/tp+03bYQKQ0y9ciqzPjSL0uWlQccTkRygojDBuTttv22j8XuNNH2/iVhLjKIFRcz7wjymv2M6eZX6WgoROU5FYYLq3NJJ0/eaaLyvkZ6dPYQKQ1RfVs05f30OlZdUYiFNHovIYCoKE4THnban2jj8k8O0/KSFrs1dEILKSyqpu6mOmjfWECnXn1tETk6vEuNYrC3GkZ8foeUnLRx55AjRw1EsYlRcVMGM62ZQ+5ZaCs4pCDqmiIwjKgrjSLwnTvvv2jn238c4+vhR2n7bhkedSFWE6tdWU/36aqr+vIpIhf6sInJm9OqRw3oP9dLxbAetv2ml9X9aaXu6De9zMChdU8qsD82i+i+qKX9FuX7lTERGhYpCDnB3evf20v77djp+35H+v+9gX7JBGMrWljHrhllU/EkFFRdW6MdsRCQrVBTGWLwrTtfWLrrqu+jY2JEuArGWWLJBCEqWllD5mkrK1pRRuqaU0tWlREr1pxKR7NMrTRbEOmL07O6hZ2cP3S92J/9t76Z7Wzc9DT3gyXaWZ5SsKKH2jbWUrimlbE0ZJStKCBfra6lFJBgqCqfB3Yl3xOnd13v8397Bl2PHYicQtYhNAAAMsUlEQVTcL1wepnhhMeWvKGf6O6dTvLiY4iXFFC8oJlSguQARyR1ZLQpmtg74EhAGvuXutw24vQC4B3gJ0AK81d13ZzNTpkRfgtixGNEjUWJHUv+3xIi2RIkejhJtjtLX1EdfYx/Rxih9jX0kuhOD1pM3NY+CWQUUnldIxUUVyct1hRSdW0TR/CIiVRF906iIjAtZKwpmFgbuAF4D7AOeMbOH3H1zRrO/Bo66+3wzuwL4Z+Ct2coEsP+/WtjxV/V4Wxx6fPiGEQhX5WE1EYqmF1DxymLypuWRmBKmuzpEzbwSqhaUUDCjYNC7/dauKC82twNx5hUB3TGa2nsojITpicWZWlbI/qNdbGg4wtq5VSydWTFsjNauKE3tPUwtK6SiOC99PXNdFcUnTjrvbeni93uOUFlcwKrZU06438D2I1ne1h1la2MbMyqKyIuETtg2kG6Xebn/PoumlTO7unjE/ctcNnA7yTE1aksLhu37SLdxJu0yb8/s66kyHH88GPNqS0eU+Wz7kq37j4WgMwa9/aBlc0/hAmCHu+8EMLP7gcuAzKJwGXBz6vKPgK+ambn7SV6tz9wvtjbxuUc2csEM6D7XiZUYsZIQXUXOfo/SWeh0FjkdRY6VhehLdFFTWkBxfjc3vW4pALc+vBm6gOfhprlLuaigaNA2PvIfG2npTJ45VFoQoTAvjLvT0tlHTWkBrV19dEWP73FcvnoGt7919ZB5b334+HC94fwZPPjcATp7Y+l1FeeHuel1S7lo0VQAvvL4Nr7wX9tJpEawsjiPd15Yx4PPHUivp7/9wPUPtbylo5f23jjuTsKhJD9EX9ypKS3A3TEzivPDdPbG0pf779P/TRo3XDyfD1yy8JT9yxzjzD66O72xBB29MeKJ5Danlg3u+3B/86H6eLrtMm/P7OvJ1tl/v8zHQ01pAbddvuKkmc+2L9m6/1gIOmPQ288FlqXXX8zszcA6d39X6vrVwMvc/f0ZbV5ItdmXuv5iqs3h4da7du1a37Bhw2nnae2Kctkdv2L/sW5CZrg7fXEnknrhig0xDHkhIxQyZlUWEU84YORHjIJImN5YnHjCeeC9F57wDne4bYRCIRzHE57elhn0D/8jH3jlCXsMrV1RLv/6rwmHktvr6oux50gXM6cUcaC1B3cnZMlsAA+890LauqP86eefJJaqN0ZyTjsSMuZWF1OcH0nnvvuaC7jm7t+l1z/UcsPY2th+wrr6x8Xs+PXzakrZebgDgLlVJexoTl4uzAulxg2e+NCrTthjGNi/3licvlgCMMIh2Hu0O1l0gIQ7sQTkhSEWT243PxxidtXxvg/3zn7gNgb+zUbSLvP2SCjE9qbkmCyYWkYskRhynUM9Hkj1ZeaUYn58/dCZhzPSvmTr/mMh6IxBbz/bzOxZd197qnbjYpbTzK4zsw1mtqG5ufmM1tHU3pN8l4kRMjt+jD/jxW2I7eJAyIxY3IknEhREku8OCyJh3JPrHdk2nEgodOK2PFkYADY0HBmUt387pDLgJPcAnPS6ksUn2X5rY1v/5pLrtfRm0i9K/bn72w7sT+by7mj8pOPSf7k7GsdI9rezL5bO4P05M9Y7XP8KIuH0GPf3ddB4pTqULFB+Qt+HMtQ2hmp/qnaZt0fjiXRfo/HEsOvsv1/m46H/MRFPJIbNPJyR9iVb9x8LQWcMevu5IptFYT8wO+P6rNSyIduYWQSoIDnhfAJ3v9Pd17r72tra2jMKM7WskHDIcJyEO+k9JE+/dg6S+U41EjbCoRC9seQLZW8sjhnp48un3oYRSyRO3FbGnsLauVWD8vZvh1QGjOQhGSO9roR7OseiaeX9m0uu19ObSd4/I3d/24H9yVxelDf0qbH949J/uSgvnNwLcqckP5LOYP05M9Y7XP96Y/H0GPf3ddB4pTqU3GezE/o+lKG2MVT7U7XLvD0vHEr3NS8cGnad/ffLfDz0PybCodCwmYcz0r5k6/5jIeiMQW8/V2SzKDwDLDCzc80sH7gCeGhAm4eAd6Quvxl4IlvzCRXFedyyfhnVJfnJJyfJ4+215YVUlw7+0riS/BAYVBbnA3DL+mXcsn4p8YTTmTq2fdPrlp6wW3mybVSm2tWUFVKclxz2/p5evnrGoMnmiuI8bnrd8e1B8th8OGRUFiXX1Z+tP8fs6mI+eMmC9LH8/u1/8JL5ACfknl1dfML6h1oejSeYUhRJHUo6cVyqSgqoLsmnprSAaDyRvpxwT98nFvd07oGTzQP7F094eoyBdB+rSgqoKS1IjV/ysFU4ZFSVnNj34f7mQ/VxYPtTtcu8vScaT/e1Jxofdp1DPR4SnpyLuWX98JmHM9K+ZOv+YyHojEFvP1dkbU4BwMxeC3yR5Cmp33H3fzSzW4EN7v6QmRUC9wKrgSPAFf0T08M50zmFfgPPBAHSZ7rsaungaGeUBVNLB51lM/DMmJOdmXCybejso5P3L3OZzj46875k6/5jIeiMQW8/W0Y6p5DVopANZ1sUREQmowk10SwiImNDRUFERNJUFEREJE1FQURE0lQUREQkTUVBRETSVBRERCRt3H1OwcyagYazXE0NMOyX7gVIuU6Pcp2+XM2mXKfnTHLNdfdTfk/QuCsKo8HMNozkQxxjTblOj3KdvlzNplynJ5u5dPhIRETSVBRERCRtshaFO4MOMAzlOj3KdfpyNZtynZ6s5ZqUcwoiIjK0ybqnICIiQ5iwRcHMvmNmTanfgR7qdjOzL5vZDjPbaGZrciTXq8ys1cyeS/27aYxyzTazJ81ss5ltMrMPDtFmzMdshLnGfMzMrNDMfmdmf0zlumWINgVm9v3UeD1tZnU5kusaM2vOGK93ZTtXxrbDZvYHM3t4iNvGfLxGmCvI8dptZs+ntjvoNwOy8pz01E8ETrR/wJ8Aa4AXhrn9tcBPSf5a5cuBp3Mk16uAhwMYr3OANanLZcA2YGnQYzbCXGM+ZqkxKE1dzgOeBl4+oM37gG+kLl8BfD9Hcl0DfHWsH2OpbX8IuG+ov1cQ4zXCXEGO126g5iS3j/pzcsLuKbj7L0n+mttwLgPu8aSngClmdk4O5AqEux9099+nLrcD9cDMAc3GfMxGmGvMpcagI3U1L/Vv4ATdZcC/pS7/CLjEzIb7SfCxzBUIM5sF/AXwrWGajPl4jTBXLhv15+SELQojMBPYm3F9HznwYpPyitTu/0/NbNlYbzy1276a5LvMTIGO2UlyQQBjljrk8BzQBDzm7sOOl7vHgFagOgdyAbwpdbjhR2Y2O9uZUr4I/D2QGOb2QMZrBLkgmPGCZEH/uZk9a2bXDXH7qD8nJ3NRyFW/J/lx9FXAV4AHx3LjZlYK/AfwN+7eNpbbPplT5ApkzNw97u7nA7OAC8xs+Vhs91RGkOsnQJ27rwQe4/i786wxs9cBTe7+bLa3dTpGmGvMxyvDK919DXApcL2Z/Um2NziZi8J+ILPiz0otC5S7t/Xv/rv7I0CemdWMxbbNLI/kC+/33P2BIZoEMmanyhXkmKW2eQx4Elg34Kb0eJlZBKgAWoLO5e4t7t6buvot4CVjEOdCYL2Z7QbuBy42s+8OaBPEeJ0yV0Dj1b/t/an/m4D/BC4Y0GTUn5OTuSg8BLw9NXv/cqDV3Q8GHcrMpvcfRzWzC0j+jbL+QpLa5reBene/fZhmYz5mI8kVxJiZWa2ZTUldLgJeA2wZ0Owh4B2py28GnvDU7GCQuQYcc15Pcp4mq9z9Y+4+y93rSE4iP+HuVw1oNubjNZJcQYxXarslZlbWfxn4M2DgWYuj/pyMnM2dc5mZ/TvJs1JqzGwf8EmSk264+zeAR0jO3O8AuoBrcyTXm4H3mlkM6AauyPYTI+VC4Grg+dTxaICPA3MysgUxZiPJFcSYnQP8m5mFSRahH7j7w2Z2K7DB3R8iWczuNbMdJE8uuCLLmUaa6wYzWw/EUrmuGYNcQ8qB8RpJrqDGaxrwn6n3OxHgPnf/mZn9H8jec1KfaBYRkbTJfPhIREQGUFEQEZE0FQUREUlTURARkTQVBRERSVNREDkLlvwW111mVpW6Xpm6XhdsMpEzo6IgchbcfS/wdeC21KLbgDvdfXdgoUTOgj6nIHKWUl/D8SzwHeDdwPnuHg02lciZmbCfaBYZK+4eNbMPAz8D/kwFQcYzHT4SGR2XAgeBnPimVJEzpaIgcpbM7HySXzz3cuBvx+LHmkSyRUVB5Cykvp316yR/52EP8DngX4JNJXLmVBREzs67gT3u/ljq+teAJWZ2UYCZRM6Yzj4SEZE07SmIiEiaioKIiKSpKIiISJqKgoiIpKkoiIhImoqCiIikqSiIiEiaioKIiKT9fxqX68YzFWbQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(5)\n",
    "n_sample = 100\n",
    "a = 3.0\n",
    "b = -9.0\n",
    "\n",
    "x = uniform(1, 5, size=n_sample)\n",
    "x = np.sort(x)\n",
    "\n",
    "q = logistic(a * x + b)\n",
    "y = binomial(n=1, p=q)\n",
    "\n",
    "exog, endog = sm.add_constant(x), y\n",
    "\n",
    "# Logistic regression\n",
    "mod = sm.GLM(endog, exog, family=sm.families.Binomial(link=sm.families.links.logit()))\n",
    "res = mod.fit()\n",
    "display(res.summary())\n",
    "\n",
    "y_pred = res.predict(exog)\n",
    "\n",
    "idx = x.argsort()\n",
    "x_ord, y_pred_ord = x[idx], y_pred[idx]\n",
    "plt.plot(x_ord, y_pred_ord, color='m')\n",
    "plt.scatter(x, y,  s=20, alpha=0.8)\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "# plt.savefig('logistic_regression.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## custom GLM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(5)\n",
    "n_sample = 100\n",
    "a = 0.6\n",
    "b = -0.4\n",
    "sd = 0.8\n",
    "\n",
    "x = uniform(1, 5, size=n_sample)\n",
    "mu = np.exp(a * x + b)\n",
    "y = normal(mu, sd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Y')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuQnFd55/Hv0z13aTySpZFiXWwZY7TYmAjvYGANNgVxCrBLZh1CnISANhAVhF1MoEJwqNhl1+7GqQRXyGVNqWwiCJdAQIBWGIIqJPLaxdqMjNbGMkI2NrYlWzOSkDSa0Vy6+9k/unvcandPX6bfS/f7+1S55vZK79ELc573POc555i7IyIiyZWKugEiIhItBQIRkYRTIBARSTgFAhGRhFMgEBFJOAUCEZGEUyAQEUk4BQIRkYRTIBARSbiuqBtQj5UrV/qGDRuiboaISFvZu3fvUXcfrnVdWwSCDRs2MDo6GnUzRETaipn9op7rlBoSEUk4BQIRkYRTIBARSTgFAhGRhFMgEBFJOAUCEZEInJya4+CRCU5OzUXdlPYoHxUR6SR7Doxx+67981/fct0lXL1xVWTt0YhARCREJ6fmuH3XftIpo7crTSbr3Lpzf6QjAwUCEZEQjU1MAzCbyfHk+GleODXNoRNTfOfRw5G1SYFARCREqwb7yObg0IkzZ31/230/P2tUEOYcggKBiEiIhga62XrVhQC4gwNrhvpJmc2PFvYcGOOGux7gg1/cyw13PcCeA2OBtkmBQEQkZNdetoa1y/o5b1kfLx9eSm93CrP8aKF0DmGgp4t0yrh9V7BzCAoEIiIhGxro5rbNl9KVMqbnsmRzzi3XXcLQQPf8qKC3Kz3/0f3FuYUgqHxURCQCV29cxY71yxmbmGbVYB9DA91AflQAMJPJ0tuVZiaTnR8tBEUjAhGRCJycmntJEID8aOGW6y4hm3MmZzJnjRaCohGBiMgiVOvQF1JrQVm10UJQFAhERJrUzArh8gVlM5kst+/az471y18yMgg6ABQpNSQi0oRmq3uimAyuJbBAYGafM7MxM/tJyff+0sx+amaPmNk3zWxZUPcXEQlStQ79yfHTZy0EK18YVjoZXPwY9GRwLUGmhrYDfwd8oeR7u4Gb3T1jZn8B3Az8SYBtEBEJRKXqnjNzWT72tX2kUwbAOzet4Vv7Xtw6opg6uuW6S7h9134mZzKYEfhkcC3m7sH95WYbgF3u/qoKP/vPwLvc/Xdr/T0jIyOuw+tFJG6KcwTukHOYnstwTn83vV1ppmYzPHN8ivPPHWCgp4uZTH69wI4PXcnQQHdTk8yNMrO97j5S67ooJ4t/H/hqhPcXEVmU0uqeiek5/vjrj8ynilJm4IWP5FNHkzMZxiam5yeCoxwFlIokEJjZp4AM8KUFrtkKbAU4//zzQ2qZiEhjih16cQ6gmCrKuYOR/0g85gKqCb1qyMy2ANcBv+sL5KXcfZu7j7j7yPDwcGjtExFpRvlCMICPvOXlAA0vDAv79LJQRwRm9jbgE8DV7j4V5r1FRIJWaSHYe99wYUNzAVGcXhZk+ehXgB8CG83sOTN7P/kqokFgt5ntM7PPBnV/EZEoDA10c/HqwflOv/zrhUSx8ygEOCJw99+u8O17grqfiEi7q7Q2oXSCOShaWSwiEhNRLTZTIBARiYkodh4FbTonIhIrYe88CgoEIiKxE/ZiM6WGREQSToFARCThFAhERBJOgUBEJEJhbydRiSaLRUQCVm3L6Si2k6hEgUBEJEDVOvt6zy4Og1JDIiIBWWjvoDidXaxAICISkIU6+zidXaxAICISkIU6+6i2k6hEcwQiIgEpdvbVDqqPYjuJShQIRKSjhHEofCNqdfZxOLtYgUBEOkZcyjHLxaGzX4jmCESkI0R1ulcnUCAQkY4Qp3LMdqNAICIdIQ7lmHHYLqIZmiMQkY5Qq0InaHGdn6iHAoGIdIyoyjHjtF1EMwJLDZnZ58xszMx+UvK9c81st5kdLHxcHtT9RSSZhga6uXj1YKgdcLvPTwQ5R7AdeFvZ9z4J/Ku7Xwz8a+FrEZG2Fof5icUILBC4+33A8bJvXw98vvD554F3BnV/EZEwvf+NFzKbiX67iGaEPUew2t2fL3z+ArA65PuLiLTU2ZPEzgfe9DKuvWxN2wQBiLB81N0d8Go/N7OtZjZqZqPj4+MhtkxEkmKx5Z7li9h6ulLcc/9TLW5l8MIeERwxs/Pc/XkzOw8Yq3ahu28DtgGMjIxUDRgiIo3sL1S89uCRCT69+2fz32+m3LPSJPHkTIaxiem2GhGEHQh2Au8D7ih8/HbI9xeRDtNI/X7x2mzOOXTiDMNLe1ixtK/pcs/SSeJi2Wg7TRIXBVk++hXgh8BGM3vOzN5PPgBcY2YHgV8rfC0i0pRG9hcqvbY7ne/6jk7Okc150+WecTpTYDECGxG4+29X+dFbg7qniCRLI6mZ0muzOSdlRs6duWyOTA7MoK8rzcEjE/R1pZnOZOtKNcXlTIHF0MpiEWlbjaRm+rrSzGZyQIaBni5WLu1hbGKG2UyOrrTxzk1r2LL9ISZnMhybnGXl0l4GetJ1zR3EfZvpWrTpnIi0ldJKn3pTM3sOjLFl+0PMZXM8c3yK8YkZ+rrTfOa3XsPd7xth+5Yr+Na+wwAcL6SVjk/OAiRiK2uNCESkbVSbGF4oNVM6NzA82MeS3i5mMzm2b7mC9SsGADh4ZAKAlBk4dKVSZD2fPsrmvO2qgBqlEYGItIWFJoYX2l+ofB5hoKeL7nSK6cJ2EPBiKinnDgaZXA4rfN2OVUCNUiAQkbbQ7MZu9ewDVEwxASzvzweT5QM9QLhbWUdFqSERaQvN1uzXe05BaYqpkaqhTqBAICJtYTEHz9Rb4tnu1T/NUiAQkbaxmJr9pHby9VAgEJG2og699TRZLCKScAoEIiIJp0AgIrG02LMCpH6aIxCR2Glka2lZPI0IRCRWGtlaWlpDgUBEYqXZFcTSPAUCEYmVeraEkNZSIBCRWOmUU7/aiSaLRSR2mllB3MgB9nI2BQIRiaVGVhCrymhxlBoSkbamKqPF04hARNpGafoH8pVEE9MZoL4D7KUyBQIRaQul6Z/TMxmy2RxL+/Id/UwmRzplDZ1TIC+KJDVkZn9kZo+Z2U/M7Ctmpv/FRKSq0vTPbCbHkVMzHJ2c4/CJaWYzWdyd2cziq4ySuq1F6CMCM1sLfAS4xN3PmNnXgBuB7WG3RUTaQ3ExWVcqxQunJgEwwHGOTs5x3jl93PnuVzPY19101VCSJ5yjmizuAvrNrAsYAA5H1A4RaQPFNM/kbAbDsML30ynD3ZnJZBle2lf1APtakj7hHHogcPdDwF8BzwDPAyfd/fvl15nZVjMbNbPR8fHxsJspIi22mLRLcZFZcRSQTuWDwFwmx1zWcWDL9ofYc2CsqbYlfVuL0AOBmS0HrgcuBNYAS8zsPeXXufs2dx9x95Hh4eGwmykiLbTnwBg33PUAH/ziXm6464GmOuyrN67i2x9+I7dtvpR1y5ewZlkfZsavnNPLqsG+Rb3FJ31biyhSQ78GPOXu4+4+B+wA/lME7RCRELQy7TI00M3vvO4Cvv3hK7nluktYt7yf4UJnvZi3+KRvaxFF+egzwOvNbAA4A7wVGI2gHSISgkppl8XW+Q8NdPPaDStIp4yZTLYlZaPNbGvRKaKYI3gQ+DrwMPBooQ3bwm6HiIQjqLRLEG/xQwPdTU84tzNz96jbUNPIyIiPjmrQINKuiqWZ7mDW2tJMbTZXnZntdfeRWtdpZbGIBK6RtEujHXsjm9NJZQoEIhKKejrsJC/qipJ2HxWRWEj6oq4oKRCISCwkfVFXlBQIRCQW+rrSzGZyTM3mt5VO2qKuKGmOQEQiV5wbmMvmeOHUNCuW9LKkN52oRV1RUiAQkUiVzg0MD/axpLeL2UyO7VuuYP2KgaiblwhKDYlIU1q1d3/53MBATxfd6RTThQVoEjyNCERiLo4LplpZ5lm68lgnjEVDIwKRGGvFrp2t1uoyz6Rv+BYHGhGIxFRph1t8U7591352rF8eaScZxCZySd7wLQ4UCERiKogOtxWCSuVoq4joKDUkElNxPSxFqZzOoxGBSEwVO9zbd+1nciYzv2tnHDpcpXI6iwKBSIzFucNVKqdzKBCIxNxiO9w4lp9KvFQNBGZ2L/CH7v50eM0RkVbSts5Sj4Umi/8B+L6ZfcrM9Boh0ma0rbPUq+qIwN3/2cy+C/wZMGpm/wjkSn5+ZwjtE5EmxbX8VOKn1hzBLDAJ9AKDlAQCEYk3bd0g9VpojuBtwJ3ATuByd58KrVUismhxLj+VeFloRPAp4Dfd/bFW39TMlgF3A68CHPh9d/9hq+8j0kmaqf6Jc/mpxMdCcwRvCvC+nwG+5+7vMrMeQJuOiyxgMdU/qveXWkLfYsLMhoCrgHsA3H3W3U+E3Q6ROKhnT39V/0jQolhQdiEwDvyDmf0qsBe4yd0nI2iLSGTqfctX9Y8ELYpN57qAy4G73P015KuSPll+kZltNbNRMxsdHx8Pu40igWrkLb8Vm8+16jQx6UxRBILngOfc/cHC118nHxjO4u7b3H3E3UeGh4dDbaBI0Cq95bu/+P1Si93ts/Rwm+v//gG+/OAvFBDkLKGnhtz9BTN71sw2uvsB4K3A/lp/TqSTNFrj32z1T+nIYzaT49CJM9y68zG23fdzbtt8qbabECC68wj+G/AlM3sE2AT8z4jaIRKJZt7yhwa6uXj1YEPzAsURRlcqxfMnp0mZYRgOmnCWeZHsPuru+4CRKO4t0oggd+4Mo8a/OMKYnMm8+E2DJT1dTM9lNeEsgLahFqkqjJ07g67xL448bt25n5w7AGuG+snkctpuQubpqEqRCjqpdv/qjav49oev5LbNl7J2WT/plOl4STmLRgQiFXRa7f7QQDe/87oLuPayNdpuQl5CIwKRCuJ6cPxi6KQyqUYjApEK4rJzZ6s6b51UJgtRIBCpIuqdO5vpvCsFjtL5juKahdt37WfH+uUaGQigQCCyoKh27qy38y7t+Pc9+8uKgaPT5juk9RQIRCJULfVTT+ddOmLI5mB6LsM5/d0vCRw6qUxq0WSxSERK9wC64a4H2HNgbP5ntSary8tb3Z1jk7N0pfK/0qV7F5WuYj51Zo4zs1k+fs0rNBqQeQoEIhGotU6h1hYU5SOGJb35wf3kbH4FcXnguHrjKj5+zUZy7qRTxqd3/+yswCPJptSQSABqVfvUk/pZaLK6PN2TyeVYubQXg4pVTien5vj07gP096Q1YSwvoUAg0mL1VPvUm7evNlldqbz1jhsuY1OVwKEJY1mIAoFIE6q98ddb7dOKdQrVRgyV/g5NGMtCFAhEGrTQG38jb96tWKdQb3lrXBbISTwpEIg0oNYbf6Nv3vV25K1YYRz1AjmJLwUCkQbUeuMP4s27ldtDRLVATuJNgUCkAfW88bfyzVvbQ0gYtI5AEuvk1BwHj0w0dMZAvUdMNnOsZCWNHHIv0iyNCCSRFpNuCTPXrmofCYNGBJI4rTh9rFVv/PXcp9FD7kUapRGBtEy7HHzSDourSp+lqn0kaAoE0hLtdPBJ3NMt1Z6lAoAEJbLUkJmlzezHZrYrqjZIa7TbQe9hpVuamYxut2cpnSHKEcFNwOPAORG2QVqgHVIt5UrTLX1daaYzWU5OzbWsvc2OkNrxWUr7i2REYGbrgGuBu6O4v7RWux70PjTQzeETZ9iy/aGKZwI0azFv9e36LKW9RZUa+mvgE0Cu2gVmttXMRs1sdHx8PLyWScPatbIlqDTM2MQ02ZzjDtmcN1T7367PUtpb6KkhM7sOGHP3vWb25mrXufs2YBvAyMiIh9Q8aVI7VrYElYY5eOQ0h06cAcDMWLmkm/6errrf6tvxWUp7i2KO4Epgs5m9A+gDzjGzL7r7eyJoi7RQu+1jE0T1UPEAmFWDvRw9PUvOnfHTs3zmtxp7q2+3ZyntLfTUkLvf7O7r3H0DcCPwAwUBiUIQaZjiKOPcJb1cNLyUDSuWsHZZPxevHmxVs0VaTusIEqZdFn2FpdVpmJceIQnplGmyV2It0kDg7v8O/HuUbUiSdlr0FWbAamUaRgfASDvSiCAh2mk747gGrHqDkyZ7pd0oEAQkbimYdlmoFNeA1Whw0mSvtBMFggDE8Y027vvrFMUxYMU1OIm0irahbrG47hUT94VKxX15+goBIE4ra3U4jHQ6jQhaLI5vtEVxzV2Xj6DeuWkN39p3uKHJ1iBTcY2OpuKWFhSpRYGgxeKegolb7rpS2uVb+w6zfcsVTGeydXWmQafiGqkEimNaUKQWBYIWU/lgY6qNoKYz2aqLsErfuIFQ8vf1jKY0lyDtSoEgAHFNwcRRoyOo8jfu97/xQiCcVFz5aKo8BRTntKDIQhQIAhK3FEwt5W/Zpfv0BxnMGhlBVXrj3nbfU4CHnoqrlALatH45EN+0oEg1CgRyVqc2NZvFPb/Z67HJWVYu7WWgJx1orrveQ2IqvXFnshk+8KaXcc/9T4WWiquaAvrQlUoLSltSIEi40k6tK5XimeNTZ/38+OQsAz39gee6hwa62ffsLxecaK2WRrr2sjVce9ma0FJxT46fZjaTY7Avf5/SFJDSgtKOtI4g4UrfsueyOVJm+R84dKVSOJAyC7xuvp71FwuthRga6Obi1YPzHe/JqTke/sUvefgXx1u6hmPPgTE+9rV9vHBqmoNjE0xMz70kBVTeFpG404gg5oKuSS99y+5Op8gV0kIYZHL5wJBzDzzXXe9Eaz1v3HsOjPHJHY9y9PQMACuW9PAXv/HqRae2isGqpyvFmqF+Dp88w6ETZ1i7bIDbNisFJO1LgSDGwqhJL52szWSzrFzai3v+mMXjU7MsH+iZv3eQHV0j1UMLTcSfnJrj1p37OXp6Zn50c2xyllt3Psa3P7y41FZpsOrtSrOkt4tTZ+a4892v5vILzm367xWJmgJBQdxWgy62Jr2Rf0/5WzaEVzVU1Kr1F/nzgvMjmWIgcIdM1hddxvnSswZy9HanuGhYh85Ie1MgIJ6rQRdTk97Mv6f8LTuKYLhQ2qfewLZqsI90qiTFBThOV3rxh8NosaB0qsQHgoXevIHIRgnNblXR7qtbK6V9GglsQwPd3Lb5kpfMEdy2+dKW/PtVFSSdKPGBoNqb93cePcw99z81f13Yo4Rm3z7jvrq10RRcM4Ht6o2r+N5NV/Hk+GnAuWi4tRU87bZYUKSWxAeCSm/eOXe23fcUPV3RvlXXu9CqVJw3vWsmZdVsYBsa6ObyC5a3oNUinS/x6wgq1aZvveplpFPx2H9+aKCbwyfOsGX7Q3zwi3u54a4H2HNg7Kxrinv5F4NEkOcOlN6r0T/XzDkNpYGt+DEugU2kUyR+RACVq2buuf+put6qy1Mdra4+qpUaqfaWHUQeezGT6ot5s9cErUiwQg8EZrYe+AKwGnBgm7t/Jux2lCvP+9bT+VQ7UKX071jsvMJCHSgwf//iwq/SINHKznKxk9CLSVlpglYkWFGMCDLAx939YTMbBPaa2W5331/rD4apVudT3jFOzWb4mx88wfnnDjDQ09WyeYWFOtCxiWkmZzIcn5rLh1SD5f3dgUwMVwpIp87M8aOnj/HaDStq3m+xb/aaoBUJTuiBwN2fB54vfD5hZo8Da4FYBQJYuPMp7xhTZuDML2IqT300mzJaqAM9dWaOY5OzQH5foEwux/Gp2flzf1upPCAdOz3N+OlZ/vt3HiedMj5+zUYuXr10wX+f3uxF4inSOQIz2wC8BngwynY0o7xjzLmDMb+QqfTN/d5HDvPn3/0p6ZSRTlnDKaNqHeh0Jr8lxPHJWbLupMxYPtDDdGFitZVKA9KpM3OMn55l1WAvg33dHJ+c4aav/pi1y/pr/vv0Zi8SP5FVDZnZUuAbwEfd/VSFn281s1EzGx0fHw+/gTWUV+cAfOQtLwc4q1rngSeOctNX87tVHj45zZnZTF3VMpXuV76j5arBPgZ60qxb3s8F5w6wbnk/S3rTLa+oKVYKbVq/nB0fupI/u+6VrF3Wz7lLesnmnKOn86OS7nSq7mogEYmPSEYEZtZNPgh8yd13VLrG3bcB2wBGRka80jWNanVFT6U39fe+4cKzqo+u+9v/A+RTNzl3jk7OcV463ZI8fulbejbngVTUVKoUeu2GFaRTxkwmizuF3UltPhDEaQGbiNQWRdWQAfcAj7v7nWHdN6j9hCrt0VP8+uCRCdIpm6/oSZmRyeXIubfsrT3IvHs9J3Flc/kYvXJJ93xwaGWdf9w2AxTpRFGMCK4Efg941Mz2Fb73p+5+bxA3Ozk1x5PjE9y6c3/oK4XzG6AZK5f2cPT0LJlcDoCb3/4f2mLLg4VKV0sD0MEjE3x6989aXucfx80ARTpRFFVD9wMWxr2KHclsJscLp6ZZu6x/fi/5MNIXpamb84b6yOacm9/+St7x6vMCu2cr1ar9Lz0Z7MqXD4e6kE5EWqdjVxaXdiSDfd28cGqaQyfOcPGqLjK5XGjbFLRzyWQjtf+tHpXEffM8kU7SsYGgvCMpHi04MT1HT1cq1G0Kgi6ZDDKPHlUgi/PmeSKdpmMDQXlH0tudYu2yfu589yYuGl4a+ltlvZ11o516WMdZhv28tMeQSHjMvSWVmYEaGRnx0dHRhv9csZN0Z74jiWKysd7OutFO/eTUHDfc9cBZefRsztnxoSs7psNU1ZBI88xsr7uP1LquY0cEEI/8fL2Tns1MjiYhj66VyCLB6/jzCCqtyA1Tpc660tkG9V5XSnv1i0grdHwgiFq9nXUznXqtQ2iaPURGRJKlo1NDcVDvpGezk6PV0l9ajCUi9eroyeIw1ZrUDKpqqNrf0emTyCJSmyaLQ1TP23e9k56tmBxNwiSyiLSO5ggWqdlD2YOkSWQRaYQCwSI1U+0TtFqTyCIipZQaWqS4boUQhzUUItIeEjsiaFVpZZzfvqNeQyEi7SGRI4JWl1bq7VtE2lniRgRBTe520tu3FqKJJEviRgRjE9Nkc/ljI7M5V2llGS1EE0mexI0IDh45zaETZ3j62CRPjJ/m2OnpWEzuxkEcS2FFJHgdHwhK0xwnp+b49O4DrBrsJWWGuzN+epaPX/MKjQaIZymsiASvo1ND5WmO97/xQgDOXdLLUH8Pc9kcc9kcF68ejKqJsRLXUlgRCVbHjggqpTm23fdzsrl8R5dOGWaQTpk6uoI4l8KKSHAiGRGY2duAzwBp4G53v6PV96iU5shknQ+86ULuuf8pHX9YhUphRZIn9EBgZmng74FrgOeAH5nZTnffv/CfbEy1NMe1l63h2svWqKNbgE4FE0mWKFJDVwBPuPvP3X0W+Cfg+lbfZKE0RyfV/IuILFYUqaG1wLMlXz8HvC6IG9VKc5ycmuPJ8QnAuGh4qQKDiCRSbKuGzGwrsBXg/PPPb/rvqZbm2HNgjD/5xiMcm5wFYOXSXu644TItnhKRxIkiNXQIWF/y9brC987i7tvcfcTdR4aHh1vagJNTc9y68zGOTc6SMiNlxtHTM9y6U4unRCR5oggEPwIuNrMLzawHuBHYGWYDittMGDYfCMyMbC6nxVMikjihp4bcPWNm/xX4F/Llo59z98fCbMOqwT7SKcNxcj7fLtKplNYUiEjiRLKgzN3vdfdXuPtF7v4/wr7/0EA3t22+lBVLesi5k3Nn5dJebtusNQUikjyxnSwO2tUbV/EvH71aVUMikniJDQSQHxlcfsG5UTdDRCRSHbvXkIiI1EeBQEQk4RQIREQSToFARCThFAhERBLO3D3qNtRkZuPAL5r84yuBoy1sTivFtW1xbRfEt21xbRfEt21xbRfEt22NtusCd6+5R09bBILFMLNRdx+Juh2VxLVtcW0XxLdtcW0XxLdtcW0XxLdtQbVLqSERkYRTIBARSbgkBIJtUTdgAXFtW1zbBfFtW1zbBfFtW1zbBfFtWyDt6vg5AhERWVgSRgQiIrKAjggEZvY5Mxszs59U+bmZ2d+Y2RNm9oiZXR6jtr3ZzE6a2b7Cf7eE1K71ZvZvZrbfzB4zs5sqXBPJc6uzbaE/NzPrM7OHzOz/Fdp1W4Vres3sq4Vn9qCZbQi6XQ20bYuZjZc8sw+E0bbCvdNm9mMz21XhZ5E8szraFeXzetrMHi3cd7TCz1v7u+nubf8fcBVwOfCTKj9/B/BdwIDXAw/GqG1vBnZF8MzOAy4vfD4I/Ay4JA7Prc62hf7cCs9haeHzbuBB4PVl1/wh8NnC5zcCX41R27YAfxf2/9cK9/4Y8OVK/5tF9czqaFeUz+tpYOUCP2/p72ZHjAjc/T7g+AKXXA98wfP+L7DMzM6LSdsi4e7Pu/vDhc8ngMeBtWWXRfLc6mxb6ArP4XThy+7Cf+WTbNcDny98/nXgrWZmMWlbJMxsHXAtcHeVSyJ5ZnW0K85a+rvZEYGgDmuBZ0u+fo4YdCwl3lAY0n/XzC4N++aFofhryL9Flor8uS3QNojguRVSCfuAMWC3u1d9Zu6eAU4CK2LSNoDfKKQSvm5m68NoF/DXwCeAXJWfR/XMarULonlekA/i3zezvWa2tcLPW/q7mZRAEGcPk18G/qvA3wLfCvPmZrYU+AbwUXc/Fea9a6nRtkiem7tn3X0TsA64wsxeFcZ961FH2/43sMHdXw3s5sW38MCY2XXAmLvvDfpejaizXaE/rxJvdPfLgbcDHzazq4K8WVICwSGgNJqvK3wvcu5+qjikd/d7gW4zWxnGvc2sm3xH+yV331HhksieW622RfncCvc8Afwb8LayH80/MzPrAoaAY2G1a6G2ufsxd58pfHk38B9DaM6VwGYzexr4J+AtZvbFsmuieGY12xXR8yre+1Dh4xjwTeCKskta+ruZlECwE3hvYab99cBJd38+6kYBmNmvFPOhZnYF+f9NAu84Cve8B3jc3e+sclkkz62etkXx3Mxs2MyWFT7vB64Bflp22U7gfYXP3wX8wAuze1G3rSyHvJn83Eug3P1md1/n7hvITwT/wN3fU3ZZ6M+snnZF8bwK911iZoPFz4FfB8qrDlv6u9kRZxab2VfIV5GsNLPngFvJT5bh7p8F7iU/y/4EMAXPgC7CAAABfklEQVT8lxi17V3Ah8wsA5wBbgyj4yD/RvR7wKOFvDLAnwLnl7QtqudWT9uieG7nAZ83szT5wPM1d99lZrcDo+6+k3wA+0cze4J8kcCNAbepkbZ9xMw2A5lC27aE1LaXiMkzq9WuqJ7XauCbhfecLuDL7v49M/sgBPO7qZXFIiIJl5TUkIiIVKFAICKScAoEIiIJp0AgIpJwCgQiIgmnQCDSIMvvjvqUmZ1b+Hp54esN0bZMpDkKBCINcvdngbuAOwrfugPY5u5PR9YokUXQOgKRJhS2wNgLfA74A2CTu89F2yqR5nTEymKRsLn7nJn9MfA94NcVBKSdKTUk0ry3A88DsdmBVKQZCgQiTTCzTeQ3dns98EdhHXQkEgQFApEGFXY9vYv8OQnPAH8J/FW0rRJpngKBSOP+AHjG3XcXvv5fwCvN7OoI2yTSNFUNiYgknEYEIiIJp0AgIpJwCgQiIgmnQCAiknAKBCIiCadAICKScAoEIiIJp0AgIpJw/x/0VxKr+kkb3QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y,  s=20, alpha=0.8)\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "# plt.savefig('custom_GLM_data.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "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>   100</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    98</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>        <td>Gaussian</td>     <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>          <td>log</td>       <th>  Scale:             </th> <td> 0.69261</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -122.53</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Tue, 24 Sep 2019</td> <th>  Deviance:          </th> <td>  67.876</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>20:35:38</td>     <th>  Pearson chi2:      </th>  <td>  67.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>6</td>        <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</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>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -0.3115</td> <td>    0.075</td> <td>   -4.130</td> <td> 0.000</td> <td>   -0.459</td> <td>   -0.164</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.5789</td> <td>    0.018</td> <td>   32.984</td> <td> 0.000</td> <td>    0.545</td> <td>    0.613</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:                  100\n",
       "Model:                            GLM   Df Residuals:                       98\n",
       "Model Family:                Gaussian   Df Model:                            1\n",
       "Link Function:                    log   Scale:                         0.69261\n",
       "Method:                          IRLS   Log-Likelihood:                -122.53\n",
       "Date:                Tue, 24 Sep 2019   Deviance:                       67.876\n",
       "Time:                        20:35:38   Pearson chi2:                     67.9\n",
       "No. Iterations:                     6                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -0.3115      0.075     -4.130      0.000      -0.459      -0.164\n",
       "x1             0.5789      0.018     32.984      0.000       0.545       0.613\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Y')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8nGWd///XNTOZSSbnpmnpKU3pkR6gYEVYVmBl2R8oC4roorKKp36tK3hgPeCqCL/1q98V/C66LlpPCAIiWLSColVOq6LQYqVN2rSUng9J2pwzOczMfX3/mElIwySZmc7MPcm8n49HHznMnbkv7pD7c12f63Ndt7HWIiIihcvjdgNERMRdCgQiIgVOgUBEpMApEIiIFDgFAhGRAqdAICJS4BQIREQKnAKBiEiBUyAQESlwPrcbkIzp06fb+vp6t5shIjKpbNmy5bi1tnai4yZFIKivr2fz5s1uN0NEZFIxxuxP5jilhkRECpwCgYhIgVMgEBEpcAoEIiIFToFARKTAKRCIiBQ4BQIRkQKnQCAikmOdoTC7m7vpDIXdbgowSRaUiYhMFU83tXDbo43DX3/hiuVctHSGiy3SiEBEJGc6Q2Fue7QRr8cQ8HmJRC23bGx0fWSgQCAikiMt3f0ADEYc9rT2cKyrn8MdIR7bdmT4mNCuELs+sotITyRn7VIgEBHJkRnlxUQdONzRd9L31z/zMp2hME7EYfu7Gjl6XzPtx/pz1i4FAhGRHKkMFrH2wgUAWAsWmF1ZgscYWrr7eeYTOwht7uHhK6Jc+/PNPN3UkpN2KRCIiOTQm1bNZk5VCbOqillUW0agyIMxwNY+nG+2sv0c2HeuD6/HcNujuZk/UCAQEcmhymARt165Ap/H0B+OEnUsn794Kcc+tIeuSnjq7V4AAj4v1r4yr5BNKh8VEcmxi5bOYMO8alq6+5lRXszRG/Zw7MAAP/sQdPkdAngZiEQxJjavkG0aEYiI5FhnKDwcBAZ/1cGx7x+j7jN1vPfG5UQdS+9AhKhj+cIVy6kMFmW9PRoRiIikYeTNPJWb9cgFZeUd8LFveChfU079F+s5vchz0kghF0EAFAhERFKW7urgkxaUeTxc9WCUwZDD3O8uxlMUS9BUBotyFgCGKDUkIpKCkTfzoD+16p6hid+Az8trnrIs2A2/vgq6TjPZbva4shYIjDHfN8a0GGO2j/jeV40xO40xLxpjHjHGVGXr/CIi2TDyZj700VrY09pz0kZyiTaWG5r4rd4b4fW/cGhaBVtfl5sJ4fFkMzV0N/BfwD0jvrcJuNlaGzHG/B/gZuDTWWyDiEhGDd20ByJRAr5YdU9fOMonfrIVryfWs3/z6tn8bOsr20YMpY4qg0V8/uKlnLisiZ4y2Pg2+MI/5mZCeDxZGxFYa58B2kZ97zfW2qENNP4EzM3W+UVEsqEyWMQXrnilumcwYrHW4vd5CPpjfeuvP/ESwKtSR9ZaZtzRTk0bLPzBEu6/6QLXdx4FdyeL3wc8ONaLxpi1wFqAurq6XLVJRGRCI9cBdPeH+eTDLw6nijzGgI1/JJY66h2I0NLdT+/93bQ80EL9bfXUv2W2m/8JJ3FlstgY829ABLhvrGOsteuttWustWtqa2tz1zgRkSRUBotYPLOchbXlQCxVBOBYCyb+Mf59Y6D8QJSXbniJqkuqmP/Z+a61O5GcBwJjzPXAFcC7rI1fKRGRSWp0qgjgxjcsAhheGPb5NyzjwHt24a30svy+5Rhv4ioht55cltPUkDHmMuBTwEXW2lAuzy0iki2jt4yoDBbx7vMXDH995MMv0bwzxFmbzsI/05/wPdx8clk2y0cfAJ4FlhpjDhlj3k+siqgc2GSM2WqM+Va2zi8ikktDqaKhCqChr/seOkHzD5uZ/7n5VF9SnfBnT2VtQiZkbURgrX1Hgm9/L1vnExHJN707etn14V1UXlRJ/S31Yx6XaG3C0ARzLkpLtbJYRCQLor1RGq5pwFvmZfn9Y88LwMlrE4Y+5mrnUVAgEBHJOGstTWubCO0Isfz+5QRmB8Y9fvSEcy53HgVtOiciknGHv3mYlvtbWPDvC8acFxgt0YRzrigQiIhkUOcfO9nz8T3U/GMNdTenthjWjZ1HQakhEZGMGWwepOFtDQTqAiy7ZxnG4+6uosnSiEBEJAOciEPjtY1E2iKc/ezZFFW5u5FcKjQiEBHJgL0376XjqQ6WfHsJ5avLJzzerVXEiWhEICJyiloebOHg7QeZ/eHZnPbu0yY83s1VxIloRCAicgp6Xuxh5/t2UvE3FSz6v4tOei1Rr9/tVcSJaEQgIpKm8Ikw29+8HV+VjxU/XYHH/0rfeqxev9uriBPRiEBEJA1OxKHh7Q0MHB5gxU9XEDjtlUVj4/X63V5FnIgCgYhIGvbctIeOJ2KTw5XnVZ702ljPNR7q9bu5ijgRpYZERFJ0+K7DHP76YeZ+bC6zrp/1qtcTPdd4ZK/fzVXEiWhEICKSgrZft7H7ht3UXFHDwtsXJjwmmV7/6G2r3aQRgYhMGZ2hcFZ72b0NvTS8vYHSFaWccf8Z4+4omm+9/vEoEIjIlJDt2vzBlkG2XbENT4mHVb9Yha984tunW3sHpUqpIRGZ9LJdmx/tj7L9zdsZPDbIqo2rKK5zr8InGxQIRGTSG69K51RZa2l6fxNdz3ax7N5lVJxb8apj8mm7iHQoNSQik95EVTqnYv//vz/2bIEvLWDGNa9ONeXbdhHp0IhARCa9bNXmN/+4mX237GPme2YmfLZAPm4XkY6sjQiMMd8HrgBarLUr49+bBjwI1AP7gLdba9uz1QYRKRyZrtLpfLaTndfvpPL1lSz99lKMeXWFUD5uF5GObI4I7gYuG/W9zwC/s9YuBn4X/1pEJCMyVZvft6+P7VdtJzA3wIoNK/AEEt8q83G7iHRkLRBYa58B2kZ9+yrgh/HPfwi8OVvnFxFJR6QzwrYrtmHDllWPrsI/3T/msZXBIm66dCl9g1G6+8N5sV1EOnI9RzDTWns0/vkxYGaOzy8iMiYnHNtIrq+pjxUPr6B0Wem4xz/d1MIdm5rwegxRx3LTpUsm3UQxuDhZbK21gB3rdWPMWmPMZmPM5tbW1hy2TEQmk3RLN0f/nLWWpg800f6bdpZ8ewnVl1RP+PNDE8XlxUWU+L3csWnXpJsohtyXjzYbY2ZZa48aY2YBLWMdaK1dD6wHWLNmzZgBQ0QKV6qlm0NbUOxu7uaOTbtO+rl59/TQfE8z9bfVM+t9r95IbrSpMlEMuQ8EG4H3AF+Jf/x5js8vIlPEyB750NqB2x5tZMO86oQ34qGgEXUshzv6qC3zU1NWzEAkymM3N/LGR2DW2lnM/9z8pM6fzbULuZa11JAx5gHgWWCpMeaQMeb9xALApcaY3cDfx78WEUlZKquJRwaNIm/stne8Nza5u+pFw+U/g+J/qGTxNxcnLBNNJB+fK5CurI0IrLXvGOOlS7J1ThEpHKn0yEcGjahj8RiDYy1zdkR5471wYD6svHMe3YNRGIwmvRZhMu0wOh5tMSEik9JQj/y2RxvpHYhgDGP2yIt9XgYjDhAh6PcxvcxPya4w19xnaZ0Od18XxfPY9vj7GIL+2Cgjme0iJssOo+NRIBCRSWP08waS6ZEPzQ2Eow7HuvqpKQ1Q1+ll3c/8eGo93Hf9IEXTivB5PBxoCwGweEY5EccZd85hKlEgEJFJYawKofF65CPnBmrLiykN+Cg57nDjfX5MkaX6J0sIPdtA0OelPxzFYMBAOOpQXDR5q4BSpU3nRCTvpbu52+gJ5ZoBL9d/z0O0O8qZvz6TOaurgNg8Q5HXg8ViraXI65nUVUCpUiAQkbyX7vMGTppQDlmu/u8olR2waMNyyleXn1T50x+OUlPqZ3pZgP5wdFJXAaVKqSERyXvp1uwP3ei//NNGrlofpfYYeL9Vx5xLaoaPGT3PAEz6KqBUKRCISN5LpUJotL+dV8O//7qSnkOdLPjRMuZfe1rC9x/5XoUSAIYoEIjIpJBOzf7Qs4Z7/qeTM+49g5nXap/LRBQIRGTSSKVm3xlwaHhrA+2/bWfZD5Yx810KAmNRIBCRKccZjG0n3fbLNpZ8ZwmnvefV6SB5haqGRGRKccIOje9o5MTGEyz+78XM/sBst5uU9xQIRCTvpPuMASfisOOfd3B8w3EW3bmIOevmZKmFU4tSQyKSV1J9xsAQJ+Kw8907aX2wlYW3L2TujXOz2cwpRSMCEckb6a4gdiIOO67bQcsDLZz+ldOZd9O8HLV4alAgEJG8kc4KYifssOOdO2h9sJXT/+N06j5dl5O2TiUKBCKSN0auIB76ON4KYmfQofGfGml9KJYOqvukgkA6FAhEJG+k8tSvaF+U7W/ZzvFHYhPDSgelT5PFIpJXkllBHOmJsP2q7XQ82cGSby+h9Lpadjd3F9T+QJmkQCAieWe8FcSRzggvvulFup7tYtk9y2h6rYfb7vrD8OvJVhnJK5QaEpFJY/D4IFvfsJXu57pZ/uBySq6uSavKSE6mEYGITAqtL/XQ9MbtRA8OULy+nkOv9UFrD3BylVGhPFUsk1wJBMaYjwMfACywDXivtXb8J0yISMF66vGDdLxjD8V9ljuvGWDXrgY8u6E6WITX48HrMSk9p0BOlvPUkDFmDnAjsMZauxLwAtfmuh0iMjkc/UM7vdfsoSgMt79rgB1zoziOxRhoD4WJOg6DkYmrjMaT7pYWU4VbqSEfUGKMCQNB4IhL7RCRPNb+RDu7rtxG2A8P/C/Y73EwsSUGsQfNxz/edtVyZlWWpFU1lO6WFlNJzkcE1trDwO3AAeAo0Gmt/U2u2yEi+a31p628ePmLBOqKufsjcHx67KZv469HHYfBqMOJ3gG+8PNGjnT0pTUS0GSzO6mhauAqYAEwGyg1xlyX4Li1xpjNxpjNra2tuW6miGRQqqmXw986TMPbGyhfU85rfn82H3/XcmBoTsDgMRBxwGNgbnUQvy+9G3g6W1pMRW6khv4e2GutbQUwxmwA/gb40ciDrLXrgfUAa9assaPfREQmh1RSL9Za9n5+Lwe+dIBpb5rGigdX4C31ctG0VxaZFfu8vHCgnf94fCeVQT9eTyxFlE610MgtLQp5stmNdQQHgPOMMUFjjAEuAXa40A4RybJUUi9O2KHpfU0c+NIBZn1gFit/thJvqXf49cpgEYtnljOvJsjFS2dQ7PcScRxg4j2JxpLKlhZTWc5HBNbaPxtjHgZeACLAX4j3/EVkakmUeknUc490RWi4poH2Te3Mv2U+9bfUE+snJjZ0A7/t0UZ6ByIYQ9o38GS2tJjqXKkastbeAtzixrlFJHeSSb307etj2xXb6GvqY+n3lzLrvbOSeu9M3sDH29KiEGiLCRHJmolSL53PdvLC615g4NAAZz5+JsF/mp7SpPJQuqiQb+KZoC0mRCSrxuq5Nz/QzM737iQwN8Dqp1ez2fRq8ziXaEQgIlk3sudurWXfrfvY8c4dVJxbwTl/OodInV/1/C5SIBCRnIn2R9nxrh3s++I+Zr5nJmdtOgv/dL/q+V2m1JCI5MRg8yDb37Kdrme7WPDlBdR9um64MqjY52Uw4gARgn5fwdbzu0WBQESyrmd7D9uu2Ea4JcyKh1dQ+9ba4deGFpyFow7HuvqpKQ1QGvAWZD2/WxQIRCSrWh5qYed7d+Kr8LH6mdVUrKkYfm3kgrPa8mJKAz4GIw53X38u82qCLra6sCgQiEhWOBGHlz/zMofuOETF+RWseGgFgTmBk44ZPTcQ9PuwNkJ/JJrz9hYyBQKRPNcZCufdqteJ2jTYPEjDPzXQ+XQncz4yh4V3LMTjf3Vtivb6yQ8KBCJ5LB/3yp+oTZ3PdtJwTQOR9gjL7l3GadedNuZ7ZXKrCEmfsTb/N/Zcs2aN3bx5s9vNEMmpzlCYq+/6w0mPYYw6lg3rLnDtRjlemypKfBy56wgvfewlAvMCrNywkrKzypJ+33wb9UwFxpgt1to1Ex2nEYFInkp2w7Z8aNOx5hBHbjlC873N1FxRw7J7llFUnXwbC32vH7cpEIjkqXzMnydqU3UbdFy1i77tvdTfWs/8z83HeMbeOVTyj1YWi+SpfNwrf3SbFjRYbvimh/ChAVY9tor6L9QrCExCGhGI5LF83Cv/oqUz+OmcKvZ8cS/d3z9K2VlBVvx0BSWnl7jdNEmTAoFInjvV/HmmJ2L7D/Wz/71NdP+2nZnvmcmSu5bgLfFO/IOSt8YMBMaYXwIfttbuy11zRCSTMl1+2vJgC7s+tAsn7LDkO0uY9f5Z4z5JTCaH8eYIfgD8xhjzb8YY98ejIpKSVJ4XPJFwR5jG6xppvLaRkqUlrNm6htkfmK0gMEWMOSKw1j5kjPkV8HlgszHmXsAZ8frXctA+EUlTpspP259qZ+e7dzJwZID62+qpu7kOj091JlPJRHMEg0AvEADKGREIRCS/nWr5qTPgsPdzezl4x0FKFpdwzrPnUPHaiol/UCad8eYILgO+BmwEzrHWhnLWKhE5ZaeyfUPPth52XLeD3hd7mb1uNgu/uhBvqSaEp6rxRgT/BrzNWtuQ6ZMaY6qA7wIrAQu8z1r7bKbPIzKVpFP9k2r5qXUsh/7zEC/f/DK+ah+rHltFzRtrMtF8yWPjzRG8PovnvRN43Fp7jTHGD2jjcZFxnEr1T7Llp/0H+9n5np10PNlBzVU1LP3OUvy1/rTbLJNHzmd8jDGVwIXA9wCstYPW2o5ct0PEbZ2hMLubuyes4slk9c9Ymu9v5vlVz9P1XBdLv7uUlY+sVBAoIG4sKFsAtAI/MMacBWwBPmqt7XWhLSKuSKWHn83N58LtYXZ/eDctP26h4m8qOOOeMyhZqBXChcaNGjAfcA5wl7X2bGJVSZ8ZfZAxZq0xZrMxZnNra2uu2yiSNan28EdW/wx9TKX6Z6yRR/sT7Ww+czOtD7ey4N8XsPrp1QoCBcqNEcEh4JC19s/xrx8mQSCw1q4H1kPseQS5a55IdqXawz+V6p/RI4+bLl3Kosog4dubaf2vI5QsLeHsZ88+6TnCUnhyHgistceMMQeNMUuttU3AJUDjRD8nMlWkU9+fzuZzI0ceAZ+Xtt4Bvvr1rXzosWJmHjPwzmrWfGcl3qDKQgudW5vO3QDcF68Yehl4r0vtEBlXNp6clW4PP9XN504aeQw6nP+Yw2V/LKa/1PDjD8Lupe1swKESBYJC50ogsNZuBSZ8fJqIm7L5vOBcbC89NMKo3R3h8p9AbXMRf1oZ4flr/QxWeLAuP+1M8oe2oRZJYHRaZSAS5bZHG9kwrzqjI4Ns3oSDA/D556vhwXY6qix3/tMAzWf5qCnz5MXTziR/aOcokQQSTeha+8r385m1luYHmnlu2XPwUDu1H5nNGVvO5oOfPZMSvy9vnnYm+UMjApEE8vF5wcno29PHrnW7aN/UTvlryznz8TMpP7scgNNCpZxWWQJYFtaWKwjIMAUCkQROpWQzU1KZqI72Rzn41YMc+N8HMEWGRd9YxJx1czDe2PMCsjnfIZOfAoHIGNx8XnAqN+4TvzrB7ht207+nn9q31bLoPxcRmB0Yfj0X8x0yuWmOQGQclcEiFs/MbRol2ZXHfXv7+MuVL7LtjduwHli4cQX+b9TRX3Xyn/Vknu+Q3NCIQMQlY6V+Jlp5HO2Nsv/L+znw1YMMWsvvL4cnzu8lsu2vlO6K/UmPHEFM1vkOyR0FAhEXjJf6GevGXVsWoPmBZvZ8cg+DhwfZdjY89Y/QN81Lc0vsuVHTSgNEHOek1M/I+Y6uvjCOtdx8+TKlhWSYUkMiOTZR6mfoxh117HCp5+frF7D3sgZ2vHMH/pl+pm9cys/eBYPTfYSjDh5jMBjCUSdh6ueipTO46dKlONbi9Rju2LSLp5ta3LoEkmc0IhDJsWQ2nRuaqD7a1MXg7cdo+9e9hGqLWLJ+CbPeN4uugQjsamIgEqXI68GxsX0Zi7yJF4t1hsLcsamJEr9XE8byKhoRiGTBeA+dSWZb6Uh3hLYvH6LlggbaHz7OvE/N43W7X8fsD87GeM1Jo4b+cJTpZQFqSv30h6MJF4tpwljGoxGBSIZNVPo53hoFJ+xw9DtH2XfrPsItYWZcO4MF/3sBJQte/ZyA0eWtwJilrpowlvEYa/N/q/81a9bYzZs3u90MkWFjVfx0hsJcfdcfTqrZjzqWDesueNXNeeR7VBT7aH2olb2f20vfS31UXljJwq8upOLczD0nYChAWctw8NGisqnNGLPFWjvhBp8aEYikaLwefyoPnRmq6Gn/XTtbPr2Hni09lK4qZdVjq5h2+TSMMRltt5sL5CS/KRCIpGCiVbqppGC6nu9i77/tpX1TO4G6AMt+uIyZ75o5vC3E6PNm4gae7R1PZXJSIBBJwUQ9/mT2KOr+Szf7btnHiV+cwFfjY+HXFjJ73Wy8xYkfEKN9giTbFAhEUpBMj3+sFEzPiz3s++I+jj9yHF+1jwVfWsCcG+bgKx/7z1D7BEkuKBBIwUon3ZLsrqQjUzC9Db3su3UfrQ+14q30Un9rPXM/Ohdf5cR/fqk+6F4kHQoEUpBOJd2S7KRrqCnEvlv30fLjFrxlXuZ/fj5zPz6Xourkb+Aq+5Rc0IIyyZjxFlHlk2R39xzPeLuS9mzvofEdjTx3xnMc33icuk/Xcd7e81hw24KkgsDI65houwk9WUwyTSMCyYjJNKGZrXRL15+72P/l/Zz4+Qm8ZV7mfXIe826ah3+GP+n3GOs6quxTssm1QGCM8QKbgcPW2ivcaoecusk2oZnJdIu1lvbftnPgywfoeLIDX7WP+bfMZ+6NcwkVw/7ufmaETFLXYaLrmI/XUqYGN0cEHwV2AJlbOimumGwTmqMnfB1rWXvh6Sm9h41aWh9p5cBXDtCzpQf/bD8L71jIrLWz8JX50hohTbbrKFOHK3MExpi5wJuA77pxfsmsZDZRyzcXLZ3BhnUX8IHXLwAM3/v9Xq6+6w8Tbs0cDUU5/K3DPLfsORrf1ki0M8qS9Us47+XzmPeJefjKfGnPQUzG6yhTg1uTxf8JfApwXDq/ZNBkntD83u/34vdNfMMebB5k7xf38qf5f2L3ut34qn0sf2g55+48l9kfnI0n8MqfUro7fU7m6yiTW85TQ8aYK4AWa+0WY8zF4xy3FlgLUFdXl6PWSbom44RmMqmYnm09HPq/h2i+rxk7aKm5soZ5/zqPyr+tHHMvoBnlxUQdS1d/mFK/j4jjJN2zn4zXUSY/N+YILgCuNMa8ESgGKowxP7LWXjfyIGvtemA9xHYfzX0zJVWTbUJzzEdClgRo3dDK4W8cpuOpDjxBD7M+MIu5H51LcElwwvfderCdgYjD8Z4BAGpK/fyft56Z0qK1yXQdZfLLeSCw1t4M3AwQHxH86+ggIJILoyeNgyH4164Z7Fi1hYEDAwTmBzj9P05n1vtnUTQtuRvz0PxAebGP6qA/vvrYsHpedZb/a0TSp3UEBShTO1lOBRctncGP2v0c/MZB+ja0Y/tbKPm7KhbduYjp/zg94U6g4xmdbqooKVLlj+Q9VwOBtfYp4Ck321BoJtPCr2wGrGhflNaHWznyrSN0/bELT4mH0949kzk3zKFsZVna76stIWQy0oiggEymhV/ZCli9Db0cWX+E5nuaiXREKFlcwsKvLeS0609LaQ+gsSS7KZ1IPlEgyLJ8SsNMlgVLmQ5Y0VCU1odaObI+1vs3fkPt1bXM+l+zqLqoKukngSX7u1Tlj0w2CgRZlG9pmMmStshUwOp5sYej3znKsXuPEe2MUrKkhIW3L2Tmu2fir01+/x9I/Xepyh+ZTBQIsiQf0zD5nLYY2ds+lYAVbgvT8kALR39wlJ4tPbHe/zW1zF47m8oLx679n6ht+fa7FMkkBYIsydc0TD6mLRL1tlMJWE7Y4dAvWjl69zH6f92BHbSUrS5j0Z2LmPHOGfinp9b7Hy2V32U+pQJFkqVAkCX5nIbJp7TFmL3tdRewYd0FY95UrbX0bO2h+Z5mDt17FE5ECQVh27nwuptOZ82bM7caPdnfZb6lAkWSpUCQJfmchskn4/W2Ez34pX9/P80PNHP03mb6G0PgN+xYatl+FRxa6aUPh9/seZmzQ7Mydq2T+V0qfSSTmQJBFuVjGibfJNPbHmwdpPXhVlrub6Hz950AHKyHbVfD1pVRugIOtfHjA2QnBZfodzkyDZSvqUCRZCgQZFk+pWEmMvLGBrHeerHPS38kmrVANlZvOzgIxx46RsuPW2jb1AZRCC4PMuuLdXx24ADdtbGed3Qwwom2fkoDPoJ+X1ZTcCN/l6PTQDddugTIz1SgyEQUCAQ4+cYWGoxibWyfvxO9g0wvCxD0e7OW8x7qbR873Ev0t1103XiYPz61AztoKa4vpu6Tdcx4xwxKV5XyUksPHT86QDDe8w76fdSUBhiMOFibmxRcojTQHZt2cdOlS7ljU5NSgTLpKBDISTc2n8fDgbbQSa+39Q4S9JdkJecdPhHm+MbjHN9wnBO/boOwpbMSGs+D8z+ygNddU3dSyWeiVFJpwMvd15+b1ZHLSHtaexiMOJQXx84zlAZaPLNs3AlukXylQCAn5bf7w1E8xuBYCxZ8Xg9Ra/EYQ9SxGcl5Dxwb4PgjsZt/+5PtEAV/XYA/nW9pWg3HT/cy4Dj89sBezu6bfdL5xkolzas5eXvozlCYPa09gGVh7asnndP1dFMLt2xs4FhXP8e6+plTVYLf5xlOA02mVKDIEAWCSSDbtekje9lFXk8sCAAYiDjOcGA4lQe8hxpDHN94nBMbT9D15y6wULKkhLpP1VH71lqOzrFsuu8Fgv7Y/5IBz9iTrRNNwj/d1MJnNmx71fMATjWtNTRy8vs8zK4s4UhnH4c7+phTFeTWK5UGkslLgSDP5aI2fWQvOxKNMr0sgLUWa6EtNEh10D987mRvdk7YofP3nZzYeILjG4/T/3Js1FH+2nLqb6un9i21BJcHh9M+TvzxkMlOto7V8+4MhbllYyPHewbwxN/7RO8gt2xs4Of/cmpprZEjp4DPS2nAR1dfmK+9/UzOmT8t7fcVcZsCwSj5tDL0VGvTU/lvGd3LhtSrhsIdYdp/086RUSQYAAAQa0lEQVTxnx+n7ZdtRDoimICh+u+rqft0HTVX1BCYHUj4s5lad9HS3U80PooZCgTWQiR66mmt0fMTEcchUORhYW152u8pkg8UCEbIt5Whp1Kbns5/y+he9kTnsNbSu72Xtl+1ceKxE3T+oROiUFRbxPS3TKfmyhqmXToNb6l33PcZMl7KJ9mgNqO8GK9nRHoLsFh8XnPKpZxaJChTlQJB3Fi977unl+WsGmW0dLepyOYq13BHmI7fdXDiVydoe7yNwcODAJStLov1+t9YQ8V5FSk/2WtIopRPKkGtMljErVcuf9Ucwa1XrsjI70+LBGUqUiCIS9T7bunq553f/RN+nwfI/Qgh3R5oJle5OhGH7ue7ad/UTtuv22ITvVHwVnipvrSaaZdNo+byGgJzEqd8EkklZZVOULto6Qwe/+iFWakagsm1SFAkGQoEcaN736HBCG2hQeqmBYdXrLqxd8zIHuhQvr4zFJ4wPQLprXK11tK3q4/237XT/tt22p9oJ9oZBQPlrymn7jN1TPv/plFxXgWeIk/K/z2ppqzSDWqVwSLOma8HxoskQ4EgbnTvOxx14itq4+WMLu4dUxksYuvB9nFvoKN72amMJPr399P+RDsdT3bQ8WQHA4diKZXA/AAz3jaD6kurqb6kmqKaV/bXMeEolSkGgnR69/m8i6vIVKFAMMLo3vf1dz+X1A0oUaojk9VHE91Ax+plj5XLHjg8QPuTr9z4+/fGet1F04uouriK6r+vpuqSKkoWlpy0qvdUJ9PT6d1rglYk+3IeCIwx84B7gJmABdZba+/MdTvGMjL/m8wNKNHNEcho9dF4N9CR5xpa+DUySFSU+Ojf28/RZ47T+UwnHc900L8n9nO+ah9VF1Ux92Nzqfq7KkpXlGI8iSd5MzEBnW7vXhO0ItnlxoggAtxkrX3BGFMObDHGbLLWNk70g7k20Q0o0c3xlo0NgMHvy1zFzng30JbufnoHIrSFwmDBWDijq4iX7zxA4MUBOv6nY7iyxzfNR9WFVcz5lzlU/V0VZWeWjXnjH22sYLSntZvy4qKkbtCn0rvXBK1I9uQ8EFhrjwJH4593G2N2AHOAvAsEMP4NKNHNMTQQxRhLwOcf/t5Q772rL0xTcxdLZ1a8am+cidow1g20o6WfGdsjnH/Iy6LDXhYc9hAcMHRyEP8cP5V/W0nV66uovKiS0uVj9/gnkigYhQYjfOInL+KNTxXcdOlSFs8sGzcoqHcvkn+MHbHwJucnN6YeeAZYaa3tGvXaWmAtQF1d3Wv279+f8/ZNpDMU5uq7/nDSiGAw4jB6RBB1LJetOI31//Py8M/e+IZF3HDJkpTO19EzyJG/dlK8bYDwC710PdtFb0MvWHCwHK21vDTP4ehCDx/91Jkse01NWg9rH8tQGsxacKxlIOJQXuwj4PPS1jtAS/cAc6pK8HqM64vxRASMMVustWsmPM6tQGCMKQOeBr5krd0w3rFr1qyxmzdvzk3DUjTy5jjUUwdO+t66i07n5ke2A+DzeIg4DgBPfOLicUcGg82DdD3fRfdz3XT9OfYx0hEBwFvppfL8SvxrSrm95SBH5kE4aIZX1G5Yd0FGe9tDk99DJazd/RE++fBfCfp9RB3LntYeHGuprynFGIg6NuNtEJHUJBsIXKkaMsYUAT8F7psoCGRDJit6xkp1jPze8/tOALEgMPQxHHVoau4aDgThE2G6X+ime0s3PVt66Hq+i4H9sTJOPFC6spTat9VS8boKKs6vILgsOJzm+eem8ljgcWxWqmoSTYivnher0R+IRIdHCMYYirwevB6jxzSKTCJuVA0Z4HvADmvt13J13qGb/+7mHu7Y1DT8/UykMBLNI4z83tKZFUBsS2efx0Ogx2HJMQ8zHuhh+47t9GzpoX9f//DPFp9eTMV5FVTcWEH5a8spP6d83P16spl3H7NaaN0Fw/MWUSc2CpleWoTXYzJa659PmwCKTFVujAguAP4Z2GaM2Rr/3mettb/M1gmHerRRx3K4o48Z5QGmlQayvlrYWsvA4QFK/hLi34/M5PAf26lr9jC9MzYyaOMQxQuLKT+3nNnrZlP+mnLKzimjqDr1tmSrqma80tWRAWh3czd3bNqV0Vr/fNsEUGSqcqNq6PdA5mYwJzCyRzu0LfHxnkEqS/wZXS0c7Y8S2hGi56899L7YS89fe+jZ2kOkLZbTn2lg7oJSBv7Gz7Rzq5h34bTYTb8qv3u5E9X+DwWgxTPLuWBRbc4W0YlI5kz5lcUje7RRJ5bHdqwlHHWIOKScwrBRS9+ePnq399K7rTf2cXsvod0hiMaO8ZR4KF1VSu1bayk7q4yys8soPbMUX1n2Lne2Uiip1P5nclSSyY3zRGR8Uz4QjO7RTi8torVnkHDUGS5zTHRjsY6lf38/ocYQvQ0jbvg7Qjj9saofDJQsLIlN5F5TS+mZpZSdWUbJopK0t2GG1G/q2U6huFH7rz2GRHLH1XUEyTrV8tHRJZ43XbqExTPLmVFeTLnPS99LfYR2hAjtDNG7I3azD+0I4fQ5w+/hn+OndGXp8L+yVWUEzwjiDSb30BVI7gaf6k090VqGqVK6mag0V3MEIsnL6/LRXLtwSS0PXv1ajv21k5LDUew9IfqaTrCzKUTfy33DKR2AwLwAweVBZl80m+DyIKXLSwmeEUxrAnekZG7w6eTFp3IKRauQRXKjIALBkW8fYfe63cNfm4AhuDhI2Vll1L69luCy4PC/bOTxk73Bp3NTn+opFO0xJJJ9BREIqi6uYtHXFxFcGqRkSQnF84pPKYefqmRv8Onc1CeazFUdvohMpCACQemyUkqXlbp2/mRv8OnuzjlWCkV1+CKSjIKYLM4HqUx8ZqIXP5UnkUUkOZosdsF4N/BUJj4zkRefypPIIpJZCgQZkkwaJpcTn1N9EllEMie1p49LQiOrgoJ+H16P4bZHG+kMhV1r09B8Q9Sx9A5EiDpWz/oVkYQ0IsiAfE3DqA5fRJJR8IEgExOz+ZyGUR2+iEykoANBpsorT+Wh7CIibivYQJDpbY6VhhGRyapgA0FLdz9Rx+IxhqhjM5LXnyppGK1GFiksBRsIdjf3cLijDwBjDNNLiyjx+/Iir+8mrUYWKTwFWT7aGQpzx6YmZpQH8BiDtZbWnkFuunRJQfeA87EMVkSyryBGBKNTHUPlntNKA1SW+AlHHcJRh8Uzy11uqbvytQxWRLJrygeCRKmO1fOqgVfKPSMOeD2m4NNC+VwGKyLZM6VTQ2OlOgCtuk1Aq5FFCpMrIwJjzGXAnYAX+K619ivZOM94qQ6Veyam6yJSeHIeCIwxXuCbwKXAIeB5Y8xGa23j+D+ZuolSHVOl3DPTdF1ECosbqaFzgZestS9baweBHwNXZeNEyaQ6OkNhdjd3c/BEiN3N3aqQEZGC40ZqaA5wcMTXh4DXjT7IGLMWWAtQV1eX9snGS3UMTST3DkQ40TvI9LIAQb9XtfMiUlDydrLYWrveWrvGWrumtrb2lN6rMljE4pnlrxoJDE0ct8VHAW29gwCqnReRguJGIDgMzBvx9dz493JqaCLZYwxY8Hk82PjX1r7yuojIVOdGIHgeWGyMWWCM8QPXAhtz3YihCWPHWjAQcRxM/GvVzotIIcl5ILDWRoCPAL8GdgA/sdY25LodQxPJANUlsZRRddAPaAtpESksrqwjsNb+EvilG+ceaeREcrHPS38kqtp5ESk4U36LiYmoZl5ECl3eVg2JiEhuKBCIiBQ4BQIRkQKnQCAiUuAUCERECpwCgYhIgVMgEBEpcMZa63YbJmSMaQX2n+LbTAeOZ6A5maZ2pUbtSk2+tgvyt21TqV3zrbUT7to5KQJBJhhjNltr17jdjtHUrtSoXanJ13ZB/ratENul1JCISIFTIBARKXCFFAjWu92AMahdqVG7UpOv7YL8bVvBtatg5ghERCSxQhoRiIhIAlMqEBhjvm+MaTHGbB/jdWOM+box5iVjzIvGmHPypF0XG2M6jTFb4/++kKN2zTPGPGmMaTTGNBhjPprgmJxfsyTblfNrZowpNsY8Z4z5a7xdtyY4JmCMeTB+vf5sjKnPk3Zdb4xpHXG9PpDtdo04t9cY8xdjzKMJXsv59UqyXW5er33GmG3x825O8Hrm/yattVPmH3AhcA6wfYzX3wj8CjDAecCf86RdFwOPunC9ZgHnxD8vB3YBy92+Zkm2K+fXLH4NyuKfFwF/Bs4bdcyHgW/FP78WeDBP2nU98F+5/n8sfu5PAPcn+n25cb2SbJeb12sfMH2c1zP+NzmlRgTW2meAtnEOuQq4x8b8CagyxszKg3a5wlp71Fr7QvzzbmKPDp0z6rCcX7Mk25Vz8WvQE/+yKP5v9CTbVcAP458/DFxijDF50C5XGGPmAm8CvjvGITm/Xkm2K59l/G9ySgWCJMwBDo74+hB5cIOJOz8+tP+VMWZFrk8eH5KfTaw3OZKr12ycdoEL1yyeTtgKtACbrLVjXi8bez53J1CTB+0CeGs8lfCwMWZettsU95/ApwBnjNdduV5JtAvcuV4QC+K/McZsMcasTfB6xv8mCy0Q5KsXiC0FPwv4BvCzXJ7cGFMG/BT4mLW2K5fnHs8E7XLlmllro9ba1cBc4FxjzMpcnHciSbTrF0C9tfZMYBOv9MKzxhhzBdBird2S7XOlIsl25fx6jfC31tpzgMuBfzHGXJjtExZaIDgMjIzsc+Pfc5W1tmtoaG+t/SVQZIyZnotzG2OKiN1s77PWbkhwiCvXbKJ2uXnN4ufsAJ4ELhv10vD1Msb4gErghNvtstaesNYOxL/8LvCaHDTnAuBKY8w+4MfAG4wxPxp1jBvXa8J2uXS9hs59OP6xBXgEOHfUIRn/myy0QLAReHd81v08oNNae9TtRhljThvKixpjziX2e8n6zSN+zu8BO6y1XxvjsJxfs2Ta5cY1M8bUGmOq4p+XAJcCO0cdthF4T/zza4AnbHyGz812jcohX0ls3iWrrLU3W2vnWmvriU0EP2GtvW7UYTm/Xsm0y43rFT9vqTGmfOhz4B+A0dWGGf+b9J3KD+cbY8wDxKpJphtjDgG3EJs4w1r7LeCXxGbcXwJCwHvzpF3XAOuMMRGgD7g2238McRcA/wxsi+eXAT4L1I1omxvXLJl2uXHNZgE/NMZ4iQWen1hrHzXG3AZsttZuJBbA7jXGvESsQODaLLcp2XbdaIy5EojE23V9DtqVUB5cr2Ta5db1mgk8Eu/j+ID7rbWPG2M+BNn7m9TKYhGRAldoqSERERlFgUBEpMApEIiIFDgFAhGRAqdAICJS4BQIRFJkYruj7jXGTIt/XR3/ut7dlomkR4FAJEXW2oPAXcBX4t/6CrDeWrvPtUaJnAKtIxBJQ3wLjC3A94EPAquttWF3WyWSnim1slgkV6y1YWPMJ4HHgX9QEJDJTKkhkfRdDhwF8mIHUpF0KRCIpMEYs5rY5m7nAR/PxQOORLJFgUAkRfFdT+8i9pyEA8BXgdvdbZVI+hQIRFL3QeCAtXZT/Ov/Bs4wxlzkYptE0qaqIRGRAqcRgYhIgVMgEBEpcAoEIiIFToFARKTAKRCIiBQ4BQIRkQKnQCAiUuAUCERECtz/Aw/oFgRODtheAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "exog = sm.add_constant(x)\n",
    "endog = y\n",
    "\n",
    "# Normal distribution + log link function\n",
    "mod = sm.GLM(endog, exog, family=sm.families.Gaussian(sm.families.links.log()))\n",
    "res = mod.fit()\n",
    "display(res.summary())\n",
    "\n",
    "y_pred = res.predict(exog)\n",
    "\n",
    "idx = x.argsort()\n",
    "x_ord, y_pred_ord = x[idx], y_pred[idx]\n",
    "\n",
    "plt.plot(x_ord, y_pred_ord, color='m')\n",
    "plt.scatter(x, y,  s=20, alpha=0.8)\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "# plt.savefig('custom_GLM_fit.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}