02_linear_model_for_classification.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Methods for Classification"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * Since our predictor $G(x)$ takes values in a discrete set $G$, we can always divide the input space into a collection of regions labeled according to the classification.\n",
    "\n",
    " * The boundaries of these regions can be rough or smooth, depending on the prediction\n",
    "function. For the current model  these decision boundaries are *linear*.\n",
    "\n",
    "* Suppose there are $K$ classes, labeled $1, 2, ..., K$, and the fitted linear model for the $k^{th}$ indicator response variable is $\\hat{f}_k(x) = \\hat{β}_{k0} + \\hat{β}_k^T x.$, That mean we compute $\\hat{\\beta}_0$ and $\\hat{\\beta}^T$ for each class. Here $x = (x_1, ..., x_d)$.\n",
    "\n",
    "* The decision boundary between class $k$ and $l$ is that set of points for which $\\hat{f}_k(x) = \\hat{f}_l(x)$, that is, the set ${x : (\\hat{β}_{k0} − \\hat{β}_{l0}) + (\\hat{β}_k − \\hat{β}_l)^T x = 0}$, an affine set or *hyperplane*.\n",
    "\n",
    "* Above approach is a member of *discriminant functions* $δ_k(x)$, which classify $x$ to the class with the largest value for its discriminant function."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The other approach, is that model the posterior probabilities $Pr(G = k|X = x)$. Like a Logit, $P(G = k| X =x) = \\frac{e^z}{1 + e^z}$, where $z = B_{k0} + B_{k1}x$. We can compute *log-oods* $log(\\frac{Pr(G = k|X = x)}{1 -Pr(G = k|X = x)}) = z = B_{k0} + B_{k1}^{T}x $. The decision boundary is defined by $\\{x|β_{k0} + β_{k1}^T x = 0 \\}$, this since, we usually fix like a threshold $p=50\\%$, so this result in $z=0$ like a threshold\n",
    "\n",
    "* The decision in both is linear *hyperplane*.\n",
    "\n",
    "* The methods that explicitly look for \"separating hyperplanes.\" are two. The first is the wellknown *perceptron* model of Rosenblatt (1958), with an algorithm that finds a separating hyperplane in the training data, if one exists. \n",
    "* The second method, due to Vapnik (1996), finds an optimally separating hyperplane if one exists, else finds a hyperplane that minimizes some measure of overlap in the training data.\n",
    "\n",
    "<div align = \"center\">\n",
    "  <img src = \"assets/lin_model_class/Captura_lin_for_class_01.PNG\" />  \n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use only $X_1, X_2, ..., X_p$ like a features, the boundary decision is the left but if we add other combination of features like $X_1*X_2, ... ,$ or $X_1^2,X_2^2, ..$  the boundary decision is like right."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Discriminant Analysis"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $k = 1, 2, ..., K$ be the labels to classes and Suppose $f(x|G=k) = f_k(x)$ is the density function conditioned to $G=k$ and let $π_k$ be the prior probability of class $k$, with $\\sum_{k=1}^{K} π_k = 1$.\n",
    "\n",
    "Bayes Theorem  gives us:\n",
    "\n",
    "$$Pr(G = k|X = x) = \\frac{f_k(x)π_k} {f(x)} = \\frac{f_k(x)π_k} {\\sum_{l=1}^{K}f_l(x)\\pi_l} $$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a valid, since $\\sum_{k=1}^{K} Pr(G=k|X=x) = 1$, we only need to put $\\frac{f_k(x)\\pi_k}{\\sum f_l(x)\\pi_l}$ instead $Pr(G=k|X=x)$.\n",
    "\n",
    "As we can see we use density function $f_k$ instead probability function $Pr(X=x|G=k)$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that we model each class density as multivariate Gaussian:\n",
    "\n",
    "$$f_k(x) = \\frac{1 }{(2π)^{p/2}|Σ_k|^{1/2}} e^{−1/2(x−µ_k)^T Σ_k^{-1}(x−µ_k)}$$\n",
    "Linear discriminant analysis (LDA) arises in the special case when we assume that the classes have a common covariance matrix $Σ_k = Σ_{∀k}$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In comparing two classes $k$ and $ℓ$ it is sufficient to look at the log-ratio.\n",
    "\n",
    "$$log \\frac{Pr(G = k|X = x)}{Pr(G = ℓ|X = x)} = log{f_k(x)} + log{π_k} - (log{f_ℓ(x)} + log{π_ℓ})$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the hyperplane $log{f_k(x)} + log{π_k} -  (log{f_ℓ(x)} + log{π_ℓ}) = 0$. This means $log{f_k(x)} + log{π_k} = log{f_ℓ(x)} + log{π_ℓ}$. \n",
    "\n",
    "As we can see we can use one of them to compute the hyperplane. So we choose $log{f_k(x)} + log{π_k}$, this is the linear discriminant function $δ_k(x)$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$δ_k(x) = x^TΣ^{−1}µ_k − \\frac{1}{2}µ_k^T Σ^{−1} µ_k + log π_k$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An equivalent description of the decision rule, with $G(x) = \\arg\\max_k{δ_k}(x)$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function $δ_k(x)$ is linear in $x$, so all the decision boundaries are linear."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div align = \"center\">\n",
    "    <img src = \"assets/lin_model_class/Captura_lin_for_class_02.PNG\" />\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice we do not know the parameters of the Gaussian distributions, and will need to estimate them using our training data:\n",
    "    \n",
    "  * $\\hat{π}_k = N_k/N$, where $N_k$ is the number of class-k observations\n",
    "  * $\\hat{µ}_k = \\sum_{g_i=k}^{N} x_i/N_k$\n",
    "  * $\\hat{Σ} = \\sum_{k = 1}^{K} \\sum_{g_i = k}(x_i − \\hat{µ}_k)(x_i − \\hat{µ}_k)^T /(N − K)$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $δ_k(x)>δ_l(x)$ The LDA rule classifies to class $k$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the $Σ_k$ are not assumed to be equal, then the convenient cancellations in do not occur. \n",
    "In particular the pieces quadratic in $x$ remain. We then get *quadratic discriminant functions (QDA)*"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$δ_k(x) = −\\frac{1}{2} log|Σ_k| − \\frac{1}{2}(x − µ_k)^T Σ_{k}^{−1}(x − µ_k) + log π_k$\n",
    "\n",
    "The decision boundary between each pair of classes $k$ and $l$ is described by a quadratic equation $\\{x : δ_k(x) = δ_ℓ(x)\\}$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div align = \"center\">\n",
    "    <img src = \"assets/lin_model_class/Captura_lin_for_class_03.PNG\" />\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The left plot shows the quadratic decision boundaries obtained using LDA\n",
    "in the five-dimensional space $X_1, X_2, X_1X_2, X_{1}^{2}, X_2^2$.\n",
    "\n",
    "* The The right plot shows the quadratic decision boundaries found by QDA in 2-dimentional space $X_1, X_2$\n",
    "* The differences are generallysmall but QDA is the preferred approach"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regularized Discriminant Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot  as plt\n",
    "from sklearn import model_selection, metrics\n",
    "from sklearn.preprocessing import PolynomialFeatures, StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "data_drime = pd.read_stata(r'data_dta/CRIME1.dta')\n",
    "\n",
    "# Generate the variable arr86, where a person is labeled equal to 0, if he has not committed a crime in 1986\n",
    "# otherwise 1.\n",
    "data_drime['arr86'] = data_drime.narr86.where(data_drime.narr86==0, 1)\n",
    "\n",
    "# Selecting the features and target\n",
    "features= ['pcnv', 'avgsen', 'tottime', 'ptime86', 'qemp86']\n",
    "target = 'arr86'\n",
    "\n",
    "X = np.array(data_drime[features])\n",
    "\n",
    "# Generate polynomial and interaction features\n",
    "poly = PolynomialFeatures(2)\n",
    "X_poly = poly.fit_transform(X)\n",
    "\n",
    "# PCA\n",
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=X_poly.shape[1])\n",
    "X_trans = pca.fit_transform(X_poly)\n",
    "condition = pca.explained_variance_ratio_.cumsum() <= 0.997\n",
    "X_pca = X_trans[:,condition]\n",
    "\n",
    "y = np.array(data_drime[target])\n",
    "\n",
    "# Split the data to train and test\n",
    "X_train, X_test, y_train, y_test = model_selection.train_test_split(X_pca, y, random_state=123, test_size=0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearDiscriminantAnalysis()"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "\n",
    "clf = LinearDiscriminantAnalysis()\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "# svd: Singular value decomposition (default). \n",
    "# Does not compute the covariance matrix, \n",
    "# therefore this solver is recommended for data with a large number of features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "from assets.lin_model_class import utils_metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy 0.7273\n",
      "Initial AUC value is 0.4515\n",
      "\n",
      "Probability Conditioned on Predicts\n",
      "> P(not committed a crime|low risk)   72.97%\n",
      "> P(not committed a crime|high risk)  60.00%\n",
      "> P(committed a crime|low risk)       27.03%\n",
      "> P(committed a crime|high risk)      40.00%\n",
      "\n",
      "Probability Conditioned on Outcomes\n",
      "> P(low risk|not committed a crime)   99.40%\n",
      "> P(high risk|not committed a crime)  0.60%\n",
      "> P(low risk|committed a crime)       98.92%\n",
      "> P(high risk|committed a crime)      1.08%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "utils_metrics.assessment(y_test, X_test, clf)\n",
    "utils_metrics.prob_cond_pred(y_test, X_test, clf)\n",
    "utils_metrics.prob_cond_out(y_test, X_test, clf)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Cross - Validation for LDA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy 0.7273\n",
      "Initial AUC value is 0.4515\n",
      "\n",
      "Probability Conditioned on Predicts\n",
      "> P(not committed a crime|low risk)   72.97%\n",
      "> P(not committed a crime|high risk)  60.00%\n",
      "> P(committed a crime|low risk)       27.03%\n",
      "> P(committed a crime|high risk)      40.00%\n",
      "\n",
      "Probability Conditioned on Outcomes\n",
      "> P(low risk|not committed a crime)   99.40%\n",
      "> P(high risk|not committed a crime)  0.60%\n",
      "> P(low risk|committed a crime)       98.92%\n",
      "> P(high risk|committed a crime)      1.08%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fit_to_grid_cv = model_selection.GridSearchCV(\n",
    "\n",
    "    estimator=LinearDiscriminantAnalysis(\n",
    "        solver='lsqr'),\n",
    "    scoring= 'roc_auc',\n",
    "    cv=5,\n",
    "    param_grid={\n",
    "        'shrinkage': np.linspace(0,1, 10)},\n",
    "    verbose=False,\n",
    "    return_train_score=True\n",
    ")\n",
    "\n",
    "fit_to_grid_cv.fit(X_train, y_train)\n",
    "\n",
    "clf = LinearDiscriminantAnalysis(solver='lsqr', **fit_to_grid_cv.best_params_)\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "utils_metrics.assessment(y_test, X_test, clf)\n",
    "utils_metrics.prob_cond_pred(y_test, X_test, clf)\n",
    "utils_metrics.prob_cond_out(y_test, X_test, clf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array(data_drime[features])\n",
    "\n",
    "# standrization\n",
    "scaler = StandardScaler()\n",
    "X_trans = scaler.fit_transform(X)\n",
    "\n",
    "y = np.array(data_drime[target])\n",
    "\n",
    "# Split the data to train and test\n",
    "X_train, X_test, y_train, y_test = model_selection .train_test_split(X_trans, y, random_state=123, test_size=0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QuadraticDiscriminantAnalysis()"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis\n",
    "import numpy as np\n",
    "\n",
    "clf = QuadraticDiscriminantAnalysis()\n",
    "clf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy 0.7273\n",
      "Initial AUC value is 0.4515\n",
      "\n",
      "Probability Conditioned on Predicts\n",
      "> P(not committed a crime|low risk)   72.97%\n",
      "> P(not committed a crime|high risk)  60.00%\n",
      "> P(committed a crime|low risk)       27.03%\n",
      "> P(committed a crime|high risk)      40.00%\n",
      "\n",
      "Probability Conditioned on Outcomes\n",
      "> P(low risk|not committed a crime)   99.40%\n",
      "> P(high risk|not committed a crime)  0.60%\n",
      "> P(low risk|committed a crime)       98.92%\n",
      "> P(high risk|committed a crime)      1.08%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "utils_metrics.assessment(y_test, X_test, clf)\n",
    "utils_metrics.prob_cond_pred(y_test, X_test, clf)\n",
    "utils_metrics.prob_cond_out(y_test, X_test, clf)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The logistic regression model arises from the desire to model the posterior probabilities of the $K$ classes via linear functions in $B$, while at the same time ensuring that they sum to one and remain in $[0, 1]$. The model has\n",
    "the form:\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Pr(G = k|X = x) = \\frac{\\exp(β_{k0} + β_{k}^Tx)}{1 + \\sum_{ℓ=1}^{K-1} \\exp(β_{ℓ0} + β_{ℓ}^T x)} , k = 1,..., K - 1$$\n",
    "\n",
    "$$Pr(G = K|X = x) = \\frac{1}{1 + \\sum_{ℓ=1}^{K-1} \\exp(β_{ℓ0} + β_{ℓ}^T x)}$$\n",
    "\n",
    "* The model can be specified in terms of $K − 1$ log-odds or logit transformations\n",
    "* Although the model uses the last class as the denominator in the odds-ratios, the choice of denominator is arbitrary."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$log \\frac{Pr(G = 1|X = x)}{Pr(G = K|X = x)} = β_{10} + β_{1}^T x$$\n",
    "\n",
    "$$log \\frac{Pr(G = 2|X = x)}{Pr(G = K|X = x)} = β_{20} + β_{2}^T x$$\n",
    "\n",
    "$$ ... $$\n",
    "\n",
    "$$log \\frac{Pr(G = K - 1|X = x)}{Pr(G = K|X = x)} = β_{(K - 1)0} + β_{(K- 1)}^T x$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To emphasize the dependence on the entire parameter set $θ = \\{β_{10}, β_{1}^T , ... , β_{(K−1)0}, β_{K−1}^T \\}$, we denote the probabilities $Pr(G = k|X = x) = p_k(x; θ)$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting Logistic Regression Models"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $K = 2$ the model is of binary response."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Pr(G = 1|X = x) = \\frac{\\exp(β_{10} + β_{1}^Tx)}{1 + \\exp(β_{10} + β_{1}^T x)}$$\n",
    "\n",
    "$$Pr(G = 2|X = x) = \\frac{1}{1 + \\exp(β_{10} + β_{1}^T x)}$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can note only we need one of these equations, since, if we find, $Pr(G=2|X=x)$, automatically we can compute $Pr(G=1|X=x)$.\n",
    "\n",
    "\n",
    "The ecuation, $Pr(G = 1|X = x) = \\frac{1}{1 + \\exp(-(β_{10} + β_{1}^T x))}$ incorpore the *sigmoid function*, $\\sigma(z)=1/( 1 + \\exp(-z)$ , with $z = β_{10} + β_{1}^T x$. This take real values and map it to range [0, 1], so it create the probability"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Sigmoid Function\n",
    "> \n",
    "> $$\\sigma(z)=1/( 1 + \\exp(-z))$$\n",
    ">\n",
    "> The sigmoid function (named because it looks like an s) is also called the *logistic function*, and gives logistic regression its name. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAu0AAAGsCAYAAABpUpkzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDvElEQVR4nO3deXhU5eH28XsyyUz2AbJBIEDYl8gWFAGpggqiPytqFaQV95aqKFCXUvuqpba0VVvcQKm7IuK+tLigraDFjZCIbLKTkJUkkMk6k5k57x+B1JSATEhyJjPfz3XNNTNPzmHuXCMzt4fnPMdiGIYhAAAAAAErzOwAAAAAAI6P0g4AAAAEOEo7AAAAEOAo7QAAAECAo7QDAAAAAY7SDgAAAAQ4SjsAAAAQ4MLNDnAifD6fCgoKFBcXJ4vFYnYcAAAA4KQZhqHKykqlpqYqLOz4x9I7RGkvKChQWlqa2TEAAACAVpeXl6cePXocd5sOUdrj4uIkNfxC8fHxJqcBAAAATp7T6VRaWlpj1z2eDlHaj0yJiY+Pp7QDAAAgqJzI9G9ORAUAAAACHKUdAAAACHCUdgAAACDAUdoBAACAAEdpBwAAAAIcpR0AAAAIcJR2AAAAIMBR2gEAAIAAR2kHAAAAApzfpX3t2rW68MILlZqaKovForfeeusH91mzZo0yMzMVGRmpPn366PHHH29JVgAAACAk+V3aq6urNXz4cD366KMntP2ePXt0/vnna8KECcrOztZvfvMb3XLLLXr99df9DgsAAACEonB/d5g6daqmTp16wts//vjj6tmzpxYvXixJGjx4sNavX68HHnhAl156abP7uFwuuVyuxudOp9PfmAAAAEDQaPM57Z9//rkmT57cZGzKlClav3696uvrm91n0aJFcjgcjbe0tLS2jgkAAAAErDYv7UVFRUpJSWkylpKSIo/Ho9LS0mb3WbBggSoqKhpveXl5bR0TAAAAQcowDLk8XlXU1Kuook57Squ1pcCprH0H9Z+dpdqUX2F2xB/k9/SYlrBYLE2eG4bR7PgRdrtddru9zXMBAAAg8Hi8PlW5PKqs86ja7VFVnUdVroZbtcujKpdXNS6Pqt1eVbsatqlxeRvu3V7VuL2qPfy41u1VTb1XXp9xzNebmtFVS3+W2Y6/of/avLR37dpVRUVFTcZKSkoUHh6uhISEtn55AAAAtDPDMFTt9upQjVuHauobbrVuOWs9ctbVq6K2Xs7aejnrPKqorVdlXb0q6zyN9zVub5tls4ZZFB1hVaTNqqiIhltXR2SbvV5rafPSPnbsWL377rtNxj788EONHj1aERERbf3yAAAAOEkuj1fl1W6VVblVWuVSWZVbB2vcKq/+3n11vcpr3I1F3XOcI9snyhYepjh7uGIjwxVja7iPtYcrxh6uWLtV0baGxzE2q6KP3NvCFW2zKtpmVdTh51ERRx5bFWHtmJcp8ru0V1VVaefOnY3P9+zZo5ycHHXp0kU9e/bUggULlJ+fr+eff16SNHv2bD366KOaP3++brjhBn3++ed66qmntGLFitb7LQAAAOAXn89QeY1bRRV1OlDl0gGnq+G+0qWSyjodqHSp9HBJr6zztOg1bOFh6hwdoU5RNjmiIhQfFXH4PrzhPrJhLC4yXHGR4YqPPPI4QrH2cNnCO2bBbgt+l/b169dr4sSJjc/nz58vSbrqqqv07LPPqrCwULm5uY0/T09P16pVqzRv3jw99thjSk1N1cMPP3zM5R4BAABwcrw+QyWVdSo4VKfCiloVHKpVYUWdip11KqqoU7GzoZjXe0/8aHh4mEUJsTYlxNiVEGtTlxibOkcfvo+xqUu0TZ1jItQlxqZOUTZ1io5QZIS1DX/L0GIxjpwVGsCcTqccDocqKioUHx9vdhwAAABTebw+FVbUaf/BWuUdrNH+8hrtP1ir/QdrlX+oVsXOuhOanmKxSAkxdqXE25UUZ1dSrF3J8Q33SXGRSoi1KTG24Xl8VPgxFxFBy/jTcdtl9RgAAAD4x+XxKq+8RntLa7S3rFq55TXaW1ajfWXV2n+w9riroUgNJ1x2jY9UaqdIpXaKUldHpLrGN9xSDj9OirN32DneoYbSDgAAYBLDMFRa5dauA1XafaD68H2Vdh2o1v6DNTpeL7dZw9S9c5R6dI5Sj87Rh+8bHqd2ilRyXKSsYRwZDxaUdgAAgHZQXu3W9uJK7Siu1HfFldpeXKXtxZU6VNP8FeIlKdYerl4J0eqdEKOeCdHqnRCtXgkx6pUQrZS4SIVRykMGpR0AAKAVeX2G9pY1XHFza6FTWwqd2lLgVEmlq9ntLRapR+co9U2KVZ/EWPVNjmm8T4q1M48ckijtAAAALeb1Gdp1oEob91fo2/2HtDG/QtsKK1Vb3/zFgXp0jtLAlDj1T4nTgJRYDUiJU7/kWFZZwQ+itAMAAJyggkO12pB7UNm5h/Tt/gptKqho9uqdURFWDewapyGp8RrSLV6Du8VrUNc4xdipXmgZ/ssBAABohtvj0+aCCmXtayjpG3IPqrCi7qjtom1WZXR3aFh3h07p4VBGd4d6J8RwEihaFaUdAABAUl29V9m5h/TlnjJ9ubtcG3IPyuXxNdnGGmbR4G5xGpnWWcPTOml4D4f6JMVS0NHmKO0AACAkuT0+bcg9qHU7S/XF7nLl5B2S29u0pHeOjlBmr84a2bOzRvXsrOFpDkXbqE9of/xXBwAAQoJhGNpRUqVPd5Tqsx0H9OWe8qPmoyfH2TWmT4LGpHfR6X26qG9SLKu3ICBQ2gEAQNCqrKvXZztK9fG2Eq3dfuCoZRcTY20a3y9RY/skaEyfBPVOiKakIyBR2gEAQFDZU1qtj7cW69/fleirPeWq9/73sqL28DCdlt5FE/on6ox+SRrUNY4LFKFDoLQDAIAOzTAMbdxfofc3F+mDTUXaXVrd5Od9EmM0aVCyJg5KVmavzqyJjg6J0g4AADocr8/Q+r3ljUW94HtLMUZYLRqTnqCJg5I1aVCy0hNjTEwKtA5KOwAA6BB8PkMbcg/q7ZwCvbepUKVV7safRdusmjgoWVMzuurMAUmKi4wwMSnQ+ijtAAAgoG0rcurtnAK9k1Og/EO1jePxkeE6Z0iKpmZ004T+iUx7QVCjtAMAgIBT4qzT6xvy9XZOvrYVVTaOx9ismpLRVT8enqrx/RIVYQ0zMSXQfijtAAAgILg9Pv1rW4leXZ+nT7YfkNfXsOpLhNWiswYma9qI7po0KFlRNo6oI/RQ2gEAgKl2FFfqlfV5emNDvsqq/ztPPbNXZ/0ks4fOz+gmRzRz1BHaKO0AAKDduT0+fbC5SC98vk9f7S1vHE+MtevSzO66LDNN/ZJjTUwIBBZKOwAAaDdFFXV66atcrfgqVwcOX53UGmbR2YOSdfnoNJ05MIl56kAzKO0AAKBNGYahr/aU67nP9+qDzcWNc9WT4uyaeVpPzRzTUynxkSanBAIbpR0AALQJj9enDzYXa9naXfpmf0Xj+Gm9u+jKsb00ZWhX2cI5qg6cCEo7AABoVTVuj15dv19PfrZbeeUN66rbwsN06ajumjW2twZ3izc5IdDxUNoBAECrKK9265n/7NHzn+9TRW29JKlzdISuHNtbs8b2UmKs3eSEQMdFaQcAACeltMqlv3+6Wy98vk81bq8kqVdCtK4/I10/yUxjXXWgFVDaAQBAi5RU1mnZmt168ct9qqv3SZIyusfrprP6afLQrrKGWUxOCAQPSjsAAPBLibNOS9fs0ktf5srlaSjrw3o4dOvZ/TVpULIsFso60Noo7QAA4IRU1NZr2dpdevqzvaqtb5gGMyKtk249p7/OGpBEWQfaEKUdAAAcV129V89/vldLPtmlQzUNJ5iOSOuk+ecO0IT+iZR1oB1Q2gEAQLM8Xp/e2JCvv320XYUVdZKkfsmxun3KQE0ekkJZB9oRpR0AABxl7fYD+v0/tmhHSZUkqZsjUvPOGaBLRnVXuJULIgHtjdIOAAAa7Smt1h/+uUUfbS2RJDmiInTTxL6aNba3IiNYuhEwC6UdAACosq5ej/57p57+bI/qvYbCwyy6cmwvzT17gBzREWbHA0IepR0AgBDm8xl6bcN+/eX971Ra5ZIkTeifqHsuHKJ+yXEmpwNwBKUdAIAQta3Iqd+88a025B6SJPVOiNb/+78hrLUOBCBKOwAAIabW7dVDH+/Qk5/ulsdnKMZm1S1n99c149NlC+ckUyAQUdoBAAghn3xXot++tUn7D9ZKkqYMTdG9Px6qbo4ok5MBOB5KOwAAIaCksk4L392if2wslCSlOiL1u4sydO6QFJOTATgRlHYAAIKYYRh6Kydf97y9Wc46j8Is0rXj0zXv3AGKsVMDgI6Cv60AAASpkso63fXmJq3eUixJyugerz9dMkwZ3R0mJwPgL0o7AABBxjAMvbuxUHe/vUmHauoVYbXo1rP76xdn9lUEVzMFOiRKOwAAQaSsyqX/9/Ymrfq2SJI0pFu8Hrx8uAZ3izc5GYCTQWkHACBIfLSlWHe+vlFl1W6Fh1l086R+umliP46uA0GA0g4AQAdXV+/VolVb9dzn+yRJg7rG6YHLhjN3HQgilHYAADqw7cWVmvNStr4rrpQkXXdGuu44b6Ds4VaTkwFoTZR2AAA6IMMw9OKXubrvH1vk8viUGGvTA5cN11kDk82OBqANUNoBAOhgDla7defrG/Xh4aUcfzQgSQ9eNlxJcXaTkwFoK5R2AAA6kJy8Q7rxxSwVVNQpwmrRnecN0rXj0xUWZjE7GoA2RGkHAKADODIdZuG7m1XvNdQ7IVqPzhzFyaZAiKC0AwAQ4GrcHt315ia9mZ0vSZoyNEX3XzZc8ZERJicD0F4o7QAABLDdB6r0yxc36LviSlnDLLrzvIG6YUIfWSxMhwFCCaUdAIAA9f6mQt326kZVuTxKjLXr0ZkjdXqfBLNjATABpR0AgADj8xla/PEOPfzxDknSab276NGZI5UcH2lyMgBmobQDABBAatwe/eqVb/TepiJJ0rXj07Xg/EGKsIaZnAyAmSjtAAAEiPxDtbrhufXaUuhUhNWiP1x8ii4fnWZ2LAABgNIOAEAAyNpXrl+8kKXSKrcSY216/GeZGt27i9mxAAQISjsAACZ7LWu/fvPGt3J7fRrcLV5/n5WpHp2jzY4FIIBQ2gEAMInPZ+j+D7/T0k92SWpYf/2vl49QjJ2vZwBN8akAAIAJXB6v7nhto97OKZAkzZnUT/POGaCwMNZfB3C0Fp2KvmTJEqWnpysyMlKZmZn69NNPj7v98uXLNXz4cEVHR6tbt2665pprVFZW1qLAAAB0dBW19brq6a/0dk6BwsMsuv8nw/SryQMp7ACOye/SvnLlSs2dO1d33XWXsrOzNWHCBE2dOlW5ubnNbv/ZZ59p1qxZuu6667R582a9+uqr+vrrr3X99defdHgAADqa/EO1uuzxdfpid7li7eF6+upTdRkrxAD4ARbDMAx/dhgzZoxGjRqlpUuXNo4NHjxY06ZN06JFi47a/oEHHtDSpUu1a9euxrFHHnlEf/nLX5SXl9fsa7hcLrlcrsbnTqdTaWlpqqioUHx8vD9xAQAIGJsLKnTts1+r2OlScpxdz1xzqoamOsyOBcAkTqdTDofjhDquX0fa3W63srKyNHny5CbjkydP1rp165rdZ9y4cdq/f79WrVolwzBUXFys1157TRdccMExX2fRokVyOByNt7Q0jkAAADq2T3cc0PQnvlCx06UBKbF686bxFHYAJ8yv0l5aWiqv16uUlJQm4ykpKSoqKmp2n3Hjxmn58uWaPn26bDabunbtqk6dOumRRx455ussWLBAFRUVjbdjHZEHAKAj+MfGAl377Neqcnl0ep8uenX2OHXvFGV2LAAdSItORLVYmp4oYxjGUWNHbNmyRbfccovuvvtuZWVl6f3339eePXs0e/bsY/75drtd8fHxTW4AAHREL3+VqzkrslXvNXTBsG567trT5IiKMDsWgA7GryUfExMTZbVajzqqXlJSctTR9yMWLVqk8ePH6/bbb5ckDRs2TDExMZowYYLuu+8+devWrYXRAQAIbMvW7tIfV22TJF1xWk/dNy1DVlaIAdACfh1pt9lsyszM1OrVq5uMr169WuPGjWt2n5qaGoWFNX0Zq9UqqeEIPQAAwcYwDN3/wbbGwj77zL7648UUdgAt5/fFlebPn68rr7xSo0eP1tixY7Vs2TLl5uY2TndZsGCB8vPz9fzzz0uSLrzwQt1www1aunSppkyZosLCQs2dO1ennXaaUlNTW/e3AQDAZD6foXve2awXvtgnSbrjvIG68ax+JqcC0NH5XdqnT5+usrIyLVy4UIWFhcrIyNCqVavUq1cvSVJhYWGTNduvvvpqVVZW6tFHH9WvfvUrderUSZMmTdKf//zn1vstAAAIAPVen25/9Ru9lVMgi0X6/UUZ+tnpvcyOBSAI+L1Ouxn8WcMSAAAz1Ht9umVFtt7bVKTwMIsevHy4LhrR3exYAAKYPx3X7yPtAACgKbfHp5tf2qAPtxTLZg3T0p+N0tmDm1+gAQBagtIOAMBJcHm8uml5tj7aWixbeJieuDJTEwcmmx0LQJChtAMA0EIuj1e/fHGD/rWtRPbwMC2bNVpnDkgyOxaAIERpBwCgBerqvZr9YpY++e6A7OFheuqqU3VG/0SzYwEIUpR2AAD8VFfv1c9fyNLa7QcUGRGmp686VeP6UdgBtB1KOwAAfnB5/lvYoyKsevrqUzW2b4LZsQAEOUo7AAAnqN7r080vZTcW9meuOVWn96GwA2h7YWYHAACgI/D6DM1bmaPVWxpWiXnyqtEUdgDthtIOAMAP8PkM3fn6Rv1jY6EirBY98bNMjWcOO4B2RGkHAOA4DMPQPe9s1mtZ+xVmkR6eMVITB7EOO4D2RWkHAOAYDMPQove26YUv9slikR68fLimntLN7FgAQhClHQCAY1j80Q4tW7tbkvTHi0/RxSN7mJwIQKiitAMA0IynPtujhz7eIUm658IhuuK0niYnAhDKKO0AAPyPt7Lz9ft/bJEk3TZ5gK4Zn25yIgChjtIOAMD3/Pu7Et326jeSpGvG99ZNE/uZnAgAKO0AADTakHtQN764QR6foYtGpOr/XTBEFovF7FgAQGkHAECSdhRX6tpnv1ZtvVdnDkjS/T8ZrrAwCjuAwEBpBwCEvPxDtZr19Fc6VFOvEWmdtPRno2QL5ysSQODgEwkAENLKq92a9dSXKqyoU7/kWD1z9amKtoWbHQsAmqC0AwBCVl29V9c/97V2HahWN0eknr/2NHWOsZkdCwCOQmkHAIQkn8/QvJU52pB7SPGR4Xr+2tOU2inK7FgA0CxKOwAgJP1x1Va9t6lINmuY/j5rtPqnxJkdCQCOidIOAAg5z/xnj578bI8k6f7LhmlMnwSTEwHA8VHaAQAh5YPNRVp4+Gqnd5w3UBeN6G5yIgD4YZR2AEDIyM49qFtfzpZhSFec1lO/PLOv2ZEA4IRQ2gEAIWFfWbWuf2696up9mjgwSb+/aChXOwXQYVDaAQBBr6KmXtc887XKqt0amhqvR2eOUriVr0AAHQefWACAoFbv9enGl7K0u7RaqY5IPX31qYqxc/EkAB0LpR0AELQMw9C972zWf3aWKdpm1ZNXnaqU+EizYwGA3yjtAICg9dy6vVr+Za4sFumhGSM1JDXe7EgA0CKUdgBAUPrku5LGpR1/fd4gnTskxeREANBylHYAQNDZUVypOS9ly2dIl2X20M9/1MfsSABwUijtAICgUl7t1nXPrVely6PT0rvoDxefwtKOADo8SjsAIGi4PT7NfiFLueU16tklWo//LFO2cL7qAHR8fJIBAIKCYRi6551N+mpvueLs4XrqqtHqEmMzOxYAtApKOwAgKLz4Za5WfJWnMIv0yMyR6p8SZ3YkAGg1lHYAQIf35e4y/e6dzZKkO84bpLMGJpucCABaF6UdANCh5R+q1Y3LN8jjM/Tj4an6BSvFAAhClHYAQIdV6/bqFy+sV1m1W0NT4/XnS4exUgyAoERpBwB0SIZh6NdvbNSmfKe6xNj0xJWZirJZzY4FAG2C0g4A6JD+/uluvZ1ToPAwi5b8dJR6dI42OxIAtBlKOwCgw1m7/YD+9N42SdLdFw7R6X0STE4EAG2L0g4A6FDyyms0Z0W2fIY0fXSarjy9l9mRAKDNUdoBAB1GXb1Xs1/MUkVtvYanddLCaUM58RRASKC0AwA6BMMw9Nu3NmlzQcOJp0t/Okr2cE48BRAaKO0AgA7hpa9y9VrWfoVZpEevGKnUTlFmRwKAdkNpBwAEvOzcg7r3e1c8Hdcv0eREANC+KO0AgIBWVuXSjcs3qN5raMrQFK54CiAkUdoBAAHL4/VpzopsFVbUqU9SjB64bDgnngIISZR2AEDAeuDD7Vq3q0zRNque+Fmm4iIjzI4EAKagtAMAAtIHm4v0+JpdkqS//GSY+qfEmZwIAMxDaQcABJzcshrd9uo3kqTrzkjX/w1LNTkRAJiL0g4ACCh19V79cnmWKus8GtWzk349dZDZkQDAdJR2AEBAWfiPLdpc4FTn6Ag9OnOUIqx8VQEAn4QAgIDxVna+XvoyVxaLtHgGF1ACgCMo7QCAgLCjuFIL3vhWkjRnUn+dOSDJ5EQAEDgo7QAA01W7PPrl8g2qrffqjH6JuvXs/mZHAoCAQmkHAJjKMAzd9ea32llSpZR4uxbPGCFrGBdQAoDvo7QDAEy14qs8vZVTIGuYRY/OHKXEWLvZkQAg4FDaAQCm2VLg1L3vbpYk3T5loE7t3cXkRAAQmFpU2pcsWaL09HRFRkYqMzNTn3766XG3d7lcuuuuu9SrVy/Z7Xb17dtXTz/9dIsCAwCCQ7XLo5tf2iC3x6dJg5L18wl9zI4EAAEr3N8dVq5cqblz52rJkiUaP368nnjiCU2dOlVbtmxRz549m93n8ssvV3FxsZ566in169dPJSUl8ng8Jx0eANAxGYah3761SbtLq9U1PlIPXDZcYcxjB4BjshiGYfizw5gxYzRq1CgtXbq0cWzw4MGaNm2aFi1adNT277//vmbMmKHdu3erS5cT+2dPl8sll8vV+NzpdCotLU0VFRWKj4/3Jy4AIAC98nWe7nh9o6xhFr3889OZFgMgJDmdTjkcjhPquH5Nj3G73crKytLkyZObjE+ePFnr1q1rdp933nlHo0eP1l/+8hd1795dAwYM0G233aba2tpjvs6iRYvkcDgab2lpaf7EBAAEsO3Flbr7nU2SpPnnDqCwA8AJ8Gt6TGlpqbxer1JSUpqMp6SkqKioqNl9du/erc8++0yRkZF68803VVpaqhtvvFHl5eXHnNe+YMECzZ8/v/H5kSPtAICOrcbt0U3LN6iu3qcJ/RP1yzP7mh0JADoEv+e0S5LF0nTeoWEYR40d4fP5ZLFYtHz5cjkcDknSX//6V/3kJz/RY489pqiooy9RbbfbZbez5BcABJt739msHSVVSo6z62/TRzCPHQBOkF/TYxITE2W1Wo86ql5SUnLU0fcjunXrpu7duzcWdqlhDrxhGNq/f38LIgMAOqI3s/frlfX7FWaRHpoxkvXYAcAPfpV2m82mzMxMrV69usn46tWrNW7cuGb3GT9+vAoKClRVVdU4tn37doWFhalHjx4tiAwA6Gh2H6jSXW82zGO/9ewBGts3weREANCx+L1O+/z58/Xkk0/q6aef1tatWzVv3jzl5uZq9uzZkhrmo8+aNatx+5kzZyohIUHXXHONtmzZorVr1+r222/Xtdde2+zUGABAcHF5vJqzIls1bq/G9knQzZP6mR0JADocv+e0T58+XWVlZVq4cKEKCwuVkZGhVatWqVevXpKkwsJC5ebmNm4fGxur1atXa86cORo9erQSEhJ0+eWX67777mu93wIAELD+9N42bS5wqkuMTYtnjJCVeewA4De/12k3gz9rWAIAAsdHW4p1/fPrJUlPXz1akwY1f/4TAISiNlunHQCAE1VUUafbX/tGknTt+HQKOwCcBEo7AKDVeX2Gbn05Wwdr6jU0NV53Th1odiQA6NAo7QCAVvfov3bqyz3lirFZ9ejMUbKHW82OBAAdGqUdANCqvtpTroc+3i5J+v20DKUnxpicCAA6Pko7AKDVHKpx69aXs+UzpEtGdtclo7geBwC0Bko7AKBVGIahO1/fqMKKOvVOiNbCaRlmRwKAoEFpBwC0ipe+ytUHm4sVYbXokStGKdbu96VAAADHQGkHAJy0HcWV+v0/tkiS7pgySKf0cJicCACCC6UdAHBS6uq9mrMiW3X1Pk3on6jrzkg3OxIABB1KOwDgpPzpvW3aVlSpxFibHrx8uMLCLGZHAoCgQ2kHALTYx1uL9ey6vZKk+y8bruS4SHMDAUCQorQDAFqkxFmn21/bKEm6dny6Jg5MNjkRAAQvSjsAwG8+n6H5r3yj8mq3hnSL151TB5odCQCCGqUdAOC3v3+6W5/tLFVUhFUPXzFS9nCr2ZEAIKhR2gEAftm4/5Du/+A7SdI9Fw5Rv+RYkxMBQPCjtAMATli1y6NbX86Rx2doakZXTT81zexIABASKO0AgBO28N0t2lNarW6OSC265BRZLCzvCADtgdIOADghq74t1Mr1ebJYpL9ePkKdom1mRwKAkEFpBwD8oIJDtfr16w3LO/7yzL4a2zfB5EQAEFoo7QCA4/L6DM1dmSNnnUfDezg079wBZkcCgJBDaQcAHNfja3bpqz3lirFZ9dCMkYqw8tUBAO2NT14AwDFl5x7UX1dvlyT97qIM9U6MMTkRAIQmSjsAoFlVh5d39PoM/d+wbrp0VHezIwFAyKK0AwCadffbm5RbXqPunaL0h4tZ3hEAzERpBwAc5Z1vCvTGhnyFWaTFM0bIERVhdiQACGmUdgBAE/sP1uiuN7+VJN08sZ9O7d3F5EQAAEo7AKCRx+vTvJU5qqzzaGTPTrrl7P5mRwIAiNIOAPieJZ/s0td7DyrWHq6Hpo9UOMs7AkBA4NMYACBJytp3UA99vEOStPCioeqZEG1yIgDAEZR2AIAq6+o1d2W2vD5DPx6eqotHsrwjAAQSSjsAQHe/vVl55bXq3ilK912cwfKOABBgKO0AEOLezsnXm9kNyzs+NGOE4iNZ3hEAAg2lHQBCWF55jX775iZJ0s2T+ms0yzsCQECitANAiPJ4fZq7MkeVLo8ye3XWLZP6mR0JAHAMlHYACFGP/GunsvYdVJw9XIunj2B5RwAIYHxCA0AIWr+3XI/8q2F5x/suzlBaF5Z3BIBARmkHgBBTUVuvW1/Okc+QLhnZXReNYHlHAAh0lHYACCGGYei3b21S/qFa9ewSrd9dNNTsSACAE0BpB4AQ8saGfL37TYGsYRYtnjFCcSzvCAAdAqUdAELEvrJq3f12w/KO887pr1E9O5ucCABwoijtABAC6r0+3fJyjqrdXp2W3kW/PIvlHQGgI6G0A0AI+Nvq7fom75DiI8P1t+kjZA2zmB0JAOAHSjsABLl1O0u1dM0uSdKfLh2m7p2iTE4EAPAXpR0Aglh5tVvzXsmRYUgzTk3T+ad0MzsSAKAFKO0AEKQMw9Adr21UsdOlvkkxuvvCIWZHAgC0EKUdAILUi1/s00dbi2WzhunhK0Yq2hZudiQAQAtR2gEgCH1XVKn7/rlVknTn1EEamuowOREA4GRQ2gEgyNTVezVnxQa5PD6dNTBJ147vbXYkAMBJorQDQJD5wz+3antxlRJj7XrgsuGyWFjeEQA6Oko7AASRDzcX6YUv9kmS/nr5cCXG2k1OBABoDZR2AAgShRW1uuP1jZKkGyak60cDkkxOBABoLZR2AAgCHq9Pt67I0aGaeg3r4dDtUwaZHQkA0Ioo7QAQBB751059tbdcMTarHp4xUrZwPt4BIJjwqQ4AHdwXu8v0yL92SJL+cPEp6p0YY3IiAEBro7QDQAd2sNqtuS/nyGdIl47qoWkju5sdCQDQBijtANBBGYah21/7RkXOOvVJjNHCi4aaHQkA0EYo7QDQQT23bq8+2loimzVMj8wcqRh7uNmRAABthNIOAB3Q5oIK/XHVNknSgvMHaWiqw+REAIC2RGkHgA6myuXRnJey5fb6dM7gZF09rrfZkQAAbaxFpX3JkiVKT09XZGSkMjMz9emnn57Qfv/5z38UHh6uESNGtORlASDkGYah3775rXaXVqtrfKT+8pPhslgsZscCALQxv0v7ypUrNXfuXN11113Kzs7WhAkTNHXqVOXm5h53v4qKCs2aNUtnn312i8MCQKh7df1+vZVTIGuYRQ9fMVJdYmxmRwIAtAO/S/tf//pXXXfddbr++us1ePBgLV68WGlpaVq6dOlx9/vFL36hmTNnauzYsS0OCwChbHtxpe5+Z5Mkaf65A3RaeheTEwEA2otfpd3tdisrK0uTJ09uMj558mStW7fumPs988wz2rVrl+65554Teh2XyyWn09nkBgChrMbt0U3LN6iu3qcJ/RP1yzP7mh0JANCO/CrtpaWl8nq9SklJaTKekpKioqKiZvfZsWOHfv3rX2v58uUKDz+x5cgWLVokh8PReEtLS/MnJgAEnXve3qwdJVVKjrPrb9NHKCyMeewAEEpadCLq/570ZBhGsydCeb1ezZw5U7/73e80YMCAE/7zFyxYoIqKisZbXl5eS2ICQFB4Y8N+vZq1X2EW6aEZI5UYazc7EgCgnfl1JY7ExERZrdajjqqXlJQcdfRdkiorK7V+/XplZ2fr5ptvliT5fD4ZhqHw8HB9+OGHmjRp0lH72e122e18KQHAzpIq/fathnnst5zdX2P7JpicCABgBr+OtNtsNmVmZmr16tVNxlevXq1x48YdtX18fLy+/fZb5eTkNN5mz56tgQMHKicnR2PGjDm59AAQxOrqvbr5pQ2qcXs1rm+C5kzqb3YkAIBJ/L7m9fz583XllVdq9OjRGjt2rJYtW6bc3FzNnj1bUsPUlvz8fD3//PMKCwtTRkZGk/2Tk5MVGRl51DgAoKl739msbUWVSoy1afGMEbIyjx0AQpbfpX369OkqKyvTwoULVVhYqIyMDK1atUq9evWSJBUWFv7gmu0AgON7PWu/Xv46TxaLtHj6SCXHRZodCQBgIothGIbZIX6I0+mUw+FQRUWF4uPjzY4DAG3qu6JKXfTYZ6qr92neOQN06zlMiwGAYORPx23R6jEAgLZR5fLol8uzGtdjnzOpn9mRAAABgNIOAAHCMAwteONb7T5Qra7xkVrMeuwAgMMo7QAQIF78Yp/e/aZA4WEWPfbTkUpgPXYAwGGUdgAIABv3H9Lv/7FVkvTrqYOU2auLyYkAAIGE0g4AJquoqdeNyzfI7fVpytAUXXdGutmRAAABhtIOACby+QzNfyVH+w/WqmeXaP3lJ8NlsTCPHQDQFKUdAEz06L936uNtJbKFh2nJT0fJERVhdiQAQACitAOAST75rkR/+2i7JOm+aRnK6O4wOREAIFBR2gHABHnlNbr15RwZhjRzTE9dPjrN7EgAgABGaQeAdlZX79XsF7NUUVuv4WmddM+FQ8yOBAAIcJR2AGhHhmHot29t0uYCp7rE2LT0p6NkD7eaHQsAEOAo7QDQjl76KlevZe1XmEV69IqRSu0UZXYkAEAHQGkHgHaSnXtQ976zWZJ0x3mDNK5fosmJAAAdBaUdANrBgUqXbly+QfVeQ+cN7apf/KiP2ZEAAB0IpR0A2pjb49ONy7NUWFGnPkkxuv+yYVxACQDgF0o7ALSx3727WV/vPag4e7j+Pmu04iK5gBIAwD+UdgBoQy99mavlX+bKYpEeumKE+ibFmh0JANABUdoBoI2s31uue97ZJEm6bfJATRqUYnIiAEBHRWkHgDZQWFGr2S82nHh6wSnddONZfc2OBADowCjtANDK6uq9+sULWSqtcmlQ1zhOPAUAnDRKOwC0IsMw9Js3vtXG/RXqHB2hv88arWhbuNmxAAAdHKUdAFrRU5/t0RvZ+bKGWfTozFFK6xJtdiQAQBCgtANAK/l4a7H+sGqrJOmu8wdrPFc8BQC0Eko7ALSCbUVO3bIiW4YhXXFaT10zvrfZkQAAQYTSDgAnqbTKpeueXa9qt1dj+yRo4UVDOfEUANCqKO0AcBKOrBSTf6hW6YkxWvqzUYqw8tEKAGhdfLMAQAsZhqEFb3yrrH0HFR8ZrievGq1O0TazYwEAghClHQBaaMknu/Tm4ZVilvw0U32TYs2OBAAIUpR2AGiB9zcV6v4PvpMk/e7HQ3VGf1aKAQC0HUo7APgpO/eg5q7MkSRdPa63fnZ6L3MDAQCCHqUdAPywr6xa1z+3XnX1Pk0cmKTfXjDY7EgAgBBAaQeAE3Sw2q1rnvlaZdVuZXSP16MzRymclWIAAO2AbxsAOAF19V7d8Px67S6tVvdOUXr6qlMVYw83OxYAIERQ2gHgB/h8hn716jdav++g4iLD9cw1pyo5PtLsWACAEEJpB4Af8Of3t+mfGwsVYbXoiSszNSAlzuxIAIAQQ2kHgON44fO9emLtbknSX34yTOP6srQjAKD9UdoB4Bg+2Fyke97ZLEm6bfIAXTyyh8mJAAChitIOAM34YneZ5qzIls+QZpyappsm9jM7EgAghFHaAeB/bC6o0A3PrZfb49O5Q1J037QMWSwWs2MBAEIYpR0AvmdfWbWuevprVbo8Oi29ix65YiRrsQMATMc3EQAcVlJZp1lPf6XSKpcGd4vXk1eNVmSE1exYAABQ2gFAkpx19br66a+1r6xGPbtE67lrT1V8ZITZsQAAkERpB4CGq50+t15bCp1KjLXrhetOU3IcF08CAAQOSjuAkFbv9WnOimx9uadcsfZwPXvNqeqVEGN2LAAAmqC0AwhZXp+heStztHpLsWzhYfr7rNHK6O4wOxYAAEehtAMIST6foTtf36h/bCxUhNWiJ36WqbF9E8yOBQBAsyjtAEKOYRi6553Nei1rv6xhFj08Y6QmDko2OxYAAMdEaQcQUgzD0KL3tumFL/bJYpEevGy4pp7SzexYAAAcF6UdQEhZ/NEOLVu7W5L0x4tP0bSR3U1OBADAD6O0AwgZj6/ZpYc+3iFJuufCIbritJ4mJwIA4MRQ2gGEhCfW7NKf3tsmSbrjvIG6Zny6yYkAADhx4WYHAIC2tuSTnfrL+99Jkm49u79uPKufyYkAAPAPpR1AUHv0Xzv0wIfbJUnzzx2gW87ub3IiAAD8R2kHELQe+miH/vZRQ2G/fcpA3TSRI+wAgI6J0g4g6BiGocUf7Wg86fTO8wbpl2f1NTkVAAAtR2kHEFQMw9BfV2/XI//aKUlaMHWQfnEmhR0A0LFR2gEEDcMw9Kf3tumJw+uw//aCwbp+Qh+TUwEAcPIo7QCCgtdn6LdvfasVX+VJku7+vyG69gyWdQQABAdKO4AOz+3xad4rOfrnxkKFWaQ/XTJMl5+aZnYsAABaDaUdQIdW6/Zq9otZWrP9gCKsFj00Y6TOP6Wb2bEAAGhVLboi6pIlS5Senq7IyEhlZmbq008/Pea2b7zxhs4991wlJSUpPj5eY8eO1QcffNDiwABwREVtva586kut2X5AURFWPXnVqRR2AEBQ8ru0r1y5UnPnztVdd92l7OxsTZgwQVOnTlVubm6z269du1bnnnuuVq1apaysLE2cOFEXXnihsrOzTzo8gNBVWuXSFcu+0Pp9BxUXGa4Xrz9NZw5IMjsWAABtwmIYhuHPDmPGjNGoUaO0dOnSxrHBgwdr2rRpWrRo0Qn9GUOHDtX06dN19913n9D2TqdTDodDFRUVio+P9ycugCC0r6xaVz/ztfaUVisx1qbnrx2jIal8NgAAOhZ/Oq5fR9rdbreysrI0efLkJuOTJ0/WunXrTujP8Pl8qqysVJcuXY65jcvlktPpbHIDAEnKyTukS5as057SanXvFKVXfjGWwg4ACHp+lfbS0lJ5vV6lpKQ0GU9JSVFRUdEJ/RkPPvigqqurdfnllx9zm0WLFsnhcDTe0tJYBQKA9NGWYs1Y9rnKqt0amhqvN28cpz5JsWbHAgCgzbXoRFSLxdLkuWEYR401Z8WKFbr33nu1cuVKJScnH3O7BQsWqKKiovGWl5fXkpgAgsgLX+zTz19Yr7p6n84ckKSVvxir5PhIs2MBANAu/FryMTExUVar9aij6iUlJUcdff9fK1eu1HXXXadXX31V55xzznG3tdvtstvt/kQDEKR8PkP3f/idln6yS5I0fXSa7rs4QxHWFh1zAACgQ/LrW89msykzM1OrV69uMr569WqNGzfumPutWLFCV199tV566SVdcMEFLUsKIOS4PF7NfyWnsbDPO2eA/nTpKRR2AEDI8fviSvPnz9eVV16p0aNHa+zYsVq2bJlyc3M1e/ZsSQ1TW/Lz8/X8889Laijss2bN0kMPPaTTTz+98Sh9VFSUHA5HK/4qAILJgUqXZr+Ypax9BxUeZtGiS07RZaM5vwUAEJr8Lu3Tp09XWVmZFi5cqMLCQmVkZGjVqlXq1auXJKmwsLDJmu1PPPGEPB6PbrrpJt10002N41dddZWeffbZk/8NAASdTfkV+vnz61VQUae4yHAt+ekoTejPGuwAgNDl9zrtZmCddiB0vPdtoea/8o1q673qkxijv181Wn1ZIQYAEIT86bh+H2kHgLZgGIYe/nin/vbRdknShP6JevSKUXJER5icDAAA81HaAZiu1u3Vba9+o39+WyhJunZ8un5z/iCFc8IpAACSKO0ATLa3tFq/XL5BWwudirBa9PuLMjTjtJ5mxwIAIKBQ2gGY5sPNRfrVq9+oss6jhBiblvx0lMb0STA7FgAAAYfSDqDdebw+Pbh6e+P665m9OuuxmaPU1cEVTgEAaA6lHUC7Kq1y6ZYV2Vq3q0ySdM343vrN+YO5YBIAAMdBaQfQbrL2levG5RtU7HQp2mbVny8dpguHp5odCwCAgEdpB9DmvD5Dj6/Zpb+t3i6Pz1DfpBg9cWWm+iXHmR0NAIAOgdIOoE0VVdRp3socfb67YTrMhcNTteiSUxRr5+MHAIATxbcmgDbz4eYi3fH6Rh2qqVe0zarf/XiofpLZQxaLxexoAAB0KJR2AK2urt6r+/65RS9+kStJyuger4dnjFSfpFiTkwEA0DFR2gG0qi0FTs1dma3txVWSpJ//qI9umzxQtnBWhwEAoKUo7QBaRb3Xp6Wf7NLDH++Qx2coMdauv14+XD8akGR2NAAAOjxKO4CT9l1RpX71ao425TslSVOGpugPF5+ixFi7yckAAAgOlHYALebx+vTE2t166KMdcnt9ckRFaOFFQ/Xj4amcbAoAQCuitANokZ0llfrVqxv1Td4hSdI5g5P1x4tPUXJ8pLnBAAAIQpR2AH6pq/dqySe79Pgnu+T2+hQXGa57LxyqS0Z15+g6AABthNIO4ISt21mqu97apD2l1ZKkSYMajq53dXB0HQCAtkRpB/CDyqpc+sOqrXpjQ74kKTnOrnt/PFRTM7pydB0AgHZAaQdwTD6fodey9uuP723VoZp6WSzSlaf30m1TBio+MsLseAAAhAxKO4BmZece1O/e3aKcwyeaDu4Wr0WXnKIRaZ1MzQUAQCiitANoothZpz+/v61xKkyMzapbz+mva8enK9zKVU0BADADpR2ApIZVYZ76bI8e+/dO1bi9kqTLMnvo9ikDWcYRAACTUdqBEGcYht7bVKRF721VXnmtJGlkz06698KhGs5UGAAAAgKlHQhh/9lZqj+/v00b91dIklLi7VowdbAuGsEVTQEACCSUdiAEbcqv0J/f36ZPd5RKapi3fv2EPvr5j/ooxs7HAgAAgYZvZyCE7C2t1gMffqd/bCyUJEVYLfrpmF66eVI/JcbaTU4HAACOhdIOhIA9pdV67N879WZ2vrw+QxaLNG1Ed807Z4B6JkSbHQ8AAPwASjsQxHaWVOmxf+/U2zn58hkNYxMHJun2KYM0JDXe3HAAAOCEUdqBILSjuFKP/Gun3t1YIONwWT97ULLmnN2fiyMBANABUdqBILJ+b7mWrd2t1VuLG8v6uUNSdMuk/jqlh8PccAAAoMUo7UAH5/UZWr2lSMvW7taG3EON4+cN7ao5Z/fT0FTKOgAAHR2lHeigat1evbZhv576dLf2ltVIkmzWMF0yqruun5CufslxJicEAACthdIOdDC5ZTV68ct9emV9ng7V1EuSHFERuvL0Xpo1rpeS4yJNTggAAFobpR3oAHw+Q2u2H9Dzn+/VJ9sPNM5XT+sSpevGp+vyU9MUbeOvMwAAwYpveSCAlVW59PqG/Xrxi1zlltc0jv9oQJJmnd5LEwclyxpmMTEhAABoD5R2IMB4vD6t3XFAr3y9Xx9tLZbn8ALr8ZHhumx0mn52ei+lJ8aYnBIAALQnSjsQIPaUVuvV9Xl6fcN+FTtdjePDezh0xWk9ddGI7oqyWU1MCAAAzEJpB0xUWuXSqm8L9XZOgbL2HWwc7xJj08Uju+uy0T00qCtXLgUAINRR2oF2VuXy6MPNRXo7p0Cf7SyV9/D0lzCLdOaAJF0+Ok1nD06RLTzM5KQAACBQUNqBdlDl8ujf20r0/qYifbS1WC6Pr/Fnw3s49OMR3fV/w7opJZ7lGgEAwNEo7UAbqaip1+qtxXp/U5HW7jgg9/eKenpijC4akaofD09Vn6RYE1MCAICOgNIOtKK88hr9a1uJPtparM93lTWu/CJJvROidV5GN51/Sled0t0hi4WlGgEAwImhtAMnweP1aUPuIX28rVj/3lai7cVVTX4+MCVO52V01XkZXTWoaxxFHQAAtAilHfBT/qFafbbjgD7dUapPd5Sqora+8WfWMIsye3XWpEHJmjwkhakvAACgVVDagR/grKvXF7vK9NnOUn22o1S7S6ub/LxTdITOGpCkSYNTdGb/JDmiI0xKCgAAghWlHfgfFbX1+npPub7cU6Yv95RrU36Fvjc1XdYwi4b3cOiM/kn6Uf9EjezZWdYwpr0AAIC2Q2lHyCt21ilr30F9vbdcX+4u19Yipwyj6TbpiTE6o1+izuifqLF9ExQfydF0AADQfijtCCn1Xp+2FjqVte+gNuQe0oZ9B5V/qPao7fokxmhMny4ak56g09K7KLVTlAlpAQAAGlDaEbS8PkO7D1Rp4/4Kbdx/SBvzK7SlwNnkwkZSw5VIB3WN16henTQmPUFj0rsomYscAQCAAEJpR1BwebzaUVylLYVObSlwakuhU5vzK1Tt9h61rSMqQqN6dlJmr84a1bOzhqV1UqydvwoAACBw0VTQoRiGoYKKOm0vrtT2okptK6rU1kKndpZUNbmQ0RHRNqsyUh0a1sOhU3o4NKxHJ/XqEq0wThwFAAAdCKUdAcnrM5R/sFa7Squ0q6RKO4qrtL2kUjuKq1Tl8jS7T6foCA3pFq/B3eI1pFu8hvVwqE9SLCu7AACADo/SDtMYhqHSKrdyy6u1t7RGe0qrtbu0SrtKqrWnrFru/5l7fkR4mEV9kmLUPyVOA1PiNDS1oah3c0RyxVEAABCUKO1oU3X1XhUcqlXewVrtP1ijvPLaxpK+r6y62TnnR9isYeqdGK0+ibEakBLbUNK7xql3Qoxs4WHt+FsAAACYi9KOFjMMQ4dq6lVQUauCQ3UqrKhV/qGGxwWHGkp6sdN13D/DYpFSHVHqlRCtXgkx6psUo75JseqbFKvunaOY2gIAACBKO46h2uXRgUqXDlS5VOJ0qchZp2JnnYoq6hofFzvrVFff/BSW74uKsCqtS5TSOkerR+co9UyIUe/DJT2tS5Ts4dZ2+I0AAAA6Lkp7iPD6DDlr61VW7VZ5tVtlVS6VHr4vq3KrrNql0kq3SirrdKDSddxpK/8rMdambo4opXaKVDdHlLp3ilK3TpGNJb1LjI255gAAACeB0t7B+HyGKl0eOWvrVVFbL2dtvQ7V1utgjVuHahrGDtW4dbCm4b68+r+Pm1kR8biiIqxKjrcrKdauFEekUuIi1dVhV0p8pLrGR6qrI1Ip8ZGKjOBIOQAAQFuitLcTwzBUV+9TtdujapdHlXUeVbkaHlcdvlXWeVRZV6+quobHzsPPGx43FPRKl0eGn+X7++Iiw5UYa1dCjE0JsTYlxNqVGNNwnxBrU3JcpJLi7EqKs3PBIQAAgADRola2ZMkS3X///SosLNTQoUO1ePFiTZgw4Zjbr1mzRvPnz9fmzZuVmpqqO+64Q7Nnz25x6PaUV16jfWU1qq33qrbeqzr34fvDz2vdXlW7PapxNzxuvK/3qMZ1+GeH7/090n08kRFhio+MkCMqQp2jbXJER6hTVIQ6RUeoU7RNjqgIdYmxqXN0QznvHG1Tp+gIRVhZdQUAAKCj8bu0r1y5UnPnztWSJUs0fvx4PfHEE5o6daq2bNminj17HrX9nj17dP755+uGG27Qiy++qP/85z+68cYblZSUpEsvvbRVfom29GrWfj388Y5W/TNjbFbFRoYrxh6uOHvD/ZHHcZHhiouM+J/7cMVHNRT0+MgIxUeFc/ImAABACLEYhn+TLcaMGaNRo0Zp6dKljWODBw/WtGnTtGjRoqO2v/POO/XOO+9o69atjWOzZ8/WN998o88///yEXtPpdMrhcKiiokLx8fH+xD1pL32Zq+c/36vICKuiIqyKsjXcR0ZYFWULOzwWrmibVdGHfxZ9+HmM/b/3MXarYmzhioqwKoxlDAEAAEKePx3XryPtbrdbWVlZ+vWvf91kfPLkyVq3bl2z+3z++eeaPHlyk7EpU6boqaeeUn19vSIiIo7ax+VyyeX67/reTqfTn5itauaYnpo55uh/QQAAAADai18TnEtLS+X1epWSktJkPCUlRUVFRc3uU1RU1Oz2Ho9HpaWlze6zaNEiORyOxltaWpo/MQEAAICg0qKzEv93zW3DMI67Dndz2zc3fsSCBQtUUVHReMvLy2tJTAAAACAo+DU9JjExUVar9aij6iUlJUcdTT+ia9euzW4fHh6uhISEZvex2+2y2+3+RAMAAACCll9H2m02mzIzM7V69eom46tXr9a4ceOa3Wfs2LFHbf/hhx9q9OjRzc5nBwAAANCU39Nj5s+fryeffFJPP/20tm7dqnnz5ik3N7dx3fUFCxZo1qxZjdvPnj1b+/bt0/z587V161Y9/fTTeuqpp3Tbbbe13m8BAAAABDG/12mfPn26ysrKtHDhQhUWFiojI0OrVq1Sr169JEmFhYXKzc1t3D49PV2rVq3SvHnz9Nhjjyk1NVUPP/xwh1ijHQAAAAgEfq/TbgYz12kHAAAA2oI/HZdr2gMAAAABjtIOAAAABDhKOwAAABDgKO0AAABAgKO0AwAAAAGO0g4AAAAEOEo7AAAAEOD8vriSGY4sJe90Ok1OAgAAALSOI932RC6b1CFKe2VlpSQpLS3N5CQAAABA66qsrJTD4TjuNh3iiqg+n08FBQWKi4uTxWIxO05IcTqdSktLU15eHlejDTG896GL9z408b6HLt578xiGocrKSqWmpios7Piz1jvEkfawsDD16NHD7BghLT4+nr/IIYr3PnTx3ocm3vfQxXtvjh86wn4EJ6ICAAAAAY7SDgAAAAQ4SjuOy26365577pHdbjc7CtoZ733o4r0PTbzvoYv3vmPoECeiAgAAAKGMI+0AAABAgKO0AwAAAAGO0g4AAAAEOEo7AAAAEOAo7QAAAECAo7TDby6XSyNGjJDFYlFOTo7ZcdDG9u7dq+uuu07p6emKiopS3759dc8998jtdpsdDW1gyZIlSk9PV2RkpDIzM/Xpp5+aHQltbNGiRTr11FMVFxen5ORkTZs2Td99953ZsdDOFi1aJIvForlz55odBcdAaYff7rjjDqWmppodA+1k27Zt8vl8euKJJ7R582b97W9/0+OPP67f/OY3ZkdDK1u5cqXmzp2ru+66S9nZ2ZowYYKmTp2q3Nxcs6OhDa1Zs0Y33XSTvvjiC61evVoej0eTJ09WdXW12dHQTr7++mstW7ZMw4YNMzsKjoN12uGX9957T/Pnz9frr7+uoUOHKjs7WyNGjDA7FtrZ/fffr6VLl2r37t1mR0ErGjNmjEaNGqWlS5c2jg0ePFjTpk3TokWLTEyG9nTgwAElJydrzZo1+tGPfmR2HLSxqqoqjRo1SkuWLNF9992nESNGaPHixWbHQjM40o4TVlxcrBtuuEEvvPCCoqOjzY4DE1VUVKhLly5mx0ArcrvdysrK0uTJk5uMT548WevWrTMpFcxQUVEhSfwdDxE33XSTLrjgAp1zzjlmR8EPCDc7ADoGwzB09dVXa/bs2Ro9erT27t1rdiSYZNeuXXrkkUf04IMPmh0Frai0tFRer1cpKSlNxlNSUlRUVGRSKrQ3wzA0f/58nXHGGcrIyDA7DtrYyy+/rKysLK1fv97sKDgBHGkPcffee68sFstxb+vXr9cjjzwip9OpBQsWmB0ZreRE3/vvKygo0HnnnafLLrtM119/vUnJ0ZYsFkuT54ZhHDWG4HXzzTdr48aNWrFihdlR0Mby8vJ06623avny5YqMjDQ7Dk4Ac9pDXGlpqUpLS4+7Te/evTVjxgy9++67Tb68vV6vrFarfvrTn+q5555r66hoZSf63h/5MC8oKNDEiRM1ZswYPfvsswoL4//5g4nb7VZ0dLReffVVXXzxxY3jt956q3JycrRmzRoT06E9zJkzR2+99ZbWrl2r9PR0s+Ogjb311lu6+OKLZbVaG8e8Xq8sFovCwsLkcrma/Azmo7TjhOTm5srpdDY+Lygo0JQpU/Taa69pzJgx6tGjh4np0Nby8/M1ceJEZWZm6sUXX+SDPEiNGTNGmZmZWrJkSePYkCFDdNFFF3EiahAzDENz5szRm2++qU8++UT9+/c3OxLaQWVlpfbt29dk7JprrtGgQYN05513Mj0qADGnHSekZ8+eTZ7HxsZKkvr27UthD3IFBQU666yz1LNnTz3wwAM6cOBA48+6du1qYjK0tvnz5+vKK6/U6NGjNXbsWC1btky5ubmaPXu22dHQhm666Sa99NJLevvttxUXF9d4DoPD4VBUVJTJ6dBW4uLijirmMTExSkhIoLAHKEo7gOP68MMPtXPnTu3cufOo/0HjH+qCy/Tp01VWVqaFCxeqsLBQGRkZWrVqlXr16mV2NLShI0t8nnXWWU3Gn3nmGV199dXtHwhAs5geAwAAAAQ4ziQDAAAAAhylHQAAAAhwlHYAAAAgwFHaAQAAgABHaQcAAAACHKUdAAAACHCUdgAAACDAUdoBAACAAEdpBwAAAAIcpR0AAAAIcJR2AAAAIMD9fxRady7IfGFhAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 900x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (5, 3))\n",
    "\n",
    "x_grid = np.linspace(-5, 5, 100)\n",
    "f = lambda x:1/(1 + np.exp(-x))\n",
    "ax.plot(x_grid, f(x_grid));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we work with two classes is usually better work with $0/1$ response $y$. \n",
    "\n",
    "Suppose that, in actually the classes are instead $1 \\rightarrow \\text{no enrolled}$, $2 \\rightarrow \\text{enrolled}$\n",
    "\n",
    "So  $ y = 1 $ if $i$ is $\\text{enrolled}$ and $y = 0$ if $i$ is $\\text{no enrolled}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let be the probabilities $Pr(G = k|X = x) = p_k(x; θ)$.\n",
    "\n",
    "Let $p(x; θ) = p_{\\text{enrolled}}(x; θ) = Pr(G = \\text{enrolled}|X = x)$, and $p_{\\text{no enrolled}}(x; θ) = 1 − p(x; θ) = Pr(G = \\text{no enrolled}|X = x)$. This means that $y$ ~ $\\text{Bern}(p(x; θ))$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$p_{k}(x;\\theta)= p(x;\\theta)^{y}\\;(1-p(x;\\theta))^{(1 - y)}$$\n",
    "\n",
    "Where $p_{k}(x;\\theta)$ is $PMF$ of $y$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pugging in likehood usig the $N$ samples\n",
    "\n",
    "$$l(\\theta) = \\sum_{i=1}^{N}log \\;p(x_i;\\theta)^{y_i}\\;(1-p(x_i;\\theta))^{(1 - y_i)} $$\n",
    "\n",
    "$$l(\\theta) = \\sum_{i=1}^{N} y_i\\; log \\;p(x_i;\\theta) + (1 - y_i)\\;log(1-p(x_i;\\theta)) $$\n",
    "\n",
    "In term of $\\log(\\text{odds}) = -(β_{10} + β_{1}^T x)$\n",
    "\n",
    "$$l(\\theta) = \\sum_{i=1}^{N} -y_i \\log(odds) + log(1 + e^{\\log(odds)} )$$\n",
    "\n",
    "and its derivative for a sample is  $\\frac{\\partial l}{\\partial log(odds)} = -y + \\frac{e^{log(odds)}}{ 1 + e^{log(odds)}}$. Recall,  $ \\frac{e^{log(odds)}}{ 1 + e^{log(odds)}} = p$, so the derivative can be rewritten like $\\frac{\\partial(l (y, \\log(odds)))}{\\partial \\log(odds)} =-y  + p$\n",
    "\n",
    "Let $\\theta = \\{β_{10} + β_{1}\\}$ be\n",
    "\n",
    "$ \\frac{\\partial \\log(\\text{odds})}{\\partial β_{1}^T}  = -x$\n",
    "\n",
    "$\\frac{\\partial Ll}{\\partial log(odds)} . \\frac{\\partial \\log(\\text{odds})}{\\partial β_{1}^T} = (-y +p)-x = x(y - p)$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the coeficients we can use the  *Rapson Newton method* on gradient of $L$ or *Gradient Decent* \n",
    "\n",
    "$$\\frac{\\partial l(\\theta)}{\\partial β_{1}^T} = \\sum_{i=1}^{N} x_i\\; (y_i - p(x_i;\\theta))  = 0$$\n",
    "\n",
    "$$\\frac{\\partial l(\\theta)}{\\partial β_{10}} = \\sum_{i=1}^{N}(y_i - p(x_i;\\theta))  = 0$$\n",
    "\n",
    "For *Rapson New Method*\n",
    "\n",
    "$${\\beta_{1}^T}_{new} = {\\beta_{1}^T}_{old} - \\frac{\\partial \\partial l(\\theta)}{\\partial \\beta_{1}^T \\partial \\beta_{1}^T}^{-1} \\frac{\\partial l(\\theta)}{\\partial \\beta_{1}^T}$$\n",
    "\n",
    "$${\\beta_{10}}_{new} = {\\beta_{10}}_{old} - \\frac{\\partial \\partial l(\\theta)}{\\partial \\beta_{10} \\partial \\beta_{10}}^{-1} \\frac{\\partial l(\\theta)}{\\partial \\beta_{10}}$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The negative of this is called *negative log likehood* or *Cross Entropy* at level of record.\n",
    "> \n",
    "> $$L(y, p)=-[y\\log(p) + (1-y)\\log(1-p)]$$\n",
    "> \n",
    "> And it is used in neural network like a loss function for binary response in classification model. A general version of this is \n",
    "> \n",
    "> $$L(\\textbf{y}, \\textbf{p})=-\\sum_{k=1}^{K}y_k\\log(p_k)$$\n",
    "> \n",
    ">An example:\n",
    "> $\\textbf{p} = [0.1, 0.1, 0.6, 0.2]$ and $\\textbf{y}=[0, 0, 1, 0]$ the cross entropy would be $0.5108$. Validated with torch CrossEntropyLoss. This package take as inputs the *logits*, internally it coonvert to probabilities and then compute the cross entropy loss"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Logistic regression models are usually fitted by maximum likelihood.\n",
    "> \n",
    "> The log likelihood can be derived by defining, for each individual, $d_{ik} = 1$ if alternative $k$ is chosen by individual $i$, and $0$ if not, for the $K - 1$ possible outcomes. Then, for each $i$, one and only one of the $d_{ik}’s$ is $1$.\n",
    ">\n",
    "> Recall the *conditional probabilities* are *probabilities* and the maximun likelihood maximize the conditional probabilities $Pr(g_i=k|x=x_i)$ given the $\\theta$ parameters\n",
    ">\n",
    "> <!-- $$L(\\theta) = p_k(x; θ)  = \\prod_{i=1}^{N} \\prod_{k=1}^{K-1}p_{g_i}(x_i;\\theta) $$ -->\n",
    "> <!-- $$L(\\theta)  = \\prod_{i=1}^{N} \\prod_{k=1}^{K} d_{ik} Pr(g_i = k|X = x_i; θ) $$ -->\n",
    "> $$l(\\theta) = \\sum_{i=1}^{N}\\sum_{k=1}^{K} d_{ik} log \\;p_{g_i}(x_i;\\theta) $$\n",
    "> where $p_{g_i}(x_i; θ) = Pr(g_i = k|X = x_i; θ)$. The probability of the $i^{th}$ record belong to class $k$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Descent"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is with *gradient descent* find optimal parameters, parameters that minimize the *loss function* on all train examples. \n",
    "\n",
    "$\\theta = (\\beta_{10}, \\beta_{1}^{T})$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\hat{\\theta}=\\arg\\min_{\\theta} \\frac{1}{N} \\sum_{i=1}^{N} L(y_i, p(\\textbf{x}_i;\\theta))$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gradient descent is a method that finds a minimum of a function by figuring out in which direction (in the space of the parameters $\\theta$) the function’s slope is rising, and moving in the opposite direction.\n",
    "\n",
    "For logistic regression, this loss function is conveniently convex. A **convex** function has at most one minimum; there are no local minima to get stuck in, so *gradient descent* starting from any point is guaranteed to find the minimum.\n",
    "\n",
    "\n",
    "> By contrast, the loss for multi-layer neural networks is non-convex, and gradient descent may\n",
    "> get stuck in local minima for neural network training and never find the global optimum."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters are updated following the following rule\n",
    "\n",
    "$$ \\theta_{t} =\\theta_{t-1} - \\alpha \\nabla L(y, p(\\textbf{x};\\theta)) $$\n",
    "\n",
    "Where $\\alpha$ is *learning rate* and $\\nabla L(y, p(\\textbf{x};\\theta))$ is the gradient for all parameters $\\theta$ that is represented by partial derivative of $L(y, p(\\textbf{x};\\theta))$ "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\nabla L(y, p(\\textbf{x};\\theta)) = \\begin{bmatrix}\n",
    "\\frac{\\partial L}{\\partial B_{10}} \\\\\n",
    "\\frac{\\partial L}{\\partial B_{11}} \\\\\n",
    "\n",
    "... \\\\\n",
    "\n",
    "\\frac{\\partial L}{\\partial B_{1D}}\n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see how its works, we can think the partial derivatives like slope. If the slope is negative for a $\\beta$, the $\\beta$ must be updated positively $+$ (in opposite direction). Why?\n",
    "\n",
    "<center>\n",
    "\n",
    "<img src = \".\\assets\\lin_model_class\\gradient_descent.png\" height = 300>\n",
    "\n",
    "</center>\n",
    "\n",
    "If the slope is negative it means that we can minimize the loss function increasing the $\\beta$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead using the log likelihood, we will use the negative log likelihood (cross-entropy) to update the parameters ($\\theta$) in the logistic model (binomial). so instead maximize, we will minimize it.\n",
    "\n",
    "\n",
    "$$ L(y, p(\\textbf{x};\\theta)) = - [y\\; log \\;p(\\textbf{x};\\theta) + (1 - y)\\;log(1-p(\\textbf{x};\\theta))] $$\n",
    "\n",
    "For the second feature\n",
    "\n",
    "$$\\frac{ \\partial L(y, p(\\textbf{x};\\theta))}{\\partial \\beta_{2}} = -(y-\\hat{y}) x_{2} $$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Stochastic Gradient Descent Algorithm\n",
    "\n",
    "Stochastic gradient descent is an online algorithm that minimizes the loss function\n",
    "by computing its gradient after each training example, and nudging q in the right\n",
    "direction (the opposite direction of the gradient). (An “online algorithm” is one that\n",
    "processes its input example by example, rather than waiting until it sees the entire\n",
    "input."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function *STOCHASTIC GRADIENT DESCENT*($ L(), f (), x, y$) returns $\\theta$\n",
    "\n",
    "where: $L$ is the loss function\n",
    "- $f$ is a function parameterized by $\\theta$\n",
    "- $\\textbf{x}$ is the set of training inputs $x_1, x_2, ..., x_N$\n",
    "- $y$ is the set of training outputs (labels) $y_1, y_2, ..., y_N$\n",
    "\n",
    "$\\theta = 0$\n",
    "\n",
    "**repeat** until done # see caption\n",
    "\n",
    "For each training tuple $(\\textbf{x}_i, y_i)$ (in random order)\n",
    "1. Optional (for reporting): # How are we doing on this tuple?\n",
    "  - Compute $\\hat{y}_i = p(\\textbf{x}_i; \\theta$) # What is our estimated output $\\hat{y}_i$?\n",
    "  - Compute the loss $L(\\hat{y}_i, y_i)$ # How far off is $\\hat{y}_i$ from the true output $y_i$?\n",
    "2. $g = \\nabla L(y, p(\\textbf{x};\\theta))$ # How should we move q to maximize loss?\n",
    "3. $\\theta_{t} =\\theta_{t-1} - \\alpha g $ # Go the other way instead\n",
    "\n",
    "**return** $\\theta$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *learning rate* $\\alpha$ is a **hyperparameter** that must be adjusted. If it’s too high,\n",
    "the learner will take steps that are too large, overshooting the minimum of the loss\n",
    "function. If it’s too low, the learner will take steps that are too small, and take too\n",
    "long to get to the minimum. \n",
    "\n",
    "It is common to start with a higher learning rate and then\n",
    "slowly decrease it, so that it is a function of the iteration k of training; the notation\n",
    "hk can be used to mean the value of the learning rate at iteration k"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stochastic gradient descent is called stochastic because it chooses a single random\n",
    "example at a time, moving the weights so as to improve performance on that single\n",
    "example.\n",
    "\n",
    "it’s common to compute the gradient over batches of training instances rather than a single instance"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In *batch training* we compute the gradient over the entire dataset.\n",
    "\n",
    "In *mini-batch training*: we train on a group of $m$ examples (perhaps 512, or 1024) that is less than the whole dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm, uniform\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [],
   "source": [
    "unif = uniform()\n",
    "e = norm()\n",
    "size = 700\n",
    "random_state = 123"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.DataFrame({'A':unif.rvs(size, random_state),\n",
    "                     'B':unif.rvs(size, 2*random_state)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['y'] = 8*data['A'] + 2*data['B'] + e.rvs(size, random_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7.821120803194643, 2.127896076708602)"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epochs = 10\n",
    "ca, cb = 0, 0\n",
    "alpha = 0.01\n",
    "\n",
    "while epochs>0:\n",
    "\n",
    "    data_random = data.sample(frac=1, random_state=epochs)\n",
    "    data_splits = np.array_split(data_random, 100)\n",
    "    \n",
    "    while data_splits:        \n",
    "        data_current = data_splits.pop()\n",
    "        y_current = data_current['y'].copy()\n",
    "        y_hat = ca*data_current['A'] + cb*data_current['B']\n",
    "        la = -(y_current - y_hat)*data_current['A']\n",
    "        lb = -(y_current - y_hat)*data_current['B']\n",
    "        ca, cb = ca - alpha*np.sum(la), cb - alpha*np.sum(lb)\n",
    "    epochs-=1\n",
    "\n",
    "ca, cb"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regularization\n",
    "\n",
    "There is a problem with learning weights that make the model perfectly match the\n",
    "training data. If a feature is perfectly predictive of the outcome because it happens\n",
    "to only occur in one class, it will be assigned a very high weight. \n",
    "\n",
    "The weights for\n",
    "features will attempt to perfectly fit details of the training set, in fact too perfectly,\n",
    "modeling noisy factors that just accidentally correlate with the class. This problem is\n",
    "called overfitting. A good model should be able to generalize well from the training data to the unseen test set, but a model that overfits will have poor generalization"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To avoid overfitting, a new **regularization** term $R(\\theta)$ is added to the objective function (*likelihood*) instead loss function. Why? we can undenstand better.\n",
    "\n",
    "\n",
    "$$ \\hat{\\theta} = \\arg \\max_{\\theta}  \\sum_{i=1}^{S}  y_i\\; log \\;p(x_i;\\theta) + (1 - y_i)\\;log(1-p(x_i;\\theta)) - R(\\theta)$$\n",
    "\n",
    "Where $S$ is the size of the *mini-batch*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This penalty is particularly useful in the $p \\gt N$ situation, or any situation where there are many correlated predictor variables"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two ways to compute this regularization term. \n",
    "\n",
    "\n",
    "- **L2** *regularization*  $R(\\theta) = \\sum_{d=1}^{D} \\beta_{d}^2$. \n",
    "\n",
    "  - *Ridge regression* is known to shrink the coefficients of correlated predictors towards each\n",
    "other, allowing them to borrow strength from each other. In the extreme case of k identical\n",
    "predictors, they each get identical coefficients with 1/kth the size that any single one would\n",
    "get if fit alone\n",
    "  \n",
    "- **L1** *regularization*  $R(\\theta) = \\sum_{d=1}^{D} |\\beta_{d}|$\n",
    "  - *Lasso*, on the other hand, is somewhat indifferent to very correlated predictors, and will tend\n",
    "to pick one and ignore the rest.\n",
    "  - If we fit $y = f(X, \\beta)$ where $f(X, \\beta) = \\beta_{1}x_1 + \\beta_{2}x_2 + \\beta_{3}x_3 + \\epsilon$ and $x_1 = x_2=x_3$ perfectly correlated one each other, then *coefficient Lasso* would fit to: $f(X, \\hat{\\beta}) = \\hat{\\beta}_{1}x_1 + 0x_2 + 0x_3 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ssl\n",
    "from io import StringIO\n",
    "import urllib\n",
    "import re\n",
    "from sklearn.linear_model import LinearRegression, Lasso"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Evaluation"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Accuracy**"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The common way to evaluate a model is using accuracy, the $%$ of correct prediction"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We usually assign $1$ to positive response that something that we need to know. Is is span? Does she like me?\n",
    "\n",
    "$$\\text{Yes}:1, \\text{No}:0$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Confusion Matrix*\n",
    "\n",
    "It is a table for visualizing how an algorithm performs with respect to the human gold labels, using two dimensions (system output and gold labels), and each cell labeling a set of possible outcomes\n",
    "\n",
    "<center>\n",
    "\n",
    "<img src = \".\\assets\\lin_model_class\\confusion_matrix.png\" height = 300>\n",
    "\n",
    "</center>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\text{Precision} = \\frac{TP}{FP + TP} $$\n",
    "\n",
    "$$\\text{Recall} = \\frac{TP}{FN + TP} $$\n",
    "\n",
    "$$\\text{Accuracy} = \\frac{TP + TN}{TP + TN + FP + FN}$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\text{F-Measure}$ is a simple combination of Precision and Recall. It is a harmony mean weighted by the importance of them. If both are important $\\alpha = 0.5$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$F = (\\frac{1}{\\alpha \\frac{1}{P} +  (1-\\alpha)\\frac{1}{R} })$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decision Boundary"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we make a decision about which class to apply to a test instance x? For\n",
    "a given x, we say yes if the probability P(y = 1jx) is more than .5, and no otherwise.\n",
    "boundary decision We call .5 the decision boundary"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\text{If} \\: \\hat{P}(g = 1|x) \\ge 0.5 \\: \\text{then} \\: 1\\: \\text{otherwise} \\: 0 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assign $1$ when  $\\hat{P}(g = 1|x) \\ge 0.5$  then $0.5$ is boundary decision.\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The decision boundary can be changed (increase the threshold) with the goal to improve *precision*, but the *recall* will suffer. Will apply false to some true cases that had a model output just below the raised threshold.\n",
    "\n",
    "A way to visualize this trade-off with *ROC* curve.\n",
    "\n",
    "Here appear the concept of *ROC* and *AUC*\n",
    "\n",
    "The *ROC* curve traces the relationship between the *false positive rate*  - *[FP/(TN + TP)]*  (on the $x$ axis) and the *true positive rate* - *[TP/(FN + TP)]* (on the $y$ axis) as the probability threshold $p$ is changed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we increase the threshold in order to improve precision [TP/(TP + FP)], then (TP + FP) the number of records with positive prediction  decrece since we are assing 1 to record with higher probability. The TP decrece as well, but the decrecing of (TP + FP) is strong growth. The result: *increasing precision*. The *false positive rate* decreace as well, since we are assing 1 to record with higher probability to be positive and we are keeping with these. The consecuence of this action is decresing in *true positive rate*. Again, this happend since we are assing 1 to record with higher positive probability, the quantity of records that over the new threshold usually are smaller than one with small threshold."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div align = \"center\">\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/1/13/Roc_curve.svg\" height = \"400\">\n",
    "</div>\n",
    "\n",
    "Source: https://en.wikipedia.org/wiki/Receiver_operating_characteristic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot  as plt\n",
    "from sklearn import model_selection, metrics\n",
    "from sklearn.preprocessing import PolynomialFeatures, StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "import re\n",
    "from io import StringIO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Binary Choice Model, the response is a *binary variable*. This variable can be make it based on a quantitative variable (using bining) or in actually can be a binary response $0$ (spam) and $1$ (no spam).\n",
    "\n",
    "In econometrics *binary variable* hides an choice, this choice is based on *Utility* of this choice "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset extracted from: https://pages.stern.nyu.edu/~wgreene/Text/Edition7/tablelist8new.htm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(r'assets\\lin_model_class\\pages.stern.nyu.edu_~wgreene_Text_Edition7_TableF5-1.txt') as table_5:\n",
    "    data = table_5.read()\n",
    "    data_info = re.search(r'([\\*].+?(?=\\*\\/))\\*\\/\\n(.+)', data, re.DOTALL).group(1)\n",
    "    data_cleaned = re.search(r'([\\*].+?(?=\\*\\/))\\*\\/\\n(.+)', data, re.DOTALL).group(2)\n",
    "    data_formatted = re.sub(r'[ ]+', r',', data_cleaned)\n",
    "    data_consol = pd.read_csv(StringIO(data_formatted), sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*\n",
      "doc LFP  \"A dummy variable = 1 if woman worked in 1975, else 0\";\n",
      "doc WHRS \"Wife's hours of work in 1975\";\n",
      "doc KL6  \"Number of children less than 6 years old in household\";\n",
      "doc K618 \"Number of children between ages 6 and 18 in household\";\n",
      "doc WA   \"Wife's age\";\n",
      "doc WE   \"Wife's educational attainment, in years\";\n",
      "doc WW   \"Wife's average hourly earnings, in 1975 dollars\";\n",
      "doc RPWG \"Wife's wage reported at the time of the 1976 interview\\\n",
      " (not the same as the 1975 estimated wage).\\\n",
      " To use the subsample with this wage, one needs to select 1975\\\n",
      " workers with LFP=1, then select only those women with non-zero RPWG.\\\n",
      " Only 325 women work in 1975 and have a non-zero RPWG in 1976.\";\n",
      "doc HHRS \"Husband's hours worked in 1975\";\n",
      "doc HA   \"Husband's age\";\n",
      "doc HE   \"Husband's educational attainment, in years\";\n",
      "doc HW   \"Husband's wage, in 1975 dollars\";\n",
      "doc FAMINC \"Family income, in 1975 dollars.\\\n",
      " This variable is used to construct the property income variable.\";\n",
      "doc MTR  \"This is the marginal tax rate facing the wife, and is taken from\\\n",
      " published federal tax tables (state and local income taxes are excluded).\\\n",
      " The taxable income on which this tax rate is calculated\\\n",
      " includes Social Security, if applicable to wife.\";\n",
      "doc WMED \"Wife's mother's educational attainment, in years\";\n",
      "doc WFED \"Wife's father's educational attainment, in years\";\n",
      "doc UN   \"Unemployment rate in county of residence, in percentage points.\\\n",
      " This taken from bracketed ranges.\";\n",
      "doc CIT  \"Dummy variable = 1 if live in large city (SMSA), else 0\";\n",
      "doc AX   \"Actual years of wife's previous labor market experience\";\n",
      "?Source:  1976 Panel Study of Income Dynamics,\n",
      "?based on data for the previous year, 1975.  \n",
      "?Of the 753 observations, the first 428 are for women with positive hours\n",
      "?worked in 1975, while the remaining 325 observations are for women who\n",
      "?did not work for pay in 1975.  A more complete discussion of the data is\n",
      "?found in Mroz [1987], Appendix 1. \n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(data_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h2o\n",
    "from h2o.estimators import H2OGeneralizedLinearEstimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Checking whether there is an H2O instance running at http://localhost:54321 ..... not found.\n",
      "Attempting to start a local H2O server...\n",
      "; OpenJDK 64-Bit Server VM JBR-11.0.13.7-1751.21-jcef (build 11.0.13+7-b1751.21, mixed mode)\n",
      "  Starting server from C:\\Users\\caima\\anaconda3\\Lib\\site-packages\\h2o\\backend\\bin\\h2o.jar\n",
      "  Ice root: C:\\Users\\caima\\AppData\\Local\\Temp\\tmpmja0rchn\n",
      "  JVM stdout: C:\\Users\\caima\\AppData\\Local\\Temp\\tmpmja0rchn\\h2o_caima_started_from_python.out\n",
      "  JVM stderr: C:\\Users\\caima\\AppData\\Local\\Temp\\tmpmja0rchn\\h2o_caima_started_from_python.err\n",
      "  Server is running at http://127.0.0.1:54321\n",
      "Connecting to H2O server at http://127.0.0.1:54321 ... successful.\n",
      "Warning: Your H2O cluster version is too old (3 years, 5 months and 27 days)! Please download and install the latest version from http://h2o.ai/download/\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div style=\"overflow:auto\"><table style=\"width:50%\"><tr><td>H2O_cluster_uptime:</td>\n",
       "<td>01 secs</td></tr>\n",
       "<tr><td>H2O_cluster_timezone:</td>\n",
       "<td>America/Bogota</td></tr>\n",
       "<tr><td>H2O_data_parsing_timezone:</td>\n",
       "<td>UTC</td></tr>\n",
       "<tr><td>H2O_cluster_version:</td>\n",
       "<td>3.30.0.4</td></tr>\n",
       "<tr><td>H2O_cluster_version_age:</td>\n",
       "<td>3 years, 5 months and 27 days !!!</td></tr>\n",
       "<tr><td>H2O_cluster_name:</td>\n",
       "<td>H2O_from_python_caima_jto53b</td></tr>\n",
       "<tr><td>H2O_cluster_total_nodes:</td>\n",
       "<td>1</td></tr>\n",
       "<tr><td>H2O_cluster_free_memory:</td>\n",
       "<td>2.926 Gb</td></tr>\n",
       "<tr><td>H2O_cluster_total_cores:</td>\n",
       "<td>0</td></tr>\n",
       "<tr><td>H2O_cluster_allowed_cores:</td>\n",
       "<td>0</td></tr>\n",
       "<tr><td>H2O_cluster_status:</td>\n",
       "<td>accepting new members, healthy</td></tr>\n",
       "<tr><td>H2O_connection_url:</td>\n",
       "<td>http://127.0.0.1:54321</td></tr>\n",
       "<tr><td>H2O_connection_proxy:</td>\n",
       "<td>{\"http\": null, \"https\": null}</td></tr>\n",
       "<tr><td>H2O_internal_security:</td>\n",
       "<td>False</td></tr>\n",
       "<tr><td>H2O_API_Extensions:</td>\n",
       "<td>Amazon S3, Algos, AutoML, Core V3, TargetEncoder, Core V4</td></tr>\n",
       "<tr><td>Python_version:</td>\n",
       "<td>3.10.9 final</td></tr></table></div>"
      ],
      "text/plain": [
       "--------------------------  ---------------------------------------------------------\n",
       "H2O_cluster_uptime:         01 secs\n",
       "H2O_cluster_timezone:       America/Bogota\n",
       "H2O_data_parsing_timezone:  UTC\n",
       "H2O_cluster_version:        3.30.0.4\n",
       "H2O_cluster_version_age:    3 years, 5 months and 27 days !!!\n",
       "H2O_cluster_name:           H2O_from_python_caima_jto53b\n",
       "H2O_cluster_total_nodes:    1\n",
       "H2O_cluster_free_memory:    2.926 Gb\n",
       "H2O_cluster_total_cores:    0\n",
       "H2O_cluster_allowed_cores:  0\n",
       "H2O_cluster_status:         accepting new members, healthy\n",
       "H2O_connection_url:         http://127.0.0.1:54321\n",
       "H2O_connection_proxy:       {\"http\": null, \"https\": null}\n",
       "H2O_internal_security:      False\n",
       "H2O_API_Extensions:         Amazon S3, Algos, AutoML, Core V3, TargetEncoder, Core V4\n",
       "Python_version:             3.10.9 final\n",
       "--------------------------  ---------------------------------------------------------"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h2o.init()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parse progress: |█████████████████████████████████████████████████████████| 100%\n"
     ]
    }
   ],
   "source": [
    "data_frame = h2o.H2OFrame(data_consol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 314,
   "metadata": {},
   "outputs": [],
   "source": [
    "response = 'LFP'\n",
    "numerical_feats = [\"WA\", \"WE\", \"HHRS\", \"HA\", \"HE\", \"HW\", \"WMED\", \"WFED\", \"UN\", \"AX\"]\n",
    "categorical_feats = [\"KL6\", \"K618\", \"CIT\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_frame[response] = data_frame[response].asfactor()\n",
    "data_frame[numerical_feats] = data_frame[numerical_feats].asnumeric()\n",
    "data_frame[categorical_feats] = data_frame[categorical_feats].asfactor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {},
   "outputs": [],
   "source": [
    "train, valid, test = data_frame.split_frame(\n",
    "    ratios = [0.6, 0.2],\n",
    "    seed = 1234\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 301,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*\n",
      "doc LFP  \"A dummy variable = 1 if woman worked in 1975, else 0\";\n",
      "doc WHRS \"Wife's hours of work in 1975\";\n",
      "doc KL6  \"Number of children less than 6 years old in household\";\n",
      "doc K618 \"Number of children between ages 6 and 18 in household\";\n",
      "doc WA   \"Wife's age\";\n",
      "doc WE   \"Wife's educational attainment, in years\";\n",
      "doc WW   \"Wife's average hourly earnings, in 1975 dollars\";\n",
      "doc RPWG \"Wife's wage reported at the time of the 1976 interview\\\n",
      " (not the same as the 1975 estimated wage).\\\n",
      " To use the subsample with this wage, one needs to select 1975\\\n",
      " workers with LFP=1, then select only those women with non-zero RPWG.\\\n",
      " Only 325 women work in 1975 and have a non-zero RPWG in 1976.\";\n",
      "doc HHRS \"Husband's hours worked in 1975\";\n",
      "doc HA   \"Husband's age\";\n",
      "doc HE   \"Husband's educational attainment, in years\";\n",
      "doc HW   \"Husband's wage, in 1975 dollars\";\n",
      "doc FAMINC \"Family income, in 1975 dollars.\\\n",
      " This variable is used to construct the property income variable.\";\n",
      "doc MTR  \"This is the marginal tax rate facing the wife, and is taken from\\\n",
      " published federal tax tables (state and local income taxes are excluded).\\\n",
      " The taxable income on which this tax rate is calculated\\\n",
      " includes Social Security, if applicable to wife.\";\n",
      "doc WMED \"Wife's mother's educational attainment, in years\";\n",
      "doc WFED \"Wife's father's educational attainment, in years\";\n",
      "doc UN   \"Unemployment rate in county of residence, in percentage points.\\\n",
      " This taken from bracketed ranges.\";\n",
      "doc CIT  \"Dummy variable = 1 if live in large city (SMSA), else 0\";\n",
      "doc AX   \"Actual years of wife's previous labor market experience\";\n",
      "?Source:  1976 Panel Study of Income Dynamics,\n",
      "?based on data for the previous year, 1975.  \n",
      "?Of the 753 observations, the first 428 are for women with positive hours\n",
      "?worked in 1975, while the remaining 325 observations are for women who\n",
      "?did not work for pay in 1975.  A more complete discussion of the data is\n",
      "?found in Mroz [1987], Appendix 1. \n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(data_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explaining Models with Partial Effects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{\\partial Pr(G=1|X=x)}{\\partial X_{d}} = Pr(G=1|X=x) ( 1 - Pr(G=1|X=x)) \\beta_{d}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Average Partial Effect* (APE) \n",
    "\n",
    "If the feature *d* is not dummy\n",
    "\n",
    "$$ APE = \\frac{1}{n} \\sum_{i=1}^{n} \\frac{\\partial Pr(G=1|X=x_i)}{\\partial X_{d}}$$\n",
    "\n",
    "\n",
    "but, if the feature *d* is dummy\n",
    "\n",
    "$$ APE = \\frac{1}{n} \\sum_{i=1}^{n} [Pr(G=1|X=x_{i,d}, x_d=1) - Pr(G=1|X = x_{i,d}, x_d=0)]$$\n",
    "\n",
    "Where $x_{i,d}$ means all features except $d$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Odds\n",
    "\n",
    "It is defined like $Odds = P(G=1|X=x)/P(G=0|X=x)$. \n",
    "\n",
    "For logistic $$P(G=1|X=x)/P(G=0|X=x) = \\exp(\\beta^{T}{x})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The effect of change a variable dummy on odds. *Odds Ratio*\n",
    "\n",
    "$$Odds \\: Ratio = \\frac{Odds(x_d=1)}{Odds(x_d=0 )} = \\exp(\\beta_{d})$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 348,
   "metadata": {},
   "outputs": [],
   "source": [
    "logit_model = H2OGeneralizedLinearEstimator(\n",
    "    family = 'binomial', \n",
    "    link = 'Logit',\n",
    "    compute_p_values= True,\n",
    "    # Specify if the intercept should be estimated or not\n",
    "    intercept = True,\n",
    "    # lambda_ controls the amount of regularization applied to the model.\n",
    "    lambda_ = 0,\n",
    "    # alpha = 0 controls the distribution between the ℓ1 [α]\n",
    "    # (LASSO) and ℓ2 [α -1] (ridge regression) penalties\n",
    "    # if α=1, Lasso, α = 0, Ridge\n",
    "    alpha = 1,\n",
    "    # Specify whether to automatically remove collinear columns during model-building.\n",
    "    #remove_collinear_columns = False,\n",
    "    solver = 'IRLSM',\n",
    "    #seed = 1234)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 349,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "glm Model Build progress: |"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "███████████████████████████████████████████████| 100%\n"
     ]
    }
   ],
   "source": [
    "logit_model.train(x=numerical_feats + categorical_feats, \n",
    "                  y = response, \n",
    "                  training_frame = data_frame,\n",
    "                  validation_frame=valid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 358,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ModelMetricsBinomialGLM: glm\n",
      "** Reported on test data. **\n",
      "\n",
      "MSE: 0.18289410470707004\n",
      "RMSE: 0.4276612031819932\n",
      "LogLoss: 0.5395201397908074\n",
      "Null degrees of freedom: 152\n",
      "Residual degrees of freedom: 130\n",
      "Null deviance: 213.61636782004203\n",
      "Residual deviance: 165.09316277598705\n",
      "AIC: 211.09316277598705\n",
      "AUC: 0.8072186110160794\n",
      "AUCPR: 0.8256290179944596\n",
      "Gini: 0.6144372220321588\n",
      "\n",
      "Confusion Matrix (Act/Pred) for max f1 @ threshold = 0.5836558570900445: \n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>Error</th>\n",
       "      <th>Rate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>53.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>0.2838</td>\n",
       "      <td>(21.0/74.0)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>18.0</td>\n",
       "      <td>61.0</td>\n",
       "      <td>0.2278</td>\n",
       "      <td>(18.0/79.0)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Total</td>\n",
       "      <td>71.0</td>\n",
       "      <td>82.0</td>\n",
       "      <td>0.2549</td>\n",
       "      <td>(39.0/153.0)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             0     1   Error           Rate\n",
       "0      0  53.0  21.0  0.2838    (21.0/74.0)\n",
       "1      1  18.0  61.0  0.2278    (18.0/79.0)\n",
       "2  Total  71.0  82.0  0.2549   (39.0/153.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Maximum Metrics: Maximum metrics at their respective thresholds\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>metric</th>\n",
       "      <th>threshold</th>\n",
       "      <th>value</th>\n",
       "      <th>idx</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>max f1</td>\n",
       "      <td>0.583656</td>\n",
       "      <td>0.757764</td>\n",
       "      <td>81.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>max f2</td>\n",
       "      <td>0.212131</td>\n",
       "      <td>0.864745</td>\n",
       "      <td>134.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>max f0point5</td>\n",
       "      <td>0.772923</td>\n",
       "      <td>0.784314</td>\n",
       "      <td>43.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>max accuracy</td>\n",
       "      <td>0.601834</td>\n",
       "      <td>0.745098</td>\n",
       "      <td>73.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>max precision</td>\n",
       "      <td>0.972718</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>max recall</td>\n",
       "      <td>0.119228</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>144.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>max specificity</td>\n",
       "      <td>0.972718</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>max absolute_mcc</td>\n",
       "      <td>0.772923</td>\n",
       "      <td>0.499336</td>\n",
       "      <td>43.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>max min_per_class_accuracy</td>\n",
       "      <td>0.588075</td>\n",
       "      <td>0.743243</td>\n",
       "      <td>77.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>max mean_per_class_accuracy</td>\n",
       "      <td>0.601834</td>\n",
       "      <td>0.745895</td>\n",
       "      <td>73.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>max tns</td>\n",
       "      <td>0.972718</td>\n",
       "      <td>74.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>max fns</td>\n",
       "      <td>0.972718</td>\n",
       "      <td>78.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>max fps</td>\n",
       "      <td>0.015524</td>\n",
       "      <td>74.000000</td>\n",
       "      <td>152.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>max tps</td>\n",
       "      <td>0.119228</td>\n",
       "      <td>79.000000</td>\n",
       "      <td>144.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>max tnr</td>\n",
       "      <td>0.972718</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>max fnr</td>\n",
       "      <td>0.972718</td>\n",
       "      <td>0.987342</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>max fpr</td>\n",
       "      <td>0.015524</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>152.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>max tpr</td>\n",
       "      <td>0.119228</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>144.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         metric  threshold      value    idx\n",
       "0                        max f1   0.583656   0.757764   81.0\n",
       "1                        max f2   0.212131   0.864745  134.0\n",
       "2                  max f0point5   0.772923   0.784314   43.0\n",
       "3                  max accuracy   0.601834   0.745098   73.0\n",
       "4                 max precision   0.972718   1.000000    0.0\n",
       "5                    max recall   0.119228   1.000000  144.0\n",
       "6               max specificity   0.972718   1.000000    0.0\n",
       "7              max absolute_mcc   0.772923   0.499336   43.0\n",
       "8    max min_per_class_accuracy   0.588075   0.743243   77.0\n",
       "9   max mean_per_class_accuracy   0.601834   0.745895   73.0\n",
       "10                      max tns   0.972718  74.000000    0.0\n",
       "11                      max fns   0.972718  78.000000    0.0\n",
       "12                      max fps   0.015524  74.000000  152.0\n",
       "13                      max tps   0.119228  79.000000  144.0\n",
       "14                      max tnr   0.972718   1.000000    0.0\n",
       "15                      max fnr   0.972718   0.987342    0.0\n",
       "16                      max fpr   0.015524   1.000000  152.0\n",
       "17                      max tpr   0.119228   1.000000  144.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Gains/Lift Table: Avg response rate: 51.63 %, avg score: 37.79 %\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>group</th>\n",
       "      <th>cumulative_data_fraction</th>\n",
       "      <th>lower_threshold</th>\n",
       "      <th>lift</th>\n",
       "      <th>cumulative_lift</th>\n",
       "      <th>response_rate</th>\n",
       "      <th>score</th>\n",
       "      <th>cumulative_response_rate</th>\n",
       "      <th>cumulative_score</th>\n",
       "      <th>capture_rate</th>\n",
       "      <th>cumulative_capture_rate</th>\n",
       "      <th>gain</th>\n",
       "      <th>cumulative_gain</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td></td>\n",
       "      <td>1</td>\n",
       "      <td>0.013072</td>\n",
       "      <td>0.639338</td>\n",
       "      <td>0.968354</td>\n",
       "      <td>0.968354</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.676793</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.676793</td>\n",
       "      <td>0.012658</td>\n",
       "      <td>0.012658</td>\n",
       "      <td>-3.164557</td>\n",
       "      <td>-3.164557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td></td>\n",
       "      <td>2</td>\n",
       "      <td>0.026144</td>\n",
       "      <td>0.594426</td>\n",
       "      <td>0.968354</td>\n",
       "      <td>0.968354</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.600844</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.638818</td>\n",
       "      <td>0.012658</td>\n",
       "      <td>0.025316</td>\n",
       "      <td>-3.164557</td>\n",
       "      <td>-3.164557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td></td>\n",
       "      <td>3</td>\n",
       "      <td>0.032680</td>\n",
       "      <td>0.586369</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.774684</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.590882</td>\n",
       "      <td>0.400000</td>\n",
       "      <td>0.629231</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.025316</td>\n",
       "      <td>-100.000000</td>\n",
       "      <td>-22.531646</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td></td>\n",
       "      <td>4</td>\n",
       "      <td>0.045752</td>\n",
       "      <td>0.568982</td>\n",
       "      <td>0.968354</td>\n",
       "      <td>0.830018</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.575920</td>\n",
       "      <td>0.428571</td>\n",
       "      <td>0.614000</td>\n",
       "      <td>0.012658</td>\n",
       "      <td>0.037975</td>\n",
       "      <td>-3.164557</td>\n",
       "      <td>-16.998192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td></td>\n",
       "      <td>5</td>\n",
       "      <td>0.052288</td>\n",
       "      <td>0.566186</td>\n",
       "      <td>1.936709</td>\n",
       "      <td>0.968354</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.568569</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.608321</td>\n",
       "      <td>0.012658</td>\n",
       "      <td>0.050633</td>\n",
       "      <td>93.670886</td>\n",
       "      <td>-3.164557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td></td>\n",
       "      <td>6</td>\n",
       "      <td>0.104575</td>\n",
       "      <td>0.541213</td>\n",
       "      <td>0.242089</td>\n",
       "      <td>0.605222</td>\n",
       "      <td>0.125000</td>\n",
       "      <td>0.555786</td>\n",
       "      <td>0.312500</td>\n",
       "      <td>0.582053</td>\n",
       "      <td>0.012658</td>\n",
       "      <td>0.063291</td>\n",
       "      <td>-75.791139</td>\n",
       "      <td>-39.477848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td></td>\n",
       "      <td>7</td>\n",
       "      <td>0.150327</td>\n",
       "      <td>0.480997</td>\n",
       "      <td>1.106691</td>\n",
       "      <td>0.757843</td>\n",
       "      <td>0.571429</td>\n",
       "      <td>0.509726</td>\n",
       "      <td>0.391304</td>\n",
       "      <td>0.560041</td>\n",
       "      <td>0.050633</td>\n",
       "      <td>0.113924</td>\n",
       "      <td>10.669078</td>\n",
       "      <td>-24.215740</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td></td>\n",
       "      <td>8</td>\n",
       "      <td>0.202614</td>\n",
       "      <td>0.430170</td>\n",
       "      <td>1.210443</td>\n",
       "      <td>0.874643</td>\n",
       "      <td>0.625000</td>\n",
       "      <td>0.455192</td>\n",
       "      <td>0.451613</td>\n",
       "      <td>0.532983</td>\n",
       "      <td>0.063291</td>\n",
       "      <td>0.177215</td>\n",
       "      <td>21.044304</td>\n",
       "      <td>-12.535729</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td></td>\n",
       "      <td>9</td>\n",
       "      <td>0.300654</td>\n",
       "      <td>0.400928</td>\n",
       "      <td>0.774684</td>\n",
       "      <td>0.842047</td>\n",
       "      <td>0.400000</td>\n",
       "      <td>0.415778</td>\n",
       "      <td>0.434783</td>\n",
       "      <td>0.494764</td>\n",
       "      <td>0.075949</td>\n",
       "      <td>0.253165</td>\n",
       "      <td>-22.531646</td>\n",
       "      <td>-15.795267</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td></td>\n",
       "      <td>10</td>\n",
       "      <td>0.398693</td>\n",
       "      <td>0.370549</td>\n",
       "      <td>0.903797</td>\n",
       "      <td>0.857232</td>\n",
       "      <td>0.466667</td>\n",
       "      <td>0.385842</td>\n",
       "      <td>0.442623</td>\n",
       "      <td>0.467980</td>\n",
       "      <td>0.088608</td>\n",
       "      <td>0.341772</td>\n",
       "      <td>-9.620253</td>\n",
       "      <td>-14.276821</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td></td>\n",
       "      <td>11</td>\n",
       "      <td>0.503268</td>\n",
       "      <td>0.351449</td>\n",
       "      <td>1.089399</td>\n",
       "      <td>0.905474</td>\n",
       "      <td>0.562500</td>\n",
       "      <td>0.361992</td>\n",
       "      <td>0.467532</td>\n",
       "      <td>0.445956</td>\n",
       "      <td>0.113924</td>\n",
       "      <td>0.455696</td>\n",
       "      <td>8.939873</td>\n",
       "      <td>-9.452573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td></td>\n",
       "      <td>12</td>\n",
       "      <td>0.601307</td>\n",
       "      <td>0.335017</td>\n",
       "      <td>1.291139</td>\n",
       "      <td>0.968354</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.342174</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.429035</td>\n",
       "      <td>0.126582</td>\n",
       "      <td>0.582278</td>\n",
       "      <td>29.113924</td>\n",
       "      <td>-3.164557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td></td>\n",
       "      <td>13</td>\n",
       "      <td>0.699346</td>\n",
       "      <td>0.321140</td>\n",
       "      <td>1.420253</td>\n",
       "      <td>1.031705</td>\n",
       "      <td>0.733333</td>\n",
       "      <td>0.327986</td>\n",
       "      <td>0.532710</td>\n",
       "      <td>0.414870</td>\n",
       "      <td>0.139241</td>\n",
       "      <td>0.721519</td>\n",
       "      <td>42.025316</td>\n",
       "      <td>3.170472</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td></td>\n",
       "      <td>14</td>\n",
       "      <td>0.797386</td>\n",
       "      <td>0.307851</td>\n",
       "      <td>1.291139</td>\n",
       "      <td>1.063602</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.313431</td>\n",
       "      <td>0.549180</td>\n",
       "      <td>0.402398</td>\n",
       "      <td>0.126582</td>\n",
       "      <td>0.848101</td>\n",
       "      <td>29.113924</td>\n",
       "      <td>6.360241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td></td>\n",
       "      <td>15</td>\n",
       "      <td>0.895425</td>\n",
       "      <td>0.287210</td>\n",
       "      <td>0.774684</td>\n",
       "      <td>1.031969</td>\n",
       "      <td>0.400000</td>\n",
       "      <td>0.299359</td>\n",
       "      <td>0.532847</td>\n",
       "      <td>0.391116</td>\n",
       "      <td>0.075949</td>\n",
       "      <td>0.924051</td>\n",
       "      <td>-22.531646</td>\n",
       "      <td>3.196896</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td></td>\n",
       "      <td>16</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.224614</td>\n",
       "      <td>0.726266</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.375000</td>\n",
       "      <td>0.265115</td>\n",
       "      <td>0.516340</td>\n",
       "      <td>0.377939</td>\n",
       "      <td>0.075949</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-27.373418</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      group  cumulative_data_fraction  lower_threshold      lift  \\\n",
       "0         1                  0.013072         0.639338  0.968354   \n",
       "1         2                  0.026144         0.594426  0.968354   \n",
       "2         3                  0.032680         0.586369  0.000000   \n",
       "3         4                  0.045752         0.568982  0.968354   \n",
       "4         5                  0.052288         0.566186  1.936709   \n",
       "5         6                  0.104575         0.541213  0.242089   \n",
       "6         7                  0.150327         0.480997  1.106691   \n",
       "7         8                  0.202614         0.430170  1.210443   \n",
       "8         9                  0.300654         0.400928  0.774684   \n",
       "9        10                  0.398693         0.370549  0.903797   \n",
       "10       11                  0.503268         0.351449  1.089399   \n",
       "11       12                  0.601307         0.335017  1.291139   \n",
       "12       13                  0.699346         0.321140  1.420253   \n",
       "13       14                  0.797386         0.307851  1.291139   \n",
       "14       15                  0.895425         0.287210  0.774684   \n",
       "15       16                  1.000000         0.224614  0.726266   \n",
       "\n",
       "    cumulative_lift  response_rate     score  cumulative_response_rate  \\\n",
       "0          0.968354       0.500000  0.676793                  0.500000   \n",
       "1          0.968354       0.500000  0.600844                  0.500000   \n",
       "2          0.774684       0.000000  0.590882                  0.400000   \n",
       "3          0.830018       0.500000  0.575920                  0.428571   \n",
       "4          0.968354       1.000000  0.568569                  0.500000   \n",
       "5          0.605222       0.125000  0.555786                  0.312500   \n",
       "6          0.757843       0.571429  0.509726                  0.391304   \n",
       "7          0.874643       0.625000  0.455192                  0.451613   \n",
       "8          0.842047       0.400000  0.415778                  0.434783   \n",
       "9          0.857232       0.466667  0.385842                  0.442623   \n",
       "10         0.905474       0.562500  0.361992                  0.467532   \n",
       "11         0.968354       0.666667  0.342174                  0.500000   \n",
       "12         1.031705       0.733333  0.327986                  0.532710   \n",
       "13         1.063602       0.666667  0.313431                  0.549180   \n",
       "14         1.031969       0.400000  0.299359                  0.532847   \n",
       "15         1.000000       0.375000  0.265115                  0.516340   \n",
       "\n",
       "    cumulative_score  capture_rate  cumulative_capture_rate        gain  \\\n",
       "0           0.676793      0.012658                 0.012658   -3.164557   \n",
       "1           0.638818      0.012658                 0.025316   -3.164557   \n",
       "2           0.629231      0.000000                 0.025316 -100.000000   \n",
       "3           0.614000      0.012658                 0.037975   -3.164557   \n",
       "4           0.608321      0.012658                 0.050633   93.670886   \n",
       "5           0.582053      0.012658                 0.063291  -75.791139   \n",
       "6           0.560041      0.050633                 0.113924   10.669078   \n",
       "7           0.532983      0.063291                 0.177215   21.044304   \n",
       "8           0.494764      0.075949                 0.253165  -22.531646   \n",
       "9           0.467980      0.088608                 0.341772   -9.620253   \n",
       "10          0.445956      0.113924                 0.455696    8.939873   \n",
       "11          0.429035      0.126582                 0.582278   29.113924   \n",
       "12          0.414870      0.139241                 0.721519   42.025316   \n",
       "13          0.402398      0.126582                 0.848101   29.113924   \n",
       "14          0.391116      0.075949                 0.924051  -22.531646   \n",
       "15          0.377939      0.075949                 1.000000  -27.373418   \n",
       "\n",
       "    cumulative_gain  \n",
       "0         -3.164557  \n",
       "1         -3.164557  \n",
       "2        -22.531646  \n",
       "3        -16.998192  \n",
       "4         -3.164557  \n",
       "5        -39.477848  \n",
       "6        -24.215740  \n",
       "7        -12.535729  \n",
       "8        -15.795267  \n",
       "9        -14.276821  \n",
       "10        -9.452573  \n",
       "11        -3.164557  \n",
       "12         3.170472  \n",
       "13         6.360241  \n",
       "14         3.196896  \n",
       "15         0.000000  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": []
     },
     "execution_count": 358,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logit_model.model_performance(test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 359,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients_table = logit_model._model_json['output']['coefficients_table']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 360,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients_table_frame = coefficients_table.as_data_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 361,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>names</th>\n",
       "      <th>coefficients</th>\n",
       "      <th>std_error</th>\n",
       "      <th>z_value</th>\n",
       "      <th>p_value</th>\n",
       "      <th>standardized_coefficients</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Intercept</td>\n",
       "      <td>2.5344</td>\n",
       "      <td>1.0200</td>\n",
       "      <td>2.4845</td>\n",
       "      <td>0.0130</td>\n",
       "      <td>0.6365</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>K618.1</td>\n",
       "      <td>-0.1613</td>\n",
       "      <td>0.2463</td>\n",
       "      <td>-0.6551</td>\n",
       "      <td>0.5124</td>\n",
       "      <td>-0.1613</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>K618.2</td>\n",
       "      <td>0.1576</td>\n",
       "      <td>0.2697</td>\n",
       "      <td>0.5844</td>\n",
       "      <td>0.5590</td>\n",
       "      <td>0.1576</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>K618.3</td>\n",
       "      <td>0.1879</td>\n",
       "      <td>0.3037</td>\n",
       "      <td>0.6186</td>\n",
       "      <td>0.5362</td>\n",
       "      <td>0.1879</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>K618.4</td>\n",
       "      <td>0.0734</td>\n",
       "      <td>0.4733</td>\n",
       "      <td>0.1550</td>\n",
       "      <td>0.8768</td>\n",
       "      <td>0.0734</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>K618.5</td>\n",
       "      <td>0.8141</td>\n",
       "      <td>0.6780</td>\n",
       "      <td>1.2009</td>\n",
       "      <td>0.2298</td>\n",
       "      <td>0.8141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>K618.6</td>\n",
       "      <td>-8.2713</td>\n",
       "      <td>167.0575</td>\n",
       "      <td>-0.0495</td>\n",
       "      <td>0.9605</td>\n",
       "      <td>-8.2713</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>K618.7</td>\n",
       "      <td>-11.2847</td>\n",
       "      <td>167.0571</td>\n",
       "      <td>-0.0676</td>\n",
       "      <td>0.9461</td>\n",
       "      <td>-11.2847</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>K618.8</td>\n",
       "      <td>8.2326</td>\n",
       "      <td>166.5827</td>\n",
       "      <td>0.0494</td>\n",
       "      <td>0.9606</td>\n",
       "      <td>8.2326</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>KL6.1</td>\n",
       "      <td>-1.5073</td>\n",
       "      <td>0.2620</td>\n",
       "      <td>-5.7523</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>-1.5073</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>KL6.2</td>\n",
       "      <td>-2.8945</td>\n",
       "      <td>0.5366</td>\n",
       "      <td>-5.3936</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>-2.8945</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>KL6.3</td>\n",
       "      <td>-10.7237</td>\n",
       "      <td>93.5426</td>\n",
       "      <td>-0.1146</td>\n",
       "      <td>0.9087</td>\n",
       "      <td>-10.7237</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>CIT.1</td>\n",
       "      <td>0.0788</td>\n",
       "      <td>0.1982</td>\n",
       "      <td>0.3974</td>\n",
       "      <td>0.6911</td>\n",
       "      <td>0.0788</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>WA</td>\n",
       "      <td>-0.0903</td>\n",
       "      <td>0.0253</td>\n",
       "      <td>-3.5633</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>-0.7288</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>WE</td>\n",
       "      <td>0.2517</td>\n",
       "      <td>0.0539</td>\n",
       "      <td>4.6705</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.5740</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>HHRS</td>\n",
       "      <td>-0.0004</td>\n",
       "      <td>0.0002</td>\n",
       "      <td>-2.7450</td>\n",
       "      <td>0.0061</td>\n",
       "      <td>-0.2617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>HA</td>\n",
       "      <td>-0.0092</td>\n",
       "      <td>0.0243</td>\n",
       "      <td>-0.3792</td>\n",
       "      <td>0.7046</td>\n",
       "      <td>-0.0741</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>HE</td>\n",
       "      <td>-0.0457</td>\n",
       "      <td>0.0396</td>\n",
       "      <td>-1.1528</td>\n",
       "      <td>0.2490</td>\n",
       "      <td>-0.1380</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>HW</td>\n",
       "      <td>-0.0494</td>\n",
       "      <td>0.0246</td>\n",
       "      <td>-2.0095</td>\n",
       "      <td>0.0445</td>\n",
       "      <td>-0.2089</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>WMED</td>\n",
       "      <td>0.0136</td>\n",
       "      <td>0.0332</td>\n",
       "      <td>0.4098</td>\n",
       "      <td>0.6819</td>\n",
       "      <td>0.0458</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>WFED</td>\n",
       "      <td>0.0013</td>\n",
       "      <td>0.0312</td>\n",
       "      <td>0.0401</td>\n",
       "      <td>0.9680</td>\n",
       "      <td>0.0045</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>UN</td>\n",
       "      <td>-0.0265</td>\n",
       "      <td>0.0286</td>\n",
       "      <td>-0.9279</td>\n",
       "      <td>0.3535</td>\n",
       "      <td>-0.0826</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>AX</td>\n",
       "      <td>0.1216</td>\n",
       "      <td>0.0139</td>\n",
       "      <td>8.7589</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.9816</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        names  coefficients  std_error  z_value  p_value  \\\n",
       "0   Intercept        2.5344     1.0200   2.4845   0.0130   \n",
       "1      K618.1       -0.1613     0.2463  -0.6551   0.5124   \n",
       "2      K618.2        0.1576     0.2697   0.5844   0.5590   \n",
       "3      K618.3        0.1879     0.3037   0.6186   0.5362   \n",
       "4      K618.4        0.0734     0.4733   0.1550   0.8768   \n",
       "5      K618.5        0.8141     0.6780   1.2009   0.2298   \n",
       "6      K618.6       -8.2713   167.0575  -0.0495   0.9605   \n",
       "7      K618.7      -11.2847   167.0571  -0.0676   0.9461   \n",
       "8      K618.8        8.2326   166.5827   0.0494   0.9606   \n",
       "9       KL6.1       -1.5073     0.2620  -5.7523   0.0000   \n",
       "10      KL6.2       -2.8945     0.5366  -5.3936   0.0000   \n",
       "11      KL6.3      -10.7237    93.5426  -0.1146   0.9087   \n",
       "12      CIT.1        0.0788     0.1982   0.3974   0.6911   \n",
       "13         WA       -0.0903     0.0253  -3.5633   0.0004   \n",
       "14         WE        0.2517     0.0539   4.6705   0.0000   \n",
       "15       HHRS       -0.0004     0.0002  -2.7450   0.0061   \n",
       "16         HA       -0.0092     0.0243  -0.3792   0.7046   \n",
       "17         HE       -0.0457     0.0396  -1.1528   0.2490   \n",
       "18         HW       -0.0494     0.0246  -2.0095   0.0445   \n",
       "19       WMED        0.0136     0.0332   0.4098   0.6819   \n",
       "20       WFED        0.0013     0.0312   0.0401   0.9680   \n",
       "21         UN       -0.0265     0.0286  -0.9279   0.3535   \n",
       "22         AX        0.1216     0.0139   8.7589   0.0000   \n",
       "\n",
       "    standardized_coefficients  \n",
       "0                      0.6365  \n",
       "1                     -0.1613  \n",
       "2                      0.1576  \n",
       "3                      0.1879  \n",
       "4                      0.0734  \n",
       "5                      0.8141  \n",
       "6                     -8.2713  \n",
       "7                    -11.2847  \n",
       "8                      8.2326  \n",
       "9                     -1.5073  \n",
       "10                    -2.8945  \n",
       "11                   -10.7237  \n",
       "12                     0.0788  \n",
       "13                    -0.7288  \n",
       "14                     0.5740  \n",
       "15                    -0.2617  \n",
       "16                    -0.0741  \n",
       "17                    -0.1380  \n",
       "18                    -0.2089  \n",
       "19                     0.0458  \n",
       "20                     0.0045  \n",
       "21                    -0.0826  \n",
       "22                     0.9816  "
      ]
     },
     "execution_count": 361,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coefficients_table_frame.round(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 334,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "GLM Model: summary\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>family</th>\n",
       "      <th>link</th>\n",
       "      <th>regularization</th>\n",
       "      <th>number_of_predictors_total</th>\n",
       "      <th>number_of_active_predictors</th>\n",
       "      <th>number_of_iterations</th>\n",
       "      <th>training_frame</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td></td>\n",
       "      <td>binomial</td>\n",
       "      <td>logit</td>\n",
       "      <td>None</td>\n",
       "      <td>22</td>\n",
       "      <td>22</td>\n",
       "      <td>9</td>\n",
       "      <td>py_10_sid_96c0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       family   link regularization  number_of_predictors_total  \\\n",
       "0    binomial  logit           None                          22   \n",
       "\n",
       "  number_of_active_predictors  number_of_iterations  training_frame  \n",
       "0                          22                     9  py_10_sid_96c0  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": []
     },
     "execution_count": 334,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logit_model._model_json[\"output\"]['model_summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "data_drime = pd.read_stata(r'data_dta/CRIME1.dta')\n",
    "\n",
    "# Generate the variable arr86, where a person is labeled equal to 0, if he has not committed a crime in 1986\n",
    "# otherwise 1.\n",
    "data_drime['arr86'] = data_drime.narr86.where(data_drime.narr86==0, 1)\n",
    "\n",
    "# Selecting the features and target\n",
    "features= ['pcnv', 'avgsen', 'tottime', 'ptime86', 'qemp86']\n",
    "target = 'arr86'\n",
    "\n",
    "X = np.array(data_drime[features])\n",
    "\n",
    "# Generate polynomial and interaction features\n",
    "poly = PolynomialFeatures(2)\n",
    "X_poly = poly.fit_transform(X)\n",
    "\n",
    "pca = PCA(n_components=X_poly.shape[1])\n",
    "X_trans = pca.fit_transform(X_poly)\n",
    "condition = pca.explained_variance_ratio_.cumsum() <= 0.997\n",
    "X_pca = X_trans[:,condition]\n",
    "\n",
    "y = np.array(data_drime[target])\n",
    "\n",
    "# Split the data to train and test\n",
    "X_train, X_test, y_train, y_test = model_selection.train_test_split(X_pca, y, random_state=123, test_size=0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGhCAYAAACzurT/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAk0UlEQVR4nO3dcXSU1Z3/8c8kQyYpOMNCdCAS06lFjUYLTGwMFMVVphtXl9iuhLqCtNY2K2pjareNaRfN2WPcPV0re9ZEaWFd2sqmXZSV01CcXVsNjbYlDepWbKlQJ8LENGHPDNo6keT+/uDH1OkkIU9MuJnwfp3znMNzn3uf5zs5YeaTO8/ccRljjAAAACzJsl0AAAA4vRFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFjltl3AaAwODurw4cM644wz5HK5bJcDAABGwRijo0ePqqCgQFlZw89/ZEQYOXz4sAoLC22XAQAAxqCrq0vz5s0b9nhGhJEzzjhD0vEH4/V6LVcDAABGIx6Pq7CwMPk6PpyMCCMn3prxer2EEQAAMszJbrHgBlYAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFa5bRcwEVyzJu7c5sjEnRsAgNMRMyMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwKoxhZGmpiYFAgHl5uYqGAyqra1t2L5r166Vy+VK2y666KIxFw0AAKYOx2GkpaVFNTU1qq+vV2dnp5YuXaqKigpFIpEh+2/YsEHRaDS5dXV1adasWbrhhhved/EAACDzuYwxxsmAsrIyLVq0SM3Nzcm24uJiVVZWqrGx8aTjt2/frk984hM6ePCgioqKRnXNeDwun8+nWCwmr9d70v6sMwIAgH2jff12NDPS39+vjo4OhUKhlPZQKKT29vZRnWPTpk26+uqrRx1EAADA1OZoBdbe3l4NDAzI7/entPv9fnV3d590fDQa1c6dO/X444+P2C+RSCiRSCT34/G4kzIBAEAGGdMNrC6XK2XfGJPWNpTHHntMM2fOVGVl5Yj9Ghsb5fP5klthYeFYygQAABnAURjJz89XdnZ22ixIT09P2mzJnzLGaPPmzVq9erVycnJG7FtXV6dYLJbcurq6nJQJAAAyiKMwkpOTo2AwqHA4nNIeDoe1ePHiEcc+++yz+s1vfqNbbrnlpNfxeDzyer0pGwAAmJocf2tvbW2tVq9erdLSUpWXl2vjxo2KRCKqrq6WdHxW49ChQ9qyZUvKuE2bNqmsrEwlJSXjUzkAAJgSHIeRqqoq9fX1qaGhQdFoVCUlJWptbU1+OiYajaatORKLxbRt2zZt2LBhfKoGAABThuN1RmxgnREAADLPhKwzAgAAMN4IIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwKoxhZGmpiYFAgHl5uYqGAyqra1txP6JREL19fUqKiqSx+PRueeeq82bN4+pYAAAMLW4nQ5oaWlRTU2NmpqatGTJEj366KOqqKjQK6+8onPOOWfIMStXrtSbb76pTZs26cMf/rB6enp07Nix9108AADIfC5jjHEyoKysTIsWLVJzc3Oyrbi4WJWVlWpsbEzr/8Mf/lCrVq3SgQMHNGvWrDEVGY/H5fP5FIvF5PV6T9rfNbbLjIo5MnHnBgBgKhnt67ejt2n6+/vV0dGhUCiU0h4KhdTe3j7kmKeeekqlpaX6p3/6J5199tk677zzdPfdd+sPf/jDsNdJJBKKx+MpGwAAmJocvU3T29urgYEB+f3+lHa/36/u7u4hxxw4cEC7d+9Wbm6unnzySfX29uq2227TkSNHhr1vpLGxUffdd5+T0gAAQIYa0w2sLpcrZd8Yk9Z2wuDgoFwul7773e/qox/9qK655ho9+OCDeuyxx4adHamrq1MsFktuXV1dYykTAABkAEczI/n5+crOzk6bBenp6UmbLTlh7ty5Ovvss+Xz+ZJtxcXFMsbojTfe0Pz589PGeDweeTweJ6UBAIAM5WhmJCcnR8FgUOFwOKU9HA5r8eLFQ45ZsmSJDh8+rLfeeivZ9utf/1pZWVmaN2/eGEoGAABTieO3aWpra/Wtb31Lmzdv1r59+3TXXXcpEomourpa0vG3WNasWZPsf+ONN2r27Nn69Kc/rVdeeUXPPfecvvSlL+kzn/mM8vLyxu+RAACAjOR4nZGqqir19fWpoaFB0WhUJSUlam1tVVFRkSQpGo0qEokk+8+YMUPhcFh33HGHSktLNXv2bK1cuVL/8A//MH6PAgAAZCzH64zYwDojAABknglZZwQAAGC8EUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVYwojTU1NCgQCys3NVTAYVFtb27B9f/zjH8vlcqVtr7766piLBgAAU4fjMNLS0qKamhrV19ers7NTS5cuVUVFhSKRyIjjfvWrXykajSa3+fPnj7loAAAwdTgOIw8++KBuueUWffazn1VxcbEeeughFRYWqrm5ecRxZ511lubMmZPcsrOzx1w0AACYOhyFkf7+fnV0dCgUCqW0h0Ihtbe3jzh24cKFmjt3rq666ir96Ec/GrFvIpFQPB5P2QAAwNTkKIz09vZqYGBAfr8/pd3v96u7u3vIMXPnztXGjRu1bds2PfHEEzr//PN11VVX6bnnnhv2Oo2NjfL5fMmtsLDQSZkAACCDuMcyyOVypewbY9LaTjj//PN1/vnnJ/fLy8vV1dWlr3/967r88suHHFNXV6fa2trkfjweJ5AAADBFOZoZyc/PV3Z2dtosSE9PT9psyUguu+wy7d+/f9jjHo9HXq83ZQMAAFOTozCSk5OjYDCocDic0h4Oh7V48eJRn6ezs1Nz5851cmkAADBFOX6bpra2VqtXr1ZpaanKy8u1ceNGRSIRVVdXSzr+FsuhQ4e0ZcsWSdJDDz2kD37wg7rooovU39+v73znO9q2bZu2bds2vo8EAABkJMdhpKqqSn19fWpoaFA0GlVJSYlaW1tVVFQkSYpGoylrjvT39+vuu+/WoUOHlJeXp4suukg/+MEPdM0114zfowAAABnLZYwxtos4mXg8Lp/Pp1gsNqr7R1yzJq4Wc2Tizg0AwFQy2tdvvpsGAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYNaYw0tTUpEAgoNzcXAWDQbW1tY1q3E9+8hO53W4tWLBgLJcFAABTkOMw0tLSopqaGtXX16uzs1NLly5VRUWFIpHIiONisZjWrFmjq666aszFAgCAqcdljDFOBpSVlWnRokVqbm5OthUXF6uyslKNjY3Djlu1apXmz5+v7Oxsbd++XXv37h31NePxuHw+n2KxmLxe70n7u2aN+tSOmSMTd24AAKaS0b5+O5oZ6e/vV0dHh0KhUEp7KBRSe3v7sOP+7d/+Ta+99prWr1/v5HIAAOA04HbSube3VwMDA/L7/Sntfr9f3d3dQ47Zv3+/vvKVr6itrU1u9+gul0gklEgkkvvxeNxJmQAAIIOM6QZWl8uVsm+MSWuTpIGBAd1444267777dN555436/I2NjfL5fMmtsLBwLGUCAIAM4CiM5OfnKzs7O20WpKenJ222RJKOHj2qPXv26Pbbb5fb7Zbb7VZDQ4NefPFFud1uPfPMM0Nep66uTrFYLLl1dXU5KRMAAGQQR2/T5OTkKBgMKhwO6/rrr0+2h8NhrVixIq2/1+vVyy+/nNLW1NSkZ555Rv/5n/+pQCAw5HU8Ho88Ho+T0gAAQIZyFEYkqba2VqtXr1ZpaanKy8u1ceNGRSIRVVdXSzo+q3Ho0CFt2bJFWVlZKikpSRl/1llnKTc3N60dAACcnhyHkaqqKvX19amhoUHRaFQlJSVqbW1VUVGRJCkajZ50zREAAIATHK8zYgPrjAAAkHkmZJ0RAACA8UYYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVo0pjDQ1NSkQCCg3N1fBYFBtbW3D9t29e7eWLFmi2bNnKy8vTxdccIG+8Y1vjLlgAAAwtbidDmhpaVFNTY2ampq0ZMkSPfroo6qoqNArr7yic845J63/9OnTdfvtt+uSSy7R9OnTtXv3bn3+85/X9OnT9bnPfW5cHgQAAMhcLmOMcTKgrKxMixYtUnNzc7KtuLhYlZWVamxsHNU5PvGJT2j69On69re/Par+8XhcPp9PsVhMXq/3pP1ds0Z12jExRybu3AAATCWjff129DZNf3+/Ojo6FAqFUtpDoZDa29tHdY7Ozk61t7friiuuGLZPIpFQPB5P2QAAwNTkKIz09vZqYGBAfr8/pd3v96u7u3vEsfPmzZPH41FpaanWrVunz372s8P2bWxslM/nS26FhYVOygQAABlkTDewulyulH1jTFrbn2pra9OePXv0yCOP6KGHHtLWrVuH7VtXV6dYLJbcurq6xlImAADIAI5uYM3Pz1d2dnbaLEhPT0/abMmfCgQCkqSLL75Yb775pu6991596lOfGrKvx+ORx+NxUhoAAMhQjmZGcnJyFAwGFQ6HU9rD4bAWL1486vMYY5RIJJxcGgAATFGOP9pbW1ur1atXq7S0VOXl5dq4caMikYiqq6slHX+L5dChQ9qyZYsk6eGHH9Y555yjCy64QNLxdUe+/vWv64477hjHhwEAADKV4zBSVVWlvr4+NTQ0KBqNqqSkRK2trSoqKpIkRaNRRSKRZP/BwUHV1dXp4MGDcrvdOvfcc/XAAw/o85///Pg9CgAAkLEcrzNiA+uMAACQeSZknREAAIDxRhgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWjSmMNDU1KRAIKDc3V8FgUG1tbcP2feKJJ7R8+XKdeeaZ8nq9Ki8v165du8ZcMAAAmFoch5GWlhbV1NSovr5enZ2dWrp0qSoqKhSJRIbs/9xzz2n58uVqbW1VR0eHrrzySl133XXq7Ox838UDAIDM5zLGGCcDysrKtGjRIjU3NyfbiouLVVlZqcbGxlGd46KLLlJVVZX+/u//flT94/G4fD6fYrGYvF7vSfu7Zo3qtGNijkzMeTOxZgAARjLa129HMyP9/f3q6OhQKBRKaQ+FQmpvbx/VOQYHB3X06FHNmjX8q28ikVA8Hk/ZAADA1OQojPT29mpgYEB+vz+l3e/3q7u7e1Tn+Od//me9/fbbWrly5bB9Ghsb5fP5klthYaGTMgEAQAYZ0w2sLpcrZd8Yk9Y2lK1bt+ree+9VS0uLzjrrrGH71dXVKRaLJbeurq6xlAkAADKA20nn/Px8ZWdnp82C9PT0pM2W/KmWlhbdcsst+v73v6+rr756xL4ej0cej8dJaQAAIEM5mhnJyclRMBhUOBxOaQ+Hw1q8ePGw47Zu3aq1a9fq8ccf11/+5V+OrVIAADAlOZoZkaTa2lqtXr1apaWlKi8v18aNGxWJRFRdXS3p+Fsshw4d0pYtWyQdDyJr1qzRhg0bdNlllyVnVfLy8uTz+cbxoQAAgEzkOIxUVVWpr69PDQ0NikajKikpUWtrq4qKiiRJ0Wg0Zc2RRx99VMeOHdO6deu0bt26ZPvNN9+sxx577P0/AgAAkNEcrzNiA+uMvD+sMwIAsGFC1hkBAAAYb4QRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVWMKI01NTQoEAsrNzVUwGFRbW9uwfaPRqG688Uadf/75ysrKUk1NzVhrBQAAU5DjMNLS0qKamhrV19ers7NTS5cuVUVFhSKRyJD9E4mEzjzzTNXX1+sjH/nI+y4YAABMLS5jjHEyoKysTIsWLVJzc3Oyrbi4WJWVlWpsbBxx7LJly7RgwQI99NBDjoqMx+Py+XyKxWLyer0n7e+a5ej0jpgjE3PeTKwZAICRjPb129HMSH9/vzo6OhQKhVLaQ6GQ2tvbx1YpAAA4rbmddO7t7dXAwID8fn9Ku9/vV3d397gVlUgklEgkkvvxeHzczg0AACaXMd3A6nK5UvaNMWlt70djY6N8Pl9yKywsHLdzAwCAycVRGMnPz1d2dnbaLEhPT0/abMn7UVdXp1gslty6urrG7dwAAGBycRRGcnJyFAwGFQ6HU9rD4bAWL148bkV5PB55vd6UDQAATE2O7hmRpNraWq1evVqlpaUqLy/Xxo0bFYlEVF1dLen4rMahQ4e0ZcuW5Ji9e/dKkt566y397ne/0969e5WTk6MLL7xwfB4FAADIWI7DSFVVlfr6+tTQ0KBoNKqSkhK1traqqKhI0vFFzv50zZGFCxcm/93R0aHHH39cRUVF+u1vf/v+qgcAABnP8TojNrDOyPvDOiMAABsmZJ0RAACA8UYYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFjltl0AMpdr1sSd2xyZuHMDACYXZkYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVfGsvTit80zAATD7MjAAAAKsIIwAAwKoxhZGmpiYFAgHl5uYqGAyqra1txP7PPvusgsGgcnNz9aEPfUiPPPLImIoFAABTj+Mw0tLSopqaGtXX16uzs1NLly5VRUWFIpHIkP0PHjyoa665RkuXLlVnZ6fuuece3Xnnndq2bdv7Lh4AAGQ+lzHGOBlQVlamRYsWqbm5OdlWXFysyspKNTY2pvX/8pe/rKeeekr79u1LtlVXV+vFF1/U888/P6prxuNx+Xw+xWIxeb3ek/bPxJsUqTkVNf9RJtYsTVzd3CgMZI7Rvn47mhnp7+9XR0eHQqFQSnsoFFJ7e/uQY55//vm0/h//+Me1Z88evfvuu04uDwAApiBHH+3t7e3VwMCA/H5/Srvf71d3d/eQY7q7u4fsf+zYMfX29mru3LlpYxKJhBKJRHI/FotJOp6wRsXRXI8zoy3BMWpOQc3vkYk1SxNW90TW7CuamPPGXp+Y80rU/F7UnGqiapZGX/eJ1+2TvQkzpnVGXC5Xyr4xJq3tZP2Haj+hsbFR9913X1p7YWGh01LHnc9nuwLnqPnUoOZTg5pPDWo+NTKxZsl53UePHpVvhEGOwkh+fr6ys7PTZkF6enrSZj9OmDNnzpD93W63Zs+ePeSYuro61dbWJvcHBwd15MgRzZ49e8TQ41Q8HldhYaG6urpGdS/KZJCJNUuZWTc1nxrUfGpQ86lBzamMMTp69KgKCgpG7OcojOTk5CgYDCocDuv6669PtofDYa1YsWLIMeXl5dqxY0dK29NPP63S0lJNmzZtyDEej0cejyelbebMmU5KdcTr9WbML80JmVizlJl1U/OpQc2nBjWfGtT8RyPNiJzg+KO9tbW1+ta3vqXNmzdr3759uuuuuxSJRFRdXS3p+KzGmjVrkv2rq6v1+uuvq7a2Vvv27dPmzZu1adMm3X333U4vDQAApiDH94xUVVWpr69PDQ0NikajKikpUWtrq4qKjt8pE41GU9YcCQQCam1t1V133aWHH35YBQUF+pd/+Rd98pOfHL9HAQAAMtaYbmC97bbbdNtttw157LHHHktru+KKK/SLX/xiLJeaUB6PR+vXr097S2gyy8Sapcysm5pPDWo+Naj51KDmsXG86BkAAMB44ovyAACAVYQRAABgFWEEAABYRRgBYB23rgGntzF9miZTvfHGG2publZ7e7u6u7vlcrnk9/u1ePFiVVdXT4rl5oHTkcfj0Ysvvqji4mLbpcCiaDSq5uZm7d69W9FoVNnZ2QoEAqqsrNTatWuVnZ1tu0RMkNPm0zS7d+9WRUWFCgsLFQqF5Pf7ZYxRT0+PwuGwurq6tHPnTi1ZssR2qY50dXVp/fr12rx5s+1SUvzhD39QR0eHZs2apQsvvDDl2DvvvKPvfe97KYvjTQb79u3TCy+8oPLycl1wwQV69dVXtWHDBiUSCd1000368z//c9sljuj//u//9O///u/av3+/5s6dq5tvvnnSBez3fs3De23YsEE33XRT8isiHnzwwVNZ1og6Ozs1c+ZMBQIBSdJ3vvMdNTc3KxKJqKioSLfffrtWrVplucrMt2fPHl199dUKBALKy8vTT3/6U/3N3/yN+vv7tWvXLhUXF2vXrl0644wzbJea0e644w6tXLlSS5cutV1KKnOaKC0tNTU1NcMer6mpMaWlpaewovGxd+9ek5WVZbuMFL/61a9MUVGRcblcJisry1xxxRXm8OHDyePd3d2TruadO3eanJwcM2vWLJObm2t27txpzjzzTHP11Vebq666yrjdbvM///M/tstMMXfuXNPb22uMMebAgQNmzpw5Zs6cOWb58uVm3rx5xufzmX379lmuMpXL5TILFiwwy5YtS9lcLpe59NJLzbJly8yVV15pu8wUCxcuNM8884wxxphvfvObJi8vz9x5552mubnZ1NTUmBkzZphNmzZZrnJ4XV1d5ujRo2nt/f395tlnn7VQ0dCWLFli7r333uT+t7/9bVNWVmaMMebIkSNmwYIF5s4777RV3ojeeusts3HjRrN27VrzF3/xF6aiosKsXbvWfPOb3zRvvfWW7fJSnHhenj9/vnnggQdMNBq1XZIxxpjTJozk5uaaV199ddjj+/btM7m5uaewotH5r//6rxG3b3zjG5Puhb2ystJce+215ne/+53Zv3+/ue6660wgEDCvv/66MWZyhpHy8nJTX19vjDFm69at5s/+7M/MPffckzx+zz33mOXLl9sqb0gul8u8+eabxhhjVq1aZZYtW2befvttY4wx77zzjrn22mvNX//1X9ssMc39999vAoFAWrBzu93ml7/8paWqRvaBD3wg+bu7cOFC8+ijj6Yc/+53v2suvPBCG6WN6PDhw+bSSy81WVlZJjs726xZsyYllEy2/4d5eXnmtddeS+4PDAyYadOmme7ubmOMMU8//bQpKCiwVd6wfvnLX5qCggIzc+ZMs2LFCvO5z33O3HrrrWbFihVm5syZ5uyzz55Uv9sul8v893//t/nCF75g8vPzzbRp08xf/dVfmR07dpiBgQFrdZ02YSQQCJjNmzcPe3zz5s0mEAicwopG50SKdblcw26T6QnFGGPOOuss89JLL6W03Xbbbeacc84xr7322qR7EjTGGK/Xa/bv32+MOf4k6Ha7TUdHR/L4yy+/bPx+v63yhvTeMDLUC/wLL7xg5s2bZ6O0Ef3sZz8z5513nvniF79o+vv7jTGTO4zMnj3b7Nmzxxhz/Hd77969Kcd/85vfmLy8PBuljWjNmjXmsssuMz//+c9NOBw2paWlJhgMmiNHjhhjjocRl8tluco/KioqMrt3707uHz582LhcLvP73//eGGPMwYMHJ+UfjMuWLTOrVq0yiUQi7VgikTCf+tSnzLJlyyxUNrT3Pm/09/eblpYW8/GPf9xkZ2ebgoICc8899ySfC0+l0yaMPPzwwyYnJ8esW7fObN++3Tz//PPmhRdeMNu3bzfr1q0zHo/HNDc32y4zTUFBgXnyySeHPd7Z2TnpXtjPOOMM88orr6S133777WbevHnmueeem3Q1vzeMGGPMjBkzUv5K++1vfzvpnghdLpfp6ekxxhz/Pfnf//3flOMHDx40Ho/HRmkndfToUbNmzRpz8cUXm5deeslMmzZt0oaRm266ydxyyy3GGGNuuOEG89WvfjXl+P33328uvvhiG6WNqKCgwPz0pz9N7r/zzjtmxYoVZsGCBaavr2/S/VHwhS98wZSUlJidO3eaZ555xlx55ZUpL+I//OEPzbnnnmuxwqHl5eWN+Lv78ssvT6qw+t4w8l6vv/66Wb9+vSkqKrLye3HahBFjjPmP//gPU1ZWZtxud3JWwe12m7KyMtPS0mK7vCFdd9115mtf+9qwx/fu3Tup/roxxphLL73UbNmyZchj69atMzNnzpxUT4LGGHPJJZeYnTt3Jvdffvll8+677yb329raJt3MmcvlMhdffLFZuHChmTFjhnniiSdSjj/77LPm7LPPtlTd6GzdutX4/X6TlZU1acPIoUOHzAc/+EFz+eWXm9raWpOXl2c+9rGPmVtvvdVcfvnlJicnx/zgBz+wXWaa6dOnm1//+tcpbe+++66prKw0l1xyiXnppZcm1f/Do0ePmpUrVyafnxcvXmwOHDiQPL5r1y7zve99z2KFQysoKDDbt28f9viTTz45qd5eGi6MnDA4OGiefvrpU1jRcafVR3urqqpUVVWld999V729vZKk/Px8TZs2zXJlw/vSl76kt99+e9jjH/7wh/WjH/3oFFZ0ctdff722bt2q1atXpx3713/9Vw0ODuqRRx6xUNnw/vZv/1YDAwPJ/ZKSkpTjO3funHSfplm/fn3K/gc+8IGU/R07dky+O+b/xKpVq/Sxj31MHR0dyW/+nmwKCgrU2dmpBx54QDt27JAxRj/72c/U1dWlJUuW6Cc/+YlKS0ttl5nmQx/6kF566SXNnz8/2eZ2u/X9739fN9xwg6699lqL1aWbMWOGWlpa9M477+jYsWOaMWNGyvFQKGSpspHdeuutuvnmm/XVr35Vy5cvl9/vl8vlUnd3t8LhsO6//37V1NTYLjOpqKhoxI9Iu1wuLV++/BRW9P+va8zp8dFeADidfPnLX9bevXu1a9eutGPHjh3TJz/5Se3YsUODg4MWqpta/vEf/1EbNmxIrl8lHV/Ib86cOaqpqdHf/d3fWa5w8iOMAMAUdOzYMf3+97+X1+sd8vjAwIDeeOONSTsjlYkOHjyo7u5uSdKcOXOSa9Pg5FgOHgCmILfbPWwQkaTDhw/rvvvuO4UVTX2BQEDl5eUqLy9PBpGuri595jOfsVzZ5MfMCACchl588UUtWrQo5V4pjD9+zqNzWt3ACgCni6eeemrE4wcOHDhFlUxt/JzHBzMjADAFZWVlyeVyjfiNyC6Xi7/Y3yd+zuODe0YAYAqaO3eutm3bpsHBwSG3X/ziF7ZLnBL4OY8PwggATEHBYHDEF8KT/TWP0eHnPD64ZwQApqBMXDAxE/FzHh/cMwIAAKzibRoAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVf8P0xX/yJcb4H0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_drime.narr86.value_counts(normalize=True).plot(kind='bar', color='#0011ff')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Avoiding Harms in Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0828626352604234"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(-1*(np.array([0.1, 0.9])*np.log(np.array([0.9, 0.1]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot  as plt\n",
    "from sklearn.linear_model import LogisticRegression, LogisticRegressionCV\n",
    "from sklearn.preprocessing import PolynomialFeatures, StandardScaler\n",
    "from sklearn import model_selection, metrics\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.metrics import f1_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "data_drime = pd.read_stata(r'data_dta/CRIME1.dta')\n",
    "\n",
    "# Generate the variable arr86, where a person is labeled equal to 0, if he has not committed a crime in 1986\n",
    "# otherwise 1.\n",
    "data_drime['arr86'] = data_drime.narr86.where(data_drime.narr86==0, 1)\n",
    "\n",
    "# Selecting the features and target\n",
    "features= ['pcnv', 'avgsen', 'tottime', 'ptime86', 'qemp86']\n",
    "target = 'arr86'\n",
    "\n",
    "X = np.array(data_drime[features])\n",
    "y = np.array(data_drime[target])\n",
    "\n",
    "X_norm = (X - X.mean(axis=0, keepdims=True))/X.std(axis=0, keepdims=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = model_selection.train_test_split(\n",
    "                X_norm, y, test_size=0.25)\n",
    "logit_model = LogisticRegression(solver='newton-cg')\n",
    "logit_model.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy 0.7243\n",
      "Initial AUC value is 0.6703\n",
      "\n",
      "Probability Conditioned on Predicts\n",
      "> P(not committed a crime|low risk)   72.32%\n",
      "> P(not committed a crime|high risk)  20.00%\n",
      "> P(committed a crime|low risk)       27.68%\n",
      "> P(committed a crime|high risk)      80.00%\n",
      "\n",
      "Probability Conditioned on Outcomes\n",
      "> P(low risk|not committed a crime)   99.59%\n",
      "> P(high risk|not committed a crime)  0.41%\n",
      "> P(low risk|committed a crime)       95.88%\n",
      "> P(high risk|committed a crime)      4.12%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "utils_metrics.assessment(y_test, X_test, logit_model)\n",
    "utils_metrics.prob_cond_pred(y_test, X_test, logit_model)\n",
    "utils_metrics.prob_cond_out(y_test, X_test, logit_model)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Validation for Logit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=5, estimator=LogisticRegression(solver='liblinear'),\n",
       "             param_grid={'C': array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),\n",
       "                         'penalty': ['l1', 'l2']},\n",
       "             return_train_score=True, scoring='roc_auc', verbose=False)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scaler = StandardScaler()\n",
    "\n",
    "fit_to_grid_cv = model_selection.GridSearchCV(\n",
    "    # estimator=Pipeline(steps=[(\"scaler\", scaler),\\\n",
    "    #     (\"logistic\", LogisticRegression( class_weight=\"balanced\",solver='liblinear'))]),\n",
    "    estimator = LogisticRegression(solver='liblinear'),#  class_weight=\"balanced\"\n",
    "    scoring= 'roc_auc',\n",
    "    cv=5,\n",
    "    param_grid={\n",
    "        'penalty': [\"l1\", \"l2\"],\n",
    "        'C': np.linspace(0.1,1, 10)}, \n",
    "    verbose=False,\n",
    "    return_train_score=True\n",
    ")\n",
    "\n",
    "fit_to_grid_cv.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy 0.7258\n",
      "Initial AUC value is 0.4437\n",
      "\n",
      "Probability Conditioned on Predicts\n",
      "> P(not committed a crime|low risk)   72.93%\n",
      "> P(not committed a crime|high risk)  66.67%\n",
      "> P(committed a crime|low risk)       27.07%\n",
      "> P(committed a crime|high risk)      33.33%\n",
      "\n",
      "Probability Conditioned on Outcomes\n",
      "> P(low risk|not committed a crime)   99.20%\n",
      "> P(high risk|not committed a crime)  0.80%\n",
      "> P(low risk|committed a crime)       98.92%\n",
      "> P(high risk|committed a crime)      1.08%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "logit_model = LogisticRegression(solver='liblinear', C=0.02, penalty=\"l2\")\n",
    "logit_model.fit(X_train, y_train)\n",
    "\n",
    "utils_metrics.assessment(y_test, X_test, logit_model)\n",
    "utils_metrics.prob_cond_pred(y_test, X_test, logit_model)\n",
    "utils_metrics.prob_cond_out(y_test, X_test, logit_model)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical Significance Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot  as plt\n",
    "from sklearn.linear_model import LogisticRegression, LogisticRegressionCV\n",
    "from sklearn.preprocessing import PolynomialFeatures, StandardScaler\n",
    "from sklearn import model_selection, metrics\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.metrics import f1_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "data_drime = pd.read_stata(r'data_dta/CRIME1.dta')\n",
    "\n",
    "# Generate the variable arr86, where a person is labeled equal to 0, if he has not committed a crime in 1986\n",
    "# otherwise 1.\n",
    "data_drime['arr86'] = data_drime.narr86.where(data_drime.narr86==0, 1)\n",
    "\n",
    "# Selecting the features and target\n",
    "features= ['pcnv', 'avgsen', 'tottime', 'ptime86', 'qemp86']\n",
    "target = 'arr86'\n",
    "\n",
    "X = np.array(data_drime[features])\n",
    "y = np.array(data_drime[target])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `old model` is defined by normalization on features and the `new model` is defined by the generation of the polynomios of the normalized features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [OLD MODEL] Normalization \n",
    "X_norm = (X - X.mean(axis=0, keepdims=True))/X.std(axis=0, keepdims=True)\n",
    "\n",
    "# [NEW MODEL] Generating polynomial and interaction features \n",
    "poly = PolynomialFeatures(2)\n",
    "X_poly = poly.fit_transform(X_norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f1_test_data_model(X:np.ndarray):\n",
    "    X_train, X_test, y_train, y_test = model_selection.train_test_split(\n",
    "                X, y, test_size=0.25)\n",
    "    logit_model = LogisticRegression(solver='newton-cg')\n",
    "    logit_model.fit(X_train, y_train)\n",
    "    f1 = f1_score(y_test, logit_model.predict(X_test))\n",
    "    return f1, y_test, X_test, logit_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "δ:  0.0861244019138756\n"
     ]
    }
   ],
   "source": [
    "f1_old, y_test_old, X_test_old, lm_old = f1_test_data_model(X_norm)\n",
    "f1_new, y_test_new, X_test_new, lm_new = f1_test_data_model(X_poly)\n",
    "δ_inint = f1_new - f1_old\n",
    "print(\"\\u03B4: \",δ_inint)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\delta_\\text{init}$ score was computed using $(y_{test}^{new}, X_{test}^{new})$ and $(y_{test}^{old}, X_{test}^{old})$\n",
    "$$\\delta_\\text{init} = F_{1}^\\text{new} - F_{1}^\\text{old} = 0.077185$$\n",
    "\n",
    "The hypothesis are $$H_0: \\delta<=0$$ $$H_1: \\delta>0$$\n",
    "\n",
    "Given this dataset we resample it using *bootstrap technique*. This process will result in a distribution of $\\delta$. With this distribution we can calculate the $\\text{p-value}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f1_score_by_idx_model(X, y, idx, model):\n",
    "    predict = model.predict(X[idx,:])\n",
    "    return f1_score(y, predict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# size of the dataset\n",
    "size = len(X_test_new)\n",
    "f1_score_boots = []\n",
    "x = 0\n",
    "while x < 1000:\n",
    "    # resampling of the dataset with replacement\n",
    "    idx = np.random.choice(range(size), size, replace=True)\n",
    "    \n",
    "    # compute the f_1 score having as an input the ressampling dataset\n",
    "    f1_score_boots_new = f1_score_by_idx_model(X_test_new, y_test_new, idx, lm_new)\n",
    "    f1_score_boots_old = f1_score_by_idx_model(X_test_old, y_test_old, idx, lm_old)\n",
    "    f1_score_boots.append({\n",
    "        'f1_score_boots_new': f1_score_boots_new,\n",
    "        'f1_score_boots_old': f1_score_boots_old\n",
    "    })\n",
    "    x += 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert to dataframe\n",
    "df_f1_score_boots = pd.DataFrame(f1_score_boots)\n",
    "# calculate the δ\n",
    "δ = df_f1_score_boots.eval('f1_score_boots_new - f1_score_boots_old')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resampling was generate under $\\delta_\\text{init}>0$, but to compute the $\\text{p-value}$ we need to keep as true the $h_0$, that is, $\\delta<=0$, but $\\delta_\\text{init}$ point out the otherwise, so we need to transform (only by localization) the $\\delta$ generade by resampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x26532e8a940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAEvCAYAAABmC5raAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbnElEQVR4nO3de5BdZZnv8e9DAkQwCjGNZggxjRVQFGmghVFEUWQECsFgqVAjh5lBAio4lAflJtKeKkqdET0zNTNouCh6uDpcZ4QRpBRqqgyYQAwgMASN0hAgBiXOCGjic/7Yq+MmdKd30r32ervz/VTt2mu/e71rPW+vneLHukZmIkmSVLKtmi5AkiRpNAYWSZJUPAOLJEkqnoFFkiQVz8AiSZKKZ2CRJEnFm1rXgiNiF+BbwGuAPwILM/MfImIGcDUwF1gBfCgzf131OQs4AVgHfDIzv7exdcycOTPnzp1b1xAkSVKXLVmy5FeZ2bNhe9R1H5aImAXMysx7ImI6sAR4P/BXwDOZ+cWIOBPYMTPPiIg9gCuB/YA/A74P7JaZ60ZaR39/fy5evLiW+iVJUvdFxJLM7N+wvbZDQpm5MjPvqaZ/CzwI7AwcBVxWzXYZrRBD1X5VZr6QmT8HltMKL5IkaQvXlXNYImIusDdwF/DqzFwJrVAD7FTNtjPwWFu3wapNkiRt4WoPLBHxcuBa4LTMXLOxWYdpe8nxqohYEBGLI2LxqlWrxqtMSZJUsNpOugWIiK1phZXLM/O6qvmpiJiVmSur81yertoHgV3aus8GnthwmZm5EFgIrXNYaitekqQu+8Mf/sDg4CDPP/9806XUbtq0acyePZutt966o/nrvEoogEuABzPzK21f3QQcD3yxer+xrf2KiPgKrZNu5wF311WfJEmlGRwcZPr06cydO5fWf0Ynp8xk9erVDA4O0tvb21GfOvewHAAcB9wXEUurtrNpBZVrIuIE4JfABwEy84GIuAb4KbAW+MTGrhCSJGmyef755yd9WAGICF71qlexKad21BZYMvM/Gf68FICDR+hzPnB+XTVJklS6yR5WhmzqOL3TrSRJKl6tJ91KkqTNNzAw0PXlnXvuuVx33XVss8029PX18Y1vfGNca9hcBhZJkrTehRdeyOLFi5k7dy7PPfdc0+Ws5yEhSZK03sc//nGOOOIIDjnkEG699damy1nPPSyT3DjvTWxsHVK3jfeu+KbWIW2Kxx9/nEWLFrFs2TJ+85vf8Pa3v5158+axxx57NF2ae1gkSVLLDTfcwAEHHMBWW23FjBkz+MAHPsDtt9/edFmAgUWSJFXWrl3L2rVr13/OTDLLuKm8h4Q0Zh52kjaPh51Umne+850cd9xxnHXWWWQm119/PZdffnnTZQEGFkmSitXtwNnX18eCBQvYf//9yUxOOukk+vr6ulrDSAwskiRpvVNPPZVTTz216TJewnNYJElS8QwskiSpeAYWSZJUPAOLJEkqnoFFkiQVz8AiSZKKZ2CRJEnFM7BIkqQXefbZZ5k/fz777rsve+65JxdffHHTJXnjOEmSSjXeN7rtdHnXXnst06dPZ8mSJQA899xz41vIZnAPiyRJepF99tmHO+64g/7+fs477zy23XbbpksysEiSpD959tln+cxnPsOyZctYtGgRP/jBD7jxxhs3a1m//vWvx60uA4skSVrv61//Ou9973t55StfydSpU3nrW9/Kk08++ZL5PvrRj76k7dxzzwXgC1/4Aqeffjof+9jHxq0uA4skSVrv3nvv5Y1vfOOLPu+555489thjnHjiiZx++uncdNNNLF++nHPOOYf58+cD8OSTT7J27VruuusurrzySnp7e1mxYgVf/vKXx6UuA4skSVpvxx135N577wXgu9/9LmvWrOFtb3sbDz30ENtssw2f/OQnmTFjBocddhjnn38+22+/PdAKNn19fey2224cdNBB7L///nzkIx/h9NNPH5e6vEpIkiSt9+lPf5oPf/jDXHXVVfT29nLdddex1VZbccghh7DLLrtwyimnsP/++6/fszJlyhQAli5dytFHH83SpUvZa6+9WLZsGXvttde41VVbYImIS4EjgKcz801V29XA7tUsOwC/ycy+iJgLPAg8XH23KDNPrqs2SZImgvG+rLkTvb293H333S9pP+OMM1i3bh1z5sxh+fLl7L777vzqV7+ip6cHgOXLlzNv3jxuueUWDjzwQB5//HEuvvhiZs6cyRve8IYx11XnHpZvAv8EfGuoITM/PDQdERcAz7bN/2hm9tVYjyRJ2kxf+tKXXtI2c+bM9eeoXHLJJQCcdtppAOy7774ceeSR47b+2gJLZt5Z7Tl5iYgI4EPAu+tavyRJmjyaOun2QOCpzHykra03Iu6NiDsi4sCROkbEgohYHBGLV61aVX+lkiSpcU0FlmOBK9s+rwTmZObewKeAKyLiFcN1zMyFmdmfmf1Dx80kSdLk1vXAEhFTgaOBq4faMvOFzFxdTS8BHgV263ZtkiSpTE3sYXkP8FBmDg41RERPREyppncF5gE/a6A2SZIalZlNl9AVmzrO2gJLRFwJ/AjYPSIGI+KE6qtjePHhIIB3AMsi4ifAvwInZ+YzddUmSVKJpk2bxurVqyd9aMlMVq9ezbRp0zruU+dVQseO0P5Xw7RdC1xbVy2SJE0Es2fPZnBwkC3hopJp06Yxe/bsjuf3TreSJBVi6623pre3t+kyiuSzhCRJUvEMLJIkqXgGFkmSVDwDiyRJKp6BRZIkFc/AIkmSimdgkSRJxTOwSJKk4hlYJElS8QwskiSpeAYWSZJUPAOLJEkqnoFFkiQVz8AiSZKKZ2CRJEnFM7BIkqTiGVgkSVLxDCySJKl4BhZJklQ8A4skSSqegUWSJBXPwCJJkopnYJEkScWbWteCI+JS4Ajg6cx8U9U2AJwIrKpmOzszb66+Ows4AVgHfDIzv1dXbZp4BgYmxzqkbhvowg+7G+uQ6tzD8k3g0GHav5qZfdVrKKzsARwDvLHq8y8RMaXG2iRJ0gRSW2DJzDuBZzqc/Sjgqsx8ITN/DiwH9qurNkmSNLE0cQ7LKRGxLCIujYgdq7adgcfa5hms2iRJkroeWC4EXgf0ASuBC6r2GGbeHG4BEbEgIhZHxOJVq1YNN4skSZpkuhpYMvOpzFyXmX8ELuJPh30GgV3aZp0NPDHCMhZmZn9m9vf09NRbsCRJKkJXA0tEzGr7OB+4v5q+CTgmIraNiF5gHnB3N2uTJEnlqvOy5iuBg4CZETEInAccFBF9tA73rABOAsjMByLiGuCnwFrgE5m5rq7aJEnSxFJbYMnMY4dpvmQj858PnF9XPZIkaeLyTreSJKl4BhZJklQ8A4skSSqegUWSJBXPwCJJkopnYJEkScUzsEiSpOIZWCRJUvEMLJIkqXgGFkmSVDwDiyRJKp6BRZIkFc/AIkmSimdgkSRJxTOwSJKk4hlYJElS8QwskiSpeAYWSZJUPAOLJEkq3tSmC9iSDQw0XYE0MQ34j0fa4riHRZIkFc/AIkmSimdgkSRJxTOwSJKk4hlYJElS8WoLLBFxaUQ8HRH3t7X9fUQ8FBHLIuL6iNihap8bEc9FxNLq9bW66pIkSRNPnXtYvgkcukHbbcCbMvPNwH8BZ7V992hm9lWvk2usS5IkTTC1BZbMvBN4ZoO2WzNzbfVxETC7rvVLkqTJo8lzWP4GuKXtc29E3BsRd0TEgSN1iogFEbE4IhavWrWq/iolSVLjGgksEXEOsBa4vGpaCczJzL2BTwFXRMQrhuubmQszsz8z+3t6erpTsCRJalTXA0tEHA8cAfxlZiZAZr6Qmaur6SXAo8Bu3a5NkiSVqavPEoqIQ4EzgHdm5u/a2nuAZzJzXUTsCswDftbN2qRuPJ7GR+BoMurGs518fpRqCywRcSVwEDAzIgaB82hdFbQtcFtEACyqrgh6B/B/ImItsA44OTOfGXbBkiRpi1NbYMnMY4dpvmSEea8Frq2rFkmSNLF5p1tJklQ8A4skSSpeR4ElIt5UdyGSJEkj6XQPy9ci4u6I+PjQ838kSZK6paPAkplvB/4S2AVYHBFXRMQhtVYmSZJU6fgclsx8BPgs1X1UgH+snrx8dF3FSZIkQefnsLw5Ir4KPAi8G3hfZr6hmv5qjfVJkiR1fB+WfwIuAs7OzOeGGjPziYj4bC2VSZIkVToNLIcDz2XmOoCI2AqYlpm/y8xv11adJEkSnZ/D8n3gZW2ft6vaJEmSatdpYJmWmf899KGa3q6ekiRJkl6s08DyPxGxz9CHiNgXeG4j80uSJI2bTs9hOQ34TkQ8UX2eBXy4lookSZI20FFgycwfR8Trgd2BAB7KzD/UWpkkSVKl0z0sAG8B5lZ99o4IMvNbtVQlSZLUpqPAEhHfBl4HLAXWVc0JGFgkSVLtOt3D0g/skZlZZzGSJEnD6fQqofuB19RZiCRJ0kg63cMyE/hpRNwNvDDUmJlH1lKVJElSm04Dy0CdRUiSJG1Mp5c13xERrwXmZeb3I2I7YEq9pUmSJLV0dA5LRJwI/Cvw9appZ+CGmmqSJEl6kU5Puv0EcACwBiAzHwF2qqsoSZKkdp0Glhcy8/dDHyJiKq37sEiSJNWu08ByR0ScDbwsIg4BvgP828Y6RMSlEfF0RNzf1jYjIm6LiEeq9x3bvjsrIpZHxMMR8d7NGYwkSZqcOg0sZwKrgPuAk4Cbgc+O0uebwKHDLOf2zJwH3F59JiL2AI4B3lj1+ZeI8KReSZIEdH6V0B+Bi6pXRzLzzoiYu0HzUcBB1fRlwA+BM6r2qzLzBeDnEbEc2A/4UafrkyRJk1enzxL6OcOcs5KZu27i+l6dmSurvisjYujE3Z2BRW3zDVZtkiRJm/QsoSHTgA8CM8axjhimbdiTeiNiAbAAYM6cOeNYgiRJKlVH57Bk5uq21+OZ+X+Bd2/G+p6KiFkA1fvTVfsgsEvbfLOBJ0aoZWFm9mdmf09Pz2aUIEmSJppObxy3T9urPyJOBqZvxvpuAo6vpo8HbmxrPyYito2IXmAecPdmLF+SJE1CnR4SuqBtei2wAvjQxjpExJW0TrCdGRGDwHnAF4FrIuIE4Je0Di2RmQ9ExDXAT6vlfyIz13U+DEmSNJl1epXQuzZ1wZl57AhfHTzC/OcD52/qeiRJ0uTX6VVCn9rY95n5lfEpR5Ik6aU25Sqht9A61wTgfcCdwGN1FCVJktSu08AyE9gnM38LEBEDwHcy86N1FSZJkjSk01vzzwF+3/b598Dcca9GkiRpGJ3uYfk2cHdEXE/rhm7zgW/VVpUkSVKbTq8SOj8ibgEOrJr+OjPvra8sSZKkP+n0kBDAdsCazPwHYLC6wZskSVLtOr3T7Xm0nqp8VtW0NfD/6ipKkiSpXad7WOYDRwL/A5CZT7B5t+aXJEnaZJ2edPv7zMyISICI2L7GmoowMNB0BdLENOA/Hkk16HQPyzUR8XVgh4g4Efg+cFF9ZUmSJP3JqHtYIiKAq4HXA2uA3YHPZeZtNdcmSZIEdBBYqkNBN2TmvoAhRZIkdV2nh4QWRcRbaq1EkiRpBJ2edPsu4OSIWEHrSqGgtfPlzXUVJkmSNGSjgSUi5mTmL4HDulSPJEnSS4y2h+UGWk9p/kVEXJuZH+hCTZIkSS8y2jks0Ta9a52FSJIkjWS0wJIjTEuSJHXNaIeE9oqINbT2tLysmoY/nXT7ilqrkyRJYpTAkplTulWIJEnSSDq9D4skSVJjDCySJKl4BhZJklQ8A4skSSpep7fmHzcRsTutpz8P2RX4HLADcCKwqmo/OzNv7m51kiSpRF0PLJn5MNAHEBFTgMeB64G/Br6amV/udk2SJKlsTR8SOhh4NDN/0XAdkiSpYE0HlmOAK9s+nxIRyyLi0ojYcbgOEbEgIhZHxOJVq1YNN4skSZpkGgssEbENcCTwnarpQuB1tA4XrQQuGK5fZi7MzP7M7O/p6elGqZIkqWFdP4elzWHAPZn5FMDQO0BEXAT8e1OFSXUZGJgc65C6baALP+xurEObr8lDQsfSdjgoIma1fTcfuL/rFUmSpCI1soclIrYDDgFOamv+u4joo/VU6BUbfCdJkrZgjQSWzPwd8KoN2o5rohZJklS+pq8SkiRJGpWBRZIkFc/AIkmSimdgkSRJxTOwSJKk4hlYJElS8QwskiSpeAYWSZJUPAOLJEkqnoFFkiQVz8AiSZKKZ2CRJEnFM7BIkqTiGVgkSVLxDCySJKl4BhZJklQ8A4skSSqegUWSJBXPwCJJkopnYJEkScUzsEiSpOIZWCRJUvEMLJIkqXgGFkmSVLypTaw0IlYAvwXWAWszsz8iZgBXA3OBFcCHMvPXTdQnSZLK0uQelndlZl9m9lefzwRuz8x5wO3VZ0mSpKIOCR0FXFZNXwa8v7lSJElSSZoKLAncGhFLImJB1fbqzFwJUL3v1FBtkiSpMI2cwwIckJlPRMROwG0R8VCnHauAswBgzpw5ddUnSZIK0sgelsx8onp/Grge2A94KiJmAVTvT4/Qd2Fm9mdmf09PT7dKliRJDep6YImI7SNi+tA08BfA/cBNwPHVbMcDN3a7NkmSVKYmDgm9Grg+IobWf0Vm/kdE/Bi4JiJOAH4JfLCB2iRJUoG6Hlgy82fAXsO0rwYO7nY9kiSpfCVd1ixJkjQsA4skSSqegUWSJBXPwCJJkopnYJEkScUzsEiSpOIZWCRJUvEMLJIkqXgGFkmSVDwDiyRJKp6BRZIkFc/AIkmSimdgkSRJxTOwSJKk4hlYJElS8QwskiSpeAYWSZJUPAOLJEkqnoFFkiQVz8AiSZKKZ2CRJEnFM7BIkqTiGVgkSVLxDCySJKl4BhZJklS8rgeWiNglIn4QEQ9GxAMR8bdV+0BEPB4RS6vX4d2uTZIklWlqA+tcC/zvzLwnIqYDSyLituq7r2bmlxuoSZo0BgYmxzqkbhvowg+7G+uYrLoeWDJzJbCymv5tRDwI7NztOiRJ0sTR6DksETEX2Bu4q2o6JSKWRcSlEbHjCH0WRMTiiFi8atWqbpUqSZIa1FhgiYiXA9cCp2XmGuBC4HVAH609MBcM1y8zF2Zmf2b29/T0dKtcSZLUoEYCS0RsTSusXJ6Z1wFk5lOZuS4z/whcBOzXRG2SJKk8TVwlFMAlwIOZ+ZW29llts80H7u92bZIkqUxNXCV0AHAccF9ELK3azgaOjYg+IIEVwEkN1CZJkgrUxFVC/wnEMF/d3O1aJEnSxOCdbiVJUvEMLJIkqXgGFkmSVDwDiyRJKp6BRZIkFc/AIkmSimdgkSRJxTOwSJKk4hlYJElS8QwskiSpeAYWSZJUPAOLJEkqXhNPa5Y0wQ0MjPzdD3940Lis46CDfjguy5FKMrCxfzwTaB1NcA+LJEkqnoFFkiQVz8AiSZKKZ2CRJEnFM7BIkqTiGVgkSVLxDCySJKl43odFUpHG634uG+O9XjQZTdZ7vbiHRZIkFc/AIkmSimdgkSRJxSsusETEoRHxcEQsj4gzm65HkiQ1r6jAEhFTgH8GDgP2AI6NiD2arUqSJDWtqMAC7Acsz8yfZebvgauAoxquSZIkNay0wLIz8Fjb58GqTZIkbcFKuw9LDNOWL5ohYgGwoPr43xHxcO1V1W8m8Kumi+gixzu5TZjx3nHHuCxmwox3nGxp44Utb8yjjvfzn/98net/7XCNpQWWQWCXts+zgSfaZ8jMhcDCbhZVt4hYnJn9TdfRLY53cnO8k9uWNl7Y8sZc6nhLOyT0Y2BeRPRGxDbAMcBNDdckSZIaVtQelsxcGxGnAN8DpgCXZuYDDZclSZIaVlRgAcjMm4Gbm66jyybVIa4OON7JzfFOblvaeGHLG3OR443MHH0uSZKkBpV2DoskSdJLGFhqEhEzIuK2iHiket9xhPmGfRRBRFwdEUur14qIWFq1z42I59q++1qXhjSqcRjzQEQ83ja2w9u+O6ua/+GIeG83xjOacRjv30fEQxGxLCKuj4gdqvZitvFoj8qIln+svl8WEfuM1rfTv1tTNnfMEbFLRPwgIh6MiAci4m/b+oz4227aGLfxioi4rxrT4rb2YrfxGLbv7m3bb2lErImI06rvJvL2fX1E/CgiXoiI0zvp29j2zUxfNbyAvwPOrKbPBL40zDxTgEeBXYFtgJ8Aewwz3wXA56rpucD9TY+vjjEDA8Dpw/TZo5pvW6C36j9lEoz3L4Cp1fSXhvqXso07+X0ChwO30LqH0p8Dd3Uw7lH/bhN0zLOAfarp6cB/jfbbbvo1lvFW360AZg6z3CK38VjHu8FyngReOwm2707AW4Dz28dQ4r9h97DU5yjgsmr6MuD9w8wz6qMIIiKADwFX1lfquBmXMY+w3Ksy84XM/DmwvFpO08Y03sy8NTPXVvMtonXfoZJ0sq2OAr6VLYuAHSJi1ih9O/m7NWWzx5yZKzPzHoDM/C3wIOXfqXss23hjSt3G4zXeg4FHM/MX9Zc8JqOONzOfzswfA3/YhL6NbF8DS31enZkrAar3nYaZp5NHERwIPJWZj7S19UbEvRFxR0QcOJ5Fj9F4jPmUajfspW27GUt9ZMN4bWOAv6H1f3VDStjGndQ+0jwb69vJ360pYxnzehExF9gbuKutebjfdtPGOt4Ebo2IJdG6C/mQUrfxuGxfWvcI2/B/Iifq9t2cvo1sXwPLGETE9yPi/mFenT6wcdRHEQDH8uJ/GCuBOZm5N/Ap4IqIeMWmV795ah7zhcDrgD5a47yggz616sY2johzgLXA5VVTo9u4vbRh2jb8u480T2PbbIzGMubWlxEvB64FTsvMNVXzSL/tpo11vAdk5j7AYcAnIuId41lcDcZj+24DHAl8p+37ibx96+hbi+LuwzKRZOZ7RvouIp4a2k1c7U58epjZNvoogoiYChwN7Nu2zheAF6rpJRHxKLAbsJguqHPMmflU27IuAv59tD5168I2Ph44Ajg4qwPCTW/jNp383UeaZ5uN9O3k79aUsYyZiNiaVli5PDOvG5phI7/tpo1pvJk59P50RFxP6zDCnZS7jcc03sphwD3t23SCb9/N6dvI9nUPS31uAo6vpo8HbhxmntEeRfAe4KHMHBxqiIieiJhSTe8KzAN+VkP9m2NMY97gOPF84P625R4TEdtGRC+tMd9dQ/2baqzjPRQ4AzgyM3831KGgbdzJozJuAv5XtPw58Gy1i3hjfTv5uzVls8dcnW92CfBgZn6lvcNGfttNG8t4t4+I6QARsT2tk8jb/82WuI3H8psesuFe74m+fTenbzPbtxtn9m6JL+BVwO3AI9X7jKr9z4Cb2+Y7nNbVBI8C52ywjG8CJ2/Q9gHgAVpnbN8DvK/psY7XmIFvA/cBy2j9g5jV9t051fwPA4c1PdZxGu9yWseIl1avr5W2jYerHTh56HdJa7fxP1ff3wf0dzDuYf9upbw2d8zA22ntMl/Wtk0PH+233fRrDOPdtfqN/qT6vU6IbTzG3/R2wGrglRsscyJv39fQ2puyBvhNNf2Kkfo2uX29060kSSqeh4QkSVLxDCySJKl4BhZJklQ8A4skSSqegUWSJBXPwCJJkopnYJEkScUzsEiSpOL9f3SkcMr+7i3ZAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 648x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (9, 5))\n",
    "δ.plot(kind='hist', ax =ax , color = 'k', alpha = 0.5, label = 'δ')\n",
    "δ_shift = δ - δ_inint\n",
    "δ_shift.plot(kind='hist', ax =ax, color = 'b', alpha = 0.5, label = r'$δ_{shift}$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value:  0.00%\n"
     ]
    }
   ],
   "source": [
    "print(\"p-value: \" ,f\"{np.mean(δ_shift>δ_inint):.2%}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $\\text{p-value}=0\\%$  we can reject the null hypothesis and conclude new is better than old"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resources:\n",
    "\n",
    "* The Elements of Statistical Learning Data Mining, Inference, and Prediction - Second Edition - Trevor Hastie - Robert Tibshirani - Jerome Friedman. Springer.\n",
    "\n",
    "* Data Page: https://hastie.su.domains/ElemStatLearn/\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multinomial Logitic Regression"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use *multinomial logistic regression*, also called *softmax regression* (in older NLP literature you will sometimes see the name *maxent classifier*)\n",
    "\n",
    "The output of $\\textbf{y}$ for each $X$ is a vector of lenght $K$, it is **one-hot vector** version of $y$, such as, $y$ can be classified in range from $1$ to $K$. It can take one and only one of these classes.\n",
    "\n",
    "If $y$ belong to the $k$ class the **one-hot vector** will be $[0, 0, ..., 1, ..., 0]$ where $y_k=1$ and set all remaining elements of $\\textbf{y}$ of zeros. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The job of the classifier is produce an estimate vector $\\hat{\\textbf{{y}}}$, it is a vector of probabilities $[\\hat{y}_{1},..., \\hat{y}_{K} ]$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This job in actually is done by *Softmax* function. This function map a vector $z=[z_1, z_2, ..., z_K]$ of values into distribution of probabilities. Each in a range $(0, 1)$ and all sum to $1$.\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a vector $\\textbf{z}$ of dimensionality $K$, the softmax is defined as\n",
    "\n",
    "$$ \\text{softmax}(\\textbf{z}) = [\\frac{\\exp(z_1)}{\\sum_{k=1}^{K}\\exp(z_k)}, \\frac{\\exp(z_2)}{\\sum_{k=1}^{K}\\exp(z_k)}, ..., \\frac{\\exp(z_K)}{\\sum_{k=1}^{K}\\exp(z_k)}]$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.array([0.6,1.1,-1.5,1.2,3.2,-1.1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def softmax(z:np.array):\n",
    "    return np.exp(z)/np.sum(np.exp(z))    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.05482541, 0.09039182, 0.00671372, 0.09989841, 0.73815494,\n",
       "       0.0100157 ])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "softmax(z)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Softmax can be applied in *Logistic Regression*. We only need to assign for each observation and the for $k^{th}$ category:\n",
    "\n",
    "$$z_{k} = β_{k0} + β_{k}^T\\textbf{x} $$\n",
    "\n",
    "\n",
    "$$ \\hat{y_k} =  P(y = k|\\textbf{x}) = \\frac{\\exp(β_{k0} + β_{k}^T\\textbf{x})}{\\sum_{k=1}^{K} \\exp{(β_{k0} + β_{k}^T\\textbf{x})}}  $$\n",
    "\n",
    "Where $ \\hat{y_1} + \\hat{y_2} + ... + \\hat{y_K} = 1$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Features in Multinomial Logistic"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there is a $\\beta_0$ for each class, we need to evaluate the effect of this on eah class\n",
    "\n",
    "We suppose $\\beta_{2} = 4$ for the $k=1$, $\\beta_{2} = -1$ for the $k=2$ and $\\beta_{2} = 5.3$ for the $k=3$. Here $2$ means the second features, so we can say the second feature has a different effects on each class. e.g. the second feature increase the probability to belong to class $1$ (and to class $3$, but with more power), while the decrease the probability to belong to class $2$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross Entropy \n",
    "\n",
    "Loss Function of Sofmax "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In *Maximum Likehood* the likehood (or log likehood) is maximized.\n",
    "\n",
    "The requirement for a function is consired `loss function` is that this is minimized, so whe can multiply by `-`. So we get `negative log likehood` (called cross entropy) and this is the loss function for classification multiple categories."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log likelihood can be derived by defining, for each individual, $d_{k} = 1$ if alternative $k^{th}$ is chosen by individual, and $0$ if not, for the $K - 1$ possible outcomes. Then, for each individual one and only one of the $d_{k}’s$ is $1$.\n",
    "\n",
    "Recall the *conditional probabilities* are *probabilities* and the maximun likelihood maximize the conditional probabilities $p_{k} = P(y=k|\\textbf{x};\\theta)$ given the $\\theta$ parameters. The probability of the individual choices to class $k$\n",
    "\n",
    "$$L(\\theta) = -\\sum_{k=1}^{K} d_{k} log \\;p_{k} $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}