{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "tp4-000",
   "metadata": {},
   "source": [
    "# Practical 4 — Logistic regression\n",
    "\n",
    "In this practical, we analyse and implement a logistic regression for classification (yes the name is not perfect, but you should understand after the lesson)\n",
    "\n",
    "**Objectives**\n",
    "- connect a linear score, a sigmoid probability and a class decision;\n",
    "- compute binary logistic loss for a few observations;\n",
    "- fit a regularized logistic model in a scaling pipeline;\n",
    "- choose `C` and a decision threshold using development data only;\n",
    "- interpret a ROC curve and report one untouched test evaluation.\n",
    "\n",
    "Some exercises are inspired by *Machine Learning with PyTorch and Scikit-Learn*, Chapter 3, \"Modeling class probabilities via logistic regression\" of Raschka S."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-001",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.datasets import load_breast_cancer\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import (\n",
    "    accuracy_score,\n",
    "    confusion_matrix,\n",
    "    f1_score,\n",
    "    precision_score,\n",
    "    recall_score,\n",
    "    roc_auc_score,\n",
    "    roc_curve,\n",
    ")\n",
    "from sklearn.model_selection import StratifiedKFold, cross_val_score, train_test_split\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "RANDOM_STATE = 42\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-002",
   "metadata": {},
   "source": [
    "## 1. Data and protocol\n",
    "\n",
    "We use the [Wisconsin breast cancer dataset](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_breast_cancer.html). Each row contains 30 measurements computed from an image of a cell nucleus. In the original encoding, `0` means malignant and `1` means benign. **Here we recode malignant as class 1**, so recall means the fraction of malignant cases detected.\n",
    "\n",
    "This is a teaching dataset and a predictive exercise, not a clinical decision system. Split the data before fitting any preprocessing. Keep the test set untouched until Section 6. \n",
    "\n",
    "Documentation: [`train_test_split`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = load_breast_cancer()\n",
    "X = dataset.data\n",
    "y = (dataset.target == 0).astype(int)  # 1 = malignant\n",
    "feature_names = np.asarray(dataset.feature_names)\n",
    "\n",
    "X_dev, X_test, y_dev, y_test = train_test_split(\n",
    "    X, y, test_size=0.20, stratify=y, random_state=RANDOM_STATE\n",
    ")\n",
    "X_train, X_valid, y_train, y_valid = train_test_split(\n",
    "    X_dev, y_dev, test_size=0.25, stratify=y_dev, random_state=RANDOM_STATE\n",
    ")\n",
    "\n",
    "for name, target in ((\"train\", y_train), (\"valid\", y_valid), (\"test\", y_test)):\n",
    "    print(f\"{name:>5}: n={len(target):3d}, malignant prevalence={target.mean():.1%}\")\n",
    "print(\"X shape:\", X.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-004",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: inspect the range of at least three features on X_train.\n",
    "# Are their scales similar? What does y_train.mean() represent?\n",
    "# Why use stratify in both splits?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-006",
   "metadata": {},
   "source": [
    "## 2. Score, sigmoid and loss\n",
    "\n",
    "For a linear score $z=\\mathbf w^\\top\\mathbf x+b$, logistic regression assumes\n",
    "\n",
    "$$\n",
    "\\operatorname{logit}(p)=\\log\\frac{p}{1-p}=z,\n",
    "\\qquad\n",
    "p=\\sigma(z)=\\frac{1}{1+e^{-z}}.\n",
    "$$\n",
    "\n",
    "The logit maps probabilities to real scores; the sigmoid maps scores back to probabilities. For the three scores below, compute $p$, then verify that $\\operatorname{logit}(p)$ recovers $z$. What happens at $z=0$?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-007",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "scores = np.array([-2.0, 0.0, 2.0])\n",
    "\n",
    "# TODO: write sigmoid(z) using np.exp.\n",
    "# TODO: compute probabilities = sigmoid(scores).\n",
    "# TODO: recover scores with np.log(probabilities / (1 - probabilities)).\n",
    "# Compare both score arrays with np.allclose.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-009",
   "metadata": {},
   "source": [
    "For a binary label $y\\in\\{0,1\\}$, the loss for a predicted probability $p=P(Y=1\\mid\\mathbf x)$ is\n",
    "\n",
    "$$\n",
    "\\ell(y,p)=-y\\log p-(1-y)\\log(1-p).\n",
    "$$\n",
    "\n",
    "Evaluate it for `(y, p) = (1, 0.9), (1, 0.1), (0, 0.9), (0, 0.1)`. Which predictions receive the largest penalty? This is the negative log-likelihood of the Bernoulli model discussed in CM03.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-010",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "labels_small = np.array([1, 1, 0, 0])\n",
    "probs_small = np.array([0.9, 0.1, 0.9, 0.1])\n",
    "\n",
    "# TODO: implement binary_log_loss(y_true, p) and compute all four losses.\n",
    "# For these values, no clipping is required. Keep the function vectorized.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-012",
   "metadata": {},
   "source": [
    "## 3. A complete logistic pipeline\n",
    "The 30 features have different scales. Keep [`StandardScaler`](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html) inside a [`Pipeline`](https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html), followed by [`LogisticRegression`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html). In cross-validation, each fold must fit its scaler using only that fold's training observations.\n",
    "\n",
    "The classifier's `decision_function` returns the linear score $z$; `predict_proba` returns the corresponding probability. Check this relationship numerically. We use the default $L_2$ regularization and `C=1.0` for the initial model. Identify the role of $C$ according to $\\lambda$ in the lesson.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-013",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: build a Pipeline with steps named \"scale\" and \"classifier\".\n",
    "# Use StandardScaler() and LogisticRegression(C=1.0, max_iter=2000).\n",
    "# Fit on X_train, y_train only.\n",
    "# Compute z_valid and p_valid for X_valid.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-015",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: compare model.predict(X_valid) with (p_valid >= 0.5).astype(int).\n",
    "# Print accuracy, precision, recall, F1 and ROC-AUC on validation data.\n",
    "# Which metric uses probabilities instead of thresholded labels?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-017",
   "metadata": {},
   "source": [
    "## 4. Choose the regularization strength\n",
    "\n",
    "`C` is the inverse of regularization strength: a smaller value shrinks coefficients more. Compare a short logarithmic list of `C` values with five-fold [`StratifiedKFold`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.StratifiedKFold.html) on **training data only**. Use [`cross_val_score`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_score.html) with `scoring=\"roc_auc\"`. We will use the validation set to derive a possibly better threshold.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-018",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "C_values = np.logspace(-3, 2, 8)\n",
    "cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)\n",
    "\n",
    "# TODO: for each C, create a fresh pipeline, compute five CV ROC-AUC scores,\n",
    "# and store their mean and standard deviation. Select the best C.\n",
    "# Keep the validation and test sets out of this loop.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-020",
   "metadata": {},
   "source": [
    "Fit a new pipeline on the whole training subset with `best_C`. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-021",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: fit selected_model using best_C on X_train, y_train."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-023",
   "metadata": {},
   "source": [
    "## 5. From probabilities to a ROC operating point\n",
    "Compute the selected model's validation probabilities. Plot the [`roc_curve`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_curve.html) and report [`roc_auc_score`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_auc_score.html). The diagonal corresponds to random ranking, with an AUC of $0.5$. Because the sigmoid is increasing, the ROC is the same whether we pass `decision_function` scores or class 1 probabilities.\n",
    "\n",
    "Each threshold gives one point $(\\mathrm{FPR},\\mathrm{TPR})$ on this curve. Locate the operating point for $\\tau=0.5$ and compare it with $\\tau=0.8$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-024",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: compute p_valid and z_valid from selected_model.\n",
    "# Use roc_curve(y_valid, p_valid) to obtain fpr and tpr.\n",
    "# Plot the ROC and the random-ranking diagonal.\n",
    "# Verify ROC-AUC is the same for z_valid and p_valid.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-026",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: for thresholds 0.5 and 0.8, compute predictions and report precision, recall, FPR and F1 for each threshold.\n",
    "# Explain the direction of movement on the ROC curve.\n",
    "# Verify that tau=0.8 corresponds to z >= log(0.8 / 0.2).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-028",
   "metadata": {},
   "source": [
    "To make threshold selection explicit, suppose a false negative costs 5 and a false positive costs 1. On the validation subset only, evaluate candidate thresholds in `np.linspace(0.1, 0.9, 17)` and choose the one with the smallest total cost $5\\,\\mathrm{FN}+\\mathrm{FP}$. This illustrative cost ratio is **not** a clinical recommendation. Compare the chosen point with $\\tau=0.5$ on the ROC plot.\n",
    "\n",
    "If several thresholds have the same minimum cost, the code below chooses the first one. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-029",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "candidate_thresholds = np.linspace(0.1, 0.9, 17)\n",
    "\n",
    "# TODO: compute validation cost 5 * FN + FP for each candidate threshold.\n",
    "# Store the costs in an array and select the first minimum with np.argmin.\n",
    "# Do not inspect the test set.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-031",
   "metadata": {},
   "source": [
    "## 6. One final test evaluation\n",
    "\n",
    "Freeze the fitted model and the selected threshold. Evaluate them **once** on `X_test, y_test`. Report accuracy, precision, recall, F1, ROC-AUC and the confusion matrix. Compare the chosen threshold with the default `0.5` on the same test probabilities. The test results are an estimate of performance for this fixed protocol; do not tune either `C` or the threshold after seeing them. If the chosen threshold performs worse on test than on validation, report that discrepancy.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-032",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: compute p_test from selected_model without refitting it.\n",
    "# Compare tau=0.5 and chosen_threshold on the test set.\n",
    "# Report the requested metrics and confusion matrices.\n",
    "# State whether the cost-based choice changed the FP/FN trade-off.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-034",
   "metadata": {},
   "source": [
    "### Takeaways\n",
    "\n",
    "1. What remains fixed when the decision threshold changes?\n",
    "2. Why can ROC-AUC remain constant while recall changes?\n",
    "3. Which choices used only training data, which used validation data, and what was reserved for test?\n",
    "4. Why do the reported metrics not establish clinical usefulness?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp4-036",
   "metadata": {
    "tags": [
     "optional"
    ]
   },
   "source": [
    "## Optional — Gradient of the logistic loss\n",
    "\n",
    "Raschka's derivation and the bonus CM03 slides show that, for one observation,\n",
    "\n",
    "$$\n",
    "\\frac{\\partial\\ell_i}{\\partial z_i}=p_i-y_i.\n",
    "$$\n",
    "\n",
    "Therefore, without regularization,\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf w}J=\\frac{1}{n}X^\\top(\\mathbf p-\\mathbf y),\n",
    "\\qquad\n",
    "\\frac{\\partial J}{\\partial b}=\\frac{1}{n}\\sum_i(p_i-y_i).\n",
    "$$\n",
    "\n",
    "On a tiny synthetic design matrix, compute these gradients at $\\mathbf w=0,b=0$ and make one gradient-descent step. This illustrates the **shape** of the update, not scikit-learn's numerical solver.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp4-037",
   "metadata": {
    "tags": [
     "exercise",
     "optional"
    ]
   },
   "outputs": [],
   "source": [
    "X_tiny = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [-1.0, -1.0]])\n",
    "y_tiny = np.array([1, 1, 1, 0])\n",
    "w_tiny = np.zeros(2)\n",
    "b_tiny = 0.0\n",
    "learning_rate = 0.5\n",
    "\n",
    "# TODO: compute z, p, grad_w and grad_b.\n",
    "# TODO: update w_tiny and b_tiny once.\n",
    "# Compare mean binary_log_loss before and after the update.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cedd0ae4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "insa-iml-practicals (3.10.20)",
   "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.10.20"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}