{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "tp3-00",
   "metadata": {},
   "source": [
    "# Practical 3 — Linear regression and regularization\n",
    "\n",
    "In this practical, you will connect the equations of ordinary least squares to their numerical implementation, then compare OLS, Ridge and Lasso on a scientific regression problem.\n",
    "\n",
    "By the end of the session, you should be able to:\n",
    "\n",
    "- solve a small OLS problem from the normal equations;\n",
    "- build OLS, Ridge and Lasso estimators with scikit-learn;\n",
    "- explain why scaling matters for regularized linear models;\n",
    "- select a regularization strength without using the test set;\n",
    "- read Ridge and Lasso coefficient paths;\n",
    "- compare regression models against a simple baseline;\n",
    "- distinguish predictive coefficients from causal effects.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-01",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.base import clone\n",
    "from sklearn.datasets import load_diabetes\n",
    "from sklearn.dummy import DummyRegressor\n",
    "from sklearn.linear_model import Lasso, LinearRegression, Ridge\n",
    "from sklearn.metrics import (\n",
    "    mean_absolute_error,\n",
    "    root_mean_squared_error,\n",
    "    r2_score,\n",
    ")\n",
    "from sklearn.model_selection import (\n",
    "    GridSearchCV,\n",
    "    KFold,\n",
    "    cross_val_score,\n",
    "    train_test_split,\n",
    ")\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "RANDOM_STATE = 42\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-02",
   "metadata": {},
   "source": [
    "## 1. OLS from the normal equations\n",
    "\n",
    "A laboratory instrument measures a signal $x_i$ and must predict a reference concentration $y_i$. We use the affine model\n",
    "\n",
    "$$\n",
    "\\widehat y_i = w x_i + b.\n",
    "$$\n",
    "\n",
    "To estimate both parameters with the matrix formulation, augment the design matrix with a column of ones:\n",
    "\n",
    "$$\n",
    "X_{\\mathrm{aug}}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "x_1 & 1\\\\\n",
    "\\vdots & \\vdots\\\\\n",
    "x_n & 1\n",
    "\\end{bmatrix},\n",
    "\\qquad\n",
    "\\boldsymbol\\theta=\n",
    "\\begin{bmatrix}w\\\\b\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "What are the dimensions of $X_{\\mathrm{aug}}$, $X_{\\mathrm{aug}}^\\top X_{\\mathrm{aug}}$ and $\\boldsymbol\\theta$ for the six observations below?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-03",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_calibration = np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0])\n",
    "y_calibration = np.array([1.5, 2.0, 2.9, 3.3, 4.2, 4.6])\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.scatter(x_calibration, y_calibration, color=\"tab:blue\")\n",
    "plt.xlabel(\"Measured signal\")\n",
    "plt.ylabel(\"Reference concentration\")\n",
    "plt.title(\"Calibration observations\")\n",
    "plt.grid(alpha=0.25)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-04",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: build X_aug with the measurement column followed by a column of ones.\n",
    "# Expected shape: (6, 2)\n",
    "\n",
    "# TODO: compute A = X_aug.T @ X_aug and c = X_aug.T @ y_calibration.\n",
    "# What are the shapes of A and c?\n",
    "\n",
    "# TODO: solve A @ theta = c with np.linalg.solve.\n",
    "# Extract w_manual and b_manual from theta.\n",
    "\n",
    "# TODO: compute y_manual and residuals_manual.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-06",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# Read the LinearRegression documentation, then fit the same affine model.\n",
    "# TODO: reshape x_calibration as a two-dimensional design matrix.\n",
    "# TODO: instantiate and fit LinearRegression.\n",
    "# TODO: compare coef_, intercept_ and predictions with the manual solution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-08",
   "metadata": {},
   "source": [
    "## 2. Scientific dataset and independent test set\n",
    "\n",
    "We now use scikit-learn's [diabetes dataset](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_diabetes.html). Each row describes one patient through ten baseline clinical measurements. The continuous target is a quantitative measure of disease progression one year after baseline.\n",
    "\n",
    "We load the data with `scaled=False` to retain the original feature units. This differs from the default `scaled=True`, which returns centered and scaled features (see the documentation).\n",
    "\n",
    "*Disclaimer* : This is a **predictive exercise**. The fitted coefficients do not establish causal effects and must not be interpreted as medical recommendations.\n",
    "\n",
    "1. Is this a regression or classification problem?\n",
    "2. What is one observation?\n",
    "3. Which data must remain untouched during model selection?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-09",
   "metadata": {},
   "outputs": [],
   "source": [
    "diabetes = load_diabetes(scaled=False)\n",
    "X = diabetes.data\n",
    "y = diabetes.target\n",
    "feature_names = np.asarray(diabetes.feature_names)\n",
    "\n",
    "X_dev, X_test, y_dev, y_test = train_test_split(\n",
    "    X,\n",
    "    y,\n",
    "    test_size=0.25,\n",
    "    random_state=RANDOM_STATE,\n",
    ")\n",
    "\n",
    "print(\"Development shape:\", X_dev.shape, y_dev.shape)\n",
    "print(\"Test shape:\", X_test.shape, y_test.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-10",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# Use development data only in this section.\n",
    "# TODO: print the mean, standard deviation, minimum and maximum\n",
    "# of each feature with its name.\n",
    "\n",
    "# TODO: draw a boxplot of the development features.\n",
    "# TODO: draw a histogram of y_dev.\n",
    "\n",
    "# Why would coefficient penalties be unfair across variables\n",
    "# expressed on very different scales?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-12",
   "metadata": {},
   "source": [
    "## 3. OLS reference model and residuals\n",
    "\n",
    "Feature scaling does not change OLS predictions when the model includes an intercept, although it changes the numerical coefficients. We nevertheless use the same preprocessing structure for all three models so that their coefficients are expressed on comparable standardized features.\n",
    "\n",
    "Documentation: [`Pipeline`](https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html), [`StandardScaler`](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html), and [`LinearRegression`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-13",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: create a Pipeline containing StandardScaler and LinearRegression.\n",
    "# Fit it on X_dev, y_dev only.\n",
    "# Compute development predictions and report MAE, RMSE and R².\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-15",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: compute the development residuals y_dev - y_dev_ols.\n",
    "# Plot them against the fitted values and add a horizontal line at zero.\n",
    "# Look for systematic structure, changes in spread and isolated residuals.\n",
    "# Why is this not a final estimate of generalization performance?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-17",
   "metadata": {},
   "source": [
    "## 4. Ridge and Lasso pipelines\n",
    "\n",
    "- **Ridge** shrinks coefficients and is often stable when predictors are correlated.\n",
    "- **Lasso** also shrinks coefficients and may set some of them exactly to zero.\n",
    "- Both penalties act on coefficient magnitudes, so the features must be scaled consistently.\n",
    "\n",
    "Documentation: [`Ridge`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html) and [`Lasso`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-18",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: build ridge_pipeline and lasso_pipeline.\n",
    "# Each pipeline must contain a fresh StandardScaler followed by a regressor.\n",
    "# Start with alpha=1.0 and use a sufficiently large max_iter for Lasso. \n",
    "# What does alpha represent?\n",
    "\n",
    "# Fit one temporary pipeline on development data.\n",
    "# Check that its scaler.mean_ equals X_dev.mean(axis=0) (not statistics estimated from the complete dataset). \n",
    "\n",
    "# Why must the scaler remain inside the object passed to cross-validation?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-20",
   "metadata": {},
   "source": [
    "## 5. Select the regularization strengths\n",
    "\n",
    "Use a five-fold shuffled `KFold` on the development set. RMSE is the primary selection metric and the final test set remains untouched.\n",
    "\n",
    "We explore logarithmic grids because useful regularization strengths may span several orders of magnitude:\n",
    "\n",
    "```python\n",
    "ridge_alphas = np.logspace(-4, 4, 41)\n",
    "lasso_alphas = np.logspace(-4, 1, 41)\n",
    "```\n",
    "\n",
    "scikit-learn maximizes scores. The search therefore uses `scoring=\"neg_root_mean_squared_error\"` (RTFM). Negate the reported score to recover RMSE.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-21",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "ridge_alphas = np.logspace(-4, 4, 41)\n",
    "lasso_alphas = np.logspace(-4, 1, 41)\n",
    "\n",
    "# TODO: instantiate KFold with five folds, shuffling and RANDOM_STATE.\n",
    "# TODO: create one GridSearchCV for each pipeline.\n",
    "# Address alpha with the parameter name regressor__alpha.\n",
    "# TODO: fit both searches on development data only.\n",
    "# Print each selected alpha and its positive cross-validated RMSE.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-23",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: evaluate a mean DummyRegressor and the OLS pipeline with\n",
    "# cross_val_score, using the same cv object and scoring metric.\n",
    "# Compare their mean RMSE with the best Ridge and Lasso validation RMSE.\n",
    "# Why must all models use the same folds?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-25",
   "metadata": {},
   "source": [
    "## 6. Coefficient paths and interpretation\n",
    "\n",
    "Track the standardized coefficients as the regularization strength changes. Every model below is fitted on development data. The paths help inspect model behaviour; they are not a test evaluation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-26",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "ridge_path_alphas = np.logspace(-4, 4, 50)\n",
    "lasso_path_alphas = np.logspace(-4, 2, 50)\n",
    "\n",
    "# TODO: for every Ridge alpha, fit a cloned pipeline on development data\n",
    "# and store regressor.coef_. Repeat for Lasso.\n",
    "\n",
    "# TODO: produce two plots with logarithmic alpha axes.\n",
    "# Use feature_names in the legends and label both axes.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-28",
   "metadata": {},
   "source": [
    "### Interpret the paths\n",
    "\n",
    "1. What happens to Ridge coefficients as alpha increases?\n",
    "2. Which Lasso coefficients become exactly zero?\n",
    "3. Are all coefficient trajectories necessarily monotonic?\n",
    "4. Why does a zero Lasso coefficient not prove that the corresponding variable has no scientific relevance?\n",
    "5. Why are coefficients easier to compare after standardization?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-30",
   "metadata": {},
   "source": [
    "## 7. Freeze the models and evaluate once\n",
    "\n",
    "Before revealing the test results, write down:\n",
    "\n",
    "- the OLS, selected Ridge and selected Lasso learning procedures;\n",
    "- the selected alpha values;\n",
    "- RMSE as the primary metric and the mean-prediction baseline;\n",
    "- one limitation of this experiment.\n",
    "\n",
    "Then evaluate exactly once on `X_test`, `y_test`:\n",
    "\n",
    "- the mean-prediction baseline;\n",
    "- OLS fitted on all development data;\n",
    "- `ridge_search.best_estimator_`;\n",
    "- `lasso_search.best_estimator_`.\n",
    "\n",
    "Report MAE, RMSE and $R^2$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-31",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "def regression_metrics(y_true, y_pred):\n",
    "    # TODO: return MAE, RMSE and R² in a dictionary.\n",
    "    return {}\n",
    "\n",
    "\n",
    "# TODO: fit the baseline and OLS on all development data.\n",
    "# Ridge and Lasso were refitted on all development data by GridSearchCV.\n",
    "# Evaluate the four frozen models on the test set exactly once.\n",
    "# Print a compact aligned comparison without pandas.\n",
    "\n",
    "# TODO: report the number of nonzero coefficients in the selected Lasso.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-33",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "source": [
    "### Final interpretation\n",
    "\n",
    "Complete these statements in your own words.\n",
    "\n",
    "1. Ridge differs from OLS because ...\n",
    "2. Lasso differs from Ridge because ...\n",
    "3. The selected alpha was chosen without test leakage because ...\n",
    "4. The coefficient values do not prove causal effects because ...\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp3-35",
   "metadata": {
    "tags": [
     "optional"
    ]
   },
   "source": [
    "---\n",
    "\n",
    "The remainder of this practical is optional and can be completed as bonus work.\n",
    "\n",
    "---\n",
    "\n",
    "## 8. Optional — Coefficient stability under resampling\n",
    "\n",
    "Repeatedly resample the development observations, refit OLS, the selected Ridge model and the selected Lasso model, then compare the variation of their standardized coefficients.\n",
    "\n",
    "Before running the experiment, predict which model should have the most stable coefficients.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp3-36",
   "metadata": {
    "tags": [
     "exercise",
     "optional"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO (optional): use a fixed NumPy generator and 100 bootstrap samples.\n",
    "# For each sample, clone and fit ols_pipeline, ridge_search.best_estimator_\n",
    "# and lasso_search.best_estimator_.\n",
    "# Store the standardized coefficient vectors.\n",
    "# Compare their feature-wise standard deviations.\n",
    "\n",
    "# Do not evaluate the final test set inside this loop.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (cours-ml)",
   "language": "python",
   "name": "cours-ml"
  },
  "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
}