Chapter 8 - General Linear Model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "8b336d96-b179-4978-ad97-aad1cee6eac3",
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'arvis'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "Input \u001b[1;32mIn [36]\u001b[0m, in \u001b[0;36m<cell line: 10>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      8\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mstatsmodels\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mapi\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01msm\u001b[39;00m\n\u001b[0;32m      9\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpymc3\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpm\u001b[39;00m\n\u001b[1;32m---> 10\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01marvis\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01maz\u001b[39;00m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'arvis'"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.special\n",
    "#plt.style.use('fivethirtyeight')\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import pymc3 as pm\n",
    "import arvis as az"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d844190b-3d5b-41b1-8f1f-430edb3a5f88",
   "metadata": {
    "tags": []
   },
   "source": [
    "# <font face=\"gotham\" color=\"orange\"> Generalized linear models </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22630629-eff9-4518-9951-9a3142aef855",
   "metadata": {},
   "source": [
    "The **generalized linear model** (GLM) is simply an umbrella term for all linear models which are not simple/multiple linear regression models, in econometrics we have seen logit, probit, tobit and etc., they are all generalized linear models.\n",
    "\n",
    "As a reminder, the linear regression model's predicted values take the form\n",
    "$$\n",
    "\\hat{Y}=X\\beta\n",
    "$$\n",
    "To formulate a generalized linear model, we impose an arbitrary function onto the linear combination, i.e.\n",
    "$$\n",
    "\\tilde{Y} = f(X\\beta)\n",
    "$$\n",
    "where $f$ is called the **inverse link function**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b0f1b75-c672-4499-822a-73014d634973",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"orange\"> Logistic Model with a Single Variable</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36b8a6d0-1454-4fda-a71c-d3eb4ae91e37",
   "metadata": {},
   "source": [
    "The most famous type of GLM is **logistic regression** which is also called _logit regression_ in the context of econometrics. Here we provide a brief refresh of logit regression.\n",
    "\n",
    "A simple logit model hypothesizes a sigmoid shape function \n",
    "$$\n",
    "P_i = \\frac{1}{1+e^{-(\\beta_1+\\beta_2X_i)}}\n",
    "$$\n",
    "where the inverse link function is $\\frac{1}{1+e^{-z}}$, and $z$ represents the predicted value of a linear model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "850faec3-a87a-412a-83c7-0573f830a15a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAE9CAYAAAD5xEDLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsnElEQVR4nO3deZwU9Z3G8c+3ey6G4T6G+5RDRFRAEE+IF2qiMcZE470xxCSabDabRHPubrLJJrq5dk2UGKNGI5tEk2BCxAu8ECKgooAoDNdwMyMw90x3f/ePbnQkoHP0dHXXPO/Xq1/dXVXT/fwc2mfq6Cpzd0RERCQcIkEHEBERkfRRsYuIiISIil1ERCREVOwiIiIhomIXEREJERW7iIhIiOQFHSAd+vbt6yNGjEjb69XU1NC1a9e0vV6QNJbsE5ZxgMaSrcIylrCMA9I/lhUrVux1936HmxeKYh8xYgTLly9P2+stXryYmTNnpu31gqSxZJ+wjAM0lmwVlrGEZRyQ/rGY2eYjzdOmeBERkRBRsYuIiISIil1ERCREVOwiIiIhomIXEREJERW7iIhIiGS02M3sbjPbbWavHWG+mdnPzGy9ma0ys8mZzCciIpLrMr3Gfg8w+z3mnweMSd3mAL/IQCYREZHQyGixu/szQOV7LHIRcJ8nLQV6mtnAzKQTERHJfdl25rnBwNZmz8tT03YEE0dERLKNuxNPOLGE0xRPEE/4Ozd3YvF3Hjefl3An4bzzOLWMO2/PS7jj7iQS70yD5P3B5TyVwR2cgz+fmpZc/O3pyWVhy44YMzP038fcPUNvlXpDsxHAX9x94mHm/RX4vrs/l3r+JPAVd19xmGXnkNxcT2lp6ZR58+alLWN1dTUlJSVpe70gaSzZJyzjAI0lW2V6LLGEUxeD+phTH4e6mFMfS05rjDsNcWhMOE1xaIhDUyI1Le40JaApkXyN2NuPk8s0xRLEiSRL3CHuEE8k73NN3yLntpnp+53MmjVrhbtPPdy8bFtjLweGNns+BNh+uAXdfS4wF2Dq1KmeznPw6vzE2SksYwnLOEBjyVZtHYu7U9UQY/eBBnZX1bOnqoH9dU3sq23irdpG9tc2sa/u3Y+r62M0xhMtfo+CaISi/AhF+VG6FEQpzItQUBChKBqhIC9CQV6UgmiEwrwIlRW7GTZoEHlRIz8aIT9q5EUj5EeS93lRIz8SIRox8qJGxIy8iBE95JYXMcyMqBmRCEQsOT1iB28QTS0TMd6ebkbyxjs/Z6Smvf04eX9wef5hfvLnly19IWP/vrKt2OcDN5rZPGA6sN/dtRleRCQN6hrjbKmsZXNFDVsqa9m2r47dVQ3sOdDArqp6dh9ooK4pftifLSnMo2dxfvLWpYDBPbvQo0s+3YryKSmMUlKYR9fCPLoVJe9LUrcuBVG65EcpSt2iEWtx3uQfKJPSNfxA9Shs+bjbK6PFbmYPAjOBvmZWDnwbyAdw9zuABcD5wHqgFrguk/lERHJdPOFsr04w/5XtbNpbw+aKWrZUJu93VzW8a9nigiil3Yvo162QSUN60r9bIaXdC+nfrYj+3QvpV1JIz+ICehbnkx/VaU9yRUaL3d0vf5/5DnwuQ3FERHJaTUOM13ceYM32A6zZUcWaHQdYt/MA9U0J4CUABnQvYlifYs4Y24/hfYoZ1qcrw3oXM7x3MT2L8zHL3JqkZEa2bYoXEZEj2HWgnhc2VLBkw15e3PQWmypqOHj8c48u+RwzqDtXTB9OZP82Lp41nVH9ulKUHw02tGScil1EJEvtrW5gaVkFL2xI3sr21gDJEj9xRG8uPmEwEwZ2Z8Kg7gzsUfT22vfixbuZMKh7kNElQCp2EZEs4e6s3VHFI6u289Ta3azbVQUkD1ybNrI3l08bxozRfTh6YPdWHYQmnYuKXUQkYBv31jD/5e08smo763dXE40YJ43qzUUnjGPGqD4cO7gHeTp4TVpIxS4iEoDt++r4y6rtPPLKDl7dth+AaSN7890PT+S8iQPoU1IYcELJVSp2EZEMcXeeeXMvc5/ZwPPrKwCYNKQHXz//aD543EAG9ugScEIJAxW7iEgHi8UTLHhtJ3cs3sCaHQcY0L2IL541lguPH8TIvl2Djicho2IXEekg9U1xfr+inF8+U8aWylpG9+vKDz86iQ8fP5iCPO0zl46hYhcRSbP9dU3cv3Qzv35+I3urGzl+aE++fsHRnH10KREdzS4dTMUuIpIm8YRz3wub+NFjb1DVEOOMsf34zMzRTB/ZW2d4k4xRsYuIpMFr2/bztT++yqry/Zw2pi83nzeeYwb1CDqWdEIqdhGRdqhuiPGjx97gniUb6d21kJ9dfgIfmjRQa+gSGBW7iEgbLVy9k3+bv5qdB+q5YvowvnzueHp0yQ86lnRyKnYRkVbatq+Ob/95NU+s3cX4Ad24/YrJTB7WK+hYIoCKXUSkVf6wopxv/fk13OFr54/nulNG6lrlklVU7CIiLRCLJ/j+317nV89tZMaoPtx66SSG9CoOOpbIP1Cxi4i8j/21Tdz44EqefXMv1548gm9ccLQuyiJZS8UuIvIe1u+u4lP3raD8rVp+cMmxfPzEYUFHEnlPKnYRkSNY9PpuPv/gSxTmR3jwUycxdUTvoCOJvC8Vu4jIIdydO58p4wePvs6Egd2Ze/VUBvfUldckN6jYRUSaqW+Kc/NDq/jTy9u5YNJAbvvocXQpiAYdS6TFVOwiIim1jTGuvftF/r6pki+fO47PzhytM8hJzlGxi4gATQnn079ZwfLNlfzs8hO48LhBQUcSaRMVu4h0erF4gjteaWDFrlp+eMkklbrkNH0RU0Q6tUTC+cofVrFiV5xvf2gCHztxaNCRRNpFxS4inZa78635r/HwS9v4yJh8rjtlZNCRRNpNm+JFpFNyd/7r0de5f+kWPn3GKE4q2hl0JJG00Bq7iHRKP1+8gTufLuOK6cO4efZ4Hf0uoaFiF5FO557nN3LrwnVcfMJgvnPRRJW6hIqKXUQ6ld8v38q/PbKGcyaUcutHJxGJqNQlXFTsItJpLN9Uyc0Pv8ppY/ryP584QVdok1DSv2oR6RQqaxq56cGXGNKrC7dfMZnCPJ0mVsJJR8WLSOglEs6XfvcyFdWNPPzZk+lelB90JJEOozV2EQm9uc+WsWjdHr75waOZOLhH0HFEOpSKXURCbfmmSm5duI4Ljh3IlScNDzqOSIdTsYtIaDXfr/79S47V19qkU9A+dhEJJe1Xl85Ka+wiEkrary6dlYpdREJH+9WlM1Oxi0ioaL+6dHbaxy4ioaH96iIBrLGb2WwzW2dm683s5sPM72Fmj5jZK2a22syuy3RGEclND764hUXr9vAN7VeXTiyjxW5mUeB24DxgAnC5mU04ZLHPAWvc/ThgJvDfZlaQyZwiknv2Vjfwg7+9zoxRfbhK+9WlE8v0Gvs0YL27l7l7IzAPuOiQZRzoZskdYyVAJRDLbEwRyTXfW7CWuqY43/nwMdqvLp1apot9MLC12fPy1LTm/hc4GtgOvAp8wd0TmYknIrloaVkFD6/cxqdOG8VR/bsFHUckUObumXszs0uBc939+tTzq4Bp7n5Ts2U+CpwC/AswGngcOM7dDxzyWnOAOQClpaVT5s2bl7ac1dXVlJSUpO31gqSxZJ+wjAOyYyyxhPOtJXU0xuE/T+1CYbRta+vZMJZ0CctYwjIOSP9YZs2atcLdpx52prtn7AbMABY2e34LcMshy/wVOK3Z86dIlv8RX3fKlCmeTosWLUrr6wVJY8k+YRmHe3aM5ReL1/vwr/7Fn1izs12vkw1jSZewjCUs43BP/1iA5X6ETsz0pvgXgTFmNjJ1QNxlwPxDltkCnAlgZqXAOKAsoylFJCeUv1XLT594k7MnlHLm0aVBxxHJChn9Hru7x8zsRmAhEAXudvfVZnZDav4dwHeAe8zsVcCAr7r73kzmFJHc8B+PrAHg2x869Ms1Ip1Xxk9Q4+4LgAWHTLuj2ePtwDmZziUiueXJtbt4bM0uvjp7PEN6FQcdRyRr6JSyIpJz6hrjfHv+asb0L+GTp44MOo5IVtEpZUUk5/zvojcpf6uO/5tzEgV5Wj8RaU6fCBHJKet3VzP3mTI+Mnkw00f1CTqOSNZRsYtIznB3vvmn1+iSH+Vr5x8ddByRrKRiF5GcMf+V7bxQVsFXZo+nb0lh0HFEspKKXURyQmMswW2PrWPi4O5cPm1Y0HFEspaKXURywu+Wb2VrZR3/es44ohFd5EXkSFTsIpL16pvi/M9TbzJ1eC/OGNsv6DgiWU3FLiJZ7/6lm9l1oIEvnTNOl2QVeR8qdhHJajUNMX6xeAOnHtWXGaP19TaR96NiF5Gsds+STVTUNPKlc8YGHUUkJ6jYRSRr7a9r4s6nN3DW0f05YVivoOOI5AQVu4hkrbueLeNAfYwvnq21dZGWUrGLSFaqqG7g7uc2csGkgRwzqEfQcURyhopdRLLSnc+UUdcU54tnaW1dpDVU7CKSdXYdqOfeJZu4+IQhHNW/JOg4IjlFxS4iWef2ReuJJ5wvnDkm6CgiOUfFLiJZZWtlLQ/+fQsfO3Eow/oUBx1HJOeo2EUkq/zPU29iZtz0gaOCjiKSk1TsIpI1yvZU89DKbVw5fTgDe3QJOo5ITlKxi0jW+MkTb1IQjfDZWaODjiKSs1TsIpIVtlTU8pdV27n65OH0LSkMOo5IzlKxi0hWuOu5MqIR459OGRl0FJGcpmIXkcC9VdPI75Zv5cPHD6a0e1HQcURymopdRAL3m6WbqW9KMOf0UUFHEcl5KnYRCVR9U5x7l2ziA+P7M6a0W9BxRHKeil1EAvXQynIqahq1ti6SJip2EQlMPOHc9exGJg3pwfSRvYOOIxIKKnYRCcwTa3excW8Nc04fhZkFHUckFFTsIhKYuc+UMbR3F2YfMyDoKCKhoWIXkUCs2FzJis1vcf2po8iL6n9FIumiT5OIBOLOp8voWZzPpVOHBB1FJFRU7CKScWV7qnl87S6uOmk4xQV5QccRCRUVu4hk3C+f3Uh+NMLVM0YEHUUkdNpU7GbW1cyi6Q4jIuG3t7qBh1aWc8nkIfTrpou9iKRbi4rdzCJm9gkz+6uZ7QZeB3aY2Wozu9XMxnRsTBEJi/uWbKIpnuBTp+liLyIdoaVr7IuA0cAtwAB3H+ru/YHTgKXAf5nZlR2UUURCorYxxn1LN3P20aWM6lcSdByRUGrpUStnuXvToRPdvRJ4CHjIzPLTmkxEQuf3y8vZV9vEp8/Q6WNFOsr7rrGb2UKgoNnzw15T8XDFLyJyUDzh3PVcGZOH9WTKcJ0+VqSjtGRTfByYa2bFqec/6sA8IhJST6zdxdbKOl3sRaSDtaTYNwHfJFnuXYB2ndDZzGab2TozW29mNx9hmZlm9nLq4Lyn2/N+IpId7nthE4N7duGso0uDjiISai3Zx27uXmZm/wbcBbR5k3vqK3K3A2cD5cCLZjbf3dc0W6Yn8HNgtrtvMbP+bX0/EckO63dX8fz6Cr4ye5xOHyvSwVryCfsGgLuvB75DspTbahqw3t3L3L0RmAdcdMgynwAedvctqffd3Y73E5EscN8LmynIi/DxqUODjiISeu9b7O5eAWBmfd39dWB8O95vMLC12fPy1LTmxgK9zGyxma0ws6vb8X4iErCq+iYeWlHOhyYNok+JTkgj0tHM3Vu2YHKT+YXtejOzS4Fz3f361POrgGnuflOzZf4XmAqcCXQBXgAucPc3DnmtOcAcgNLS0inz5s1rT7R3qa6upqQkHN+x1ViyT1jGAS0byxObm7h/bSPfnlHEyB7Ze8LKzvZ7yQVhGQekfyyzZs1a4e5TDzevNVdfaNdBcynlQPNtcUOA7YdZZq+71wA1ZvYMcBzwrmJ397nAXICpU6f6zJkz0xAvafHixaTz9YKksWSfsIwD3n8siYTzHz9+muOHFnPdRadkLlgbdKbfS64Iyzggs2NpzVEsLVu1f28vAmPMbKSZFQCXAfMPWebPwGlmlpf6it10YG0a3ltEMuz5DXsp21PDNScPDzqKSKeR0TV2d4+Z2Y3AQiAK3O3uq83shtT8O9x9rZk9CqwCEsBd7v5ae99bRDLv3iWb6dO1gPOPHRh0FJFOozXFfks63tDdFwALDpl2xyHPbwVuTcf7iUgwtlbW8uTru/jczKMozMvefesiYdPiTfFaaxaR1rh/2WYiZnxi+rCgo4h0Kq06U4SZjTGzu83s9o4KJCK5r74pzv+9uJVzJpQyqGeXoOOIdCqtPQXUb4Dfk7xcK2Y20czuS3sqEclp81/Zzr7aJq6eMSLoKCKdTmuLPeLufyN5YZiDm+cnpj2ViOQsd+feJZsYV9qNk0bpKm4imdbaYt9uZiNJffXNzIzkSWRERABYuWUfq7cf4OqTh5P8X4SIZFJrjooH+GeSF4IZYGbXAbMBHVQnIm+774VNdCvK48PHH3q2aBHJhFatsbv7JpJl/nlgFPA0cFX6Y4lILtpdVc+CV3dw6ZShdC1s7XqDiKRDqz957h4D/pC6iYi87cFlW2mKO1fN0JnmRIKiCyOLSFo0xRM8sGwzZ4ztx8i+XYOOI9JptanYzaw912QXkRBauHonu6sadF54kYC1dY39B2lNISI57zcvbGZIry6cMbZ/0FFEOjVtiheRdntzVxXLNlZyxfThRCP6iptIkFp88JyZ/Zrk99cNGGZmdx+c5+7/1AHZRCRHPLBsCwXRCB+bOiToKCKdXmuOir+n2eNTgXvTG0VEclFtY4yHVpRz3rED6FNSGHQckU6vxcXu7k8ffGxmVc2fi0jn9eeXt1PVEOOqk3TQnEg2aOs+9sa0phCRnOTu3L90M+MHdGPK8F5BxxER2ljs7n5SuoOISO4p259g9fYDXHGSzgsvki10VLyItNlTW2J0LYhy8Qk6L7xItlCxi0ibvFXTyLKdMS6ePJgSnRdeJGuo2EWkTf6wopxYAq7UQXMiWaXdxZ66fKuIdCKJhPPAss2M7RVh/IDuQccRkWbSscb+72l4DRHJIc+t38umilpmDc0POoqIHKJFO8bMbNWRZgGl6YsjIrng/qWb6dO1gKkDokFHEZFDtPSIl1LgXOCtQ6YbsCStiUQkq+3YX8cTa3cx5/TR5Ed2Bh1HRA7R0mL/C1Di7i8fOsPMFqczkIhktweXbcGBK6YPY8MqFbtItmlpsV/v7n64Ge7+CQAzsyMtIyLh0BRPMO/Frcwc24+hvYvZEHQgEfkHLT14bpGZ3WRmw5pPNLMCM/uAmd0LXJP+eCKSTR5fs4vdVQ36iptIFmvpGvts4J+AB81sJLAPKAKiwGPAjw+3mV5EwuX+pZsZ3LMLM8f1DzqKiBxBi4rd3euBnwM/N7N8oC9Q5+77OjCbiGSR9burWbKhgi+fO45oROeFF8lWLdoUb2YzLHWFB3dvcvcdKnWRzuWBZZvJjxofP3Fo0FFE5D20dB/7NcAKM5tnZtea2YCODCUi2aWmIcYfVpRz3sSB9C0pDDqOiLyHlm6KvwHAzMYD5wH3mFkPYBHwKPC8u8c7LKWIBOpPL2+jqj7GNSePCDqKiLyPVp1S1t1fd/cfu/ts4APAc8ClwLKOCCciwXN37luymYmDuzN5WM+g44jI+2jzueLdvc7dFwAr3X1qGjOJSBZZtrGSdbuquHrGCFKH2ohIFtNFYETkPd27ZBM9i/O58LhBQUcRkRbQRWBE5Ii276vjsTW7uP60kRTl64IvIrlAF4ERkSP67bItJNy5crrONCeSK3QRGBE5rIZYnAf/voUzx5cytHdx0HFEpIVa+nW3T77HvE+kL46IZIsFr+6goqaRa07W2rpILknHwXMiEkL3LtnMqH5dOWV036CjiEgrqNhF5B+8snUfL2/dx9UnDSei88KL5JSMF7uZzTazdWa23sxufo/lTjSzuJl9NJP5RATue2EzXQuiXDJlSNBRRKSVWl3sZvaB5vet/NkocDvJ09JOAC43swlHWO4HwMLWvoeItE9FdQOPrNrORyYPoVtRftBxRKSV2rLGftsh960xDVjv7mXu3gjMAy46zHI3AQ8Bu9vwHiLSDv+3fCuNsQRXz9BBcyK5qD2b4tuy420wsLXZ8/LUtHde1GwwcDFwR9ujiUhbxOIJHli6hZNH92FMabeg44hIG7T0e+zpcrg/BvyQ5z8Bvuru8fc6L7WZzQHmAJSWlrJ48eI0RYTq6uq0vl6QNJbsk83jWLErxrZ9DXxkZKJFGbN5LK2lsWSfsIwDMjwWd2/VjeRFXwBeasPPzgAWNnt+C3DLIctsBDalbtUkN8d/+L1ed8qUKZ5OixYtSuvrBUljyT7ZPI7L577gJ3//SW+KxVu0fDaPpbU0luwTlnG4p38swHI/Qidmeo39RWCMmY0EtgGXAe86wY27jzz42MzuAf7i7n/KYEaRTunNXVUs2VDBV2aPIy+qb8KK5Kq2FHt16r6qtT/o7jEzu5Hk0e5R4G53X21mN6Tma7+6SEDue2EzBXkRPj51aNBRRKQdWl3s7n568/s2/PwCYMEh0w5b6O5+bVveQ0RaZ39dEw+vLOeDkwbSp6Qw6Dgi0g7a3iYiPLBsMzWNca4/dVTQUUSknVpU7GZ2jZntNbNKM7vPzPQ9GJGQaIjF+fXzmzhtTF8mDOoedBwRaaeWrrF/EzgbGA9sBr7XYYlEJKP+/PJ29lQ18KnTtLYuEgYt3cd+wN1fSj3+ppkt66hAIpI5iYTzy2fKOHpgd04bo6u4iYRBS9fYB5rZHDM7zcz6ATqBtEgILH5jN2/urmbO6SN5rxNCiUjuaOka+7eBScAVwLFAiZktAF4BVrn7gx2UT0Q60J1PlzGoRxEfnDQo6CgikiYtKnZ3n9v8uZkNIVn0xwLnAyp2kRzzytZ9LNtYyTcuOJp8nZBGJDTadOY5dy8neQGXBe+3rIhkp7nPlNGtKI/Lpg0LOoqIpJH+TBfphLZU1PK313ZwxfThlBRm+szSItKRVOwindCvnisjGjGuO2VE0FFEJM1U7CKdzFs1jfxueTkXHT+Y0u5FQccRkTRTsYt0Mr9Zupm6pjhzTtcJaUTCSMUu0onUN8W5d8kmZo3rx9hSnRlaJIxU7CKdyEMry6moaWTO6aODjiIiHUTFLtJJJBLOXc9uZNKQHpw0qnfQcUSkg6jYRTqJx9fuYuPeGuacPkqnjxUJMRW7SCcx95kyhvTqwuxjBgQdRUQ6kIpdpBN4cVMlKza/xfWnjiRPp48VCTV9wkVCzt25beE6+nUr5OMn6vSxImGnYhcJuefXV7BsYyU3zjqKLgXRoOOISAdTsYuEmLtz22PrGNSjiMumDQ06johkgIpdJMSeXLubl7fu4/NnjqEwT2vrIp2Bil0kpBIJ578ff4MRfYq5ZMqQoOOISIao2EVCasFrO1i74wD/fNZY8nUkvEinoU+7SAjF4gl+9PgbjC0t4UPHDQo6johkkIpdJIT+9PJ2yvbU8C9njyUa0VnmRDoTFbtIyDTGEvz0yTeYOLg75+oscyKdjopdJGR+t3wrWyvr+NI543ROeJFOSMUuEiL1TXH+56k3mTK8FzPH9gs6jogEQMUuEiIPLNvCrgMN/KvW1kU6LRW7SEjUNMT4+aL1nHJUH2aM7hN0HBEJiIpdJCTuWbKJippGvnTOuKCjiEiAVOwiIbC/rok7n97AmeP7M3lYr6DjiEiAVOwiIXDH0xs4UB/jX84ZG3QUEQmYil0kx63fXcVdz5bxkcmDOWZQj6DjiEjAVOwiOczd+cafXqO4II+vnX900HFEJAuo2EVy2J9f3s7Sskq+MnscfUsKg44jIllAxS6So/bXNfHdv67huKE9ufzEYUHHEZEsoWIXyVG3LVxHZU0j//nhiUR0oRcRSVGxi+SgVeX7uH/ZZq6eMYKJg3XAnIi8I+PFbmazzWydma03s5sPM/8KM1uVui0xs+MynVEkm8UTztf/+Bp9Swr19TYR+QcZLXYziwK3A+cBE4DLzWzCIYttBM5w90nAd4C5mcwoku1+u2wzr27bzzcuOJruRflBxxGRLJPpNfZpwHp3L3P3RmAecFHzBdx9ibu/lXq6FBiS4YwiWWt3VT0/XLiOU47qw4XHDQo6johkoUwX+2Bga7Pn5alpR/JJ4G8dmkgkh3zvr2tpaErwHxdN1NXbROSwzN0z92ZmlwLnuvv1qedXAdPc/abDLDsL+DlwqrtXHGb+HGAOQGlp6ZR58+alLWd1dTUlJSVpe70gaSzZp63jWFsR5wcv1vOh0flcMqagA5K1Xlh+J6CxZKOwjAPSP5ZZs2atcPeph53p7hm7ATOAhc2e3wLccpjlJgEbgLEted0pU6Z4Oi1atCitrxckjSX7tGUcDU1x/8Bti/zUHzzpdY2x9Idqo7D8Ttw1lmwUlnG4p38swHI/QidmelP8i8AYMxtpZgXAZcD85guY2TDgYeAqd38jw/lEstIvny1jw54a/uPCiRTlR4OOIyJZLC+Tb+buMTO7EVgIRIG73X21md2Qmn8H8C2gD/Dz1D7EmB9pc4NIJ7CqfB8/feJNzps4gFnj+wcdR0SyXEaLHcDdFwALDpl2R7PH1wPXZzqXSDbaX9fE5367kr4lBXzv4mODjiMiOSDjxS4iLePu3PzQKnbsq+f/Pj2DXl2z44A5EcluOqWsSJa674XN/O21nXxl9jimDO8VdBwRyREqdpEstKp8H//517WcOb4/1586Kug4IpJDVOwiWab5fvXbLj1OV24TkVbRPnaRLHJwv/r2ffX87tMnab+6iLSa1thFsshvlqb2q587jinDewcdR0RykIpdJEu8Wr6f7/5lLR8Y359Pnab96iLSNip2kSxwoD65X71PSQH/rf3qItIO2scuErCD+9W37avTfnURaTetsYsEyN353oK1LHhV+9VFJD1U7CIB+tmT6/nlsxu5ZsZw5pyu/eoi0n4qdpGA3PVsGT9+4g0umTyEb3/oGFIXPRIRaRcVu0gA5v19C9/961rOmziAH1xyrA6WE5G00cFzIhm2bEeMO1a9yhlj+/GTy44nL6q/r0UkffR/FJEMenLtLuauauDE4b2548opFOZFg44kIiGjYhfJkCUb9vKZB1YyrFuEX107lS4FKnURST9tihfJgJVb3uL6e5czok8xnz8mQbei/KAjiUhIaY1dpIO9vHUf1979d/p1K+T+T06npEAHyolIx1Gxi3SgP75UzsfufIHuXfK5/5PT6d+9KOhIIhJy2hQv0gHiCeeHC1/nzqfLmD6yN7+4cgq9dapYEckAFbtImh2ob+ILD77EonV7uOqk4XzrQxPI11faRCRDVOwiaVS2p5rr71vOlopa/vPiiVwxfXjQkUSkk1Gxi6TJM2/s4cbfriQvGuGB66czfVSfoCOJSCekYhdpJ3fnV89t5HsL1jK2tBu/vHoqQ3sXBx1LRDopFbtIO+yuquff56/hr6/u4LyJA7jt0uPoWqiPlYgER/8HEmmDRMJ54O9b+OGjr9MQS/Dlc8fxmTNG62IuIhI4FbtIK63dcYCv/fFVXtqyj1OO6sN3P3wsI/t2DTqWiAigYhdpsdrGGD954k1+9dxGenbJ5ycfP56Ljh+k66iLSFZRsYu0wJNrd/GtP69m2746Lp82lK/OHk/PYp1wRkSyj4pd5D28Wr6fnz31Jo+v2cXY0hJ+f8MMThzRO+hYIiJHpGIXOYS78/z6Cn7x9HqeX19Bt8I8vnzuOD512igK8nQGORHJbip2kZR4wnn0tZ3c8fQGXt22n37dCrnlvPF8YvowXWZVRHKGil06vfqmOA+v3MbcZzawqaKWkX278v2PHMvFJwymKD8adDwRkVZRsUun5O6sKt/PI69s508vb2dvdQOThvTgF1dM5pxjBhDV99FFJEep2KVTeWNXFY+8sp1HXtnOpopaCqIRzhjXj2tPHsHJo/voq2sikvNU7BJ6WytrmZ8q89d3VhExOHl0Xz478yjOPWYAPYq1/1xEwkPFLqGzr7aRpWWVLC2rYMmGvbyxqxqAKcN78e8XHsP5xw6kX7fCgFOKiHQMFbvkvKr6Jl7cVMmS9RW8UFbBmh0HcIcu+VGmjujFJZOHcMGkgQzppSuuiUj4qdglpxyob2Lt9gOs2XGANan713dWEU84BXkRpgzrxRfPGsuM0X04bkhPfe9cRDodFbtkpdrGGFsqa9m0t5bXdyZL/KWNtex59LG3l+nTtYAJg7rz2ZmjmTG6D5OH9dLX00Sk01OxSyCa4gn2VDWwY389Wypr2FxRy5bKWrZU1LK5spY9VQ1vL2sGI/t2ZWSPCNedcRQTBnZnwsDu9OtWqKPYRUQOkfFiN7PZwE+BKHCXu//XIfMtNf98oBa41t1XZjqntI67U9sYZ19dE2/VNLK/rom3ahuprGlk14F6dh9oYFdVA7sP1LO7qoHKmsZ3/bwZDOhexLDexcwa14/hfboyrHcxw/sUc1T/EooL8li8eDEzZx4V0AhFRHJDRovdzKLA7cDZQDnwopnNd/c1zRY7DxiTuk0HfpG6lzRxd2IJpzGWoK4pTl1jnIZYnLrG5PP6pvjb99UNMWoaYlTXx6huiCcfp241DTH21zWxr66J/bVNNMYTh32/aMToV1JIafdChvQqZvLwXpR2K6J/90IGdC9iaO9ihvTqos3oIiJpkOk19mnAencvAzCzecBFQPNivwi4z90dWGpmPc1soLvvyETATXtrWLkrRsPqnbgDOO7gkLpv/tzfnpZIJKclUjMT7iQOzvPksonEwWmkHjvxZtPjqWkJd+IJiCcSb9/HUvNi8eTPxBPJx03x5LymeIJY3IklEjTGnVg8QVM8wf6qWqIvPEljPEFDLEFjLEFjPJEaW+uUFObRtTBK18I8uhXm0bUwj6P6l9CzOJ+exQX07JJPr+ICehTnJx93LaBXcQF9uhYQ0ZncREQyItPFPhjY2ux5Of+4Nn64ZQYDGSn2R1fv5GcvNcBLKzLxdv/ADKJmRMyIRoy8iBFJ3UcPueVFjPxohLyokReJkJ+6L8pPTs+PGvusjiGD+lKQF6EwL0JBXoSCaOTt510KohTlJ29dUrei/Mjb07oVJQu8OD+qchYRyQGZLvbDNcOh644tWQYzmwPMASgtLWXx4sXtDgcwoCHBV493iou7vB3m4AFadjCcvfPY7N33pB5H3rWMJacdMj/S/LG9+71a7uDm7/hh51b3ilFS8tbhf9SBhtQtJQZUp27Zprq6Om2/5yCFZRygsWSrsIwlLOOAzI4l08VeDgxt9nwIsL0Ny+Duc4G5AFOnTvWZM2emLWTyIK30vV6QNJbsE5ZxgMaSrcIylrCMAzI7lkyfveNFYIyZjTSzAuAyYP4hy8wHrrakk4D9mdq/LiIikusyusbu7jEzuxFYSPLrbne7+2ozuyE1/w5gAcmvuq0n+XW36zKZUUREJJdl/Hvs7r6AZHk3n3ZHs8cOfC7TuURERMJAJ9IWEREJERW7iIhIiKjYRUREQkTFLiIiEiIqdhERkRBRsYuIiISIil1ERCREzNtyma8sY2Z7gM1pfMm+wN40vl6QNJbsE5ZxgMaSrcIylrCMA9I/luHu3u9wM0JR7OlmZsvdfWrQOdJBY8k+YRkHaCzZKixjCcs4ILNj0aZ4ERGREFGxi4iIhIiK/fDmBh0gjTSW7BOWcYDGkq3CMpawjAMyOBbtYxcREQkRrbGLiIiEiIr9CMzseDNbamYvm9lyM5sWdKb2MLObzGydma02sx8Gnac9zOxfzczNrG/QWdrKzG41s9fNbJWZ/dHMegadqbXMbHbq39R6M7s56DxtYWZDzWyRma1NfTa+EHSm9jKzqJm9ZGZ/CTpLe5hZTzP7Q+pzstbMZgSdqa3M7Iupf1+vmdmDZlbUke+nYj+yHwL/7u7HA99KPc9JZjYLuAiY5O7HALcFHKnNzGwocDawJegs7fQ4MNHdJwFvALcEnKdVzCwK3A6cB0wALjezCcGmapMY8CV3Pxo4Cfhcjo6juS8Aa4MOkQY/BR519/HAceTomMxsMPB5YKq7TwSiwGUd+Z4q9iNzoHvqcQ9ge4BZ2uszwH+5ewOAu+8OOE97/Bj4CsnfT85y98fcPZZ6uhQYEmSeNpgGrHf3MndvBOaR/OMxp7j7DndfmXpcRbI8Bgebqu3MbAhwAXBX0Fnaw8y6A6cDvwJw90Z33xdoqPbJA7qYWR5QTAf3iYr9yP4ZuNXMtpJcw82pNapDjAVOM7NlZva0mZ0YdKC2MLMLgW3u/krQWdLsn4C/BR2ilQYDW5s9LyeHCxHAzEYAJwDLAo7SHj8h+YdvIuAc7TUK2AP8OrVb4S4z6xp0qLZw920kO2QLsAPY7+6PdeR75nXki2c7M3sCGHCYWV8HzgS+6O4PmdnHSP7leFYm87XG+4wlD+hFclPjicDvzGyUZ+FXIt5nHF8DzslsorZ7r7G4+59Ty3yd5ObgBzKZLQ3sMNOy7t9TS5lZCfAQ8M/ufiDoPG1hZh8Edrv7CjObGXCc9soDJgM3ufsyM/spcDPwzWBjtZ6Z9SK5NWsksA/4vZld6e73d9R7dupid/cjFrWZ3UdyXxXA78nyTVvvM5bPAA+nivzvZpYged7iPZnK11JHGoeZHUvyg/GKmUFy0/VKM5vm7jszGLHF3ut3AmBm1wAfBM7Mxj+y3kc5MLTZ8yHk6O4qM8snWeoPuPvDQedph1OAC83sfKAI6G5m97v7lQHnaotyoNzdD249+QPJYs9FZwEb3X0PgJk9DJwMdFixa1P8kW0Hzkg9/gDwZoBZ2utPJMeAmY0FCsixCyu4+6vu3t/dR7j7CJIf/MnZWurvx8xmA18FLnT32qDztMGLwBgzG2lmBSQPBpofcKZWs+Rfib8C1rr7j4LO0x7ufou7D0l9Pi4DnsrRUif1ud5qZuNSk84E1gQYqT22ACeZWXHq39uZdPCBgJ16jf19fAr4aepgh3pgTsB52uNu4G4zew1oBK7JwTXEsPlfoBB4PLUFYqm73xBspJZz95iZ3QgsJHmU793uvjrgWG1xCnAV8KqZvZya9jV3XxBcJEm5CXgg9YdjGXBdwHnaJLUr4Q/ASpK73V6ig89CpzPPiYiIhIg2xYuIiISIil1ERCREVOwiIiIhomIXEREJERW7iIhIiKjYRUREQkTFLiIiEiIqdhFpNTO7wcxeTt02mtmioDOJSJJOUCMibZY6z/pTwA/d/ZGg84iI1thFpH1+SvKc5Cp1kSyhc8WLSJuY2bXAcODGgKOISDPaFC8irWZmU4B7gdPc/a2g84jIO7QpXkTa4kagN7AodQDdXUEHEpEkrbGLiIiEiNbYRUREQkTFLiIiEiIqdhERkRBRsYuIiISIil1ERCREVOwiIiIhomIXEREJERW7iIhIiPw/T4JSui2fY2QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = np.linspace(-8, 8)\n",
    "fig, ax = plt.subplots(figsize = (8, 5))\n",
    "ax.plot(z, 1 / (1 + np.exp(-z)))\n",
    "ax.set_xlabel('z')\n",
    "ax.set_ylabel(r'$P = 1/(1+e^{-Z})$')\n",
    "ax.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f299756-04be-41ed-841c-c766549c8264",
   "metadata": {},
   "source": [
    "No matter what value we receive from $z$ (infinite domain), the final result will be transformed into a $[0, 1]$ range. The "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc9392a-1690-4f50-adf5-c301c90b1402",
   "metadata": {},
   "source": [
    "We usually denote $Z_i = \\beta_1+\\beta_2X_i$, ranges from $(-\\infty, +\\infty) $. This model solves the issues of probabilistic range and linear model. \n",
    "\n",
    "In order to estimate the logit model, we are seeking ways to linearize it. It might look confusing at first sight, here is the steps\n",
    "$$\n",
    "P_{i}=\\frac{1}{1+e^{-Z_{i}}}=\\frac{e^{Z_i}}{e^{Z_i}}\\frac{1}{1+\\frac{1}{e^{Z_i}}}=\\frac{e^{Z_i}}{1+e^{Z_i}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edaf88fd-f65f-4e7f-b7f6-4d0eb7772b7b",
   "metadata": {},
   "source": [
    "And $1-P_i$\n",
    "$$\n",
    "1-P_{i}=\\frac{1}{1+e^{Z_{i}}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "945b4e3a-86cf-4fde-8724-0a40adb0d238",
   "metadata": {},
   "source": [
    "Combine them, we have the odds ratios\n",
    "$$\n",
    "\\frac{P_{i}}{1-P_{i}}=\\frac{1+e^{Z_{i}}}{1+e^{-Z_{i}}}=e^{Z_{i}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9813b328-5ee6-4310-b69c-aec133cc1c48",
   "metadata": {},
   "source": [
    "Here's the interesting part, take natural log, we get a linearized model and we commonly call the log odds ratios the **logit**.\n",
    "$$\n",
    "\\underbrace{\\ln{\\bigg(\\frac{P_{i}}{1-P_{i}}\\bigg)}}_{\\text{logit}}= Z_i =\\beta_1+\\beta_2X_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a079095-8e4c-42e5-9793-0f29716bd0ea",
   "metadata": {},
   "source": [
    "To estimate the model, some technical procedures have to be carried out. Logit wouldn't have any meaningful interpretation, if we simply plug in the data $Y_i$. Because we don't observe $P_i$, the results are nonsensical.\n",
    "$$\n",
    "\\ln \\left(\\frac{1}{0}\\right) \\quad \\text{if a family own a house}\\\\\n",
    "\\ln \\left(\\frac{0}{1}\\right) \\quad \\text{if a family does not own a house}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f919ae45-a0b4-4aee-a805-3609c44933cf",
   "metadata": {},
   "source": [
    "One way to circumvent this technical issue, the data can be grouped to compute\n",
    "$$\n",
    "\\hat{P}_{i}=\\frac{n_{i}}{N_{i}}\n",
    "$$\n",
    "where $N_i$ is number of families in a certain income level, e.g. $[20000, 29999]$, and $n_i$ is the number of family owning a house in the that level."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd2c1e16-05de-4441-b1bd-f9f8a0bf3d33",
   "metadata": {},
   "source": [
    "However, this is not a preferred method, since we have nonlinear tools here. Recall that owning a house is following a Bernoulli distribution whose probability mass function is\n",
    "$$\n",
    "f_i(Y_i)= P_i^{Y_i}(1-P_i)^{1-Y_i}\n",
    "$$\n",
    "The joint distribution is the product of Bernoulli PMF due to their independence\n",
    "$$\n",
    "f\\left(Y_{1}, Y_{2}, \\ldots, Y_{n}\\right)=\\prod_{i=1}^{n} f_{i}\\left(Y_{i}\\right)=\\prod_{i=1}^{n} P_{i}^{Y_{i}}\\left(1-P_{i}\\right)^{1-Y_{i}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a54810d-170b-4a52-a167-ab7f36c0880d",
   "metadata": {},
   "source": [
    "In order to obtain log likelihood function, we take natural log\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\ln f\\left(Y_{1}, Y_{2}, \\ldots, Y_{n}\\right) &=\\sum_{i=1}^{n}\\left[Y_{i} \\ln P_{i}+\\left(1-Y_{i}\\right) \\ln \\left(1-P_{i}\\right)\\right] \\\\\n",
    "&=\\sum_{i=1}^{n}\\left[Y_{i} \\ln P_{i}-Y_{i} \\ln \\left(1-P_{i}\\right)+\\ln \\left(1-P_{i}\\right)\\right] \\\\\n",
    "&=\\sum_{i=1}^{n}\\left[Y_{i} \\ln \\left(\\frac{P_{i}}{1-P_{i}}\\right)\\right]+\\sum_{i=1}^{n} \\ln \\left(1-P_{i}\\right)\n",
    "\\end{aligned}\n",
    "$$\n",
    "Replace \n",
    "$$\n",
    "\\ln{\\bigg(\\frac{P_{i}}{1-P_{i}}\\bigg)}=\\beta_1+\\beta_2X_i\\\\\n",
    "1-P_{i}=\\frac{1}{1+e^{\\beta_1+\\beta_2X_i}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc20683-c09d-4bc0-a4c0-011cdaf86bd7",
   "metadata": {},
   "source": [
    "Finally we have log likelihood function of $\\beta$\n",
    "$$\n",
    "\\ln f\\left(Y_{1}, Y_{2}, \\ldots, Y_{n}\\right)=\\sum_{i=1}^{n} Y_{i}\\left(\\beta_{1}+\\beta_{2} X_{i}\\right)-\\sum_{i=1}^{n} \\ln \\left[1+e^{\\left(\\beta_{1}+\\beta_{2} X_{i}\\right)}\\right]\n",
    "$$\n",
    "Take partial derivative w.r.t to $\\beta_2$\n",
    "$$\n",
    "\\frac{\\partial LF}{\\partial \\beta_2} = n\\beta_2 -\\sum_{i=1}^n\\frac{e^{\\beta_1+\\beta_2X_i}X_i}{1+e^{\\beta_1+\\beta_2X_i}}=0\n",
    "$$\n",
    "And stop here, apparently there won't be analytical solutions, i.e. numerical solutions are necessary. And this is exactly what Python is doing for us below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "236bebec-4523-49a1-96d0-790582a98397",
   "metadata": {},
   "source": [
    "### <font face=\"gotham\" color=\"orange\"> Apartment Price Example</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e8bd7c9-002e-4ad7-8f69-b2caf214f113",
   "metadata": {},
   "source": [
    "Here we come back to the apartment owning example. The data is an excerpt of survey of citizens in Manchester, UK."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e9603d2f-3642-4387-bcc9-003230fd1dad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.257124\n",
      "         Iterations 8\n"
     ]
    }
   ],
   "source": [
    "df = pd.read_excel('Basic_Econometrics_practice_data.xlsx', sheet_name = 'HouseOwn_Quali_Resp'); df.head()\n",
    "\n",
    "Y = df[['Own_House']]\n",
    "X = df[['Annual_Family_Income']]\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "log_reg = sm.Logit(Y, X).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a0cede6-4022-4a54-9758-f18c255a0e53",
   "metadata": {},
   "source": [
    "Optimization of likelihood function is an iterative process and global optimization has been reached as the notice shows. \n",
    "\n",
    "Print the estimation results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9b570f8e-c175-4a87-8422-40009d4b518c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:              Own_House   No. Observations:                   38\n",
      "Model:                          Logit   Df Residuals:                       36\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Tue, 06 Sep 2022   Pseudo R-squ.:                  0.6261\n",
      "Time:                        14:04:45   Log-Likelihood:                -9.7707\n",
      "converged:                       True   LL-Null:                       -26.129\n",
      "Covariance Type:            nonrobust   LLR p-value:                 1.067e-08\n",
      "========================================================================================\n",
      "                           coef    std err          z      P>|z|      [0.025      0.975]\n",
      "----------------------------------------------------------------------------------------\n",
      "const                   -6.7617      2.189     -3.088      0.002     -11.053      -2.471\n",
      "Annual_Family_Income     0.1678      0.053      3.176      0.001       0.064       0.271\n",
      "========================================================================================\n"
     ]
    }
   ],
   "source": [
    "print(log_reg.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85c0e497-6cb6-44e0-afbd-4e6b2a212689",
   "metadata": {},
   "source": [
    "Here is the estimated model\n",
    "$$\n",
    "\\hat{P_i} = \\frac{1}{1+e^{6.7617-0.1678 X_i}}\n",
    "$$\n",
    "And the fitted curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d01de24f-057b-4b25-ac2c-d0cb5a1731e9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAE9CAYAAAAbGFuyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzDklEQVR4nO3deZxU1Zn/8c/T1d0sgoJsIosgEhBFRVFwb3ADN4xm0bhEjTrOxMQsY9RJfpPNmTgxZpJJNMYkmsUoajSKBtxtcUMRRVmUgCiyCYJsDQ29Pb8/bkFVN93Vt7uWW9X1fb9e9apz7r1166nT3fX0vffcc8zdERERkcJTEnUAIiIi0j5K4iIiIgVKSVxERKRAKYmLiIgUKCVxERGRAqUkLiIiUqBKow6grXr37u1DhgzJ2P62bt3KHnvskbH9FSu1Y/rUhulTG6ZPbZi+bLThnDlz1rl7n6bLCy6JDxkyhDfeeCNj+6usrKSioiJj+ytWasf0qQ3TpzZMn9owfdloQzNb1txynU4XEREpUEriIiIiBUpJXEREpEApiYuIiBQoJXEREZECpSQuIiJSoJTERUREClTW7hM3s7uAM4G17n5wM+sN+CVwOrANuNTd38xWPCLt8chbK7nlyUWs2ljNvj26cN1pIzhnzIC8jqmtMTfdfsLIPjz/3ieh69edNgIg5Xt+75F53PfacurdiZlxwbhB3HTO6HbH3NJnOH/QFr5783NZ/zmFibe1z5yJ98h0zG15fXO/B1H/bTSV7b/ffPh+yOZgL38Efg38uYX1k4Hh8cc44DfxZ5G88MhbK7nx4XlU19YDsHJjNTc+PA8gsi+r1mJqa8zNbX/PrI92rQ9Tv+5vb4NDbYM3+57fe2Reo9fUu++q33TO6LTbudHrB2X/5xQm3tY+cybeI9Mxt/X1TX8PWtyfO3gDNMSfdy5Lfm5U9qT10Kjg3vwymmzvzoz5q7n1H+9SU1dPb6B2fRU/f/BTOleNZNLB+7ThvZsKVjy1YA23PbmIkrp6BgGs38TtD6yl66YR7LllA3z6Cey92wBrGZe1JO7uM81sSIpNpgB/dncHZplZDzPr7+6rsxWTSFvc8uSiXV9aO1XX1nPLk4siS+KtxdTWmJvbvq1q63f/tkt+z/teW97s6+57bTk3nTM67XbO9c8pzPu19plb8/MnFlK+Yyt7ew2dvY5O1NO5to5n/76ac2IHQW3N7o+aHVBbC/W1UF8ff9RBfT2d31nBT3fUEqOBMhoo9QZKaWCPPzwEA7s32vbIqi3wzD1BMmtoAHfGb9zGsw0NlOCUuGM4JTix+LPhlDjEbneIWdJrG9Jv8HaaHH/s5k+Z2f+p8cdu7ok/v/kMfP+2zLxZCube4r8b6e88SOKPt3A6/XHgZnd/KV5/Frje3XcbU9XMrgKuAujXr98RU6dOzViMVVVVdOvWLWP7K1YdsR3nrdzU4rrRA/bK+PuFacPWYmprzKm2z4QwMaXbzsmv79cF1lS37fVtFSbe+Ss20Kmuhs612xs9utRuZ1BnKKupJlazg9K6GmK1O4jV1lBaW0OstoZYXQ2x+rqMxy25taVnX+accXnG9jdhwoQ57j626fIok/g/gJ80SeLfcfc5qfY5duxY19jp+acjtuOxNz/Hyo3Vuy0f0KMLL98wMePvF6YNW4uprTG3tH0m7HzPYTdOp76Z75mYGe//5PS02zn59d8eXcet80rb9Pq2mvDfT1Hy6RoGNmyhX0MV/Rq20q9hK/vFtnN8H4MN66nftIFYy+dji4+VQIkBBmbxZU2eG5WTt4NGhea2N0taHzx9uq2W+gYHLHGW3YyYQe/unVt576bLdvtArNxYHd9/fN/x59KSEvbuDF32HwFf/2ELr287M2s2iUc5AcoKCC4lxA0EVkUUi8hurjttRKPrgABdymK7OnLlY0xtjbm57duqLGaNrok3fc8Lxg1qdP10pwvGDWpXzGE+Q9o/p7paWLMSVn4Ia1fB2tWw7mP4ZDXPb1jX8uu2Bk+x9r8zECSbLZSx1cqopowdFqOmpJQBfXrQu2c3KO8EZWVQ1gnKypMeZVBaDrESiJXGHzHmrNjCfXNWsq0e6iih1kooKS3jsuOHcdyIfSAW27Xt63PmcNS48UECiyffpxau5eYnF7G11mkwwzEaMBog/lxCp7IYP5hyMGeNGZj02pIUiTC7Zja5jg/B78VPzh2dkcsss1Psv8emxTk7qIkyiU8DrjGzqQQd2jbperjkk51/6FH3Pm1LTG2NubntM907fec14JZ6aqfbzsmvhy0MaOvPafNG+GARLF8aJO2VH8KaFcF14jRUl3ZmbUM5G0q6sLGkMz369OKwAwfDHt2h257QtRt06gKdu0DnrtCl666ylXfiubmrdmuTQ9v5u3cEsHzk7j2pj2tmf9ve/wj6D2q07NQT+7Ntz96t9k4/K496p2f77zfV/isrF2fkPcLI2ul0M7sPqAB6A2uA7wNlAO5+R/wWs18DkwhuMbusuevhTel0en5SO6ZPbZi+Vtuwvh6WLYYlC2DpIvhwEaxb07Y3sRLo1Qd67wM9e8NevaBnL+jRC3rsHTzvtTeUlqX1WaKi38P0ZWkq0tyeTnf3C1pZ78BXs/X+IiI0NARH1u/Ohffehn/Og+3bwr22dz/YdwjsMxD69Ie+/aF3f+jVF0qjPIkpkqDfRBHpWGpr4N23YO4seHsWbNqQevuychh8AAwZDgOGwMChsO/g4BS3SJ5TEheRwldXC/PfYNSLj8CDv4Ad21vetmdvGHko7H8gDB0RJG0dWUuB0m+uiBQm96BD2ivPwOwXYOsW+ja3Xbe94MDDgsQ98rDgtHhEPaZFMk1JXEQKS80OeP0FeP6xoJNac/YZCIeOh8OOhmEjoSTdm75E8pOSuIgUhk0b4Om/w4szYOuW3dfv3Zdl/Yex3xe+HFzbFikCSuIikt82rIMnHoSZM4JOa8lKy+CoE+HYU2H4wXwwcyb7KYFLEVESF5H8tGkDPP5XePHJoONast77QMUZcNxpwcApIkVKSVxE8suO7cFp8xkPwI4mY6oPPgDOvAAOG6/r3CIoiYtIvnCHN1+GqXcEp9CT7T8SzvwSjD5SPctFkiiJi0j01q+Fe2+Dt19rvHzfwfD5K+HgsUreIs1QEheR6LgHt4o9dFfjAVq694BzLobjJgUzbIlIs5TERSQamz6Fu38O85tMaHTi6XDe5cEsXyKSkpK4iOTe3Fnwx/+Fqk2JZfsOhouvheEHRReXSIFREheR3Gmoh0f+DNPvb7z81PPgs18OJiMRkdCUxEUkN6o2w503w8I3E8t69obLvw0HjokuLpECpiQuItm38kP41fdh3ZrEsoOOgCuv12AtImlQEheR7Fr0Dtz2I9hWlVh25gVw9kUasEUkTUriIpI9r1fCXbcmhk3t1AWuuA7GHBNpWCIdhZK4iGTHc9Pg3tsT9b16wtd/BPsNjy4mkQ5GSVxEMu+ph+GBOxP1/oPh2h8FE5eISMYoiYtIZs14EB76Q6I+7ED42o+gW/foYhLpoJTERSRzZjwQDKG60/CDgyPwzl2ji0mkA1MSF5HMeOEfjRP4iEOCa+CdOkcXk0gHVxJ1ACLSAbzxItzz60R95KFK4CI5oCQuIulZ+Cb87n+CGckAhnwGrvm+ErhIDiiJi0j7rfgAbvsx1NcF9X0GwrU/1jVwkRxREheR9tm8MRhKdUd1UO/ZG77539B9r0jDEikmSuIi0na1NXD7j2D92qDeqQt84ybo1TfauESKjJK4iLSNO/zlV7BkYVA3g6tugAFDIg1LpBgpiYtI2zz/GLzydKL+ua/AoeOii0ekiCmJi0h4HyyC+5OGUz3mZDj1vOjiESlySuIiEs7WLfDb/070RB88DC7+enA6XUQioSQuIq1zh7t/DuvWBPUuXeHq70JZebRxiRS5lMOumtlIYAowAHBgFTDN3d/NQWwiki+efRTmvpqoX/ot6LtvdPGICJDiSNzMrgemAga8DsyOl+8zsxtyE56IRG7VMvhb0qxkJ02BI46LLh4R2SXVkfhXgIPcvTZ5oZn9HFgA3JzNwEQkD9TVwu9vCZ4BBg2Dz18RbUwiskuqa+INQHPny/rH14lIR/f4ffDRkqBcWgZXXBc8i0heSHUk/g3gWTNbDCyPLxsMHABck+W4RCRqS9+D6VMT9c9eqgFdRPJMi0nc3Z8ws88ARxF0bDNgBTDb3etzFJ+IRKG2Bu76GTTET7p9ZjSc8tloYxKR3aTsne7uDcCsHMUiIvli+v3w8Yqg3LkrXP5tKNEdqSL5Jqt/lWY2ycwWmdmS5nq0m9leZvaYmb1tZgvM7LJsxiMiIaxaFiTxnc67DHrvE108ItKirCVxM4sBtwGTgVHABWY2qslmXwUWuvuhQAVwq5lp9AiRqDQ0wJ9/mRiVbdgoOPGMaGMSkRZl80j8KGCJuy919xqCe86nNNnGge5mZkA34FOgLosxiUgqL85IzE4WK4VLvq7T6CJ5LNRfp5k9kPwc0gASvdoh6BQ3oMk2vwYOJBgJbh5wbfw6vIjk2qZP4W93JeqTP6/e6CJ5LmXHtiQHxJ+Ht2Hfzc2K4E3qpwFzgYnAMOBpM3vR3Tc32pHZVcBVAP369aOysrINYaRWVVWV0f0VK7Vj+qJuwxGv/oP+1VsB2Na9J290H0BDgf1Mo27DjkBtmL5ctmHYJN4eK4BBSfWBBEfcyS4DbnZ3B5aY2QfASIJhXndx9zuBOwHGjh3rFRUVGQuysrKSTO6vWKkd0xdpGy5dBPfM21XtesW/c8JBR0QTSxr0e5g+tWH6ctmG2bzYNRsYbmZD453VzgemNdnmI+AkADPrB4wAlmYxJhFpqqEBpv4mUT/saCjABC5SjLJ2JO7udWZ2DfAkEAPucvcFZnZ1fP0dwI+BP5rZPILT79e7+7psxSQizXj12WB0NgiGVP3CldHGIyKhhU3izV3fbpW7TwemN1l2R1J5FXBqe/YtIhlQvRUeTurMdup5mmJUpICEPZ1+S5NnEekIZjwImzYE5R694PQvRhuPiLRJqCTu7vcmP4tIB7BhHTzz90T9vMuhc5fo4hGRNtMoDiLFato9ULMjKA8+AMZNiDYeEWkzJXGRYrRqGbz0VKL+ucs1MptIAdJfrUgxevhu2Dk44kGHw6jDo41HRNql1d7pZvatZhZvAua4+9yMRyQi2bV4PsyNzzBsBud9Jdp4RKTdwhyJjwWuJhj3fADB8KcVwO/M7DvZC01EMs49OArfadwEGDwsunhEJC1h7hPvBRzu7lUAZvZ94G/ACcAc4KfZC09EMmrhW7B4QVCOlcI5l0Qbj4ikJcyR+GCgJqleC+zn7tXAjqxEJSKZ5w7T/pKoH3ca9N4nunhEJG1hjsTvBWaZ2aPx+lnAfWa2B7Awa5GJSGYtmAPvvxuUS8vgjPOjjUdE0tZqEnf3H5vZDOBYguFXr3b3N+KrL8xmcCKSIe7wyJ8T9RMmw959ootHRDIi7NjpbxFMI1oKYGaD3f2jrEUlIpn1zuvw4T+DcmmZhlcV6SDC3GL2NeD7wBqgnuBo3IFDshuaiGSEezA6204VZwbjpItIwQtzJH4tMMLd12c7GBHJggVzYNnioFzeCSZ/Ptp4RCRjwvROX04wuIuIFKLp9yfKx0+CvfaOLhYRyagwR+JLgUoz+wdJt5S5+8+zFpWIZMaShfDPeUE5FgvmCxeRDiNMEv8o/iiPP0SkUCQfhY+bCL36RheLiGRcmFvMfpiLQEQkw5YvhXdeC8pmMPkL0cYjIhnXYhI3s1+4+zfM7DGC3uiNuPvZWY1MRNIz44FE+fBjof+g6GIRkaxIdSS+c3zGn+UiEBHJoDWrYPbMRH2y7gsX6YhaTOLuPif+/ELuwhGRjHjywcR84aMOhyHDo41HRLIizGAvxwI/APaLb2+Au/v+2Q1NRNpl0wZ45ZlEXaOziXRYYXqn/wH4JsG0o/XZDUdE0lb5ONTVBuWhI2CEBlcU6ajCJPFN7j4j65GISPpqdsDzjyfqp54X9EwXkQ4pTBJ/3sxuAR6m8WAvb2YtKhFpn1nPQVV8gMW9+wa90kWkwwqTxMfFn8cmLXNgYubDEZF2a2iAp/+eqJ88JRilTUQ6rDCDvUzIRSAikqYFc2B1fIbgzl3huEnRxiMiWRemd/r7wCzgRWCmuy/MelQi0nbJR+HHnwZd94guFhHJiTCzmI0Cfgv0An5mZkvN7O+tvEZEcmnFB7Aw3k3FSuCkKdHGIyI5ESaJ1wO18ecGYA2wNptBiUgbJR+FH34M9N4nulhEJGfCdGzbDMwDfg78zt3XZzckEWmTzRvhtecT9VPPjSwUEcmtMEfiFwAzgX8DpprZD83spOyGJSKhzZyRGNxl/5EwbFS08YhIzoTpnf4o8KiZjQQmA98AvgN0yW5oItKq+np44R+J+kRNLihSTFo9Ejezh+I91H8JdAMuAXpmOzARCWHuq7BhXVDuvhcccVy08YhIToW5Jn4z8Ka7a9x0kXyTPMTqCZOhrDy6WEQk51ImcTPrC5wJXG9mDiwEbnf3NbkITkRSWLUM3psblK0ETjwj0nBEJPdaPJ0en4J0drz6Z+CeePm1+DoRidLzjyXKY46GvftEF4uIRCLVkfitwDnu/lbSskfjA738lsSY6iKSa9Vb4ZVnE3V1aBMpSqk6tu3ZJIED4O5zge5Zi0hEWvfqs7CjOijvO1hzhosUqVRJ3Mxst17oZrZ3K68TkWxyh+eSTqVPOEtzhosUqVTJ+H+Bp8zsRDPrHn9UADPi61plZpPMbJGZLTGzG1rYpsLM5prZAjN7oa0fQKTovDcXPl4elDt3haM19pJIsWrxmri732lmq4AfAwcRzCG+ELjJ3R9r6XU7mVkMuA04BVgBzDazacmzoJlZD+B2YJK7fxTvDS8iqSQfhR9zcpDIRaQopbzFzN0fBx5PtU0KRwFL3H0pgJlNBaYQ/COw05eAh939o/j7aWIVkVQ2rIO5sxL1CWdFF4uIRC6b17YHAMuT6iviy5J9BuhpZpVmNsfMLsliPCKF7+WnwBuC8ohDoP+gaOMRkUiFGbGtvZrraePNvP8RwEkEY7G/amaz3P2fjXZkdhVwFUC/fv2orKzMWJBVVVUZ3V+xUjumr9U2dGf804/QOV5d2HsIa9Xmjej3MH1qw/Tlsg1bTeJmNtTdP2htWTNWAMmHCQOBVc1ss87dtwJbzWwmcCjQKIm7+53AnQBjx471ioqK1sIOrbKykkzur1ipHdPXahvOfwO2bg7K3fZk1EVXMErDrDai38P0qQ3Tl8s2DHM6/aFmlv0txOtmA8PNbKiZlQPnA9OabPMocLyZlZpZV4IBZN4NsW+R4jNzRqJ89EkaJ11EWj4Sj089ehCwl5mdm7RqT9h1Rq9F7l5nZtcATwIx4C53X2BmV8fX3+Hu75rZE8A7QAPwe3ef3/6PI9JBbfoU3k7q0Hb85OhiEZG8kep0+giCyU96AMldYLcAV4bZubtPB6Y3WXZHk/otwC1h9idStF5+Opg7HGD4QcEobSJS9FLdJ/4owVjpR7v7qzmMSUSSNTTAi08k6joKF5G4MNfEl5vZ381srZmtMbOHzGxg1iMTkcCit+GT1UG5azcYe3y08YhI3giTxO8m6JC2L8F93o/Fl4lILsxMOgofPxHKO0UXi4jklTBJvK+73+3udfHHHwFNXCySC1s2wpsvJ+on6FS6iCSESeKfmNlFZhaLPy4C1mc7MBEhmDO8vi4o7z8SBg6NNh4RySthkvjlwBeAj4HVwOfiy0Qkm9zhxaR7w3UULiJNtDpiW3xykrNzEIuIJFs8Hz5eEZQ7d4UjT4w2HhHJO6kGe/nPFK9zd/9xFuIRkZ1eSBpiYdwE6NTqGEsiUmRSHYlvbWbZHsBXgF4E84yLSDZUbYE5LyXqJ0yKLhYRyVupBnu5dWfZzLoD1wKXAVOBW1t6nYhkwKxnoa42KO83PHiIiDSR8pq4me0NfAu4EPgTcLi7b8hFYCJFy73xZCc6CheRFqS6Jn4LcC7BFKCj3b0qZ1GJFLOl78KqZUG5vBMcVRFpOCKSv1LdYvZtglHavgesMrPN8ccWM9ucm/BEitALSUfhR1VAlz0iC0VE8luqa+Jh7iEXkUzathXemJmo695wEUlBiVokn7z2HNTsCMoDh8LQEdHGIyJ5TUlcJF/s1qFtMphFF4+I5D0lcZE80X39x7B8aVApKw9mLBMRSUFJXCRP9F8yN1EZe0Iwd7iISAptTuJm9oyZzTCzM7MRkEhR2r6Nfh8uTNTVoU1EQmh1ApRmXAL0B8ZnOBaR4vX6C8R2jtDWfzAcMCraeESkILQ5ibv7KmAVMCfz4YgUqZlJk52oQ5uIhNSua+JmNqP1rUQklI+WwIeLg3JpGRx9UrTxiEjBSDXs6uEtrQIOy0o0IsVo5hOJ8uHHQrc9o4tFRApKqtPps4EXCJJ2Uz2yEo1IsdmxPRjgZSd1aBORNkiVxN8F/sXdFzddYWbLsxeSSBF5YyZUbwNgW/eedB1xSMQBiUghSXVN/Acp1n8t86GIFKGkEdpWH3CYOrSJSJukmgDlbynWPZKVaESKycoP4f13g3KslI+HHcywSAMSkUKjEdtEopI8TvqYo6ntrClHRaRtlMRFolCzA159NlE/Xh3aRKTtlMRFojDnRdhWFZT79IcDD4s0HBEpTKGSuJmNTH4WkTQln0o/fhKU6P9pEWm7sN8c9zZ5FpH2WrUMFi8IyrEYHHtKtPGISMFq67//uv9FJF3JR+GHHQ177R1dLCJS0HQOTySXanbAK88k6hqhTUTSoCQukktzXkp0aOu9Dxw4Jtp4RKSgtTWJe1aiECkWyafST1CHNhFJT9hvEGvyLCJttWoZLJ4flGMxOPbUaOMRkYIXNokf3+RZRNoqecrRQ8erQ5uIpC1UEnf3quRnEWmj2hp45elE/cTTo4tFRDoMXZATyYVGHdr6qUObiGSEkrhILjQaoW2yOrSJSEa0+k1iZteGWdbCayeZ2SIzW2JmN6TY7kgzqzezz4XZr0hBWfUR/HNeUNYIbSKSQWEOB77czLJLW3uRmcWA24DJwCjgAjMb1cJ2/wM8GSIWkcLzYtJR+KHjoUev6GIRkQ6ltKUVZnYB8CVgqJlNS1rVHVgfYt9HAUvcfWl8f1OBKcDCJtt9DXgIOLINcYsUhtoajdAmIlnTYhIHXgFWA72BW5OWbwHeCbHvAcDypPoKYFzyBmY2APgsMBElcemI3pgJW7cE5d79YNTh0cYjIh1Ki0nc3ZcBy4Cj27nv5gaGaTri2y+A69293qzlcWTM7CrgKoB+/fpRWVnZzpB2V1VVldH9FSu1Y/PGPHEve8XLSweO4KOZM1vcVm2YPrVh+tSG6ctlG6Y6EgfAzMYDvwIOBMqBGLDV3fds5aUrgEFJ9YHAqibbjAWmxhN4b+B0M6tz90eSN3L3O4E7AcaOHesVFRWthR1aZWUlmdxfsVI7NmPZYlgX/5UvLWP/S/6N/ffs0eLmasP0qQ3TpzZMXy7bsNUkDvwaOB94kCDpXgIcEOJ1s4HhZjYUWBnfx5eSN3D3oTvLZvZH4PGmCVykYD3/WKJ8xHGQIoGLiLRHmCSOuy8xs5i71wN3m9krIV5TZ2bXEPQ6jwF3ufsCM7s6vv6OdAIXyWtVW+C1ykR94lmRhSIiHVeYJL7NzMqBuWb2U4LObnuE2bm7TwemN1nWbPJ290vD7FOkILz8VNAzHWDwMNj/wGjjEZEOKcx94hcTHElfA2wluM59XjaDEiloDQ1Q+XiiPuEsSNFxU0SkvVo9Eo/3UgeoBn6Y3XBEOoAFc+CT1UG5azc4qiLScESk40o12Ms8dr8lbBd3PyQrEYkUuuQObceeCp06RxeLiHRoqY7EP09w9C0iYX3yMcybnahXnBldLCLS4aVK4ve6++Fm9hd3vzhnEYkUssrHweMnsA4eC/32jTYeEenQUiXxcjP7MnCMmZ3bdKW7P5y9sEQKUM0OeClpHp8Juq1MRLIrVRK/GrgQ6AE0/TZyQElcJNms5xLjpPfqC6PHRhuPiHR4qcZOfwl4yczecPc/5DAmkcLjDs88kqhPnAIlscjCEZHi0Op94krgIiEsfAtWxe/G7NQZjj8t2nhEpCiEGexFRFrz9N8T5WNPDe4PFxHJMiVxkXSt+gjmx28rM4OTz4k0HBEpHqEmQDGznsBwYNeoFe7e8sTIIsXk2UcT5UPHQ1/dViYiuRFmPvErgGsJ5gOfC4wHXgUmZjUykUJQtRlefSZR11G4iORQmNPp1wJHAsvcfQIwBvgkq1GJFIqZ04P7wwEGDYMRGo1YRHInTBLf7u7bAcysk7u/B4zIblgiBaCuDp5LGif9lM9qtjIRyakw18RXmFkP4BHgaTPbAKzKZlAiBeH1Sti4Pijv2ROOPCHScESk+ISZivSz8eIPzOx5YC/giaxGJZLvGhrgiQcS9YlnQ1l5dPGISFEK1Tt9J3d/IVuBiBSUd14Lbi0D6NwVJmqcdBHJPd0nLtJW7jD9/kT9xNM1uIuIREJJXKStFs+Hpe8F5dKyoEObiEgElMRF2mp60rXwY06GHr2ii0VEilqYwV62EEw9mmwT8AbwbXdfmo3ARPLS8qWNh1g97XPRxiMiRS1Mx7afE9xSdi9gwPnAPsAi4C6gIlvBieSdGUlH4UccB/0GRBeLiBS9MKfTJ7n7b919i7tvdvc7gdPd/X6gZ5bjE8kfa1fB7KQpAyZ9IbpYREQIl8QbzOwLZlYSfyR/czU9zS7ScU2fCt4QlA8cA0OGRxuPiBS9MEn8QuBiYC2wJl6+yMy6ANdkMTaR/PHJanglaaKTs74UXSwiInFhRmxbCrQ0ksVLmQ1HJE/9475glDaAkYfBZ0ZHGo6ICITrnd4HuBIYkry9u1+evbBE8kjTo/CzL4wuFhGRJGF6pz8KvAg8A9RnNxyRPKSjcBHJU2GSeFd3vz7rkYjkIx2Fi0geC9Ox7XEzOz3rkYjko2n36ChcRPJWmCR+LUEirzazzWa2xcw2ZzswkcgtXwqznkvUz74oulhERJoRpnd691wEIpJ3/v7HYMYygEOOgs8cHGk4IiJNtZjEzWyku79nZoc3t97d38xeWCIR++d8eOf1oGwG514WbTwiIs1IdST+LeAq4NZm1jkwMSsRiUTNHR76Q6I+fiIMHBpdPCIiLWgxibv7VfHnCbkLRyQPzJ0F778blEvLYMrF0cYjItKCMLeYYWbHsPtgL3/OUkwi0amrg4fuStQrzoDe+0QXj4hICmFGbPsLMAyYS2KwFweUxKXjqXwcPl4elDt3hdPPjzYeEZEUwhyJjwVGubtmLJOObcum4L7wnc68APbsEVk4IiKtCXOf+HxA5xOl45v2F9hWFZT79IeTpkQbj4hIK8Ik8d7AQjN70sym7XyE2bmZTTKzRWa2xMxuaGb9hWb2Tvzxipkd2tYPIJIRKz6AyumJ+hevgrLy6OIREQkhzOn0H7Rnx2YWA24DTgFWALPNbJq7L0za7APgRHffYGaTgTuBce15P5F2c4f7fwseH171wDFw6PhoYxIRCSFMEh8GvOjui9u476OAJfH5yDGzqcAUYFcSd/dXkrafBQxs43uIpG/2C/Du3KBsJXD+vwQDvIiI5LkwSXwIcJGZ7QfMIZiW9EV3n9vK6wYAy5PqK0h9lP0VYEaIeEQyZ1sVTP1toj7xLBgwJLJwRETawsJ2OjezLsCVwL8DA9w91sr2nwdOc/cr4vWLgaPc/WvNbDsBuB04zt3XN7P+KoLR4+jXr98RU6dODRVzGFVVVXTr1i1j+ytWhdqOw197kgGL3wJgR5duvH7WldSXd4oklkJtw3yiNkyf2jB92WjDCRMmzHH3sU2Xh7lP/HvAsUA34C2CJP5iiPdcAQxKqg8EVjWz/0OA3wOTm0vgAO5+J8H1csaOHesVFRUh3j6cyspKMrm/YlWQ7fj+QrjnrV3VTpd+g+OPOC6ycAqyDfOM2jB9asP05bINw5xOPxeoA/4BvADMcvftIV43GxhuZkOBlcD5wJeSNzCzwcDDwMXu/s+2BC6Slro6+PP/JeqHjIPDj40uHhGRdggzFenhZtYdOI6gp/nvzGyNu6c8ZHH3OjO7BngSiAF3ufsCM7s6vv4O4D+BXsDtFnQkqmvudIFIxs24H1Z+GJTLO8GF/6bObCJScMKcTj8YOB44kWD0tuWEO52Ou08HpjdZdkdS+QrgijbEK5K+ZYvh8XsT9SmXQK9+0cUjItJOYU6n/w/BafT/A2a7e212QxLJotoauOtWqI9PAzBsFJxyTqQhiYi0V5gkfh5wAMGkJzFASVwK12N/bXwa/fJvQ0nKGy1ERPJWi8Oumlmpmf2U4PT5n4B7gOVm9lMzK8tVgCIZ8/67MOPBRP28y6HfgOjiERFJU6qx028B9gb2d/cj3H0MwehtPYCf5SA2kczZVgW/uzkxtOrIQ2HCWdHGJCKSplRJ/EzgSnffsnOBu28G/hU4PduBiWSMO/zpF7BuTVDvsgdc+k0oCTP/j4hI/kr1LebNzSHu7vUE18dFCkPl4zDnpUT90m9Cb82uKyKFL1USX2hmlzRdaGYXAe9lLySRDProfbj/zkR9wlkQ4ahsIiKZlKp3+leBh83scoKJTxw4EugCfDYHsYmkp2oz/OYmqIvfUDFoGHzhymhjEhHJoBaTuLuvBMaZ2UTgIMCAGe7+bK6CE2m3+nr47X/DJ6uDeqcu8C83Qll5tHGJiGRQmGFXnwOey0EsIpnzwJ2JOcIBrrgO9tF09SLSsah7rnQ8Lz0Jzz6aqE+5GMYcE108IiJZoiQuHcuCOfCXpNnJjjgOzrggunhERLJISVw6jg8Xw+03JcZFHzAELvu27gcXkQ5L327SMaxdBb/8f7CjOqj37A3X/hg6d4k2LhGRLFISl8K3cT3873dhy8ag3rUbfPO/YO8+kYYlIpJtSuJS2Dauh1u+k7iVrKwcvvZD2He/aOMSEckBJXEpXBvXw8+uhzUrg3osFtwLPvygaOMSEckRJXEpTBvXw89ugI9XBPVYDK66EQ47Otq4RERyqNXBXkTyzsfLg2vg69cG9ZKSIIFrTHQRKTJK4lJYlr4H//efwbjoECTwf/kPJXARKUpK4lI43n4tGA+9ZkdQL+8E//o9GH1ktHGJiERESVzynztMvx8e+VNQBui2F3z9R7D/iGhjExGJkJK45Lft1XD3rTDnpcSy3v3gGzfBPoOii0tEJA8oiUv+WvEB3HkzrFqWWPaZ0fCv34XuPSILS0QkXyiJS/5xh+cfgwd+B3W1ieUTz4YvXAWl+rUVEQElcck3G9cHs5C9/VpiWXkn+NJX4bhTo4tLRCQPKYlLfmhoCOYBf/D3UL01sXzQ/nDlDbDv4OhiExHJU0riEr2VH8K9t8OidxovP/kcOO/yYDx0ERHZjZK4RKdqMzz6F3jhH8GR+E5994VLroWRh0YXm4hIAVASl9zbsR2eewxm3A/bqhLLS0rg1PPg7IuC6+AiIpKSkrjkzo7twVH3jAcTc3/vNPIw+OJVwTVwEREJRUlcsm/Tp8EtY5XToWpT43V9+ge3jR02HsyiiU9EpEApiUt2uMOyJfDcNHi9svH93gB794EzLoBjT4HSskhCFBEpdErikllbNsKs54PbxVZ+uPv6Xn3htM/B8ZPU61xEJE1K4pK+rVvo9/48mPc8zH8D6ut332bYgXDKuTDmGIjFch+jiEgHpCQu7fPJapg/B+a+Cu/N5cDmEnd5p2Ce74ozgyQuIiIZpSQu4WzdAv+cBwvehAVzgiTekmEHwrGnwpEnQJc9chejiEiRURKX3TU0wMcr4P2F8P67wWP1RylfsrlXf/asmAyHHwf99s1RoCIixU1JvNhVb4UVH8KKpbB8aTD958oPg3u6UynvFIyoNupwGHMMb85bSEVFRQ4CFhGRnZTEOzp3qN4G6z6Gtatg7UpYE39euwo2bQi3n1gMBg2DUWOCxD3swCa9yxdmJXwREWmZknihamiAbVuDwVO2bArGId+4Hjasg43rguedj9aOqpvTvQcMGwnDRgWP/Q6ATp0z/jFERKT9sprEzWwS8EsgBvze3W9ust7i608HtgGXuvub2Ywpcu5QXwc1NVBbAzXbYXt1cLS8fWv8eVvj5+ptUL2VD5Z9TO2mjfRsqKanb6cUTz+eWCnsMxAGDg2GPB20PwzcH/bq2ezm33tkHve9tpx6d2JmXDBuEDedM3rX+kfeWsktTy5i1cZq9u3RhetOG8E5Ywa0uH7CyD48/94nWas3ff+2au3zNqe1NpDCp5+x5IusJXEziwG3AacAK4DZZjbN3ZPPu04Ghscf44DfxJ+zr6Ee3nqVvh8uhFfrg3ubG+qDBFsfr9fXxZc1V0/abufr6uqCxFxbAzU7EuVGy2rBG1qPrxlD2/lRqyllVawby0t7MnDkARxw0HDoNwD6DoBefaAk3H3b33tkHvfMSnRwq3ffVT+5R/DFduPD86iuDW43W7mxmhsfngfAOWMGNLs+eX/ZqCe/f1ul+rwtJfLW2kAKn37Gkk+yeSR+FLDE3ZcCmNlUYAqNL55OAf7s7g7MMrMeZtbf3VPcv5QhDQ3wm5sYBfBS1t8tKzZbORusM/sN2Re67Ql79oCefaBnb7799ErmV5exuqQbm63TrnHJB2ztwssTJ7br/e57bXmLy08+rSu3PLlo1xfbTtW19dzy5CLOGTOg2fXZlvz+bZXq87aUxFtrAyl8+hlLPrEgf2Zhx2afAya5+xXx+sXAOHe/Jmmbx4Gb3f2leP1Z4Hp3f6PJvq4CrgLo16/fEVOnTk0/QHcq/vo/6e+nHRqshIbSUhpKSmkoLaWurBN1ZZ2oL+tEfVk5dWXl1Jd1avT8wZZ6tpd1orqsC9Vlndle1pmG+BH06AF77fYe81Zu2m3ZTs1tH0aqfQ7dK8YHm1pO0KMH7JXy9dnWns/cnjZMp92rqqro1q1buOCkWblow2z8beUT/R6mLxttOGHChDnuPrbp8mweiTc3JVXT/xjCbIO73wncCTB27FjP2K1M773M2nXr6LtP/6D3daw0OLUci7Vcj8Xiy0oTy5Pr5Z2CXttl5U3K5VAW1EtiMUraGOoZN06nvpl/uGJmvH9hxW7Lv3vzc6zcWL3b8gE9uvC1ZrYP4yspYvjDaV2ZOr8h5Xu2FFO2tfczp/q8zbU5pNfulZWVuk0vTblow2z8beUT/R6mL5dtmM0kvgIYlFQfCKxqxzbZ86/fY2FlJX0L4Bf2gnGDGl2fTV7enOtOG9Houh1Al7IY1502IksxrG/1PZtbn23pfOa2tjlkp90lv+hnLPmkrQeEbTEbGG5mQ82sHDgfmNZkm2nAJRYYD2zKyfXwAnTTOaO5aPxgYvFr2zEzLho/uMVrs+eMGcBPzh3NgB5dMIKjhJ+cOzqta3atxdDaeza3/qLxg7NaT+czt7XNw7SBFD79jCWfZO2aOICZnQ78guAWs7vc/b/M7GoAd78jfovZr4FJBLeYXdb0enhTY8eO9TfeSLlJm+jUUWaoHdOnNkyf2jB9asP0ZaMNzSzn18Rx9+nA9CbL7kgqO/DVbMYgIiLSUWXzdLqIiIhkkZK4iIhIgVISFxERKVBK4iIiIgVKSVxERKRAKYmLiIgUKCVxERGRApXVwV6ywcw+AZZlcJe9gXUZ3F+xUjumT22YPrVh+tSG6ctGG+7n7n2aLiy4JJ5pZvZGc6PgSNuoHdOnNkyf2jB9asP05bINdTpdRESkQCmJi4iIFCgl8fg85ZI2tWP61IbpUxumT22Yvpy1YdFfExcRESlUOhIXEREpUEWdxM1skpktMrMlZnZD1PEUAjMbZGbPm9m7ZrbAzK6NL9/bzJ42s8Xx555Rx5rvzCxmZm+Z2ePxutqwDcysh5n9zczei/8+Hq02bBsz+2b873i+md1nZp3Vhq0zs7vMbK2ZzU9a1mK7mdmN8TyzyMxOy2QsRZvEzSwG3AZMBkYBF5jZqGijKgh1wLfd/UBgPPDVeLvdADzr7sOBZ+N1Se1a4N2kutqwbX4JPOHuI4FDCdpSbRiSmQ0Avg6MdfeDgRhwPmrDMP4ITGqyrNl2i38/ng8cFH/N7fH8kxFFm8SBo4Al7r7U3WuAqcCUiGPKe+6+2t3fjJe3EHxxDiBouz/FN/sTcE4kARYIMxsInAH8Pmmx2jAkM9sTOAH4A4C717j7RtSGbVUKdDGzUqArsAq1YavcfSbwaZPFLbXbFGCqu+9w9w+AJQT5JyOKOYkPAJYn1VfEl0lIZjYEGAO8BvRz99UQJHqgb4ShFYJfAN8BGpKWqQ3D2x/4BLg7fkni92a2B2rD0Nx9JfAz4CNgNbDJ3Z9CbdheLbVbVnNNMSdxa2aZuuqHZGbdgIeAb7j75qjjKSRmdiaw1t3nRB1LASsFDgd+4+5jgK3otG+bxK/ZTgGGAvsCe5jZRdFG1SFlNdcUcxJfAQxKqg8kOJUkrTCzMoIE/ld3fzi+eI2Z9Y+v7w+sjSq+AnAscLaZfUhwGWeimd2D2rAtVgAr3P21eP1vBEldbRjeycAH7v6Ju9cCDwPHoDZsr5baLau5ppiT+GxguJkNNbNygo4H0yKOKe+ZmRFch3zX3X+etGoa8OV4+cvAo7mOrVC4+43uPtDdhxD83j3n7hehNgzN3T8GlpvZiPiik4CFqA3b4iNgvJl1jf9dn0TQx0Vt2D4ttds04Hwz62RmQ4HhwOuZetOiHuzFzE4nuDYZA+5y9/+KNqL8Z2bHAS8C80hcz/0PguviDwCDCb4cPu/uTTt+SBNmVgH8u7ufaWa9UBuGZmaHEXQMLAeWApcRHJioDUMysx8CXyS46+Qt4AqgG2rDlMzsPqCCYLayNcD3gUdood3M7LvA5QTt/A13n5GxWIo5iYuIiBSyYj6dLiIiUtCUxEVERAqUkriIiEiBUhIXEREpUEriIiIiBUpJXKQImFlV1DGISOYpiYuIiBQoJXGRImJmFWZWmTQP91/jo3VhZkea2Stm9raZvW5m3ePzS99tZvPiE41MiG97qZk9YmaPmdkHZnaNmX0rvs0sM9s7vt0wM3vCzOaY2YtmNjLKzy/S0ZRGHYCI5NwYgrmNVwEvA8ea2evA/cAX3X12fKrPaoI5z3H30fEE/JSZfSa+n4Pj++pMML3i9e4+xsz+F7iEYDTEO4Gr3X2xmY0Dbgcm5uhzinR4SuIixed1d18BYGZzgSHAJmC1u88G2DkzXXyY3V/Fl71nZsuAnUn8+fic8lvMbBPwWHz5POCQ+Ex3xwAPxg/2ATpl96OJFBclcZHisyOpXE/wPWA0Pz1ic9MoNrefhqR6Q3yfJcBGdz+s3ZGKSEq6Ji4iAO8B+5rZkQDx6+GlwEzgwviyzxBM7rAozA7jR/MfmNnn4683Mzs0G8GLFCslcRHB3WsIZrP6lZm9DTxNcK37diBmZvMIrplf6u47Wt7Tbi4EvhLf5wJgSmYjFylumsVMRESkQOlIXEREpEApiYuIiBQoJXEREZECpSQuIiJSoJTERURECpSSuIiISIFSEhcRESlQSuIiIiIF6v8D8rjmyD+gGC4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.linspace(0, 100, 500)\n",
    "P_hat = 1/(1 + np.exp(-log_reg.params[0]-log_reg.params[1]*X))\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (8, 5))\n",
    "ax.scatter(df['Annual_Family_Income'], df['Own_House'])\n",
    "ax.plot(X, P_hat, color = 'tomato', lw = 3)\n",
    "ax.set_xlabel('Income')\n",
    "ax.set_ylabel('Owning a flat = 1, Not Owning = 0')\n",
    "ax.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45609362-ff0d-40ae-ac11-2071e6533616",
   "metadata": {},
   "source": [
    "To interpret estimated coefficients, we can calculate the marginal effect by taking derivative w.r.t certain $\\beta$, in this case $\\beta_2$.\n",
    "$$\n",
    "\\frac{dP_i}{d X_i} = \\frac{dP_i}{dZ_i}\\frac{dZ_i}{dX_i}=\\frac{e^{-Z_i}}{(1+e^{-Z_i})^2}\\beta_2\n",
    "$$\n",
    "Plot both accumulative and marginal effects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "edd757bf-9790-4fab-8bc2-28d2108f6148",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgEAAAE/CAYAAADffhAhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABfDUlEQVR4nO3deXwU9f348dd7d3MHEsJNAEFEEeVUQMUDLy4vlHpV61W11tbrW221tb9atVV7aT2q4lFrVcR6lSreinghoiIihyByhPvKfW8+vz8+k+xs2CSbZJPZTd7Px2Mf+XxmZmfeO7uZfe98PvMZMcaglFJKqc7H53UASimllPKGJgFKKaVUJ6VJgFJKKdVJaRKglFJKdVKaBCillFKdlCYBSimlVCelSUAbE5F1InJCC597lIisinVMsSQig0TEiEjA61jcROQ1EbkwBuuZJCJ5LXyuiMg/RWSPiCxypv1URLaJSLGIdG9tfC3ldRz1P9ut+T9R7c/53OzrdRyq9RImCRCR+c7BNMXrWNqK82W6X23dGPOBMeYAL2NqLefgXikiPepNX+K83kFtsV1jzDRjzL/aYt1uzmsocQ6KtY9fOrOPBE4E+htjxotIEvA3YLIxJtMYs6uF22xV4tVUHK71u1/TVy3ZVkNi9dkWkVtE5KlYxNQRiMhFIvK1iJSKyFYReVBEsmO9HedzszbW622Ohv4PROQJEbndq7gSTUIkAc4XxVGAAU71NhrVAt8D59ZWRGQEkNbSlcXbWQdglHNQrH38yZm+D7DOGFPi1HsDqcA3nkQZEm0c2a7XNKod4lKtICK/AO4CbgCygMOwn8G3RCTZy9hU/EqIJAC4AFgIPAGEneIVkQEi8qKI7BCRXSJyv2veZSKyQkSKRGS5iIx1pof94nZnjrWnf0XklyKyXUS2iMgMEZkuIt+KyG4R+XWk57qfH+lFiMh4EflERPKd9d5f+88pIgucxb5yfnmd7V6XiNwoIs/XW9/fReRep5wlIo85690kIreLiL+5cbj2zxUisto5+/KAiIgzzy8ifxGRnSKyFjgp8lsW5t/Y97DWhcCT9WI6SUS+FJFCEdkoIre45tVm/D8WkQ3Au04cf3Xi+F5Efu7+VSD2zNGlTvkiEfnQiXuPs/w01/ovdn1O1orIT6J4TY0SkR8DjwKHO+/nbKD29He+iLzrLDdMRN5yPlerROQs1zrSnNe4XkQKnNeQBixwradYRA6PsP0UEblHRDY7j3ucaftHiiPK1xTN5+ZK53NTJCK3icgQ5zmFIvKc6/Me8f9ERPqI/RXb3TXtELH/30lRxNjgZ9eZ39Ax4UDnM5MvIt+IyKmu5zwhIv8Q28RULCIfOXHe42xjpYiMcS3fT0RecGL+XkSubiTeLBF50ll2vYjcLCI+Z16jn9t66+kK/B64yhjzujGmyhizDjgLmwicLyKpIlImzlk5Z1vVznMRe8y4x/WaHxCRV5199amIDKm3n/eLctnJzme7wNmP74vzvxnhdTT6GWsJETnVeU/znff4wEivw/Vaar8LeojIK87zdovIB673Jur3OO4ZY+L+AawBrgQOAaqA3s50P/AVcDeQgf11c6Qz70xgEzAOEGA/YB9nngH2c63/CeB2pzwJqAb+H5AEXAbsAJ4BugAHAeXAvvWf63p+nqu+DjjBKR+Czc4DwCBgBXCta9n6cdWtC/uPXAp0db32LcBhTv1l4GFnP/QCFgE/aWB/RhPHK0A2MNB5/VOdeVcAK4EBQA7wnrN8oIFtrQNOwH7xHOjEvdF5PQYY5HqtI7CJ6UhgGzDDmTfIWfZJ5/WlOXEsB/oD3YC33XEA84FLnfJF2M/NZc72fwpsBsSZfxIwBPs5OcbZz2MjvZ8RXl/Ye1Zv3kXAh6567euojTHD2RcXO+/FWGAncJAz/wHndeQ6cR8BpNRfTwPbvhWbOPcCegIfA7dFiiPCcyPOJ7rPzVygK/b/pAJ4B9gX+8t0OXBhFP8n84CfuubdDdzXQKy3AE9F+dmNeEzA/p+vAX4NJAPHAUXAAa7/8Z3O608F3sWe3brAeV9uB95zlvUBn2OPH8nOa18LTGkg/ieB/2KPLYOAb4EfR/O5rbeeqdjj1l7vKfAvYLZTXgDMdMpvAt8B01zzTne95t3AeOf9fhp4NtLnvrFlgR5AIXCGM+8a5zVd2pJjU5Sf0ycIHc/3B0qwzXJJwC+d9zo5iu+CO4CHnOclYc9GS3Pf43h/eB5AkwHadtUqoIdTXwlc55QPx/6TR/rgvwFc08A6m0oCygC/U+/iLD/BtfznhL6g6p7ren7Eg1uEOK4FXmokrvrr+hC4wCmfCHznlHtjD7hprmXPxTkwRbGPI8VxpKv+HHCjU34XuMI1b3Kkf8T6rx+42fmnmgq8hf0nr0sCIjzvHuBupzzIWXZf1/x3cSU5zjYaSwLWuJZNd5bt08C2X6797NR/Dxr4LBUC+a7HFNd2G0sCzgY+qLe+h4HfYQ80ZdimhvrbDFtPA3F9B0x31adgmyaafL5rvvs1XR/l52Zivf+TX7nqfwXuaer/xNkvHzllP7AVGN9ArLewdxLQ0Gc34jEBe3DfCvhc02YDt5jQ//gjrnlXAStc9RFAvlOeAGyot/6bgH9G2K4f+3873DXtJ8D85n5ugfOBrQ3sozuBt5zybcC92P+/rdgv5TuxyU0ZoePsE8CjrnVMB1bW28/7NbUsNlH6xDVPsIlvxCSgqc9YFJ/TfKCS0PH8t8Bzruf4sIngpPqvw/Vaap97KzZB26/edqN+jxPhEW9tq5FcCLxpjNnp1J9xpt2N/TW63hhTHeF5A7AHwpbYZYwJOuUy5+821/wyILO5K3VOxf4NOBT7Dx3AHiij9Qz2y/1J4IdOHUK/Zra4znz6sP9sLY1jq6tcSuj19qu33vVRxv5v7C+NwdRrCnBimoA9GB2Mza5TgP/UW8y93fpxRHytLnWvxxhT6uynTGfb07BfvPtj91s68HVTL8hlrDFmTTOWr7UPMEFE8l3TAth91QN7YG7pZ7gf4e/Nemdac/Rw/29F+bmp/39Sv94niu3+F3hIbO/z/YECY8yiZsTd0Ge3oWNCP2CjMabGNW099gxMraZeV+029gH61XtP/cAHEbbbA/tZr/8+ubfb4Oe2np1ADxEJRDge9nXmA7yPfQ/HYj/jbwGPYX99r3EdZ8O2Tfh+jCSq44UxxkRqBqrVwmNk/c/pE655Yf8HxpgaEdlI+D5uyJ+xSeabzn6fZYy5k+a9x3EvrvsEiG3/PAs4RmxP163AdcAoERmF/XANlMgdxTZiT/FGUor9gNWK5sDUkJJmrOtB7JmMocaYrtjTj9LI8vX9B5gkIv2B0wklARuxvyh6GGOynUdXY8xBbRDHFuzBtNbAaJ5kjFmPPYU6HXgxwiLPYE8lDzDGZGFPw9WPydSLo7+rPoAWEHu1yQvAX7DNTNnY09HNeV9aaiPwvus9yza2E95PsQftciJ/hk2EafVtxh6sag10prVGaz+/UTHGlGN/wZ8H/AibFMVCQ8eEzcCA2vZex0DsL8aWbOP7eu9pF2PM9AjL7sSe5az/PrVku59gjwFnuCeKSAYwDdssA7ZZ6ADs8eN9Y8xyZ5snYROEWAv7P3X6Z/RvePGYf8bC/g+c7Q8gtI8b/C4wxhQZY35hjNkXOAX4PxE5nua9x3EvrpMAYAYQBIYDo53HgdiM6wJsu/cW4E4RyXA6vkx0nvsocL3YTkUiIvuJSO2HYQnwQ7Gdy6Zi24FbagkwXURyRKQP9vRVQ7pgTx0Xi8gwbBuf2zZs+1JExpgd2NPc/8R+CFc407dg2/f+KiJdRcQntkNWQ6+rqTga8xxwtYj0F5FuwI3NeO6PgeNMqLd8/Zh2G2PKRWQ89kxHU3FcIyK5Yi+B+lUz4nCrPeuwA6h2zgpMbuG6musVYH8R+ZGIJDmPcSJyoPOr9HHgb04nJL+IHO4kLTuAGhr5rGBPZ98sIj2djmD/D2jtpXSt+dw015PY0+Gn0vq4azV0TPgUm8z/0nkPJmEP+s+2YBuLgEIR+ZXYjp1+ETlYRMbVX9A52/gc8AcR6eLE8n+04PUaYwqwHQPvE5GpzusYhP3hkIeTSBljSrG/rH9G6Ev/Y2wzRFskAa8CI8R2rg44223sh1KsP2PPASeJyPFiO5b+ApssfezMX0ID3wUicrLzGREnpqDziPo9TgTxngRciG1n2WCM2Vr7AO7H/koQ7D/rfsAG7If9bABjzH+AP2B/YRZh23lznPVe4zwv31nPy62I8d/YzonrsF/EcxpZ9nrsl1sR8EiEZW8B/uX0Rj2LyJ7Btn8/U2/6BdgvtOXAHuB57GnAlsTRmEewbatfAV8Q+Vd9RMaY74wxixuYfSVwq4gUYb+wnosijjeBpcCX2F/v1dh/0qgZY4qAq53t7cHul7nNWQehKzpqH/c0Y9uTgXOwv1i2Yi/xqh0L43rsKdvPsB2v7sK2W5diP9sfOZ+VwyKs/nZgMXb/fI19r1p77XRrPjfNYoz5CJvofGFsL/dYrDPiMcEYU4lNNqZhf53/A9v3ZmULthHEHltGY8987cQmH1kNPOUqbAKyFtvn5xls8tdsxl6a+mvsWa1CbHKzETjeGFPhWvR9bPPhIle9C6GrTmLGaV44E/gTsAv7g24x9os4kph+xowxq7D9Je7DvhenAKc47zk0/l0wFNvhuBh7puUfxpj5LXiP41pt72ilEprzC/4hY8w+TS6sEoLYyxefMcY86nUsKjacJpc84DxjzHtex6Pi/0yAUhE5p+Gmi0hARHKxHfte8jouFRvOqdWxtOHZBtU+RGSKiGQ7TVm1bfwLPQ5LOTQJUIlKsG2ge7DNASuwzQgqwYnIv7CnYa91mkxUYjsce1VG7en4GcaYssafotqLNgcopZRSnZSeCVBKKaU6KU0ClFJKqU7KsxEDfT6fSUtr8Y3klFJKqYRTWlpqjDFx8wPcsyQgLS2NkpJIY8YopZRSHZOIxFWnyLjJRpRSSinVvjQJUEoppTopTQKUUkqpTkqTAKWUUqqT0iRAKaWU6qQ0CVBKKaU6KU0ClFJKqU6qyXECRORx4GRguzHm4AjzBfg7MB0oBS4yxnwR60CVao2Xv9zEn99Yxeb8Mvplp3HDlAOYMSY3rmNqbsz1lz92WE/eW7kj6voNUw4AaHSbN7/8NbM/3UjQGPwinDthALfPGNHimJu7T2Itmu019ZpjsY1Yx9yc50f6HHj9v1FfW38u4vH40F6avIGQiBwNFANPNpAETAeuwiYBE4C/G2MmNLXhjIwMo4MFqfbw8pebuOnFrymrCtZNS0vyc8cZI2L/j26CUL0FMBDoB+JvUUzNjTnS8s2V5BcwUFUTOia4t3nzy1/z1MINez3v/MMGcvuMEa3ez+36PkW5vaZecyy2EeuYm/v8+hpcnzFgaqDG+Vs7zf03rGxc8yGsYEzkadRb3hheW7aFP766gvJqG7MYQ2qSn19PG8bUg/s0Y9v12RlvfrONP7+xivLqIOLMSQ34uH7KAUwe3hsCSZDTs6GVNJuIlBpjMmK2wlaK6i6CIjIIeKWBJOBhYL4xZrZTXwVMMsZsaWydmgSo9jLxznfZlL/3IF252Wl8dONxsdlI5Xew649Q9CLU5Ntpvi6QeRp0vxFSDmpWTM2NuaHlY6F2m0NumkcwwvHCL8J3d0xv9X5ul/epmdtr6jU35eg73qJgTyGZppJUU00KQVJNNbkZfu77wUFQVbn3o7ICqqogWAXBoPOohmCQ15fmUVVRhZ8akqghYGoIUEOGH8b27xK2LDU19svaGKds2JpfSk1NDT4MPmMQDD4MfuevYPAZ8Ish1S+u59a0focnqgFD4HcPxGx18ZYExGLY4Fxgo6ue50zbKwkQkcuBywGSk5NjsGmlmra5gS/HhqY3izGQ/zBs/wWY0vB5NUVQ+BQUzoaed0DO9SASVUzNjTkmr6UBteuO9GXont7a2Nr7tUWzPVMTJNtU0M2U062mnGxTTk5NOdk1ZfDiNigpgrISKC+DslL7t7zUeZSxoKoy8sYLgHubH/PUhmZUAWuafn6faDdkgOpoF1aJLBZJgESYFvFoYYyZBcwCeyYgBttWqkn9stMi/uLrl93KG1gZA9uvhT31jua+HNsMENzhTAjCjl9C5QroMwsk0GRMzY25oeVjoXabfpEGfxU3FkO0+7nN3qcGDOqahG/3NvrXFNG7ppjeNSX0rilhH3853P467NnF6oI9+Bs6nzyvTcKKb+IDnwBSl9Du9Tes7F4OwgqRlhdxzbd/dpdWEawxgIRaCUTwC/ToktrEtutP2+sFsSm/zFm/s27nb8DnI7dbGnTr0cBzO4ZYJAF5wABXvT+wOQbrVSombphyQMR21NqOcC2248bwBCD5IOjzEKRNtPXyxbDtaihfaOsF/wRfJvS+t8mYmhtzpOWbq6E+AbXbPHfCgIjt4+dOGNCimKN5Da1+n6qrYNsm2LQOtm+G7Vtg51bYsYX39uxs+HlOS2XkHh3RMyIUkUSJJFFGEhXip9IXILdnNj26ZUJyCiQlQVIKJCW7HkkQSAa/D/wB5+Hn87wiZn++idIgVOOjSnz4AklcfNQQjjygD/j9dcvi84HPb78AnS/vN5dv5843VlFSZagRwSDUINSA89dHSpKfW047mFPG9Hc919fIF2nbWtDGfUU+a2T9uZ2gc2AskoC5wM9F5Flsx8CCpvoDKNWeag8UMe39W/AU7P5TqN7lTOj7JPhSQ9PSxsHA+bD1J1D4Lzttz32QNJQZY65qNKbmxhxp+VhfHVDbEa6hnvKt3c+tfp8K8+H7VbBxrf3S37QOtuXZ9vFWKAuksr0mmT2+NPJ9qWT37M7oAwdCRhfI7ArpmZCSBqlpkJoOael1ZUlO4d0lm/d6TaNa+Nk7BNg4bO+e7EdGub7Jx/SltGuPJq8OOCWOvvza5P+3Hdcf76K5OmA2MAnoAWwDfgckARhjHnIuEbwf21xVClxsjFnc1Ia1Y6BKWJXfwvdjwTif38xTIfd5kKTIy5sa2HwOFP3HmZAEgxZD6sh2CbdDCgZh/WpY8w2sXQXrVsHObc1bh/ige0/o0cee8s3qDt26Q3Z3yM6xf7NybO9wpWIk3joGRnV1QFvQJEAlJFMD64+A8k9tPXl/2Gcx+Ls0/ryaMthwtG0iAEgZDYMWNZw4qHA1NfaX/YolsPIr+PZr2/kuGj16Q79B0Kc/9OwLvfpCj77QvRcEYnEyVKnoxVsSoP8BSjVHweOhBECSod+cphMAAF8a9P03rBsNpgIqlsDuv0D3m9oy2sRWVQkrvoQlC+GrhVCwp/Hlk5Jh4H4waCjkDoL+g6HfQHuKXikVkZ4JUCpawd2wdn8I7rL17r+Fnrc2bx27/go7rrdlXybsuwYCvWMbZyKrroJli2Hhe/D1Iqgob3jZbj1g2CjY90AYfID90tdf9irOxduZAE0ClIrW9uth919tOWkQDP4GfM38lWmq4ftRULnc1rOvgD4PxjTMhGOM7dD38dvw2fv22vtIMrPgwNH2i3/YaHta36Me60q1VFNJgIhMxQ7F7wceNcbcWW9+o0P1i4gfWAxsMsac3FQ8mjYrFY3qLbDHNWpYr782PwEAkAD0+jPknWTr+bMg5zrbt6CzqayARe/De/+znfwi6dMfRh0Gow+HIcPsJW9KdVDOF/gDwInYy+8/E5G5xpjlrsWmAUOdxwTgQedvrWuAFUDXaLapSYBS0dh1Fxjn1HTKWMg8veXrypgG6cdD6TtADey6E/o+HpMwE0LBHnjrJfjgtci/+nN6wWHHwoRjbdu+Up3HeGCNMWYtgHPp/WmAOwk4DXsvHwMsFJFsEelrjNkiIv2Bk4A/AP8XzQY1CVCqKdXbIP+hUL3nra07DS0CPW6BDe/YesG/ocfvIGmfVoUZ9/bshNf/Awtes53+3AJJMP4YmDgZhh5sB6dRqvOJNAx//RvyNTZU/z3AL4EoeitbniUBOTk5zJ8/36vNKxW1QZn/ZFBmBQCFlcP4YnE6ML/V6x2dM5Ls5KVANZuWXsvqomtavc64FAxCwW4oLgATgMNOCc0LJEGXLNve7/PBlt2wZYF3sSrV9gIi4h5LZ5YzpD5ENwx/xGVE5GRguzHmcxGZFG0w2jFQqcbUlMF3AyHoDDHb7znoemZs1l3yJmycYsuSDvttAn92bNYdDyrK7Wn/156Dinr3BBi4H5x8Low+TNv5VafSWMdAETkcuMUYM8Wp3wRgjLnDtUzEO/cCVwM/wt76KRXbJ+BFY8z5jcWjzQFKNabw36EEILAPdGlFX4D60k+ElJFQsdTegbDgn7aTYKIzBr74CJ59yDYBuO07DE7+IYwYpz37ldrbZ8BQERkMbALOAX5Yb5mGhuq/yXngnAm4vqkEADQJUKphxsCe+0P1nGts7/5YEYFuP4etl9v6ngeg2zV2ONtEtWs7PPMAfPVp+PR+A+HMy+DgQ/XLX6kGGGOqReTnwBvYSwQfN8Z8IyJXOPMfwt6/cjr25tGlwMWt2aY2ByjVkLKFsP5wW5Z02G8z+LNiu42aEljTH2rybb3/q5A5PbbbaA/G2Ev9Xng8fICfLtkw40dw5FR7ZzulOrl4GyxIzwQo1ZD8WaFy13NinwAA+DIg6xLY8zdbL3g88ZKAgt3wz7/Zkf7cjpkOMy+xd9lTSsUlPROgVCTBAljTF4zToW2fhZBW/0qdGKlYDt8fZMuS7Jxx6N4224q1JQvhibttz/9a/QbCj66BoQd5F5dScUrPBCiVCIrmhBKAlJGQOr7ttpUy3K6/fBGYSih8Frr9rO22Fws1QXj5SZg3J3z65Jlw+oX2Zj5KqbinSYBSkRQ8FSpnXdL2ndmyLrJJAEDBE/GdBBQXwqw7YfkXoWndesAlv4ADx3gXl1Kq2bQ5QKn6qtbDd4Ocit9ev9/Wd/oL7oE1feyZAIDBX0PKwW27zZbYtA7u+x3s3BaadtAhcNmvIDOqocqV6tTirTkgga9FUqqNFM4OlTNOaJ9b/fq7QeaMUL3gX22/zeZatRTuuj48ATj5XLjmVk0AlEpQmgQo5WZMeFNA1/Pab9tZF4XKBU+BCbbftpuyaD7c/RsoLbb1lDT42f+DGRfqiH9KJTBNApRyq1gKld/YsqSF/zpvaxkngr+XLQe3QtlH7bftxrw71/YBqK6y9axu8Ms/wZgjvI1LKdVqmgQo5Vb4dKiceRr4o74ZV+tJALqc4YrlP+237Ya8+SI8849Qve9AuOlu2GeodzEppWJGkwClahkDhc+F6lnt2BRQq4vr5kTFL4Cpaf8Yar32H3jONWDSkAPhV3+FHn28i0kpFVOaBChVq+JLqF5vy74syJjc/jGkHw3+nrZcvcW7JoHXnoMXHgvVhx4M1/0BMtvxzIhSqs1pEqBUraIXQ+XMU+3ofe2tfpNAkQdNAu+/au8BUOuAkXDt7ZCa3v6xKKXalCYBStVyJwHuL+L25m4SKHq+fZsEFn8AT7nunDhsFFx9K6Sktl8MSql2o0mAUgAVK6ByhS1LmjdNAbXSjwF/D1tuzyaB5V/AI3fZvhEAg/aHn/9OEwClOjBNApQCKHopVM6YBj4PT33XbxIontv228z7Hh64DYLVtt6nP1xzmzYBKNXBaRKgFEBxnDQF1Mo8LVRu6ySgMN8OBVzh3DCpWw+47o/QpQ1unayUiiuaBChVtR7KP3cqSZB5kqfhAJB+HIjzK7zyW6hY1TbbqaqEf9wKu7bbekqa7QTYvVfbbE8pFVc0CVCq+JVQOeM48Gd7FkodXypkTAnV2+JsgDHw7/tgzXJbF4HLb4TcQbHfllIqLmkSoFTxq6Fy5snexVFf5imhcvH/Yr/+9/4HH78Vqv/gxzBqQuy3o5SKW5oEqM6tphRK3wvVM+KgKaBW5kmA2HLZR1C9M3br/n4VzHGNBnjECTB5ZuzWr5RKCJoEqM6t9F0w5bacfCAkD/Y2HrdAL0g73KnUQMm82Ky3pAge/mPoSoCBQ+BHV9vmAKVUp6JJgOrcil1frPHQIbC+zFND5Vj0CzAG/vk32LnN1tPS4YrfQJIHoyMqpTynSYDqvIwJ7w8QT00Btdz9AkreAlPVuvW9819Y8kmoftH/Qa9+rVunUiphaRKgOq/Kb6B6gy37siB9orfxRJJ8IAQG2nJNIZQtbPm6Nq+H5103BTr+NDjkyNbFp5RKaJoEqM4r7CzAZJAk72JpiAhkTg3VS15v2Xqqq+DRP9u/AAOGwJmXtj4+pVRC0yRAdV5hScB07+JoSoYrCShuYRLwymzYsMaWA0lw6Q32r1KqU9MkQHVOwT1Q9nGonjnNu1iakn48ELDlii+gemvznr92Jcx7NlQ//SIdEEgpBWgSoDqrkjeBoC2njoNAb0/DaZS/K6S5+iuUvBn9c6sq4fG/QI1zO+L9R8CJp8c2PqVUwtIkQHVO7i/SeG4KqNXSfgHz5sDWPFtOTYdLfgE+/bdXSllRHQ1EZKqIrBKRNSJyY4T5WSLyPxH5SkS+EZGLYx+qUjFiDJS8E6pnnOhdLNFy9wsoeRNMsOnnbF5vk4BaMy+GHn1iH5tSKmE1mQSIiB94AJgGDAfOFZHh9Rb7GbDcGDMKmAT8VUR09BEVn6q+h+r1tiwZkDbe23iikTIK/M4XeHCX666HDaipgSf/HhoVcMhwOCYOx0FQSnkqmjMB44E1xpi1xphK4FngtHrLGKCLiAiQCewGqmMaqVKxUuo6C5B+dHxeGlifSPhdBZtqEvjgtdDdAf0BuOBqbQZQSu0lmqNCLrDRVc9zprndDxwIbAa+Bq4xxtTEJEKlYq3k3VA5/Xjv4miuaPsFFOyG5x8P1aedqVcDKKUiiiYJiHRXEVOvPgVYAvQDRgP3i0jXvVYkcrmILBaRxdXVeqJAecAYe9OgWhnHeRdLc2WcSN2/bNmn9jLHSF58AspKbLl3Lpx0bntEp5RKQNEkAXnAAFe9P/YXv9vFwIvGWgN8DwyrvyJjzCxjzKHGmEMDgUBLY1aq5SqWQXC7Lfu727b2ROHvDqmHOpUaKJ2/9zJrV8FHrisffnil3hxIKdWgaJKAz4ChIjLY6ex3DlD/dmYbgOMBRKQ3cACwNpaBKhUT7rMA6ceCJFg7eYar+aLk7fB5NTXw7IOh+ujD4aBD2icupVRCavIIaIypBn4OvAGsAJ4zxnwjIleIyBXOYrcBR4jI18A7wK+MMTvbKmilWiysU2ACNQXUSj8hVC6tlwR88o4dHRDskMBnXdZ+cSmlEpIYU795v31kZGSYkpIST7atOilTDau727vxAey7CpL39zam5qoph9XdwJTb+pANkDTA9gG4+VIocPoJTD8HzrjIszCVUpGJSKkxJsPrOGol2LlQpVqh/PNQAhDIhaSh3sbTEr5USHPd/rd20KPX/hNKALK7w/Sz2z82pVTC0SRAdR6l9S4NlEgXviSAjHpNAnt2wtsvhabNvARS09o/LqVUwtEkQHUeYUMFJ2B/gFruJKDkHZj7b6issPWB+8GEY72JSymVcDQJUJ1DTTmUfRSqJ9IgQfWljAZfN1sOboUVL4bm/eASHRlQKRU1PVqozqHsk1BnuuT9Iam/t/G0hvjDz2QMdC7EOWgsDB/rTUxKqYSkSYDqHBL90sD63JcK7rPL9m+Y+WPv4lFKJSRNAlTnUL9TYKJzJzL9d8NhR8PAId7Fo5RKSJoEqI4vWAhli0L19EmehRIza4qgMNWWU4Iw/WBv41FKJSRNAlTHV/YBELTllNEQ6OFlNK1nDMx9CtZ3D01L+tK7eJRSCUuTANXxlXSw/gDffA7frQhPAtyvUSmloqRJgOr43J0CMxK8P4Ax8PKTtrwhJzS97BOoKfYmJqVUwtIkQHVs1TugYqlTCUDaUZ6G02pLF8G6b225MhMCw50ZVVD6oWdhKaUSkyYBqmMrfS9UThsP/i7exdJatX0Bak06GbqcGKq7r4BQSiUkEZkqIqtEZI2I3BhhvojIvc78pSIy1pmeKiKLROQrEflGRH4fzfY0CVAdW0e6NPCbz2H9altOToFpZ4a/Ju0XoFRCExE/8AAwDRgOnCsiw+stNg0Y6jwuBx50plcAxxljRgGjgakiclhT29QkQHVsHalT4Lw5ofJRUyErB9KPpu7fuOJLCO72JDSlVEyMB9YYY9YaYyqBZ4HT6i1zGvCksRYC2SLS16nXdgxKch6mqQ1qEqA6rqoNULXGliUN0g73Np7WWLMcvv3alv1+mDzTKWdB6jhnIQOl872ITikVG7nARlc9z5kW1TIi4heRJcB24C1jzKdNbTDQmmhbIycnh/nz53u1edUJ9El7nWFZtry7/CCWLvjE24BaY/tmOMr54s/sCl8vB5YDMDhzCPtk2v/1Td/+m9VFOQ2sRCkVBwIisthVn2WMmeWUI93fvP6v+QaXMcYEgdEikg28JCIHG2OWNRpMdDHH3u7du5k0aZJXm1edwebHoNAWc/rPZNKoSZ6G02Ib18JTd9qyCNw6C/oOCM0vqYaNzwCQ220luYdMav8YlVLRqjbGHNrAvDzA9c9Nf2Bzc5cxxuSLyHxgKtBoEqDNAapjMqbjdAp87blQeezE8AQAIG0iSLItV66EqvrHDKVUgvgMGCoig0UkGTgHmFtvmbnABc5VAocBBcaYLSLS0zkDgIikAScAK5vaoCYBqmOqXAXVzpehLwtSE/QWu9s2w2cLQvVpZ++9jC8N0o4I1d2XRSqlEoYxphr4OfAGsAJ4zhjzjYhcISJXOIvNA9YCa4BHgCud6X2B90RkKTaZeMsY80pT2/SsOUCpNhV2FmASiN+zUFrljf+AqbHl4WNh0NDIy6UfH+oUWPoOZJ3XLuEppWLLGDMP+0XvnvaQq2yAn0V43lJgTHO3p2cCVMfUES4NLNgDH78dqk+PcBaglvs1lrxjm0OUUqoJmgSojsfUhJ8ST9T7Bcx/BaqrbHnwAXDAyIaXTRsHvkxbrt4AVWvbPj6lVMLTJEB1PBVLoGaPLft7Q3L9AbcSQGUFvOdqzps8014Z0BBJgrSjQ3UdQlgpFQVNAlTH424KyDiu8S/PeLXwXSgusOWcXvaqgKZk6BDCSqnm0SRAdTyJfmlgTQ289VKofsJpdpTAprj7BZS+q/0ClFJN0iRAdSymEkpdl9QlYqfAbz6HLRtsOTUdjpwa3fNSRoK/uy0Hd0BFo2OEKKWUJgGqgylbBKbUlpMGQ/Jgb+NpCfdZgKOmQHpGdM8TH6QfG6prvwClVBM0CVAdS2mCXxqY9z0s/8KWxQfH17+BWBPCmgS0X4BSqnGaBKiOpcT16zcRLw10nwUYewT06NO857tfc+n7YKpjE5dSqkPSJEB1HDUlUOa6U2CinQkozIdPXeMbTD6j+etIGgoB586jNYVQ/nlMQlNKdUyaBKiOo/RDwBlcJ/kgCPT2NJxmW/BaaHCgfYfBkBaMbyASfkWE9gtQSjVCkwDVcZQmcFNAMAjvvxqqH3dqy9eV4R5CWJMApVTDNAlQHUcidwpc8gns2WnLXbLgkCNbvi73ay/7EGrKWxebUqrD0iRAdQzBPVDu9KrHB+nHeBpOs7mHCD56GiQlt3xdSQNs3wAAUw5lC1sXm1Kqw9IkQHUMpfMBZ4S81EPAn+1hMM20eT2sXGLL4oNjTmr9OjP0UkGlVNM0CVAdQ0kCDxX83v9C5TGHQ07P1q9TOwcqpaKgSYDqGNy/dhOpU2BZCXzsir01HQLd0ie5trEIgkWxWa9SqkPRJEAlvqrNULnCliUZ0o7wNp7m+OQdqCiz5X4D4YCRsVlvoCekjHIq1VD2QWzWq5TqUDQJUImv1DXATtoR4Ev3LpbmMAbedTUFHHtKbG977L5KQG8trJSKIKokQESmisgqEVkjIjc2sMwkEVkiIt+IyPuxDVOpRiTqpYErl8DWjbacmg6Hx7gZI0P7BSilGhdoagER8QMPACcCecBnIjLXGLPctUw28A9gqjFmg4j0aqN4lQpnTPiv3ETqFOg+C3DECTYRiKW0owA/EISKJVC9EwI9YrsNpVRCi+ZMwHhgjTFmrTGmEngWqH9rsx8CLxpjNgAYY7bHNkylGlC1Fqo32LIvE9LGeRtPtPbshCWu6/ePPSX22/B3hbTxoXrp/NhvQymV0KJJAnKBja56njPNbX+gm4jMF5HPReSCWAWoVKPcp7nTjgZJ8i6W5vjoTTA1tnzASOg7oG22E3ZrYW0SUEqFiyYJiNRTydSrB4BDgJOAKcBvRWT/vVYkcrmILBaRxdXVeotTFQMlCXhpYE0NfPB6qH70tLbbVth4Ado5UCkVLpokIA9w/0zpD2yOsMzrxpgSY8xOYAEwqt4yGGNmGWMONcYcGgg02R1BqcaZmvBft4nSKXD5F7DLaTHL7ApjJ7bdttIOB0mx5cpvoSqv7ballEo40SQBnwFDRWSwiCQD5wBz6y3zX+AoEQmISDowAVgR21CVqqfiGwjusGV/d0iJ0TX2bW3Ba6Hy4ce37j4BTfGlQporydAmAaWUS5NJgDGmGvg58Ab2i/05Y8w3InKFiFzhLLMCeB1YCiwCHjXGLGu7sJUCSt8OldOPs+Pux7uC3fCVq0PgUW3YFFDL3SSgtxZWSrlEdU7eGDMPmFdv2kP16n8G/hy70JRqQokrCcg4wbs4muOjtyAYtOWhB9lRAttaxnHg3KWY0nfsZZWxHJRIKZWwEuCnk1IRmEoodY1JlZ4ASUD9DoHtcRYAIPVQ8HWx5eo8qFrTPttVSsU9TQJUYir7FEyJLScNhuR9vY0nGqu+gh1bbDk9Ew49qn22KwFIOyZU1yGElVIOTQJUYnI3BSTCWQCABa6zAIcdB8kp7bdtHUJYKRWBJgEqMYXdOjgBkoCifPjio1C9LccGiKT+oEG1AxUppTo1TQJU4gkWQpmrh30ijA/w8TsQdAbI2ncY9B/cvttPORj8PW05uAsqvm7f7Sul4pImASrxlC0AnB72KWPi/6Y4xsAHrrEB2vssANjLJ9OPDdW1SUAphSYBKhEl2qWBq5fBVmekvtR0GHdM48u3FfcZE+0cqJRCkwCViMI6BSbA/QLedw2xMeFYSEn1Jg5358Cy98FUeROHUipuaBKgEkv1Fqj8xpYlGdKP9DaephQXwecfhupHT/UulqQhEHBuA1JTDOWLvYtFKRUXNAlQicU97G3aEeDL8C6WaCx8B6qdX9z7DLUPr4iEnw3QIYSV6vQ0CVCJpTSBxgcwJvxmQV6eBagVdqmg9gtQqrOLq/v5VlVVkZeXR3l5udehqDiTmppK/9xckhKpU+DaFbB5vS0np8D4SZ6GA4QnAWUfQ00Z+NK8i0c1So+JiSs1NZX+/fuTlJTkdSiNiqskIC8vjy5dujBo0CBEb3CiHMYYdu3axbaNH9C/2ull78uC1EO8Dawp77vOAoyfBGlx0HSRlAvJB0DlKjAVNhHISIDOlZ2UHhMTU+0xKy8vj8GD23lMkGaKq+aA8vJyunfvrh92FUZE6N69O4HK+aGJ6cfaMfHjVWkJLF4QqnsxNkBD0nUI4UShx8TEVHvMSoQzOHGVBAD6YVcRiQgZfBKaEO9NAZ++C5UVttx/MAw+wNt43HS8gISix8TElCjvW9wlAfHgpZdeQkRYuXKl16HsZdCgQezcubPRZf74xz+G1Y844oiYbNvv9zN69Oi6x5133gnABx98wEEHHcTo0aMpKyvjhhtu4KCDDuKGG25o9jbqx17HVJPOolA9njsF7tUhcJrtmR8vMiYBTjzln0GwwMtoVJwTEX70ox/V1aurq+nZsycnn3xyzLYxffp08vPzW/TcW265hb/85S8Rp+fm5oYds2q3ce655zJy5EjuvvtuVq5cyejRoxkzZgzfffdds7Y9f/58Pv744xbFHS80CYhg9uzZHHnkkTz77LNeh9Ii9b9IY/UhTUtLY8mSJXWPG2+8EYCnn36a66+/niVLlpCWlsbDDz/MF198wZ///OdWx16nbBF+KbLlQH9I3r+lL6PtrfsWNq615aRke8fAeOLvDimjnUoNlH3gZTQqzmVkZLBs2TLKysoAeOutt8jNzW3WOqqrqxudP2/ePLKzs1saYoOuu+66sGNWdnY2W7du5eOPP2bp0qVcd911vPzyy5x22ml8+eWXDBkypFnr1ySgAyouLuajjz7iscceC0sCgsEg119/PSNGjGDkyJHcd999AHz22WccccQRjBo1ivHjx1NUVMQTTzzBz3/+87rnnnzyycyfPx+AzMxMfvWrX3HIIYdwwgknsGjRIiZNmsS+++7L3LlzARp9vtuMGTM45JBDOOigg5g1axYAN954I2VlZYwePZrzzjuvbpsAZ599NvPmhUavu+iii3jhhRcIBoPccMMNjBs3jpEjR/Lwww9Hvb8effRRnnvuOW699VbOO+88Tj31VEpKSpgwYQJz5sxhx44dzJw5k3HjxjFu3Dg++uijuv188cUX1+3PF154IWLsdUpct+HNmBpfv6zr+8AV66FHQ3qmd7E0JEObBFT0pk2bxquvvgrYH0nnnntu3bxFixZxxBFHMGbMGI444ghWrVoF2OPYmWeeySmnnMLkyZMpLS3lrLPOYuTIkZx99tlMmDCBxYvtgFW1ZzjXrVvHgQceyGWXXcZBBx3E5MmT65KPRx55hHHjxjFq1ChmzpxJaWlpi17L5MmT2b59O6NHj+b3v/8999xzD48++ijHHmvvrfHUU08xfvx4Ro8ezU9+8hOCQXufktdff52xY8cyatQojj/+eNatW8dDDz3E3XffzejRo/ngg8RMpuO3Z9WlbXhN9aOvNzjr5ZdfZurUqey///7k5OTwxRdfMHbsWGbNmsX333/Pl19+SSAQYPfu3VRWVnL22WczZ84cxo0bR2FhIWlpjV9uVVJSwqRJk7jrrrs4/fTTufnmm3nrrbdYvnw5F154IaeeemrUL+Pxxx8nJyeHsrIyxo0bx8yZM7nzzju5//77WbJkyV7Ln3POOcyZM4fp06dTWVnJO++8w4MPPshjjz1GVlYWn332GRUVFUycOJHJkyfv1au19gu61k033cSll17Khx9+yMknn8wPfvADwCYdtdv/4Q9/yHXXXceRRx7Jhg0bmDJlCitWrOC2224jKyuLr7+2d7Pbs2cPM2fObDD2vZKAeFVeCp++F6rHU4dAt/TjYfdfbbnkTW9jUdFZ2YaJ7zDT6OxzzjmHW2+9lZNPPpmlS5dyySWX1H3pDRs2jAULFhAIBHj77bf59a9/zQsvvADAJ598wtKlS8nJyeEvf/kL3bp1Y+nSpSxbtizsWOK2evVqZs+ezSOPPMJZZ53FCy+8wPnnn88ZZ5zBZZddBsDNN9/MY489xlVXXdVo3HfffTdPPfUUAN26deO9995j7ty5nHzyyXXHGWMMmZmZXH/99axYsYI5c+bw0UcfkZSUxJVXXsnTTz/NtGnTuOyyy1iwYAGDBw9m9+7d5OTkcMUVV9Q9N1HFbxLgkdmzZ3PttdcC9oM/e/Zsxo4dy9tvv80VV1xBIGB3WU5ODl9//TV9+/Zl3LhxAHTt2rXJ9ScnJzN1qv0SGzFiBCkpKSQlJTFixAjWrVvXrFjvvfdeXnrpJQA2btzI6tWr6d69e4PLT5s2jauvvpqKigpef/11jj76aNLS0njzzTdZunQpzz//PAAFBQWsXr16rySgtjmgOd5++22WL19eVy8sLKSoqIi333477ExLt27dGl5J9Q7XELf++L6kbdH7UOH0CO47EPYb7m08DUk/2g67bCqhcjlUbYCkgV5HpeLUyJEjWbduHbNnz2b69Olh8woKCrjwwgtZvXo1IkJVVeieFCeeeCI5OTkAfPjhh1xzzTUAHHzwwYwcOTLitgYPHlyXIBxyyCF1x8Vly5Zx8803k5+fT3FxMVOmTGky7uuuu65ZX9DvvPMOn3/+ed0xvaysjF69erFw4UKOPvroumNi7WvqCDQJcNm1axfvvvsuy5YtQ0QIBoOICH/6058wxuzV2zPSNIBAIEBNTU1d3X2ZSFJSUt1zfD4fKSkpdeXadrPGnl9r/vz5vP3223zyySekp6czadKkJi9HSU1NZdKkSbzxxhvMmTOn7pSeMYb77rsvqn+q5qqpqeGTTz7Z6wxJQ/suopK3AOeXStrh4M+OaYwxtcB1s6B46xDo5suAtGOg9C1bL3kdsi/3NiYV10499VSuv/565s+fz65du+qm//a3v+XYY4/lpZdeYt26dUyaNKluXkZGaGwMYxo/21Cr9pgItjNybXPARRddxMsvv8yoUaN44oknIjaRtpYxhgsvvJA77rgjbPrcuXMTprd/c8VvEtDIKfu28vzzz3PBBReEtYkfc8wxfPjhh0yePJmHHnqISZMm1TUHDBs2jM2bN/PZZ58xbtw4ioqKSEtLY9CgQfzjH/+gpqaGTZs2sWjRoka2urdonl9QUEC3bt1IT09n5cqVLFy4sG5eUlISVVVVEUeqOuecc3j00UdZvHgxTzzxBABTpkzhwQcf5LjjjiMpKYlvv/2W3NzcsH/glpo8eTL3339/3ZUCS5YsYfTo0XXT77nnHsA2B3Tr1i1y7GFNAbFPVGJmwxpYt9qWA0lweByfsQDInBZKAopf0yQg3jVxyr6tXXLJJWRlZTFixIiwL+CCgoK6joK1x5RIjjzySJ577jmOPfZYli9fXtcUGK2ioiL69u1LVVUVTz/9dLM7J0bj+OOP57TTTuO6666jV69e7N69m6KiIg4//HB+9rOf8f3334c1B3Tp0oXCwsKYx9GetGOgy+zZszn99NPDps2cOZNnnnmGSy+9lIEDBzJy5EhGjRrFM888Q3JyMnPmzOGqq65i1KhRnHjiiZSXlzNx4kQGDx7MiBEjuP766xk7dmyz4ojm+VOnTqW6upqRI0fy29/+lsMOO6xu3uWXX87IkSP37lyH/VJesGABJ5xwAsnJyQBceumlDB8+nLFjx3LwwQfzk5/8JGJv3to+AbWP2qsDGnPvvfeyePFiRo4cyfDhw3nooYcA26a3Z88eDj74YEaNGsV7770XOXZTAyVvhFYYz/0BFriSlbETIbPp5iFPZbj6K5S+Y5sGlGpA//79607nu/3yl7/kpptuYuLEiXWd6CK58sor2bFjByNHjuSuu+5i5MiRZGVlRb392267jQkTJnDiiScybNiwqJ5T22mv9tFUk+vw4cO5/fbbmTx5MiNHjuTEE09ky5Yt9OzZk1mzZnHGGWcwatQozj77bABOOeUUXnrppZh2DBSRqSKySkTWiMheB1mx7nXmLxWRsc70ASLynoisEJFvRGTvNyvS9qI9RRNrGRkZpqSkJGzaihUrOPDAAz2JR8Wp8i9gnR0euNp0IzBsJ0gc5q4V5XD9D6HM6bF8/V0wbJS3MTXFGPhuMFQ79zcY8J4zhoCKFx3pmBgMBqmqqiI1NZXvvvuO448/nm+//bbux0hHFOn9E5FSY0zE06wi4ge+BU4E8oDPgHONMctdy0wHrgKmAxOAvxtjJohIX6CvMeYLEekCfA7McD83kvhtDlAKoDj067qEiWTFYwIAdojg2gSgdy4cELnTU1wRsU0C+fbsDCWvaRKg2kxpaSnHHnssVVVVGGN48MEHO3QC0ELjgTXGmLUAIvIscBrg/iI/DXjS2F/wC0UkW0T6GmO2AFsAjDFFIrICyK333L1oEqDim6s/QDFHEv3Jw3bmHiHwqDgfx8Atw50EvA7c5Wk4quPq0qVL3bgAqkG5wEZXPQ/7a7+pZXJxEgAAERkEjAE+bWqDniUBOTk5e/XuLC8vZ9u2bd4EpOKOX4o5stdHdd+nqzbuw+pt8z2NKaKqSug33D5EIDkH2qDnclvwSxITewXwSTVULOXjBf+hsqan12EpR1ZWFkVFRV6HoVqovLw80lUMARFxZ0OzjDGznHKkXw/12+wbXUZEMoEXgGuNMU32WvQsCdi9e3fYpSRg20+GDRvWYS/FUM1U9CJsspdKmpRD6DtgVHy2j85+ED74ry0fehQcd6m38TTXhqPr7iZ4xMhCyD7T44BUrRUrVpCZmanHxARkjCE1NZUxY8bUn1VtjDm0gaflAQNc9f7A5miXEZEkbALwtDHmxWjijKsG1tTUVHbt2hX19aSqg3P1ByjzHU1qaqqHwTSgsgI+cQ27e1ScjhDYGPdVAiWvNbycand6TExMxhh27drVkmPWZ8BQERksIsnAOcDcesvMBS5wrhI4DCgwxmwRmyk+Bqwwxvwt2g3GVZ+A/v37k5eXx44dO7wORXnOsB//I8n5AZRfPYH+A/t7G1Ikn38ApcW23LMvHDja03BaJHMa7HDu+FjyFpgqkL3HmFDtT4+JiSs1NZX+/Zt3zDLGVIvIz4E3AD/wuDHmGxG5wpn/EDAPe2XAGqAUuNh5+kTgR8DXIrLEmfZrY4xrBLO9xdUlgkrVqVgG34+wZV9XGLozPr+Y7voFrP7Gls+4GKaf7W08LWEMfLcPVDt9jQYugPSjvI1JqQ6qsUsEvRBXzQFK1Sn+X6icMSU+E4DN60MJgN8PE0/0Np6WEgkfhEmbBJTqNDQJUPGp+JVQOfMU7+JojPuywNGHQ1YC31Qk09UvoLjRs4dKqQ5EkwAVf6p3QNknTsUX3nEtXlRWwMdvh+rxesvgaKWfYO8qCFDxlb2roFKqw9MkQMWfknmE3TUw0MPTcCL6/MNQh8AefeDAvS4DSiz+LpA+KVR3n4lRSnVYmgSo+JNoTQFHTwVfB/hXyjw1VC6uf1WSUqoj6gBHLtWhmMrwuwZmnuxdLA3ZvB5WL7Nlvx8mTvY2nlhx7+vS9yCoI9Up1dFpEqDiS+n7UON8+SQNhuTh3sYTifuWwaMOS+wOgW5J+0CKc+dDUwklb3obj1KqzWkSoOJL/aaAeBsutaoSPn4rVD9munextAV384s2CSjV4WkSoOKHMeHjA8RjU0BYh8Deid8hsD53v4CSV8EEvYtFKdXmNAlQ8aNyOVR9b8u+LpB+jLfxRBJ2y+BpHaNDoFvqIRDoa8vBXa5LNZVSHVFURzARmSoiq0RkjYjc2Mhy40QkKCI/iF2IqtMo+m+onDEldN16vNi8Ab792pYTeYTAxogPMlxnYNxnZpRSHU6TSYCI+IEHgGnAcOBcEdmrt5az3F3YGx8o1XzFrjtfZp7uXRwN+cB1FmDUYZDd3btY2lIXvVRQqc4imjMB44E1xpi1xphK4FngtAjLXYW9j/H2GManOouq9VD+uVNJgsyTPA1nL1WVHWuEwMakHw+SZsuVK6FilbfxKKXaTDRJQC6w0VXPc6bVEZFc4HTgodiFpjqVopdC5YzjwZ/lXSyRLF4AJc6liz16w/Cx3sbTlnxp4TcUKnrBu1iUUm0qmiQg0jVa9e8/fA/wK2Ma70osIpeLyGIRWVxdXR1liKpTKHI1BXQ5w7s4GvKe69LFjtghsL4uM0NlTQKU6rACUSyTBwxw1fsDm+stcyjwrNhrunsA00Wk2hjzsnshY8wsYBZARkZG/URCdVbV26DsQ6cikBmptclD61fD2pW2HEiCo6Y2vnxHkHkykARUQcUXUPk9JA/2OiqlVIxF83PmM2CoiAwWkWTgHCCst5AxZrAxZpAxZhDwPHBl/QRAqQYVzyV0w6CjINDL03D28p6rh/whR0LXbM9CaTf+LMhwXf2gZwOU6pCaTAKMMdXAz7G9/lcAzxljvhGRK0TkirYOUHUC8dwUUFwEn84P1Y+L0xsatYUurit9NQlQqkOKpjkAY8w8YF69aRE7ARpjLmp9WKrTCOZDyTuhepc4uzTwozftlQEAA4fAvgd6G0976nIqbPUDQShfCFV5kNTf66iUUjHUwXs3qbhX/CpQZcuph0DSQE/DCVNTA/NdHQKPjcN7GbQlf3dIPzZUd1/BoZTqEDQJUN4q+k+oHG8DBH3zOezYYsvpmTB+kqfheEKvElCqQ9MkQHknmA8lrlH4up7lWSgRuTsETpwMKanexeKVLqdTd5Vw2QdQrWOBKdWRaBKgvFP0kr1vPdimgOSh3sbjtmMrfP1ZqD4pDu9o2B4Cve0VGwDUQNHznoajlIotTQKUd4qeDZW7nONdHJHMf8Xe2hjg4EOhdz9v4/FS17ND5cLZ3sWhlIo5TQKUN6q3h18V4P6i8VplBXzoug/WsZ3ossBIupwJ+G257EOo2uBpOEqp2NEkQHmj6HnAGWU67UhIGtDo4u1q4buh+wR07wUjDvU2Hq8FeoYPHFT4bMPLKqUSiiYByhvuL5KucdQUYAy8/XKoftxp4PN7Fk7c6PrDULnwGe/iUErFlCYBqv1V5dme5gD4wkem89ryL2HzeltOSYWjpngbT7zInAHiXB1R8RVUfONpOEqp2NAkQLW/wjmhcvrxtgd6vHjLNSDOxMl2fAAF/i6Q6eoboR0EleoQNAlQ7a/I9QUST00BmzfAMueyQBE4YYan4cSdsCaB2aGrJ5RSCUuTANW+KpZB+ee2LMnxda+Ad/4bKo86DHp14ssCI8mYBr4sW65aC+WLvI1HKdVqmgSo9lXwr1A5cwb4u3kWSpjiQvjk7VBdzwLszZcSPoxwwb+9i0UpFROaBKj2Y6rDvziyLvIslL0smGfHBwAYMAQOGOltPPGq6/mhcuEzUFPuXSxKqVbTJEC1n5I3ILjNlgN9w68991J1Nbzruk/Aiad3rrsFNkf6MZA0yJZr9kDxXE/DUUq1jiYBqv0UPBEqdz0fJOBZKGEWzYf8XbbctRuMO9rTcOKa+MLP4BQ87lkoSqnW0yRAtY/g7vBfjVkXeheLW00NvP5cqH7cqZCU7F08iSDrIuruLFjyJlRt9DIapVQraBKg2kfhs647Bo6DlIO8jafW0k/tpYEAqelwXCe/T0A0kvaB9OOcioGCJz0NRynVcpoEqPbhbgqIlw6BxsA818BFx0zXwYGilXVJqFzwTx0zQKkEpUmAanvlS6G8dhCe5PgZIGj1Mli70pYDSbZDoIpOl9NdYwZ85xoGWimVSDQJUG0v/8FQOfN08Od4F4vbPFdfgCNOgOzu3sWSaHxp0PXcUD3/Me9iUUq1mCYBqm0Fi6DwqVC920+9i8Vt49rwIYKnxNFNjBJF1sWhctFztvOnUiqhaBKg2lbh01BTbMvJB0JanFx+95rrLMAhR0LvXO9iSVSp4yBljC2bcsj/p7fxKKWaTZMA1XaMCW8KyL4iPgbh2b4ZPlsQqk89y7tYEpkIdLsyVM9/EEyNd/EopZpNkwDVdsoXQsVSW5Z0yLrA23hqzXs29GV14BgYNNTbeBJZ1x+CL9uWq76z4wYopRKGJgGq7exxnQXoei74sz0Lpc6OLfCx60ZBp/yw4WVV03zp4X0D9jzgXSxKqWbTJEC1jeAu21msVvYV3sXi9upsO0ogwLDRsP8IT8PpENydPUtehcrvvYtFKdUsmgSotpE/C4xzV77UQyHtUG/jgb3PApx6nnexdCTJQyFjslMxkP+wp+EopaKnSYCKPVMJe+4L1bv93LtY3PQsQNvJ/lmoXPAo1JR5F4tSKmqaBKjYK5wN1VtsOdA3fFAZr+hZgLaVeRIE9rHl4C4o+Je38SiloqJJgIotY2D3X0P1blfZoYK9NvcpPQvQlsQPOdeG6rv/AiboWThKqehoEqBiq/RtqPjaliUdsn/ibTxgRwdc+G6ofur53sXSkWVfGn65YNHLXkajVEISkakiskpE1ojIjRHmi4jc68xfKiJjXfMeF5HtIrIs2u1pEqBiy30WIOuS+LhPwEtPhO5yN3I87H+wp+F0WL7M8MGDdv9J7y6oVDOIiB94AJgGDAfOFZHh9RabBgx1HpcDrmuxeQKY2pxtahKgYqdiGZS84VQk/PSwV75dBksX2bIInHFx48ur1nE3/5Qv0rsLKtU844E1xpi1xphK4FngtHrLnAY8aayFQLaI9AUwxiwAmnUTD00CVOzs+nOonHk6JA/xLhawv0JfcN3d7rDjoP9g7+LpDAJ9oKtrZEj3Z0Ip1ZRcYKOrnudMa+4yUQu09ImtlZOTw/z5873avIqxVP8mJvR4qu7WAF+sP47CNfM9jYmyEug33D5EoN8g0M9cm0vzH8X4Ho8hYqDkFRZ9+E9KqzX5UsoREJHFrvosY8wspxzp5ir129SiWSb6YFr6xNbavXs3kyZN8mrzKta2XAoFTu/79OMZO+xnjS/f1qqr4ZafwlYnYT5hBpxwgqchdSp5L0LxfwEYP/h1yJ3jcUBKxY1qY0xDo6flAQNc9f7A5hYsEzVtDlCtV7U+/LrwHv/Pu1hqzX8llACkpsP0c7yNp7PpfnOoXPQclH/tXSxKJY7PgKEiMlhEkoFzgLn1lpkLXOBcJXAYUGCM2dLSDWoSoFpv111AtS2nHQ3pR3saDkUFdlyAWiefC12zPQunU0o7FDJPDdV3/d67WJRKEMaYauDnwBvACuA5Y8w3InKFiNTegGUesBZYAzwC1F2SIyKzgU+AA0QkT0R+3NQ2xXh0CU9GRoYpKSnxZNsqhqrWw9r97VDBAAPehozjvY3p6fvhvVdsuWdfuPVhSIqDAYs6m/IvYd3YUH3Ql5A62rNwlIoHIlJqjMnwOo5aUZ0JiGLwgvOcQQuWisjHIjIq9qGquLTzllACkHYEpB/naTjkfQ/z54XqZ1+uCYBXUsdA5hmh+k49G6BUvGkyCYhy8ILvgWOMMSOB24BZqI6vYjkUPBmq97yDussDvGAMzHkYjNNB8cAxMOow7+JR0ON3oXLxy1D+uWehKKX2Fs2ZgCYHLzDGfGyM2eNUF2J7K6qObsfNgPOFmzHV+74An70PK5bYsvjgnJ94m5QoSB0JXc4M1bf/SkcRVCqORJMENHdggh8Dr7UmKJUAyj6F4pdC9Z5/9C4WgNJieNZ1H/vjToHcQZ6Fo1x6/B7w23LpO1Ayr9HFlVLtJ5okIOqBCUTkWGwS8KsG5l8uIotFZHF1dXX0Uar4Ygxsvy5U73K2bf/10gv/hELnZFR2d5hxQePLq/aTciBkXx6qb78eTJV38Sil6kSTBEQ1MIGIjAQeBU4zxuyKtCJjzCxjzKHGmEMDAc/GKVKtVfQslH3iVJKg5+2ehsN3y+H9V0P1c38KaXHT+VYB9LgFfF1tuXIl5Gu3IaXiQTRJQJODF4jIQOBF4EfGmG9jH6aKGzWlsP2XoXrOtZC8n2fhUF0NT94bqo+cAGMnehePiizQC7r/JlTf+TsI5nsWjlLKajIJiHLwgv8HdAf+ISJL6o2LrDqS3X+B6jxb9vcKHxnOC6/NgU3rbDk5Bc67UjsDxqtuV0PSIFsO7oKdt3kajlJKBwtSzVG5Fr4/CEy5rfd5BLIv9S6e9avhj9dCMGjrZ14GU2Z6F49qWuFzsPlsp+KHQV/YKwiU6iQScrAgpTAGtv0slACkjIWsi72Lp6oSHv9rKAEYMhxOnOFdPCo6Xc6EtGOcShC2Xh4a10Ep1e40CVDRKfoPlLzuVAT6PATi9y6e/z0d3gxwyS/A52E8KjpS+9lxRnEs/xTyH/I2JqU6MU0CVNOCBbDtmlA9+2eQNs67eL5bAa/9J1SfeQn0bmzoChVXUoZBzk2h+o6boKrFd0JVSrWCJgGqaduvg+BWWw709faSwNJieOTO0CnkYaPg2FO8i0e1TPcbIXl/W64phO3XNL68UqpNaBKgGlf8ChT8M1TvdR/4s7yJxRj41z2wc5utp2XARdeBTz/GCceXCr1dzQBFz0PhbO/iUaqT0qOnalhwN2y5LFTvcg509bD3/fxX4PMPQ/WLroMefbyLR7VOxrGQ5brd+dafQtXGhpdXSsWcJgEqMmNg65WhZgB/H+hzv3fxbPgO5rhGmTv2FDjkSO/iUbHR625IGmzLNQWw5UK9WkCpdqRJgIqs4HEomhOq95kF/u7exFJcCA/eDtXOePMDhsBZlzX+HJUY/F2g77+pOxSVvgd77vEyIqU6FU0C1N4qlsO2q0L1rMugi0ed74JBePiPsGOLraekwU9ugqRkb+JRsZc+0XYUrLXjJij/wrt4lOpENAlQ4WpK7YhupszWk4dD73u8i+e5WbBiSah+6Q3Qp79n4ag20uN3dgAqAFMJm2baoYWVUm1KkwAVYgxsvRQqltm6pELuHPClexPPh2/AO/8N1U/7EYw5wptYVNuSZMidHbrTYNU62HwemKCnYSnV0WkSoEJ2/y38Mq3e90LKwd7E8s3n8G/X3QEPORJOOtebWFT7SN4f+j4Zqpe8ATtv8SwcpToDTQKUVfIm7HDdIjj7csj2qPPdutXwj9tD9wXIHQQX/0LHA+gMupwG3V2jCe66HYpe9C4epTo4PaoqKF8Cm34AOJdmpR0Ove5t7BltZ/tm+PtvocLpk9CtB1xzG6SmeROPan89boP0E0P1zedB2SfexaNUB6ZJQGdXtR7ypkNNka0HciH3BfCltH8s+bvg7t9AUb6tp2fCdX+AnJ7tH4vyjvih3zOQtJ+tm3LIOwUqv/U2LqU6IE0COrPgbtg4Daqdy+98XaH/a/b+AO0tfxf8+ZehSwGTkuGq30O/fdo/FuW9QA8Y8Br4e9h6cJfzWd3ubVxKdTCaBHRWwSLn19UKW5dkyH0ZUke0fyz5u+Avv4Jtm2zd77djAQw9qP1jUfEjeT/o/wqI0xRUtRY2TtFLB5WKIU0COqNgEeRNhbKPQ9P6/suO5d7e8nfBX26ErXm27vfD5TfB6MPbPxYVf9ImQL9nqTtUVSyBDSfas1hKqVbTJKCziZQA9LoXup7T/rFs3Qh3XGf/gu39f/lNek8AFa7LqdD3cUBsveJLTQSUihFNAjqT4O4ICcDfIeeqhp/TVtauhDt/AbucNl6fD37ya00AVGRZF0KfxwglAl/AhhOgequnYSmV6DQJ6Cyq1sP6IyMkAFe3fyxffWr7ABQX2npyiu0EqAmAakz2xdDnkVC94ktYfwRUrvYuJqUSnBhjPNlwRkaGKSkp8WTbnU75V5DnugoAbBNAe58BMAbmzYGX/2XLAJlZcPWtsO8B7RuLSlz5j8PWy6gb18LfA/q/CmnjPQ1LqWiISKkxJsPrOGppEtDRFc2FLeeHxgGQZNsJsL37AJSXwT//Cp9/GJrWozdcezv0GdC+sajEV/S/8BtdSRr0eRSyfuhtXEo1QZMAhyYBbcwEYefvYNcfQtN8WfYywIxJ7RtL3vcw607YvD40bf8R8NPfQJfs9o1FdRxln0LeSeGXDHa7Fnr9CSTJs7CUaowmAQ5NAtpQ9U7767/kjdC0wD7Q/3/tOw6AMfDe/+C5R6C6KjT9uFPhrMshEGi/WFTHVPkt5J0KlatC09KOhn6zIamfd3Ep1QBNAhyaBLSR4ldhy48huC00LWOyHYbV37394sjfZe8C+NWnoWnJKfDDn8GRk9svDtXxBQthy4VQ/HJomi8H+jwMXX/gWVhKRaJJgEOTgBirKYbtv4D8WeHTu/8GevzejsfeLnHUwIdvwH8ehTLX+ztgX7jsRug3sH3iUJ2LqYHdd8GO3wCuY1rXH0Hv+8Cf5VloSrlpEuDQJCBGjIHi/8K2q6F6Y2i6v48dYCVzWvvFsmkdPPMPWLU0fPoJM2DmJfZ+AEq1pZL59qxA9YbQNH8f6H03dDkbRDwLTSnQJKCOJgExUPm9/fIveSV8epeZ0PshexOW9lBcCP/9N7z/qj0TUKtXP7jgGhg2qn3iUAogWADbroLCf4dPTz8B+jwAyft7E5dSaBJQR5OAVgjugp1/hPz7wVSGpvt7QK+7oet57fOLp6Ic3v0fvDYHSotD030+mDwTTj3f9gNQyguFL8D2q8LHxyAA2ZdBj99BoLdnoanOS5MAhyYBLRAshD3327bPmsLwedmXQ887wJ/T9nFUlNtf/a/9B4ryw+cNGw1nX277ACjltWAh7Pyt/b/BdZZKMiDn/yDnmvbtMKs6PU0CHJoENEP1Dtjzd3sgqykIn5c6AXrfA2mHtX0cBbvtJX/z50FxvTh69rWX/Y0+TNtdVfwp/8J2nC2dHz5dMmwCnfN/kNTfk9BU56JJgEOTgCiUfw57HoTCZ0Ijo9VKPsD+8s+c0bZfusbA+jXw7lxYND/8en+AnJ5w0rkw8UQI6AAtKo4ZAyWvwY5fQcWyejOTIOs8yP4ppI7TRFa1GU0CHJoENCBYBEXPQ/6DUP7Z3vOT94ecX0HWBSBtONhOUT4sfM9e7rdp3d7zu/eCKT+Ao6Zqr3+VWEwQCufA7juh4uu956eMsmcHup6nlxaqmNMkwKFJgEtNBZS8bn/xF/9v71/9ACljofuN0OWMtrvmv6QIliyELz6EZYshGNx7mSEHwolnwJgjwN9OYw8o1RaMgZJ5sOsOKPto7/mSCpkn2UsLM08CX3r7x6g6HE0CHJ0+Cajeab/4S+ZB8WtQk7/3MpIMXc6CbldCahu1te/YAss+hyWfwMolkb/4k1PsbX4nnWyTAKU6mrKFkP+wPUMQKQmXDJsIZJ4EGVP0ygLVYpoEODpdElBTbk/vl8637ZJlCwkb2cwt5WDoegFkXRz7a/1LiuDbr+GbL+Cbz20S0JAhB8LEyTDuaEiLm8+sUm0nmA+FT0P+I1DxVcPLpR4CGdPskNyp48CX2m4hqsSmSYCjwycB1dttx76yj6B0AZQvAlPR8PJJg6HrufaRcnBsYqipga158N1y+G6FfWzZ0PhzBu1vf/WPPRJ66w1YVCdWscKeGSiaA5UrG15Okm0ikHYkpE20V+oEerZfnCqhaBLg6DBJgAlC1Xrbwaj8C6j40v6t3tTEE332YJHhnGJMGdm60/1lJZC3DvLWwsa19va9m9bZa/obk5xiR/QbPta283fv1fIYlOqIjLH/3yWv2qa7so+BCM1mboF+kDIGUkdDymj7N2lw+93DQ8UtTQIcCZUEmKAddaxqg711aeWq0KNqTfiofY1JGgrpR0H6cZA5tXmDlBgDZaWwcyts3wzbN8E25+/2zVCwJ7r1+P0wYAgMH2O/+IccqL37lWqOYAGUvG379JR9EH4b48ZICiTtZ6/wST7A+TsUkvaBQN+2vdpHxQ1NAhxxkQSYIAR3Q3CHPX0f3AHVW+2NeKo2uv5uosnMvz5Js5capY2DtKPsl3+gz97L1dRAaYkdfKeowI7Dn78L9uyE/J32b+2jqV/1kXTJhiHDYMhw+9hnP0jR9kulYqZ6uz07UPYhlH5k+xJE6lzYKD8EciFpICQNgMBAmxgEeoO/FwR62b/+7no2IcElZBIgIlOBvwN+4FFjzJ315oszfzpQClxkjPmisXXGLAkwNfaLu6bADqVbU2Azdfff2nnBAqjZE/rCD+4ibCjRlvL1Bt8QYBgE94PKQVDaE8or7K/38lL7t6wUykr4fv1Wqgry6VZTRjdTTqChDoLN4Q9An/7Qf7AdsnfAvtB/X8jqFnHxm1/+mtmfbiRoDH4Rzp0wgNtnjKib//KXm/jzG6vYnF9Gv+w0bphyADPG5DY4/9hhPXlv5Y42q9fffnM19XojaWofqMTXJu+xCdozhhVLoPxLKF9iE4Pg9hhE7LP3CAn0Al93O46BLxv82eDLCv31ZTvzuoCkgy/DXuJYW9ZEwjNNJQGt+b5t6rkRt9dUEiAifuBb4EQgD/gMONcYs9y1zHTgKieoCcDfjTETGltvzJKAqp3wXTt0wilLhZJ0KEiHXemwMxV2p9tHZfuMlFdGgM3+TDYGutF/2H7sd9BQ6J0LvXKhe0/wRfePffPLX/PUwr07CJ5/2EBunzGCl7/cxE0vfk1ZVejsR1qSnzvOGMGMMbkR57c19/abq6nXG0lT+0AlvnZ/j4P5ULna1aT4LVR9Z882BrfFfnuNkWR72aM7MfClAcnOvNpHSr16Aw8CTmLhb/7fiPN8gDgPpyzimuZ6NDS9oXlh645iPeK3+ydWu76RJKA137fRPDeSaBqhxgNrjDFrnQCeBU4D3Cs+DXjS2IxioYhki0hfY0wj15/FSmbrnl4egNJk+yhz/hal2kdhaqgcbLvMuVCS2SOp7DOoH2R2ha7Z0K0ndOvBL97axLKyJLb4MimUlLrOg7klaXx03HEt2t7sTzc2OP32GSP48xur9vqCL6sK8uc3VjFjTG7E+W3Nvf3maur1RtLUPlCJr93fY3+20zw4bu95NeVQnWf7HVVtgOoNzhnLbc7f7VC9zZ7JjAVTaR+xWl9HljIaBn/ZXltr8fctMCiK5+4lmiQgF3AfRfOw2UdTy+QCYUmAiFwOXA6Qm5vL/Pnzo9h80yZk98QYoTqYRjCYRnUwjerqtPB6bbk6jcqqrlRVZ1JVnYkxEXZBEpDjPBoj4nr47C10a//WLzv1DfnlBBGq8VEtPoL46hoDRuTWG6K0BgYN68mgugnuA1ZRi/fftQdXNThv/vz5nDOgCAZEmlvUxPy21rLX3NTrjaSpfaASX3y+xz7ssXxQxLlCFUm+ApJ9ewj4ighICQFfMQEpJuArcf7aul9K8PvK8FOOTyrwSzl+sWURb/qCJaKi4mI+j+3nISAii131WcaYWU65Nd+30Tx372CiCDjSdWv1P0HRLIPzQmeBbQ6YNGlSFJuPwoNX2y9cv995BOyp8Uj1JD+kBmzZ58yrXc5dT06xveaTkuuVkyHJqbdg2NwhN80jGKEJxi/Cd+dN2mv6b+58l035e3cyys1O46oIy0fjx03E0NQ2G5rf1lr6mpt6vZG0xX5X8aXTvsfGgCkHUwo1pVBTEiqbqtBZAlMJVIbX6x4V9cpBILj330jTovpbg/0KMbZsjKvuejR3et28miiXhy5dezNpxKRYvgPVxphDG5jXmu/bqL6H64smCcgjPF/uD2xuwTJt56c3t9umWuvcCQMitk+fOyHyT+sbphwQsd3yhikHtFkMTW0z0vy21prX3Nx9Dm2z31V86bTvsYi9eom05l2mrNpDa75vk6N47l58UQT1GTBURAaLSDJwDjC33jJzgQvEOgwoaJ/+AInn9hkjOP+wgfidtn2/SKMd1GaMyeWOM0aQm52GYH+ltLbjUlMxNLXNSPPPP2xgm9Zb85qbu8+j2Qcq8el7rOJQa75vo3nuXqK9RHA6cA+2K+fjxpg/iMgVAMaYh5xLFu4HpmIvWbjYGLO4ofVBnIwToJRSSrWjKC4RbPH3baTnNhlPpx4sSCmllGpH8TZYUDTNAUoppZTqgDQJUEoppTopTQKUUkqpTkqTAKWUUqqT0iRAKaWU6qQ0CVBKKaU6KU0ClFJKqU7Ks3ECRKQGiOUA9AGgOobr66x0P7ae7sPW033YeroPW68t9mGaMSZufoB7lgTEmogsbuSmDCpKuh9bT/dh6+k+bD3dh63XGfZh3GQjSimllGpfmgQopZRSnVRHSgJmeR1AB6H7sfV0H7ae7sPW033Yeh1+H3aYPgFKKaWUap6OdCZAKaWUUs3QIZIAEZkqIqtEZI2I3Oh1PIlARAaIyHsiskJEvhGRa5zpOSLyloisdv528zrWeCcifhH5UkReceq6D5tBRLJF5HkRWel8Hg/Xfdg8InKd83+8TERmi0iq7sOmicjjIrJdRJa5pjW430TkJud7ZpWITPEm6thK+CRARPzAA8A0YDhwrogM9zaqhFAN/MIYcyBwGPAzZ7/dCLxjjBkKvOPUVeOuAVa46roPm+fvwOvGmGHAKOy+1H0YJRHJBa4GDjXGHAz4gXPQfRiNJ4Cp9aZF3G/O8fEc4CDnOf9wvn8SWsInAcB4YI0xZq0xphJ4FjjN45jinjFmizHmC6dchD3w5mL33b+cxf4FzPAkwAQhIv2Bk4BHXZN1H0ZJRLoCRwOPARhjKo0x+eg+bK4AkCYiASAd2IzuwyYZYxYAu+tNbmi/nQY8a4ypMMZ8D6zBfv8ktI6QBOQCG131PGeaipKIDALGAJ8CvY0xW8AmCkAvD0NLBPcAvwRqXNN0H0ZvX2AH8E+nSeVREclA92HUjDGbgL8AG4AtQIEx5k10H7ZUQ/utQ37XdIQkQCJM00seoiQimcALwLXGmEKv40kkInIysN0Y87nXsSSwADAWeNAYMwYoQU9bN4vTZn0aMBjoB2SIyPneRtUhdcjvmo6QBOQBA1z1/thTYaoJIpKETQCeNsa86EzeJiJ9nfl9ge1exZcAJgKnisg6bDPUcSLyFLoPmyMPyDPGfOrUn8cmBboPo3cC8L0xZocxpgp4ETgC3Yct1dB+65DfNR0hCfgMGCoig0UkGdtxY67HMcU9ERFsO+wKY8zfXLPmAhc65QuB/7Z3bInCGHOTMaa/MWYQ9nP3rjHmfHQfRs0YsxXYKCIHOJOOB5aj+7A5NgCHiUi68399PLaPj+7Dlmlov80FzhGRFBEZDAwFFnkQX0x1iMGCRGQ6tm3WDzxujPmDtxHFPxE5EvgA+JpQe/avsf0CngMGYg8uZxpj6necUfWIyCTgemPMySLSHd2HUROR0diOlcnAWuBi7A8U3YdREpHfA2djr/r5ErgUyET3YaNEZDYwCegBbAN+B7xMA/tNRH4DXILdz9caY15r/6hjq0MkAUoppZRqvo7QHKCUUkqpFtAkQCmllOqkNAlQSimlOilNApRSSqlOSpMApZRSqpPSJEAppZTqpDQJUEoppTopTQKUUkqpTur/AyE0ZzirYnBZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "beta1_hat, beta2_hat = log_reg.params[0], log_reg.params[1]\n",
    "X = np.linspace(0, 100, 500)\n",
    "\n",
    "Z = beta1_hat + beta2_hat*X\n",
    "dPdX = np.exp(-Z)/(1+np.exp(-Z))**2*beta2_hat\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize = (8, 5))\n",
    "ax1.scatter(df['Annual_Family_Income'], df['Own_House'])\n",
    "ax1.plot(X, P_hat, color = 'tomato', lw = 3, label='Accumulative Effect')\n",
    "ax1.legend(loc ='center left')\n",
    "\n",
    "ax2 = ax1.twinx() \n",
    "ax2.plot(X, dPdX, lw = 3, color = 'Gold', label = 'Marginal Effect')\n",
    "ax2.legend(loc='center right')\n",
    "ax2.set_title('Accumulative and Marginal Effect of Family Income on Owning a House')\n",
    "ax2.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56ea2185-2b23-4d67-aa5e-ae1d23fd0f28",
   "metadata": {},
   "source": [
    "As family income raises, the marginal probability will reach maximum around $40000$ pounds, to summarize the effect, we calculate the marginal effect at the mean value of the independent variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a76056a7-0eac-426f-a3c9-b1f19c8783ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03297676614392625"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_mean = np.mean(df['Annual_Family_Income'])\n",
    "Z_mean_effect = beta1_hat + beta2_hat*X_mean\n",
    "dPdX = np.exp(-Z_mean_effect)/(1+np.exp(-Z_mean_effect))**2*beta2_hat; dPdX"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e192f7de-fd11-4d79-9ca4-247d34a88328",
   "metadata": {},
   "source": [
    "At the sample mean, $1000$ pounds increase in family income increases the probability of owning a house by $3.2 \\%$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "76c9ebf5-7e5b-433e-8a27-e990c39159f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Anaconda\\lib\\site-packages\\theano\\tensor\\elemwise.py:826: RuntimeWarning: invalid value encountered in true_divide\n",
      "  variables = ufunc(*ufunc_args, **ufunc_kwargs)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "E:\\Anaconda\\lib\\site-packages\\theano\\tensor\\elemwise.py:826: RuntimeWarning: divide by zero encountered in log\n",
      "  variables = ufunc(*ufunc_args, **ufunc_kwargs)\n",
      "E:\\Anaconda\\lib\\site-packages\\theano\\tensor\\elemwise.py:826: RuntimeWarning: invalid value encountered in multiply\n",
      "  variables = ufunc(*ufunc_args, **ufunc_kwargs)\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta2, beta1]\n",
      "E:\\Anaconda\\lib\\site-packages\\theano\\scalar\\basic.py:1955: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return x / y\n",
      "E:\\Anaconda\\lib\\site-packages\\theano\\tensor\\elemwise.py:826: RuntimeWarning: invalid value encountered in impl (vectorized)\n",
      "  variables = ufunc(*ufunc_args, **ufunc_kwargs)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 01:46&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 121 seconds.\n",
      "The number of effective samples is smaller than 25% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_0:\n",
    "    beta1 = pm.Normal('beta1', mu=0, sd=10)\n",
    "    beta2 = pm.Normal('beta2', mu=0, sd=10)\n",
    "    \n",
    "    mu = beta1 + beta2*df['Annual_Family_Income']   \n",
    "    theta = pm.Deterministic('theta', pm.math.sigmoid(mu))\n",
    "    bd = pm.Deterministic('bd', -beta1/beta2)\n",
    "    \n",
    "    yl = pm.Bernoulli('yl', p=theta, observed=df['Own_House'])\n",
    "\n",
    "    trace_0 = pm.sample(1000, return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5da7c2a-75ac-4496-86c8-b33094579064",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}