1. Simulations/GLM pt. 2 - Logistic Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Logistic Regression via `sklearn` and `statsmodels`**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **0.** Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model    import LogisticRegression\n",
    "from sklearn.metrics         import accuracy_score, precision_score, recall_score, f1_score\n",
    "from numpy.random            import uniform, binomial\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.** Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of logistic regression $Y$ is a categorical value ($0$ or $1$) and it is modeled as:\n",
    "\n",
    "$\\qquad\n",
    "\\begin{cases}\n",
    "        \\;Y \\sim Bernoulli(\\mu)\\\\\n",
    "        \\;\\log(\\frac{\\mu}{1-\\mu}) = \\beta_0 + \\beta_1 \\cdot X\\\\\n",
    "        \\;\\mu = \\frac{1}{1+e^{- (\\beta_0 + \\beta_1 \\cdot X)}}\n",
    "\\end{cases}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "# Parameters\n",
    "n_sample = 500\n",
    "beta_0   = -7.0\n",
    "beta_1   = 2.0\n",
    "liminf_x = 1\n",
    "limsup_x = 10\n",
    "\n",
    "# Data simulation\n",
    "x = uniform(liminf_x, limsup_x, size=n_sample)\n",
    "\n",
    "def logistic_function(x):\n",
    "    return 1 / (1 + np.exp(-x))\n",
    "\n",
    "mu = logistic_function(beta_0 + beta_1*x)\n",
    "y  = binomial(n=1, p=mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **2.** Logistic Regression with `sklearn`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True model: y ~ Bernoulli(μ)\n",
      "True logit(μ) = log(μ/(1-μ)) = -7.0 + 2.0 * x\n",
      "Regression logit(μ) = log(μ/(1-μ)) = -6.47 + 1.92 * x\n",
      "Accuracy:  0.9466666666666667\n",
      "Precision:  0.9537037037037037\n",
      "Recall:  0.9716981132075472\n",
      "F1:  0.9626168224299064\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFNCAYAAABIc7ibAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyyklEQVR4nO3deXwdZd338c8vadqmlLYsZWlpaa3IjoAVVFBQRLAoi7hQFRUVbsCFW5ClLWDZBAQFRMC7oN4CtyACYkWg8tyA+rA9FJAdFLC2tAUKXaB7mlzPH3PCSdKTNElzMsnJ5/16nVfmmpkz85tz0pxvr7nOTKSUkCRJUveqyrsASZKkvsgQJkmSlANDmCRJUg4MYZIkSTkwhEmSJOXAECZJkpQDQ5ikXETEzyPijE48b3RELI2I6nLU1VNFxJ0R8dW865DUdcLrhElal4iYBXwzpfR/euu+I+JrwC+AFUAD8C9gSkrp9vWtUZI6w54wSX3JgymlwcAw4ErgxogY1tU76Wu9dJI6xxAmqdMiYkBEXBoR8wqPSyNiQJPlp0TE/MKyb0ZEioh3F5b9d0ScW5jeNCJuj4jFEbEwIv4WEVURcR0wGvhj4RTkKRExprCdfoXnbhwRvyrsY1FE3LauulNKDcB1wAbANk2O5eKImB0RrxVOl9Z24Fiuiog7ImIZ8NGIGBERt0TEgoj4V0R8t8m29oiImRHxVmFfPynMHxgR10fEm4XX4pGI2Lyw7L6I+GZhuioiTo+If0fE6xFxbUQMLSxrfH2+WjiWNyJiSqffZEllYwiTtD6mAB8AdgXeC+wBnA4QEQcCJwIfB94N7NvGdk4CXgGGA5sDk4GUUjoSmA18OqU0OKX0oxLPvQ4YBOwIbAZcsq6iCz1VRwF1wL8Lsy8A3lM4lncDI4EzO3AsXwTOAzYEHgD+CDxR2M5+wH9GxAGFdS8DLkspDQHGATcV5n8VGAqMAjYBjiU7fdrS1wqPjwLvAgYDP2uxzt7AtoV9nxkR27f+ikjKgyFM0vr4EnB2Sun1lNIC4CzgyMKyzwO/Sik9k1JaDkxtYzt1wJbA1imlupTS31I7BqxGxJbAJ4FjU0qLCs/9SxtP+UBELAZWAhcDX04pvR4RARwDfC+ltDCl9DbwQ+CIDhzLH1JK9xd62XYGhqeUzk4prU4pvQxc3WR7dcC7I2LTlNLSlNJDTeZvArw7pVSfUno0pfRWiX19CfhJSunllNJSYBJwRGPvYMFZKaUVKaUnyMLge9t4XSTlwBAmaX2MoNiTRGF6RJNlc5osazrd0kXAi8CfI+LliDitnfsfBSxMKS1q5/oPpZSGARsB04EPF+YPJ+tNe7RwGnAxcFdhPrTvWJrO2xoY0bitwvYmk/XyAXyDrNft+cIpx08V5l8HzCAbqzYvIn4UETUl9lXqde/XZPsArzaZXk7WWyapBzGESVof88gCR6PRhXkA84Gtmiwb1dpGUkpvp5ROSim9CzgYODEi9mtc3Mb+5wAbd3RwfaH36DjgyIjYDXiD7LTfjimlYYXH0MIg/vYeS9M65wD/arKtYSmlDVNKEwr7/2dKaSLZ6dMLgZsjYoNCT95ZKaUdgA8BnwK+UmJfpV73NcBrHXkdJOXLECapvWoKA8cbH/2AG4DTI2J4RGxKNobq+sL6NwFHRcT2ETEIaPWaYBHxqYh4d+G04BKgnuwyEpAFi3eVel5KaT5wJ3BlRGwUETUR8ZH2HExKaSFwDXBm4RTi1cAlEbFZoaaRTcZwtftYCv4f8HZEnBoRtRFRHRE7RcT7C9v+ckQML+x3ceE5DRHx0YjYuTBm7S2y05MNJbZ/A/C9iBgbEYPJTp3+NqW0pj3HLqlnMIRJaq87yHqLGh9TgXOBmcCTwFPAY4V5pJTuBH4K3Et2qrFx3NOqEtveBvg/wFLgQeDKlNK9hWXnkwW9xRHx/RLPPZIsrDwPvA78ZweO6VJgQkTsApzaWGdEvFWoZ9tOHAsppXqyXqxdya5H9gZZ4BtaWOVA4JmIWEo2SP+IlNIKYAvgZrIA9hzwF7JTlC39sjD/r4XtrwS+04HjltQDeLFWSd2i8O28p4EBvb3HppKORVJ+7AmTVDYRcVjh+lsbkY19+mNvDS2VdCySegZDmKRy+g+yU4QvkY3zOi7fctZLJR2LpB7A05GSJEk5sCdMkiQpB4YwSZKkHPRb9yo9y6abbprGjBmTdxmSJEnr9Oijj76RUhpealmvC2Fjxoxh5syZeZchSZK0ThHx79aWeTpSkiQpB4YwSZKkHBjCJEmScmAIkyRJyoEhTJIkKQeGMEmSpBwYwiRJknJQtuuERcQvgU8Br6eUdiqxPIDLgAnAcuBrKaXHylVPXm57fC4XzXiBeYtXMGJYLScfsC2H7jYy77K6TXcef16vdan9Au2u5fTbnuKGh+dQnxIB9O9Xxao1DQAMq61h6sE7cuhuI5vtZ2htDXX1DSxbXf/OdkqtO3fximb7qgr44p6jOffQnUvuvzqCiXuOara86baqI6hPiZEljqm967X1Ok6d/gyLV9QBsEH/amqqq1iyoo4Rw2oZs0ktD7y8kMbb3favDgb17/fO8va87k1rjOCdbW00qIaDdtmSe59fsM73EeCsPz7DouV177zun3rv2s9tecylXp+NBtWwqq6e5XVrv99d4bbH565Va6nfpxHDavnodsO59/kFa/3ONNpoUA0r6+pZUai15bIdttyQB15aSNO7EX/5A81/19al6e8iQACt3d24rdf9tsfnMunWJ0vWOqimigE11Sxe3vrvTdPXouXv8u9mzub+lxa+s72aKljTQLPntfdvUKnf+YaU3ql7o0E1/ODTnft9aOu978y2Wr4+tz8x/526W2q5r87+bWj5Ny+Cd963MZvU8tDLi975XQHW+fep6esBvPM3oCN/p7pa2W7gHREfAZYC17YSwiYA3yELYXsCl6WU9lzXdsePH596y8Vasz8ET7GirvhBWVtTzfmf2blPBLHuPP68XutS+62pCgioqy/+22qtltNve4rrH5rd5j5qqoIv7DGKWx6d22w/67Nu44dja/tvXF7q+EodU3vXa81tj8/l5N89QV1D5/8eret1b6vGkturDkg0q6mmKmgA6tdRZ8tj7si+a6qCiz733vX+vb3t8bmcfPMTzV6Pxu239/epK7Q3iLXn38JaUqIqNRBApMSgmioO220Ev314NiklqlIiSEThPziN0xTWr0oN9K/KPonrGxqIVt/W1kNhtPIZGiQG9qtm0oTtOHCnLZvVDHDX0/M57/bnmv1+RYnI2a86OP2T23PATluUPP5SZjz9Kuff8SxrWvye9otg8oTt+MSOJbbVyvb+/MyrXHzX86xcs3agLVVvo5qAUz+5HQA//vM/WFnid21gTTUnfeI97L/D2vXc/eyrrT6vPZpu++5nX+XCu55nTX3peucPGU5VbW3ZPi8i4tGU0viSy8oVwgo7HgPc3koI+y/gvpTSDYX2C8C+KaX5bW2zN4WwvS64p+T/KkcOq+X+0z6WQ0XdqzuPP6/XurX9llKqlnGT7mj2P7nWNP4Prz3as251BC+dP6HV/TcuX9fxNR5Te9drTUdex45qb43l2i90/Phafb1SgjVrYOXK4mPVqubtwmPKjTNZumQpA9bUMaB+NQPWrGbAmjr619dRW19HzZrCvPo6qhoaqE4NVDU00C/VU93QQFVqoF9D/TvLqlMD1Q3FZU3b1amBqlRPv4ZsG1kwSu+Eo80G989qb2jIfpZ4LF1ZVwhLiUhNA1MhTKXi9qra+PCXOuqQI3/MEyO2LdvnRVshLM/bFo0E5jRpv1KYt1YIi4hjgGMARo8e3S3FdYV5rfzRbW1+penO48/rte7I9kut295g1d712rtu4zqtrds4f13H17i8vet1dvn6aG+NXb7fRcth+XJYvJhBL77A7iuXMWTVUoasWsbQlUsZsnIZQ1YtY8jKbN4Gq1cyoH41/dfUMaC+Dq7tXzpsNazdI1HKeWU+vg5Ztu5VBpe/CqlNeXw294p7R6aUpgHTIOsJy7mcdhsxrLbk/35HDKvNoZru153Hn9dr3dp+W1u3pfb2cJWjJ6ytdRuXr+v4Go+pveu1tbxcvVTtrXGdUmLD1csZvnQRw5ctYrOlCxm+bHE2vWwhw5cuykJWIVwNXbUMfrQGgLs7s7/XOl9qX9JAkCI7MZb9DFIABA1RRQoK8wrrEBBBwzvrNn8uUXo/qbUFheeWfg5UVwXDBw9oviCC199eVfK0dqltVVcFmw8ZWHrnJdZ/9a2VhVORpbe1ZVv/Hltsb+7ildQ3NLRxjK2/LtXV2Xf/1rTxH4d+VVWM2mjteuYsWtHm89qjcdvr2tbqfjVAPp/NeYawucCoJu2tCvMqxskHbFtynFLjQNBK153Hn9drXWq/rY1NKlXLxD1H5TImbOKeo9rcf+PyUsdX6pjau15rTj5g27KNCWtPjQBVDfVs+fabjFryKqMWv8bWb73GqMWvMmLx62y+9E2GL1tM7ZpVna6vy/TrBwMGwMCBaz+azJ+3MjHztRWsrOrHqn79WVVdw6p+/VlT059txw7nsVdXsDT6UVfVjzVV1dRXVVEf2c+GqGrWro8qGqqqWFNVnS2LKuqrisvqo7rZ8oYmYegz79uKUz+5A1RVZR/wpR5VVZzzp2e58ZFXmgWp7GeJIPVOYGoeAGprqjn8fSP5zUOzae/Hd6mxf61p64sCpbQ1HvKBdo6DrKkOLvpsx8YIPtTGeMCOjjd8pINjKVvuC1jneNFRJep5tJP7LbXtR1t5PVqun8dnc54hbDrw7Yi4kWxg/pJ1jQfrbRp/0fvqtyO78/jzeq1b2297a2kcsNyeb0eO33rjdn87snHddX07suX+W347sunxtfXNpvaut67XsZzfjmz8+fPfP8LQF59nuzdm8Z4F/2b0olcZ89ZrjHhrAdX1a9qss8P694eNNoJhw1jYfwNeWFXNgupa3h44mCUDBlE3eAiLagbxRs0glvavpWbQII786LZ8ZJdRa4Wqd9r92vdnewTQ8Phczi/xDbkDdxvJym76duSp7fx25Blf/CCrBnXNtyPHb71xj/92ZGu/813x7cjG9bvi25Gl/sZ19NuRjc/vyN+Glvtdn29Hlno9oPK/HXkDsC+wKVnH+g+AGoCU0s8Ll6j4GXAg2SUqjkoprXPEfW8amC8pR6tXwwsvwFNPwZNPZo+nnoJXXun8NmtrYcsts8cWW2SPxunNN4dNN4Vhw4qPga2cQpLUZ+QyMD+lNHEdyxPwrXLtX1If0tAAzz4L998PDzwAf/87PPcc1JX+n3qrNt8cxo7NHu96V3F61KgsaA0eXHIMjiR1Rq8YmC9JzdTVwcyZcM898H//Lzz4ICxZ0r7nDhgAO+4IO++cPbbbLgtaY8bAoEFlLVuSmjKESeod5syBP/0J7rgD7rsP3n573c/ZemvYZZfssfPO2c9ttmn3uCpJKif/EknqmRoa4OGHYfr0LHw99VTb62+2Gey1V/bYY48scA0d2j21SlInGMIk9RwpwUMPwU03wc03tz2IfvRo2G8/2Gcf2HvvbAyX47Uk9SKGMEn5mzMHrr0WfvUreOml0uv07w/77gsHHQQTJsC4cYYuSb2aIUxSPlavht//Pgtef/5z6ZsRb7IJHHoofPrTWa/XYG9uI6lyGMIkda/ly+Gaa+Cii0qfbhw6FD77Wfj85+GjH4Wamu6vUZK6gSFMUvdYvBiuvBIuvRQWLGi+LCLr6fr617Oer9q+cX9VSX2bIUxSeb3+eha8rrgC3nqr+bLNNoPjj4evfS27nIQk9SGGMEnlMWcOXHwxXH01rGhxP8LRo+GUU7KeL3u9JPVRhjBJXWvWLDjnnOzbjmta3BB7221h0iT44hcd6yWpzzOESeoaS5fCeefBJZfAqlXNl+22G0yZko33qq7OpTxJ6mkMYZLW3913w9FHw7//3Xz+hz8MkyfDAQd4TS9JasEQJqnzFi2Ck07KrvXV1Pvfn12CYp998qlLknoBQ5ikzrntNjjuOHj11eK8jTeGn/wEjjwSqqpyK02SegNDmKSOee01+M534He/az7/c5+Dyy+HzTfPpy5J6mUMYZLaJyX4n/+BE06AhQuL87fYIrsI62GH5VebJPVChjBJ67Z4MXzzm3DLLc3nH3UU/PjHsNFGuZQlSb2ZIUxS2x57LLuX47/+VZy39dYwbRp84hP51SVJvZwjZyW17tZbYe+9mwewY4+Fp582gEnSejKESVpbSnDhhXD44cVbDg0Zkg3Gv+oqGDw43/okqQJ4OlJSc6tXZ71dTa/9NW4c3H47bLddfnVJUoUxhEkqevPNrPfrL38pzvvIR7LTkptskl9dklSBPB0pKfOPf8AHP9g8gH31q/DnPxvAJKkMDGGS4L774AMfgH/+szjv/POzU5IDBuRWliRVMk9HSn3dH/+YXYJi9eqsXVsL112XnZaUJJWNIUzqy37/e/jCF6CuLmtvsQVMn57dgFuSVFaGMKmvuvlmmDgR1qzJ2u96F9xzT3YhVklS2TkmTOqLbrwRjjiiGMC22SYbkG8Ak6RuYwiT+pqbboIvfQnq67P2dttlA/O32irXsiSprzGESX3JXXfBl78MDQ1Ze8cdswA2YkSuZUlSX+SYMKmvuP9++MxnioPwt98e7r0Xhg/Pty5J6qPsCZP6gieegIMOKt4HcvTo7CKsBjBJyo0hTKp0s2fDgQfCkiVZe7PN4O67HQMmSTkzhEmV7K23sh6wV1/N2kOGwIwZ8J735FuXJMkQJlWsujr43Ofg6aezdk0N/OEPsOuuuZYlScoYwqRKlBJ85zvZuK9G11wD++6bW0mSpOYMYVIlmjYN/uu/iu0zz4SvfCW/eiRJazGESZXmkUfgu98ttr/4RZg6NbdyJEmlGcKkSvLGG3D44bB6ddZ+73vh6qshIt+6JElrKWsIi4gDI+KFiHgxIk4rsXx0RNwbEY9HxJMRMaGc9UgVrb4+6/WaMydrDx0Kt9wCgwblW5ckqaSyhbCIqAauAD4J7ABMjIgdWqx2OnBTSmk34AjgynLVI1W8qVOz6381uv56GDcut3IkSW0rZ0/YHsCLKaWXU0qrgRuBQ1qsk4AhhemhwLwy1iNVrttvh3PPLbanTIFPfSq/eiRJ61TOe0eOBOY0ab8C7NlinanAnyPiO8AGwMfLWI9UmV5+GY48stjef38466z86pEktUveA/MnAv+dUtoKmABcFxFr1RQRx0TEzIiYuWDBgm4vUuqxVqzIBuIvXpy1R42C3/wGqqtzLUuStG7lDGFzgVFN2lsV5jX1DeAmgJTSg8BAYNOWG0opTUspjU8pjR/uDYelopNPhr//PZvu3x9uvhk2XeufkCSpBypnCHsE2CYixkZEf7KB99NbrDMb2A8gIrYnC2F2dUntMWMGXHFFsX3ppbDHHrmVI0nqmLKFsJTSGuDbwAzgObJvQT4TEWdHxMGF1U4Cjo6IJ4AbgK+llFK5apIqxsKFcNRRxfYhh8Cxx+ZXjySpw8o5MJ+U0h3AHS3mndlk+llgr3LWIFWclOC442D+/Ky92WbZbYq8IKsk9Sp5D8yX1FE33AA33VRsX311FsQkSb2KIUzqTebMgeOPL7a/8Q04+ODW15ck9ViGMKm3aGjIxoEtWZK1x46FSy7JtyZJUqcZwqTe4mc/g//932w6Aq69FjbcMN+aJEmdZgiTeoPZs2Hy5GL71FNh773zq0eStN4MYVJPlxJ861uwbFnW3mGH7GbdkqRezRAm9XS33JLdoLvRtGkwYEB+9UiSuoQhTOrJliyB73632D72WNjLS+tJUiUwhEk92aRJxYuybrEFnH9+vvVIkrqMIUzqqR54AK66qtj+6U9h2LDcypEkdS1DmNQTrV4NxxxTbB90EHz2s/nVI0nqcoYwqSe6+GJ45plsetAguOIK7w0pSRXGECb1NC++CGefXWyfey5svXV+9UiSysIQJvUkKcF3vgOrVmXt3XfP2pKkimMIk3qS22+Hu+7KpiOya4L165dvTZKksjCEST3FqlXwve8V28ccA+97X371SJLKyhAm9RSXXAIvvZRNDxuWjQWTJFUsQ5jUE8yb1zx0nXMObLppfvVIksrOECb1BKeeWrxB9047ZbcnkiRVNEOYlLcHH4Trry+2L7vMwfiS1AcYwqQ8pQQnnVRsH344fOxj+dUjSeo2hjApT9OnZz1hAP37w0UX5VuPJKnbGMKkvKxZA5MnF9vHHQdjx+ZXjySpWxnCpLxcey08+2w2veGGMGVKvvVIkrqVIUzKw4oV8IMfFNsnnwzDh+dXjySp2xnCpDxccQW88ko2vfnmza+UL0nqEwxhUndbvBh++MNi+8wzYfDg3MqRJOXDECZ1twsvhEWLsulx4+Doo/OtR5KUC0OY1J3mzcsuxtro3HOhpia/eiRJuTGESd3prLOyQfkAu+8On/98vvVIknJjCJO6ywsvwC9+UWxfcAFU+U9QkvoqPwGk7jJlCtTXZ9P77Qf7759vPZKkXBnCpO7w8MNwyy3F9gUX5FeLJKlHMIRJ5ZYSnHZasf35z8P48fnVI0nqEQxhUrnNmAH33ZdN9+uXfSNSktTnGcKkcmpoaN4L9s1vwjbb5FePJKnHMIRJ5XTDDfDEE9n0oEHZ1fElScIQJpXP6tVwxhnF9ve+B1tumV89kqQexRAmlcu0afCvf2XTm2wCJ5+cbz2SpB7FECaVw4oVzW/SPXkyDB2aXz2SpB6nrCEsIg6MiBci4sWIOK2VdT4fEc9GxDMR8Zty1iN1m6uvhvnzs+ktt4Tjj8+3HklSj9OvXBuOiGrgCmB/4BXgkYiYnlJ6tsk62wCTgL1SSosiYrNy1SN1mxUrml+MddIkGDgwv3okST1SOXvC9gBeTCm9nFJaDdwIHNJinaOBK1JKiwBSSq+XsR6pe7TsBTv66HzrkST1SOUMYSOBOU3arxTmNfUe4D0RcX9EPBQRB5axHqn87AWTJLVT2U5HdmD/2wD7AlsBf42InVNKi5uuFBHHAMcAjB49uptLlDrAXjBJUjuVsydsLjCqSXurwrymXgGmp5TqUkr/Av5BFsqaSSlNSymNTymNHz58eNkKltaLvWCSpA4oZwh7BNgmIsZGRH/gCGB6i3VuI+sFIyI2JTs9+XIZa5LKx14wSVIHlC2EpZTWAN8GZgDPATellJ6JiLMj4uDCajOANyPiWeBe4OSU0pvlqkkqG3vBJEkdFCmlvGvokPHjx6eZM2fmXYbU3E9/CieckE1vuSW8/LIhTJJERDyaUhpfaplXzJfWl71gkqROMIRJ68uxYJKkTjCESetj5Up7wSRJnWIIk9bHr35lL5gkqVMMYVJn1dU17wU75RR7wSRJ7WYIkzrr+uth9uxsevhwe8EkSR1iCJM6o74efvjDYvvEE2GDDfKrR5LU6xjCpM646SZ48cVsetgwOP74XMuRJPU+hjCpoxoa4Lzziu0TToAhQ/KrR5LUKxnCpI76wx/gmWey6cGD4bvfzbceSVKvZAiTOiIlOPfcYvtb34KNN86vHklSr2UIkzrirrvgscey6drabEC+JEmdYAiT2islOOecYvuYY2CzzfKrR5LUqxnCpPa67z548MFsun9/+P73cy1HktS7GcKk9mo6Fuyoo2CrrfKrRZLU6xnCpPZ44AG4555suroaTj0133okSb2eIUxqj6a9YF/+Mowdm18tkqSK0GoIi4g7ImJMN9Yi9UyPPgp33plNR8CkSfnWI0mqCG31hP0K+HNETImImu4qSOpxml4d//Ofh223za8WSVLF6NfagpTS7yLiTuAMYGZEXAc0NFn+k26oT8rX00/D739fbE+enF8tkqSK0moIK1gNLAMGABvSJIRJfcIPf1icPuQQ2GWX/GqRJFWUVkNYRBwI/ASYDuyeUlrebVVJPcE//gG//W2xPWVKfrVIkipOWz1hU4DPpZSe6a5ipB7lggugodD5e8AB8P7351uPJKmitDUm7MPdWYjUo8yaBdddV2yffnpupUiSKpPXCZNK+dGPYM2abHqffWDvvfOtR5JUcQxhUkvz5sEvflFs2wsmSSoDQ5jU0sUXw+rV2fSee8J+++VbjySpIhnCpKYWLICf/7zYPuOM7Cr5kiR1MUOY1NQll8CKFdn0rrvChAm5liNJqlyGMKnRokXws58V26efbi+YJKlsDGFSo8svh7ffzqa33x4OOyzfeiRJFc0QJkEWvi69tNieMgWq/OchSSofP2UkgKuuyk5HAowbB1/4Qr71SJIqniFMWr4cfvzjYnvSJOi3rnvbS5K0fgxh0tVXw+uvZ9OjRsGRR+ZbjySpTzCEqW9btSq7RVGjU0+F/v3zq0eS1GcYwtS3/epX2W2KALbYAr7+9XzrkST1GYYw9V0rV8J55xXb3/8+1NbmV48kqU8xhKnvuuYaeOWVbHrzzeG44/KtR5LUpxjC1DetWAE//GGxfdppMGhQfvVIkvocQ5j6pmnTYP78bHrLLeE//iPfeiRJfU5ZQ1hEHBgRL0TEixFxWhvrHR4RKSLGl7MeCciuC3b++cX25MmOBZMkdbuyhbCIqAauAD4J7ABMjIgdSqy3IXAC8HC5apGaueoqeO21bHqrreCb38y3HklSn1TOnrA9gBdTSi+nlFYDNwKHlFjvHOBCYGUZa5EyS5fChRcW25Mnw8CB+dUjSeqzyhnCRgJzmrRfKcx7R0TsDoxKKf2pjHVIRVdcAQsWZNOjR3tdMElSbnIbmB8RVcBPgJPase4xETEzImYuaPwAlTrqrbeaXx3/9NNhwID86pEk9WnlDGFzgVFN2lsV5jXaENgJuC8iZgEfAKaXGpyfUpqWUhqfUho/fPjwMpasinb55bBwYTY9dix87Wu5liNJ6tvKGcIeAbaJiLER0R84ApjeuDCltCSltGlKaUxKaQzwEHBwSmlmGWtSX7VkCVx8cbF9xhlQU5NfPZKkPq9sISyltAb4NjADeA64KaX0TEScHREHl2u/UkmXXQaLF2fT48bBkUfmWo4kSf3KufGU0h3AHS3mndnKuvuWsxb1YYsWwU9+Umz/4AfQr6y/+pIkrZNXzFflu+SS7HQkwHveAxMn5luPJEkYwlTp3nwTLr202LYXTJLUQxjCVNl+/GN4++1sevvt4QtfyLceSZIKDGGqXAsWwE9/WmxPnQrV1bmVI0lSU4YwVa6LLoJly7LpnXaCz34233okSWrCEKbK9Npr8LOfFdtnnQVV/rpLknoOP5VUmc4/H1asyKZ33RUOPTTPaiRJWoshTJVn1iy46qpie+pUe8EkST2On0yqPGeeCatXZ9Mf/CAc7A0aJEk9jyFMleWJJ+D664vtCy+EiPzqkSSpFYYwVZZJkyClbPpTn4IPfzjfeiRJaoUhTJXjvvvgzjuz6YhscL4kST2UIUyVYc0aOOGEYvurX82uDSZJUg9lCFNluOoqePLJbHrQIDj77HzrkSRpHQxh6v1eew3OOKPYPv10GDUqv3okSWoHQ5h6v1NPhSVLsulttoETT8y3HkmS2sEQpt7tgQfg178uti+/HAYMyK8eSZLayRCm3qu+Hr71rWL7sMPggAPyq0eSpA4whKn3+vnP4e9/z6Zra+GSS3ItR5KkjjCEqXdasCAbgN9oyhTYeuv86pEkqYMMYeqdTjsNFi/Opt/9bvj+93MtR5KkjjKEqfd58EH45S+L7Z/+1MH4kqRexxCm3qW+Hr797WL7kEPgk5/Mrx5JkjrJEKbeZdo0eOyxbHrgQLj00lzLkSSpswxh6j3eeCMbgN9o8mQYMya3ciRJWh+GMPUekybBokXZ9LhxcPLJ+dYjSdJ6MISpd3j4YbjmmmL7ssuy05GSJPVShjD1fC2vjP/pT8NBB+VXjyRJXcAQpp7vmmvg0Uez6QEDHIwvSaoIhjD1bG++mQ3AbzRpErzrXfnVI0lSFzGEqWebPBkWLsymx46FU07Jtx5JkrqIIUw91yOPwNVXF9uXXZbdqFuSpApgCFPP1NCQDcZPKWsfdFA2IF+SpAphCFPP9ItfZD1hkA3Gv+yyfOuRJKmLGcLU88ybB6edVmyfemp2cVZJkiqIIUw9S0rw9a8XB+OPGdM8kEmSVCEMYepZrrwSZszIpiPgv//bwfiSpIpkCFPP8fzzze8HedJJsM8++dUjSVIZGcLUM9TVwZFHwooVWXvnneHcc/OtSZKkMjKEqWc45xyYOTOb7t8frr8++1akJEkVqqwhLCIOjIgXIuLFiFhrdHVEnBgRz0bEkxHxvxGxdTnrUQ/10ENw3nnF9nnnwS675FePJEndoGwhLCKqgSuATwI7ABMjYocWqz0OjE8p7QLcDPyoXPWoh3rzTZg4Mbs4K8C++8KJJ+ZakiRJ3aGcPWF7AC+mlF5OKa0GbgQOabpCSunelNLyQvMhYKsy1qOepr4evvQlmDUraw8Zkn0bssqz5JKkylfOT7uRwJwm7VcK81rzDeDOMtajnubss4uXowC49lrY2jPSkqS+oV/eBQBExJeB8UDJ6xFExDHAMQCjR4/uxspUNrffnoWwRpMnwyGHtL6+JEkVppw9YXOBUU3aWxXmNRMRHwemAAenlFaV2lBKaVpKaXxKafzw4cPLUqy60UsvZZejaPTxjzcPZJIk9QHlDGGPANtExNiI6A8cAUxvukJE7Ab8F1kAe72MtainWL4cDj8cFi/O2qNGwQ03QHV1rmVJktTdyhbCUkprgG8DM4DngJtSSs9ExNkRcXBhtYuAwcDvIuLvETG9lc2pEqQExx4LTzyRtfv3h5tvhk03zbcuSZJyUNYxYSmlO4A7Wsw7s8n0x8u5f/UwP/85XHddsX355bDHHvnVI0lSjrwWgLrHAw/ACScU20cdBUcfnV89kiTlzBCm8nv+eTj44Oz+kAC77QZXXAER+dYlSVKODGEqr7lz4YADsivjQzb+65ZboLY237okScqZIUzls3gxHHggzJ6dtQcNgj/9CcaOzbUsSZJ6AkOYymPlyuziq08/nbX79ct6wByIL0kSYAhTOTTeE/Kvfy3O++Uvs14xSZIEGMLU1VKCb38bbr21OO+ii5pfIV+SJBnC1IVSgqlTs+uBNTrxRPj+93MrSZKknqpH3MBbFaChIQtbl1xSnPelL2W9YJIkaS2GMK2/1auzi6/+5jfFeZ/4RDYOrMrOVkmSSvETUutn6dLsQqxNA9jhh8Mf/pDdG1KSJJVkCFPnLVgAH/sYzJhRnHfssfDb38LAgfnVJUlSL2AIU+fMmgV77w2PPFKcN3UqXHklVFfnVZUkSb2GY8LUcU8+mV3za/78rB2Rha9jj823LkmSehFDmDrmb3+DT38alizJ2v37Z+PBDj8837okSeplPB2p9rvtNth//2IAGzIkGw9mAJMkqcMMYWqfadOysLVqVdbefHP4y19g331zLUuSpN7KEKa2LVwIX/wi/Md/ZBdkBRg3Dh54AHbdNdfSJEnqzRwTptbNmAFf/zrMm1ect/vucMcdWU+YJEnqNHvCtLZly+D447NvQDYNYF//Otx3nwFMkqQuYE+YmnvwQfjKV+DFF4vzNtsMrr46uzK+JEnqEvaEKbN6NUyZkl2AtWkAO+wwePppA5gkSV3MnjDBU0/BkUfCE08U5w0ZApdfns2PyK82SZIqlD1hfVl9PVx0EYwf3zyAfexjWTD7ylcMYJIklYkhrK/6f/8P9toLTjklOxUJ2U23L7sM7r4bRo/Otz5JkiqcpyP7mueeg7PPhhtvbD5//Hi47jrYbrt86pIkqY+xJ6yveO657KKrO+7YPID17w9Tp2YXXzWASZLUbewJq3TPPgvnnAO//S2k1HzZYYdlY8LGjcunNkmS+jB7wirVs8/CxImw005Zz1fTADZhAjz8MNx6qwFMkqSc2BNWSRoaslsNXXFFdmuhlj1fEybAD34Ae+yRT32SJOkdhrBK8PrrcP31cNVVzS+02sjwJUlSj2MI661WrYLbb4df/xruvBPWrFl7nYMOgjPPNHxJktQDGcJ6k9Wr4Z57srFct9wCCxeuvc7QodmNto87DrbZpvtrlCRJ7WII6+mWLs3Ged16a9bz9dZbpdf70Ifga1/LLkOxwQbdWqIkSeo4Q1hP89ZbcP/98Je/ZI+ZM0ufaoTsqvZf+Ur2sNdLkqRexRCWt0WL4G9/ywLXX/8Kjz2WfcuxNWPGwGc+k13j60MfgiqvMiJJUm9kCOtu8+fDgw8We7qefHLtS0m0tPPOcOihWfh673u9qbYkSRXAEFYuy5Zltwp65pns8dRT8Pjj8NprbT8vIgta++yTPT78Ydh00+6pWZIkdRtD2PpYuRJmzYKXXy4+/vnPLHTNmrXuHi6A6mrYffdi6Np7bxg2rMyFS5KkvBnCWtPQAG+8AXPnZo9587Kfs2cXA9fcue0LWk1tsAHstlvWw7XPPtm4rg03LM8xSJKkHssQ1tSsWdklHubOzcZu1dV1flvV1fDud8OOOxYfu+2WzXMwvSRJfV5ZQ1hEHAhcBlQD16SULmixfABwLfA+4E3gCymlWeWsqS13vLSECQ8+2O7166OKeUOGM3vY5rwybAv+PXQLZg/bgtdHjmX7j76f6c+9yaLldbAaeBxqn36JgTWzWLy8jhHDajn5gG05dLeRa233tsfnctGMF5i7eMVay4bV1hABi5bXUR1BfUoE0Ngft9GgGg7aZUvufX4B8xavWGs/p9/2FDc8PIf6Jj141RFM3HMU5x66c8njbPqcttZtrLvUfkstA945zsZjGTmslo9uN7zV+tvjtsfnctYfn8le+8JrNvXgHUvWMrTweq7rPemMtl6PrtpWZ/bRlXWtj55ShyTlJVJHT6e1d8MR1cA/gP2BV4BHgIkppWebrHM8sEtK6diIOAI4LKX0hba2O378+DRz5swur/e2x+dyyu8e5+kLD6N/Q3ZdriUDNuDVDTfhtcHZ49UNs8fsYVnYmjdkOHXVNZ3eZ21NNed/ZudmHzy3PT6XSbc+xYq6+vU+ppb7mfnvhVz/0OxW1/vyB0avFa5Ov+2pks9puW6puhv3C6y1rKYqIKCuft2/f6Vep9bc9vhcTr75ibW2W1MVXPS595aspbP7Wlcdrb0eHd12a9s6/H0jueXRuR3aR1fWtT56Sh2SVG4R8WhKaXzJZWUMYR8EpqaUDii0JwGklM5vss6MwjoPRkQ/4FVgeGqjqHKFsL0uuIe5i1ew+9zneHPQUF7fYGNW9B/Y5ftpaeSwWu4/7WNr1VGO/by6ZGWzHrCWqiN46fwJzeaNm3RHyee0XLe1ukcOqwVY72Nq+Tq1pq3Xr721tHdfnamjM9tubVuNvYcd2UdX1rU+ekodklRubYWwcp6OHAnMadJ+BdiztXVSSmsiYgmwCfBG05Ui4hjgGIDRo0eXpdh5hQ+Ex0ZuX5btr2u/rbW7cj/ritulPtBbC20t57dWd1cdT3u309Z6XbGN9urK16O157T23nTmNSjX711Pr0OS8tQrRoinlKallManlMYPHz68LPsYUegl6W4t91uuOkYMq6V6HRd5LbW8tee0nN9a3SOG1XbJMbV3G22t195ayllvZ7bd2nNae2/W9Rp0VV3ro6fUIUl5KmcImwuMatLeqjCv5DqF05FDyQbod7uTD9iWmuruvRJ9bU31OwPUm9ZRW1Ndlv1M3HNUm+uVWt7ac1rOL1V3435LLaupina/3qVep9a09j7WVEWrtXR2X+uqo7XXo6u2NXHPUR3eR1fWtT56Sh2SlKdyno58BNgmIsaSha0jgC+2WGc68FXgQeCzwD1tjQcrp8bBwE2/VVdbU0VVBMtWZ4OHG7+FOKimiuV1ze/vWFMFjbMav6H4pyfnv7Otxu0NrKlu85t4je1yfDuycdsd+XZk47x1fTuyad2tfdutO74dWep9bPntyKa1lOvbke15PbpiW+O33rhD++jKutZHT6lDkvJUtoH5ABExAbiU7BIVv0wpnRcRZwMzU0rTI2IgcB2wG7AQOCKl9HJb2yzXwHxJkqSultfAfFJKdwB3tJh3ZpPplcDnylmDJElST9QrBuZLkiRVGkOYJElSDgxhkiRJOTCESZIk5cAQJkmSlANDmCRJUg4MYZIkSTko68VayyEiFgD/zruOXmxTWtwgXb2S72Nl8H2sDL6PlaFc7+PWKaWSN77udSFM6yciZrZ25V71Hr6PlcH3sTL4PlaGPN5HT0dKkiTlwBAmSZKUA0NY3zMt7wLUJXwfK4PvY2XwfawM3f4+OiZMkiQpB/aESZIk5cAQ1gdExKiIuDcino2IZyLihLxrUudFRHVEPB4Rt+ddizovIoZFxM0R8XxEPBcRH8y7JnVcRHyv8Hf16Yi4ISIG5l2T1i0ifhkRr0fE003mbRwRd0fEPws/Nyp3HYawvmENcFJKaQfgA8C3ImKHnGtS550APJd3EVpvlwF3pZS2A96L72mvExEjge8C41NKOwHVwBH5VqV2+m/gwBbzTgP+N6W0DfC/hXZZGcL6gJTS/JTSY4Xpt8n+2I/Mtyp1RkRsBRwEXJN3Leq8iBgKfAT4BUBKaXVKaXGuRamz+gG1EdEPGATMy7ketUNK6a/AwhazDwF+XZj+NXBoueswhPUxETEG2A14OOdS1DmXAqcADTnXofUzFlgA/KpwavmaiNgg76LUMSmlucDFwGxgPrAkpfTnfKvSetg8pTS/MP0qsHm5d2gI60MiYjBwC/CfKaW38q5HHRMRnwJeTyk9mnctWm/9gN2Bq1JKuwHL6IZTH+pahTFDh5CF6hHABhHx5XyrUldI2aUjyn75CENYHxERNWQB7H9SSrfmXY86ZS/g4IiYBdwIfCwirs+3JHXSK8ArKaXGHumbyUKZepePA/9KKS1IKdUBtwIfyrkmdd5rEbElQOHn6+XeoSGsD4iIIBt78lxK6Sd516POSSlNSiltlVIaQzb4956Ukv/r7oVSSq8CcyJi28Ks/YBncyxJnTMb+EBEDCr8nd0Pv2DRm00HvlqY/irwh3Lv0BDWN+wFHEnWc/L3wmNC3kVJfdx3gP+JiCeBXYEf5luOOqrQk3kz8BjwFNlnqlfP7wUi4gbgQWDbiHglIr4BXADsHxH/JOvlvKDsdXjFfEmSpO5nT5gkSVIODGGSJEk5MIRJkiTlwBAmSZKUA0OYJElSDgxhkvqsiBgVEf+KiI0L7Y0K7TE5lyapDzCESeqzUkpzgKsoXg/oAmBaSmlWbkVJ6jO8TpikPq1wS69HgV8CRwO7Fm5BI0ll1S/vAiQpTymluog4GbgL+IQBTFJ38XSkJMEngfnATnkXIqnvMIRJ6tMiYldgf+ADwPciYst8K5LUVxjCJPVZERFkA/P/M6U0G7gIuDjfqiT1FYYwSX3Z0cDslNLdhfaVwPYRsU+ONUnqI/x2pCRJUg7sCZMkScqBIUySJCkHhjBJkqQcGMIkSZJyYAiTJEnKgSFMkiQpB4YwSZKkHBjCJEmScvD/AbViESkBS7ewAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Train and test\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size = 0.3, random_state = 42)\n",
    "\n",
    "# Logistic Regression\n",
    "#   Note:   It is needed to give both the fit and predict methods 2D arrays. \n",
    "#           Right now x_train and x_test are currently only 1 dimensional.\n",
    "#           What we shoudl do to fix it is: \n",
    "x_train = x_train.reshape(-1, 1)\n",
    "\n",
    "# Create linear regression object\n",
    "model = LogisticRegression()\n",
    "\n",
    "# Train the model using the training sets\n",
    "model.fit(x_train, y_train);\n",
    "\n",
    "# Make predictions using the testing set\n",
    "y_pred = model.predict(x_test.reshape(-1, 1))\n",
    "\n",
    "# The coefficients\n",
    "print('True model: y ~ Bernoulli(μ)')\n",
    "print('True logit(μ) = log(μ/(1-μ)) =', beta_0, '+', beta_1, '* x')\n",
    "print('Regression logit(μ) = log(μ/(1-μ)) = %.2f' % model.intercept_[0], '+ %.2f' % model.coef_[0][0], '* x')\n",
    "\n",
    "# Metrics\n",
    "print('Accuracy: ', accuracy_score(y_test, y_pred))\n",
    "print('Precision: ', precision_score(y_test, y_pred))\n",
    "print('Recall: ', recall_score(y_test, y_pred))\n",
    "print('F1: ', f1_score(y_test, y_pred))\n",
    "\n",
    "# Plot outputs\n",
    "x_test_ord = x_test[x_test.argsort()]\n",
    "\n",
    "plt.figure(figsize=(10,5));\n",
    "plt.scatter(x_test, y_test);\n",
    "plt.plot(x_test_ord, logistic_function(model.intercept_[0] + model.coef_[0][0]* x_test_ord), color='red', linewidth=3);\n",
    "plt.title('Logistic Regression');\n",
    "plt.xlabel('X');\n",
    "plt.ylabel('Y');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metrics:\n",
    "\n",
    "* **Accuracy**: predicts the overall accuracy of the model $\\big( = \\frac{correct\\;predictions}{predictions} \\big)$\n",
    "\n",
    "* **Precision**: percentage of *true* positive among all the *predicted* positive $\\big( = \\frac{true\\;positive}{true\\;positive \\,+\\, false\\;positive} = \\frac{true\\;positive}{total\\;predicted\\;positive}\\big)$\n",
    "\n",
    "* **Recall**: percentage of the *predicted* positive among all the *true* positive $\\big( = \\frac{true\\;positive}{true\\;positive \\,+\\, false\\;negative} = \\frac{true\\;positive}{total\\;actual\\;positive}\\big)$\n",
    "\n",
    "* **F1 Score**: *harmonic mean* between Precision and Recall values $\\big( = 2\\cdot\\frac{precision\\;\\cdot\\;recall}{precision \\,+\\, recall} \\big)$<br>\n",
    "*F1 Score might be a better measure to use if we need to seek a balance between Precision and Recall AND there is an uneven class distribution (large number of Actual Negatives).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A good summary:\n",
    "\n",
    "<img src=\"..\\Miscellaneous\\confusion_matrix.png\" width=\"600\"/>\n",
    "\n",
    "A good note:\n",
    "\n",
    "<img src=\"..\\Miscellaneous\\accuracy_precision.png\" width=\"600\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **3.** Logistic Regression with `statsmodels`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.168052\n",
      "         Iterations 9\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                  350\n",
      "Model:                          Logit   Df Residuals:                      348\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Wed, 15 Jun 2022   Pseudo R-squ.:                  0.7126\n",
      "Time:                        16:12:45   Log-Likelihood:                -58.818\n",
      "converged:                       True   LL-Null:                       -204.64\n",
      "Covariance Type:            nonrobust   LLR p-value:                 2.189e-65\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -6.9265      0.966     -7.172      0.000      -8.819      -5.034\n",
      "x1             2.0528      0.278      7.387      0.000       1.508       2.597\n",
      "==============================================================================\n",
      "\n",
      "Possibly complete quasi-separation: A fraction 0.26 of observations can be\n",
      "perfectly predicted. This might indicate that there is complete\n",
      "quasi-separation. In this case some parameters will not be identified.\n"
     ]
    }
   ],
   "source": [
    "# Train and test\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size = 0.3, random_state = 42)\n",
    "\n",
    "# To have the intercept in the model\n",
    "x_train = sm.add_constant(x_train)\n",
    "x_test  = sm.add_constant(x_test)\n",
    "\n",
    "# Logistic Regression\n",
    "model = sm.Logit(y_train, x_train)\n",
    "model = model.fit()\n",
    "print(model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some statistical tests.\n",
    "\n",
    "* **LogLikelihood Ratio (LLR) test**<br>\n",
    "$\\quad$ `LLR p-value = 2.189e-65` <br>\n",
    "This test has as $H_0$ the hypothesis that the model should be simplified, as $H_1$ the hypothesis that the model is good as it is. It is calculated as $-2\\cdot\\log\\frac{likelihood(simple\\;model)}{likelihood(fitted\\;model)}$. If the p-value of this test is significantly small, then one can conclude that the predictor or set of predictors significantly improve the model's fit. This is analogous to the F-test used in linear regression analysis to assess the significance of prediction.\n",
    "\n",
    "* **Z-test for the intercept and the coefficients**<br>\n",
    "$\\quad$ **`Intercept`** $\\;$ `p-value = P>|t| = 0.00`<br>\n",
    "$\\quad$ **`Coefficient`** $\\;$ `p-value = P>|t| = 0.00`<br>\n",
    "The z-test helps to determine the correlation between the response and the predictor variables. The test is used in regression to test the null hypothesis $H_0$ that the slope or the coefficient is equal to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True model: y ~ Bernoulli(μ)\n",
      "True logit(μ) = log(μ/(1-μ)) = -7.0 + 2.0 * x\n",
      "Regression logit(μ) = log(μ/(1-μ)) = -6.93 + 2.05 * x\n",
      "Accuracy:  0.7066666666666667\n",
      "Precision:  0.7066666666666667\n",
      "Recall:  1.0\n",
      "F1:  0.828125\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFNCAYAAABIc7ibAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAy+ElEQVR4nO3deXwdZb3H8c8vadqmLbSlLZSWQlELCCKiZVH0CldZBGRRQfCioLIoiOxCkYuAKGiRRVkUEFFQEAErKgouIC5wpVABWWXtQkuBtkA3muW5f8xJT5qe7DmZk+Tzfr3OK88zM2fmN3NOk29njZQSkiRJ6l1VeRcgSZI0EBnCJEmScmAIkyRJyoEhTJIkKQeGMEmSpBwYwiRJknJgCJOUi4j4fkT8bxfet3FELI2I6nLUVaki4ncRcWjedUjqOeF9wiS1JyKeBw5PKf2xry47Ig4DfgisABqB54CvppR+090aJakr3BMmaSC5N6U0AhgFXA7cGBGjenohA20vnaSuMYRJ6rKIGBIRF0fEi4XXxRExpNn4r0TE/MK4wyMiRcTbCuOujYhzC+2xEfGbiFgSEYsi4q8RURUR1wEbA78uHIL8SkRMLsxnUOG960XEjwrLWBwRM9qrO6XUCFwHDAemNFuXCyJidkS8VDhcWtuJdbkiIm6PiGXALhExISJuiYiXI+K5iPhys3ltHxEzI+L1wrIuLAwfGhHXR8SrhW1xf0RsUBh3d0QcXmhXRcQZEfFCRCyMiJ9ExMjCuKbtc2hhXV6JiK92+UOWVDaGMEnd8VVgR+BdwDbA9sAZABGxB3Ai8GHgbcDObcznJGAuMA7YADgdSCmlTwOzgY+mlEaklL5d4r3XAcOArYD1gYvaK7qwp+qzQB3wQmHw+cBmhXV5GzAROLMT6/Ip4BvAOsA/gF8DDxXm8yHg+IjYvTDtJcAlKaV1gbcCNxWGHwqMBCYBY4AvkB0+bemwwmsX4C3ACODSFtO8H9i8sOwzI+LtrW8RSXkwhEnqjv8BzkkpLUwpvQycDXy6MO5A4EcppUdTSsuBs9qYTx2wIbBJSqkupfTX1IETViNiQ+AjwBdSSosL7/1LG2/ZMSKWACuBC4BDUkoLIyKAI4ETUkqLUkpvAN8EDurEuvwqpfT3wl62rYFxKaVzUkqrUkrPAlc1m18d8LaIGJtSWppSuq/Z8DHA21JKDSmlB1JKr5dY1v8AF6aUnk0pLQWmAQc17R0sODultCKl9BBZGNymje0iKQeGMEndMYHiniQK7QnNxs1pNq55u6XpwNPAnRHxbESc1sHlTwIWpZQWd3D6+1JKo4DRwG3ABwrDx5HtTXugcBhwCfD7wnDo2Lo0H7YJMKFpXoX5nU62lw/g82R73Z4oHHLcuzD8OuAOsnPVXoyIb0dETYllldrug5rNH2BBs/Zysr1lkiqIIUxSd7xIFjiabFwYBjAf2KjZuEmtzSSl9EZK6aSU0luAfYATI+JDTaPbWP4cYL3Onlxf2Hv0ReDTEbEt8ArZYb+tUkqjCq+RhZP4O7ouzeucAzzXbF6jUkrrpJT2LCz/Pymlg8kOn34LuDkihhf25J2dUtoSeB+wN/CZEssqtd3rgZc6sx0k5csQJqmjagonjje9BgE3AGdExLiIGEt2DtX1helvAj4bEW+PiGFAq/cEi4i9I+JthcOCrwENZLeRgCxYvKXU+1JK84HfAZdHxOiIqImI/+rIyqSUFgFXA2cWDiFeBVwUEesXaprY7ByuDq9LwT+BNyLi1IiojYjqiHhHRGxXmPchETGusNwlhfc0RsQuEbF14Zy118kOTzaWmP8NwAkRsWlEjCA7dPrzlFJ9R9ZdUmUwhEnqqNvJ9hY1vc4CzgVmAg8DjwAPFoaRUvod8F3gLrJDjU3nPb1ZYt5TgD8CS4F7gctTSncVxp1HFvSWRMTJJd77abKw8gSwEDi+E+t0MbBnRLwTOLWpzoh4vVDP5l1YF1JKDWR7sd5Fdj+yV8gC38jCJHsAj0bEUrKT9A9KKa0AxgM3kwWwx4G/kB2ibOmawvB7CvNfCRzbifWWVAG8WaukXlG4Ou/fwJC+vsemP62LpPy4J0xS2UTE/oX7b40mO/fp1301tPSndZFUGQxhksrpKLJDhM+Qnef1xXzL6Zb+tC6SKoCHIyVJknLgnjBJkqQcGMIkSZJyMKj9SSrL2LFj0+TJk/MuQ5IkqV0PPPDAKymlcaXG9bkQNnnyZGbOnJl3GZIkSe2KiBdaG+fhSEmSpBwYwiRJknJgCJMkScqBIUySJCkHhjBJkqQcGMIkSZJyYAiTJEnKQdnuExYR1wB7AwtTSu8oMT6AS4A9geXAYSmlB8tVT15mzJrH9Due5MUlK5gwqpZTdt+c/badmHdZvaY31z+vbV1quUCHazljxiPc8H9zaEiJAAYPquLN+kYARtXWcNY+W7HfthPXWM7I2hrqGhpZtqph9XxKTTtvyYo1llUV8KkdNubc/bYuufzqCA7eYdIa45vPqzqChpSYWGKdOjpdW9vxrNseZcmKOgCGD66mprqK11bUMWFULZPH1PKPZxfR9LjbwdXBsMGDVo/vyHZvXmMEq+c1elgNe71zQ+564uV2P0eAs3/9KIuX163e7ntvs/Z7W65zqe0zelgNb9Y1sLxu7c+7J8yYNW+tWkt9nyaMqmWXLcZx1xMvr/WdaTJ6WA0r6xpYUai15bgtN1yHfzyziOZPIz5kxzW/a+1p/l0ECKC1pxu3td1nzJrHtFsfLlnrsJoqhtRUs2R569+b5tui5Xf5FzNn8/dnFq2eX00V1Deyxvs6+juo1He+MaXVdY8eVsPXPtq170Nbn31X5tVy+/zmofmr626p5bK6+ruh5e+8CFZ/bpPH1HLfs4tXf1eAdn8/Nd8ewOrfAZ35PdXTyvYA74j4L2Ap8JNWQtiewLFkIWwH4JKU0g7tzXfq1Kmpr9ysNftF8Agr6op/KGtrqjnvY1sPiCDWm+uf17YutdyaqoCAuobiv63WajljxiNcf9/sNpdRUxV8cvtJ3PLAvDWW051pm/44trb8pvGl1q/UOnV0utbMmDWPU37xEHWNXf991N52b6vGkvOrDkisUVNNVdAINLRTZ8t17syya6qC6Qds0+3v7YxZ8zjl5ofW2B5N8+/o96knrA5iTX9rUirZ/tqMR7jh/7LvYlAcHiSiMHnWLga0KES0pmG1g6rYZ5sNuWXmHBoLw6PZ+ObzojCvmuogUqK+odk0q+dbXEZXDBlUzal7bM5uW40vDizUcuejCzjv9sfb/c7XVAfTPvL2NefRXIm/4Xc+uoDzf/c4LTNoTcBpe7YxrxLzvPPRBVxwx5OsrF870LZZd1Vw6ke2AOA7dzzJyvq1v2tDB1Vz8m6bsWuJev7w6AIuuPOpku+LDuSWoTXVnLTbZuy65Xj+8NgCvvX7J1Z/xi3NX3ccVbW1Zft7EREPpJSmlhxXrhBWWPBk4DethLAfAHenlG4o9J8Edk4pzW9rnn0phO10/p9L/q9y4qha/n7af+dQUe/qzfXPa1u3ttxSStXy1mm3r/E/udY0/Q+vIzoybXUEz5y3Z6vLbxrf3vo1rVNHp2tNZ7ZjZ3W0xnItF1pfv+rGBoavWsHwVSsYWr+KQQ31DG6sZ0JtNVd9ahtYtSp71dW1327R/+U/n2fV8pXUNNZT01Bf/FlYxqCGegY11lPT0EBVaiy8ElUpEamR6kI/Cj+rC9NEYZqq1EhVYyNVpOJ7GwvTkP2hrOpGgJF6y76f/g4PTdi8bH8v2gpheT62aCIwp1l/bmHYWiEsIo4EjgTYeOONe6W4nvBiK7/wWxve3/Tm+ue1rTsz/1LTdjRYdXS6jk7bNE1r0zYNb2/9msZ3dLquju+OjtbYGYMa6hm18g1GrXiD9Va8zugVrzNyxVLWWbWcYatWMLxuJcNXrYDHroWlSzl/1vMMX7WcYatWrh43fNUKhjSUPpwDwBXdq3H/7r1dGnDy+NvcJ54dmVK6ErgSsj1hOZfTYRNG1Zb83++EUbU5VNP7enP989rWrS23tWlb6ugernLsCWtr2qbx7a1f0zp1dLq2xpdrL1VHaxzUUM+4ZYsZ/8arbLD0VTZYuojxb7zK+ktfZczy1xm18nXWW/46o1a8wbqrlnds4bOyHx/o7kr0NxHZq1m7rjFlBwJj9USkwjTNh6c1hmf9pvk0fZNTxOrhqdTwdua15vuzWrqquirYcOTQtdb/xddWtntYu/k8JrScR4v5NTdvyUoaGksfPqyuqmLiqDbm1WKe85asoL6LpwkMqsrmUV/4bFubZqPRa/9+mLu47eWmDnwkg6qqmDS6ljmLV1DfyvYAWDWoBsjnb3OeIWweMKlZf6PCsH7jlN03L3meUtOJoP1db65/Xtu61HJbOzepVC0H7zApl3PCDt5hUpvLbxpfav1KrVNHp2vNKbtvXrZzwlbXuNtmXPiTvzBh4Ww2XfQimy6ax6aLX2T8G68wfumrjFn2Wq8fPmuIKpYNrmVZzVBW1gymrqqG+kGDGD9uXdYbORwGD85eNTWdbv/rpeX8/KEFvBnV1FUNYlV1DfXV1TQOquH9W03gz08vZmmqor6qmsaqKhqjioaoIkXQGFU0RhT6Wbtx9fCq1e1UmKZpWPNpIQsz/7PjJtk5YdH2X82zO3B+ZHtqa6r5+Hsm8rP7ZtPRM5hKnfvXmrYuFGitntbOM/pnB8+DrKkOpn+ic+cI3t/G+YDTD9iGiZ2cV2fOpWy5LKDd80U3KlHPzC4ut+W8J207kQda2R4tp8/jb3OeIew24EsRcSPZifmvtXc+WF/T9I9moF4d2Zvrn9e2bm25Ha2l6cqxjlwdOXWT9Tp8dWTTtO1dHdly+S2vjmy+fm1d2dTR6drbjj1xdeTlv3qA4c8+zbYrX+bAUSvY4lvXwlNPsd9TT7HfsmVt1tFRjREsGboOi2vXYemwdRm84fo81zCEl2IIse66bPeOSWy12URYZx0YMQJGjOCe+Su5etZCnnuzipVDhvH6oCEMW3c4b9Y3rnV15JY98L19F/B8G1fIjamwqyNbfheh61dHTt1kvYq/OrK173xPXB3ZNH1PXB1Z6ndcZ6+ObHp/Z343tFxud66OLLU9oP9fHXkDsDMwFngJ+BpQA5BS+n7hFhWXAnuQ3aLisymlds+470sn5ksqs4UL4cEHi69Zs+DZZ7s2rwhYf32YODF7TZhQ/LnBBjBmDIwdm/0cNQqqvM2ipPblcmJ+SungdsYn4JhyLV9SP5MSPPEE3H03/O1vcN99nQ9co0bBZpsVX1OmwOTJWdgaPz47pCdJvaRPnJgvaYB65hn44x+z4HX33bBgQfvvqamBzTfPXs0D12abZXux2jk3SZJ6iyFMUuWoq4N77oHf/jZ7PfVU29MPGQLbbgvvfnfxteWW2XBJqnCGMEn5qq+Hu+6Cm26CW2+FRYtan3b0aPjgB7PX+94H73pXdjWgJPVBhjBJ+bj/frj66ix4vfJK6WmGDYNddoEPfQh23hm22cYT4iX1G4YwSb2nrg5uuQW++124997S00yaBPvtB3vtle3xGtqBG0tKUh9kCJNUfgsXwpVXwhVXwIsvrj1+4kQ44AA48EDYYQf3dkkaEAxhksrnwQezvV433JA9WLq5wYPhoIPgiCOy87sMXpIGGEOYpJ7V0JCd5/Xd72b382pp/Hj44hfhqKOym6BK0gBlCJPUM1KCX/wCvva17KaqLe2wA3z5y/CJT3hFoyRhCJPUEx5+GI45Zu09XzU12Xlexx6bhTBJ0mqGMEld9/rrcNZZ2aHHhuLDxFl33Wyv19FHw4Yb5laeJFUyQ5ikzksJfv5zOPFEmD+/OHzQIDj+eJg2DdZbL7fyJKkvMIRJ6pwnnoAvfQn+9Kc1h++yC1x2Gbz97fnUJUl9jNeES+qY5cvh9NPhne9cM4CNHw8//Wk2zAAmSR3mnjBJ7XvgAfif/4EnnywOq6rKTrg/+2wYOTK/2iSpj3JPmKTWpQQXXww77rhmAHvf+7IbsV58sQFMkrrIPWGSSluxIruh6nXXFYeNGAEXXgif/7x3uJekbjKESVrb3Lmw//4wc2Zx2Pbbw89+Bm99a351SVI/4n9lJa3pH/+AqVPXDGCf/zzcc48BTJJ6kCFMUtFVV8HOO8NLL2X9QYPg0kuz4UOG5FqaJPU3Ho6UBKtWZTdZveKK4rCxY7NnQe68c15VSVK/ZgiTBrrXX4d994W77y4Oe9e7YMYM2GSTnIqSpP7Pw5HSQLZoEXzoQ2sGsE9+Ev7+dwOYJJWZIUwaqF56KTvU2PwE/G98A264AYYNy60sSRooPBwpDURz5mR7wP7zn6wfAd//Phx5ZL51SdIAYgiTBprnnssetv3CC1m/uhquvRYOOSTXsiRpoDGESQPJvHnZHrCmAFZTAz//eXZjVklSrzKESQPFK6/Arrtme8Igu+/XjBmwxx65liVJA5Un5ksDwWuvwe67w+OPZ/1Bg+Dmmw1gkpQjQ5jU361cCR/9KDz4YNaPgOuvh733zrcuSRrgDGFSf9bYCJ/9LPz1r8VhV12V3QtMkpQrQ5jUn51xBtx4Y7E/fXr2MG5JUu4MYVJ/ddVVcN55xf4xx8BJJ+VXjyRpDYYwqT/661/hi18s9vfeGy6+ODsfTJJUEQxhUn+zYAEceCA0NGT9bbfNHkU0yDvSSFIlMYRJ/Ul9fXbS/YIFWX/sWPjVr2DEiHzrkiStxRAm9Sennw733JO1I7I9YJMm5VuTJKkkQ5jUX9x6a3b1Y5Ovfx0+/OH86pEktckQJvUHTz0Fhx1W7O+1F0yblls5kqT2GcKkvm7ZMvj4x+GNN7L+5Mlw3XVQ5T9vSapkZf0tHRF7RMSTEfF0RJxWYvzGEXFXRMyKiIcjYs9y1iP1OynBF74A//531h8yBG65BUaPzrcuSVK7yhbCIqIauAz4CLAlcHBEbNlisjOAm1JK2wIHAZeXqx6pX7rmmuw5kE0uuwze/e786pEkdVg594RtDzydUno2pbQKuBHYt8U0CVi30B4JvFjGeqT+5bnn4Pjji/3Pfc5HEklSH1LOuzdOBOY0688FdmgxzVnAnRFxLDAc8FIuqSOaHsy9dGnW32wz+N738q1JktQpeZ+5ezBwbUppI2BP4LqIWKumiDgyImZGxMyXX36514uUKs53vwt/+UvWrqqCn/wEhg3LtyZJUqeUM4TNA5rfJXKjwrDmPg/cBJBSuhcYCoxtOaOU0pUppakppanjxo0rU7lSH/HEE2vefmLaNNih5U5mSVKlK2cIux+YEhGbRsRgshPvb2sxzWzgQwAR8XayEOauLqk19fVw6KGwcmXW32YbOPPMfGuSJHVJ2UJYSqke+BJwB/A42VWQj0bEORGxT2Gyk4AjIuIh4AbgsJRSKldNUp/3rW/BP/+ZtWtqssOQgwfnW5MkqUvKeWI+KaXbgdtbDDuzWfsxYKdy1iD1G//6F5x9drF/9tnwznfmVo4kqXvyPjFfUkc0NGS3n6iry/o77ginnJJvTZKkbjGESX3B978PDz6YtYcOhWuvhUFl3ZEtSSozQ5hU6RYsgNNPL/ZPPx023zy/eiRJPcIQJlW6k0+G11/P2lOmwFe+km89kqQeYQiTKtldd8FPf1rsX3ZZ9pBuSVKfZwiTKtWqVXD00cX+Jz8Ju+6aXz2SpB5lCJMq1Xe+k90dH2CddeDCC/OtR5LUowxhUiV6/nn4+teL/XPOgQkTcitHktTzDGFSJTruOFixImtvsw186Uv51iNJ6nGGMKnS/PrXcFuzx6xecYX3BJOkfsgQJlWSVavgxBOL/cMPh/e+N796JEllYwiTKsnll8PTT2ftUaPg/PNzLUeSVD6GMKlSLFqUnYDf5MwzYcyY/OqRJJWVIUyqFOecA4sXZ+23vQ2OOSbfeiRJZWUIkyrBU09ld8Nv8q1vweDB+dUjSSo7Q5hUCU49Ferrs/YHPgD7759vPZKksjOESXn7y19gxoxi/8ILISK3ciRJvcMQJuUpJTj99GL/kENg6tT86pEk9RpDmJSnO++Ef/wja9fUwLnn5luPJKnXGMKkvKQE//u/xf4RR8Amm+RXjySpVxnCpLz89rdw//1Ze8iQNQ9LSpL6PUOYlIeUspuxNvnCF2DixPzqkST1OkOYlIcZM2DWrKxdWwunnZZrOZKk3mcIk3pbYyN87WvF/tFHw/jx+dUjScqFIUzqbbfcAo88krWHD89u1CpJGnAMYVJvamiAs84q9o89FsaNy60cSVJ+DGFSb/r5z+Gxx7L2OuvAySfnW48kKTeGMKm31NfD2WcX+8cfD2PG5FaOJClfhjCpt/zsZ/DUU1l75Eg44YR865Ek5coQJvWGujo455xi/8QTYfTo/OqRJOXOECb1hp/8BJ55JmuPHp0dipQkDWiGMKncVq2Cr3+92D/lFFh33fzqkSRVBEOYVG4/+hG88ELWHjs2uy2FJGnAM4RJ5bRyJZx7brF/6qkwYkR+9UiSKoYhTCqnq6+GuXOz9gYbZI8okiQJQ5hUPitWwDe/WeyfdhoMG5ZfPZKkimIIk8rlBz+A+fOz9oQJcNRR+dYjSaoohjCpHN58E6ZPL/anTYPa2vzqkSRVHEOYVA4//Sm8+GLWHj8eDj8833okSRXHECb1tIYG+Pa3i/0TToChQ/OrR5JUkcoawiJij4h4MiKejojTWpnmwIh4LCIejYiflbMeqVfMmAFPPpm1R46EL3wh13IkSZVpULlmHBHVwGXArsBc4P6IuC2l9FizaaYA04CdUkqLI2L9ctUj9YqU4Pzzi/2jj/bu+JKkksq5J2x74OmU0rMppVXAjcC+LaY5ArgspbQYIKW0sIz1SOX35z/DzJlZe+hQOO64fOuRJFWscoawicCcZv25hWHNbQZsFhF/j4j7ImKPMtYjlV/zvWCf+1x2g1ZJkkoo2+HITix/CrAzsBFwT0RsnVJa0nyiiDgSOBJg44037uUSpQ6aORP++MesXV0NJ5+cbz2SpIpWzj1h84BJzfobFYY1Nxe4LaVUl1J6DniKLJStIaV0ZUppakpp6rhx48pWsNQtzfeCffKTsOmm+dUiSap45Qxh9wNTImLTiBgMHATc1mKaGWR7wYiIsWSHJ58tY01SeTz5JNx6a7F/WsmLgSVJWq1sISylVA98CbgDeBy4KaX0aEScExH7FCa7A3g1Ih4D7gJOSSm9Wq6apLKZPj27MhJgr71g663zrUeSVPEiNf3h6COmTp2aZjZdfSZVgrlz4S1vgbq6rP/Xv8L7359vTZKkihARD6SUppYa5x3zpe666KJiANtpJwOYJKlDDGFSdyxeDD/4QbE/bVp+tUiS+hRDmNQdV14Jy5Zl7a22gj33zLceSVKfYQiTumrVKrjkkmL/5JMhIr96JEl9iiFM6qobboD587P2hhvCwQfnW48kqU8xhEldkRJccEGx/+Uvw5Ah+dUjSepzDGFSV9x5J/z731l7+HA46qh865Ek9TmGMKkrmu8FO/xwGD06v1okSX2SIUzqrH/9q/ig7qoqOP74PKuRJPVRhjCps77znWL7gANg8uTcSpEk9V2GMKkz5syBG28s9k8+Ob9aJEl9miFM6ozvfhfq67P2Bz8IU0s+DkySpHYZwqSOeu21NR9R5F4wSVI3GMKkjrr6anjjjay9xRY+okiS1C2GMKkj6urg4ouL/ZNOyq6MlCSpi/wrInXETTfB3LlZe/314ZBD8q1HktTnGcKk9qQE06cX+8ceC0OH5lePJKlfaDWERcTtETG5F2uRKtOf/gQPPZS1a2vhi1/Mtx5JUr/Q1p6wHwF3RsRXI6KmtwqSKk7zRxR97nMwZkx+tUiS+o1BrY1IKf0iIn4H/C8wMyKuAxqbjb+wF+qT8vXww3DHHVk7Ak44Id96JEn9RqshrGAVsAwYAqxDsxAmDQgXNvu/xsc+Bm99a361SJL6lVZDWETsAVwI3Aa8O6W0vNeqkirBvHnws58V+96cVZLUg9raE/ZV4ICU0qO9VYxUUb73vez+YAA77QQ77phvPZKkfqWtc8I+0JuFSBXljTfg+98v9t0LJknqYd4nTCrlhz/MnhUJMGUKfPSj+dYjSep3DGFSS/X1cNFFxf6JJ0J1dX71SJL6JUOY1NLNN8Ps2Vl77Fj4zGfyrUeS1C8ZwqTmWj6i6JhjYNiw/OqRJPVbhjCpub/8BR58MGsPHQpHH51vPZKkfssQJjXX/BFFhx4K66+fXy2SpH7NECY1eewx+O1vs7aPKJIklZkhTGrS/BFF++wDm2+eXy2SpH7PECYBLFgA111X7HtzVklSmRnCJIBLL4VVq7L2DjtkjymSJKmMDGHSsmVw+eXF/sknZ+eESZJURoYw6Uc/gsWLs/amm8L+++dbjyRpQDCEaWBraFjzhHwfUSRJ6iWGMA1st94Kzz2XtUePhs9+Nt96JEkDhiFMA1dKcN55xf7RR8Pw4fnVI0kaUAxhGrh+/3uYNStr19bCl7+cbz2SpAHFEKaB65vfLLaPOMJHFEmSelVZQ1hE7BERT0bE0xFxWhvTfTwiUkRMLWc90mp//Sv87W9Zu6bGm7NKknpd2UJYRFQDlwEfAbYEDo6ILUtMtw5wHPB/5apFWss3vlFsf+YzMGlSfrVIkgakcu4J2x54OqX0bEppFXAjsG+J6b4OfAtYWcZapKIHHoA77sjaVVVw6qn51iNJGpDKGcImAnOa9ecWhq0WEe8GJqWUflvGOqQ1Nb8i8oADYMqU/GqRJA1YuZ2YHxFVwIXASR2Y9siImBkRM19++eXyF6f+6/HHs3uDNZk2Lb9aJEkDWjlD2Dyg+Yk2GxWGNVkHeAdwd0Q8D+wI3Fbq5PyU0pUppakppanjxo0rY8nq9771rez+YAB77w3bbJNvPZKkAaucIex+YEpEbBoRg4GDgNuaRqaUXkspjU0pTU4pTQbuA/ZJKc0sY00ayJ5/Hq6/vtg//fTcSpEkqWwhLKVUD3wJuAN4HLgppfRoRJwTEfuUa7lSq6ZPz54VCbDzzvDe9+ZajiRpYIvUdGimj5g6dWqaOdOdZeqkBQtg8mR4882sf+edsOuuuZYkSer/IuKBlFLJ+6B6x3wNDBddVAxg220HH/5wvvVIkgY8Q5j6v8WL4fLLi/3TT4eI/OqRJAlDmAaC730Pli7N2ltuCft4SqIkKX+GMPVvS5fCJZcU+9OmZXfJlyQpZ/41Uv925ZWwaFHW3nRTOOigfOuRJKnAEKb+68034YILiv1TT4VBg/KrR5KkZgxh6r+uvRbmz8/aG24Ihx6aazmSJDVnCFP/VF+fPaKoyUknwdCh+dUjSVILhjD1Tz//OTz3XNZebz046qh865EkqQVDmPqf+no455xi/7jjYMSI/OqRJKkEQ5j6n5/9DJ56KmuPHAnHHptvPZIklWAIU/9SVwdnn13sn3gijB6dXz2SJLXCEKb+5cc/hmefzdqjR8Pxx+dajiRJrTGEqf9YtQq+/vVi/5RTYN1186tHkqQ2GMLUf1x9NcyenbXHjvVcMElSRTOEqX9YvBjOPLPYP/VUr4iUJFU0Q5j6h3POgVdfzdqbbALHHJNvPZIktcMQpr7viSfg0kuL/enTobY2v3okSeoAQ5j6vhNPzG7QCvBf/wWf+ES+9UiS1AGGMPVtt98Ov/td1o6Aiy/OfkqSVOEMYeq76uqyvWBNDj8ctt02v3okSeoEQ5j6rssugyefzNrrrgvnnptvPZIkdYIhTH3Tyy/DWWcV+2eeCeuvn1s5kiR1liFMfdOZZ8Jrr2XtKVO8Maskqc8xhKnvefhhuPLKYv/CC2Hw4PzqkSSpCwxh6ltSyh7K3diY9XfbDfbaK9eSJEnqCkOY+pZf/hLuuitrV1fDRRd5SwpJUp9kCFPfsXIlnHxysX/00bDllvnVI0lSNxjC1HdcdBE891zWXm+9Na+OlCSpjzGEqW948UX4xjeK/XPOyYKYJEl9lCFMfcPpp8OyZVl7q63gqKPyrUeSpG4yhKny3X8//PjHxf7FF8OgQbmVI0lSTzCEqbKlBMcdV+zvuy98+MP51SNJUg8xhKmy3XAD3Htv1q6pgQsuyLceSZJ6iCFMlWvZMvjKV4r9E06At70tv3okSepBhjBVrm9/G+bNy9obbABf/Wq+9UiS1IMMYapML7yQhbAm3/wmrLtufvVIktTDDGGqTF/5SnaHfIB3vxsOOyzXciRJ6mmGMFWeX/4Sbrqp2L/kEqjyqypJ6l/8y6bKsnDhmjdi/fSn4f3vz68eSZLKpKwhLCL2iIgnI+LpiDitxPgTI+KxiHg4Iv4UEZuUsx5VuJSyAPbyy1l/4sRsL5gkSf1Q2UJYRFQDlwEfAbYEDo6ILVtMNguYmlJ6J3Az8G00cF13HcyYUexfcw2MHp1bOZIklVM594RtDzydUno2pbQKuBHYt/kEKaW7UkrLC937gI3KWI8q2Zw5cOyxxf4Xvwi77ZZfPZIklVk5Q9hEYE6z/tzCsNZ8HvhdGetRpWpshM99Dl5/Peu/9a0wfXq+NUmSVGYV8RTkiDgEmAp8sJXxRwJHAmy88ca9WJl6xeWXwx//mLUjsod1Dx+eb02SJJVZOfeEzQMmNetvVBi2hoj4MPBVYJ+U0pulZpRSujKlNDWlNHXcuHFlKVY5uf9+OOmkYv8rX4GddsqvHkmSekk5Q9j9wJSI2DQiBgMHAbc1nyAitgV+QBbAFpaxFlWiV16BT3wCVq3K+ttsA2efnW9NkiT1krKFsJRSPfAl4A7gceCmlNKjEXFOROxTmGw6MAL4RUT8KyJua2V26m8aGuBTn4LZs7P+yJFwyy0wZEi+dUmS1EvKek5YSul24PYWw85s1v5wOZevCnbWWfCHPxT711+fnZAvSdIA4R3z1ft+8xs499xi/4wzYO+986tHkqQcGMLUu555JnsUUZPddsv2ikmSNMAYwtR7li+Hj38clizJ+htvDD/9KVRX51qWJEl5MISpd6QERx8NDz2U9QcPhptvhrFj861LkqScGMLUO773vewmrM37222XXz2SJOXMEKbyu/FGOP74Yv+ww+CII/KqRpKkimAIU3ndeSd85jPZ4UiAHXbIHlMUkW9dkiTlzBCm8vnnP+FjH4O6uqz/9rfDb38LtbX51iVJUgUwhKk8Hn8c9twTli3L+pMmwR13wJgx+dYlSVKFMISp582ZA7vvDq++mvXHjMkOS06a1Pb7JEkaQAxh6lmvvpoFsDlzsv7w4XD77bDFFvnWJUlShTGEqecsXAh77JEdigSoqYFf/hK23z7fuiRJqkBlfYC3BpAnn8zOAXv22awfkT2Ue9dd861LkqQK5Z4wdd/f/w7ve18xgFVVwQ9+AAcemG9dkiRVMPeEqXt+8Yvsgdxvvpn1hw3Lbs760Y/mW5ckSRXOPWHqmpTgO9/J9nY1BbD114e77zaASZLUAe4JU+c1NGSPIbr00uKwzTfProJ8y1tyK0uSpL7EEKbOWb4cPvUp+NWvisPe/36YMcMbsUqS1AkejlTHLVwIu+yyZgA78ED4wx8MYJIkdZIhTB3z5JPw3vdmz4NscvLJcMMNMHRofnVJktRHGcLUtpTgmmtgxx3XvAXFpZfC9OlZW5IkdZrnhKl1zz4LRx4Jf/pTcVhtbbb3a99986tLkqR+wN0YWltDA1x4IbzjHWsGsE03zW5BYQCTJKnbDGFa08MPZ+d+nXQSrFiRDauqghNPhEce8TmQkiT1EA9HKrNyJXzjG3D++VBfXxy+9dbwwx/CdtvlV5skSf2QIUzwt7/BEUfAE08Uhw0eDGeeCaeckrUlSVKPMoQNZG+8AdOmwWWXrTl8p53g6qthiy3yqUuSpAHAc8IGolWrsltMTJmyZgAbMSLr33OPAUySpDJzT9hAsmQJXHklXHIJvPjimuP22guuuAImTcqlNEmSBhpD2EDwwgtZ8LrqKli6dM1xG20E3/42HHQQRORTnyRJA5AhrD978EG44AK46abs3l/NjR+fnXR/9NE+dkiSpBwYwvqblSuzB2z/4Adw111rj99yy+yZj5/6FAwZ0vv1SZIkwBDWP6QE990HP/kJ3Hhjdu5XS//931n42n13n/coSVIFMIT1VY2NcO+9cPPNcMstMGfO2tNUV8OBB2Z3v3/Pe3q/RkmS1CpDWF+yZEn2LMc77oDf/Abmzy893VveAocdBoceChtv3JsVSpKkDjKEVbLly2HWLPjzn7Pgdd99a59g32T0aNh//yx8vf/9XukoSVKFM4RVivp6ePRR+Oc/4f77s5///nfroQtg3LgseH3847DLLlBT03v1SpKkbjGE5SEleOaZYuC6//7sdhIrVrT9vojs3K7dd89e730vDPIjlCSpL/IveDnV1WU3Sn36afjPf+Dxx7PXQw/B4sXtvz8CNt8cdtwRdt01e40bV/66JUlS2RnCumPZMpg3D+bOLf6cPTvby/XMM1kAa+twYkuTJsH228N222Wv97wHRo4sX/2SJCk3hrCWli+HhQvh5ZeLr+b9BQuysDV3bun7cXXUeuutGbi22y67i70kSRoQyhrCImIP4BKgGrg6pXR+i/FDgJ8A7wFeBT6ZUnq+nDW15Q+/+hu77veBHp3ni+uMZfao8Tw/egKzx01izvhN+NeICTRuMplT9tiC/baduNZ7Zsyax/Q7nmTekrXPERtVW0MELF5eR3UEDSkRQCqMHz2shr3euSF3PfEyLy5ZwYRRtZyy++arl3PGjEe44f/m0JDS6nlWR3DwDpM4d7+tS65D8/e0NW1T3aWWW2ocsHo9m9Zl4qhadtliXKv1d8SMWfM4+9ePsnh53eptdtY+W5WsZWRhey5ZXtelZbVXR2vbo6fm1ZVl9GRd3VEpdUhSXiI1+2PcozOOqAaeAnYF5gL3AwenlB5rNs3RwDtTSl+IiIOA/VNKn2xrvlOnTk0zZ87s8XpnzJrHOdf/nQcvbHPxa3izehAvjRjDgnXGsGCdscxfZywvjRhTCF0bMmfkBrxZ0/qjgWprqjnvY1uv8Ydnxqx5TLv1EVbUdeIwZjualjPzhUVcf9/sVqc7ZMeN1wpXZ8x4pOR7Wk5bqu6m5QJrjaupCgioa2j/+1dqO7Vmxqx5nHLzQ2vNt6YqmH7ANiVr6eqy2qujte3R2Xm3Nq+Pv2citzwwr1PL6Mm6uqNS6pCkcouIB1JKU0uOK2MIey9wVkpp90J/GkBK6bxm09xRmObeiBgELADGpTaKKlcI2+n8PzNv8XL+ffGBvD5kOIuGjWRR7bq8Omwkrw4byaLCz1eGj2LBOmNZMGIMi4atS4ruPQJo4qha/n7af69ZR4k9YN01cVQtC15bucYesJaqI3jmvD3XGPbWabeXfE/LaVure+KoWoBur1PL7dSatrZfR2vp6LK6UkdX5t3avJr2HnZmGT1ZV3dUSh2SVG5thbByHo6cCDR/ls5cYIfWpkkp1UfEa8AY4JXmE0XEkcCRABuX6Q7wLy5ZARG84/ibevVGpy+2+EPUst+Ty2kvbpf6g95aaGs5vLW6e2p9OjqftqbriXl0VE9uj9be09pn05VtUK7vXaXXIUl56hNPck4pXZlSmppSmjquTLdomFDYS9Lbd5pfvdxW+j25nOp21q3U+Nbe03J4a3VPGFXbI+vU0Xm0NV1HaylnvV2Zd2vvae2zaW8b9FRd3VEpdUhSnsoZwuYBk5r1NyoMKzlN4XDkSLIT9HvdKbtvTk117waw2prq1SeoN6+jtqa6LMs5eIdJbU5Xanxr72k5vFTdTcstNa6mKjq8vUttp9a09jnWVEWrtXR1We3V0dr26Kl5HbzDpE4voyfr6o5KqUOS8lTOw5H3A1MiYlOysHUQ8KkW09wGHArcC3wC+HNb54OVU9PJwM2vqqutqaIqgmWrspOHm65CHFZTxfK6xjXeX1MFTYOarlD87cPzV8+raX5Da6rbvBKvqV+OqyOb5t2ZqyObhrV3dWTzulu72q03ro4s9Tm2vDqyeS3lujqyI9ujJ+Y1dZP1OrWMnqyrOyqlDknKU9lOzAeIiD2Bi8luUXFNSukbEXEOMDOldFtEDAWuA7YFFgEHpZSebWue5ToxX5IkqafldWI+KaXbgdtbDDuzWXslcEA5a5AkSapEfeLEfEmSpP7GECZJkpQDQ5gkSVIODGGSJEk5MIRJkiTlwBAmSZKUA0OYJElSDsp6s9ZyiIiXgRfyrqMPG0uLB6SrT/Jz7B/8HPsHP8f+oVyf4yYppZIPvu5zIUzdExEzW7tzr/oOP8f+wc+xf/Bz7B/y+Bw9HClJkpQDQ5gkSVIODGEDz5V5F6Ae4efYP/g59g9+jv1Dr3+OnhMmSZKUA/eESZIk5cAQNgBExKSIuCsiHouIRyPiuLxrUtdFRHVEzIqI3+Rdi7ouIkZFxM0R8UREPB4R7827JnVeRJxQ+L3674i4ISKG5l2T2hcR10TEwoj4d7Nh60XEHyLiP4Wfo8tdhyFsYKgHTkopbQnsCBwTEVvmXJO67jjg8byLULddAvw+pbQFsA1+pn1OREwEvgxMTSm9A6gGDsq3KnXQtcAeLYadBvwppTQF+FOhX1aGsAEgpTQ/pfRgof0G2S/7iflWpa6IiI2AvYCr865FXRcRI4H/An4IkFJalVJakmtR6qpBQG1EDAKGAS/mXI86IKV0D7CoxeB9gR8X2j8G9it3HYawASYiJgPbAv+XcynqmouBrwCNOdeh7tkUeBn4UeHQ8tURMTzvotQ5KaV5wAXAbGA+8FpK6c58q1I3bJBSml9oLwA2KPcCDWEDSESMAG4Bjk8pvZ53PeqciNgbWJhSeiDvWtRtg4B3A1eklLYFltELhz7UswrnDO1LFqonAMMj4pB8q1JPSNmtI8p++whD2AARETVkAeynKaVb865HXbITsE9EPA/cCPx3RFyfb0nqornA3JRS0x7pm8lCmfqWDwPPpZReTinVAbcC78u5JnXdSxGxIUDh58JyL9AQNgBERJCde/J4SunCvOtR16SUpqWUNkopTSY7+ffPKSX/190HpZQWAHMiYvPCoA8Bj+VYkrpmNrBjRAwr/J79EF5g0ZfdBhxaaB8K/KrcCzSEDQw7AZ8m23Pyr8Jrz7yLkga4Y4GfRsTDwLuAb+ZbjjqrsCfzZuBB4BGyv6nePb8PiIgbgHuBzSNibkR8Hjgf2DUi/kO2l/P8stfhHfMlSZJ6n3vCJEmScmAIkyRJyoEhTJIkKQeGMEmSpBwYwiRJknJgCJM0YEXEpIh4LiLWK/RHF/qTcy5N0gBgCJM0YKWU5gBXULwf0PnAlSml53MrStKA4X3CJA1ohUd6PQBcAxwBvKvwCBpJKqtBeRcgSXlKKdVFxCnA74HdDGCSeouHIyUJPgLMB96RdyGSBg5DmKQBLSLeBewK7AicEBEb5luRpIHCECZpwIqIIDsx//iU0mxgOnBBvlVJGigMYZIGsiOA2SmlPxT6lwNvj4gP5liTpAHCqyMlSZJy4J4wSZKkHBjCJEmScmAIkyRJyoEhTJIkKQeGMEmSpBwYwiRJknJgCJMkScqBIUySJCkH/w/QXsMRAm0JkAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Regression Line\n",
    "y_pred = model.predict(x_test)\n",
    "\n",
    "# The coefficients\n",
    "print('True model: y ~ Bernoulli(μ)')\n",
    "print('True logit(μ) = log(μ/(1-μ)) =', beta_0, '+', beta_1, '* x')\n",
    "print('Regression logit(μ) = log(μ/(1-μ)) = %.2f' % model.params[0], '+ %.2f' % model.params[1], '* x')\n",
    "\n",
    "# Metrics\n",
    "print('Accuracy: ', accuracy_score(y_test, np.ceil(y_pred)))\n",
    "print('Precision: ', precision_score(y_test, np.ceil(y_pred)))\n",
    "print('Recall: ', recall_score(y_test, np.ceil(y_pred)))\n",
    "print('F1: ', f1_score(y_test, np.ceil(y_pred)))\n",
    "\n",
    "# Consider just the variable-column\n",
    "x_test_single_col = pd.DataFrame(x_test)[1]\n",
    "\n",
    "# Sort it (and y_pred)\n",
    "idx = x_test_single_col.argsort()\n",
    "x_test_ord, y_pred_ord = x_test_single_col[idx], y_pred[idx]\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(10,5));\n",
    "plt.scatter(x_test_single_col, y_test);\n",
    "plt.plot(x_test_ord, y_pred_ord, color='red', linewidth=3);\n",
    "plt.title('Logistic Regression');\n",
    "plt.xlabel('X');\n",
    "plt.ylabel('Y');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.10 64-bit (system)",
   "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.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "9805bccb659ffcb6952b8b29440036fac091f8170d71a5a402ad9a8f2988209f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}