tutorial_notebooks/part1_regression_tuto.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 📈 Part I: Conformal Regression\n",
    "\n",
    "Let’s consider a simple regression problem on our heteroskedastic data. We want to evaluate the uncertainty associated with the prediction using various conformal prediction methods.\n",
    "\n",
    "<div align=center>\n",
    "<img src=\"assets/regression_demo.svg\" width=\"700\"> </img>\n",
    "</div>\n",
    "\n",
    "**Links**\n",
    "- [<img src=\"https://github.githubassets.com/images/icons/emoji/octocat.png\" width=20> Github](https://github.com/deel-ai/puncc)\n",
    "- [📘 Documentation](https://deel-ai.github.io/puncc/index.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⚙️ Setup  <a class=\"anchor\" id=\"cr-setup\"></a>\n",
    "\n",
    "🐾 First, ensure we have puncc installed. You can install it if needed using pip:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install puncc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import some modules and implement the function `evaluate_cp` that evaluates a conformal predictor on a new dataset by predicting intervals for the test inputs and computing the average width and coverage of these predictions, using a specified risk level $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "from deel.puncc.metrics import regression_sharpness, regression_mean_coverage\n",
    "import logging\n",
    "\n",
    "logging.getLogger(\"matplotlib.font_manager\").setLevel(logging.ERROR)\n",
    "\n",
    "sns.set_theme()\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# Set random seed for reproducibility\n",
    "np.random.seed(0)\n",
    "\n",
    "\n",
    "# Function to evaluate the sharpness and coverage of a conformal prediction model\n",
    "def evaluate_cp(X_test, y_test, model_cp, alpha):\n",
    "    \"\"\"\n",
    "    Evaluate the performance of a model using conformal prediction.\n",
    "\n",
    "    Parameters:\n",
    "    - X_test : The input features for the test set.\n",
    "    - y_test : The true labels for the test set.\n",
    "    - model_cp : The conformal prediction model.\n",
    "    - alpha : The maximum risk level for the prediction intervals.\n",
    "\n",
    "    Returns:\n",
    "    - sharpness : The average width of the prediction intervals.\n",
    "    - coverage : The average coverage of the prediction intervals.\n",
    "    \"\"\"\n",
    "    y_pred, y_pred_lower, y_pred_upper = model_cp.predict(X_test, alpha=alpha)\n",
    "    sharpness = regression_sharpness(y_pred_lower, y_pred_upper)\n",
    "    coverage = regression_mean_coverage(y_test, y_pred_lower, y_pred_upper)\n",
    "    return sharpness, coverage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 💾 Dataset <a class=\"anchor\" id=\"cr-data\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider a synthetic 1D heteroskedastic dataset, where the variance of the noise increases with the value of the input feature.\n",
    "We generate $N$ samples as follows:\n",
    "\n",
    "- Inputs $X$ are uniformly distributed on $[0, 20]$ \n",
    "- Outputs are given by $Y = (1+\\epsilon)\\cdot X, $\n",
    "\n",
    "Such that $\\epsilon \\sim {\\cal N}(\\mu=0,\\sigma=1)$ is standard gaussian noise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Generate synthetic 1D heteroskedastic data\n",
    "n_samples = 2000\n",
    "X = np.linspace(0, 20, n_samples).reshape(-1, 1)\n",
    "y = X.squeeze() * (1 + np.random.randn(n_samples))\n",
    "\n",
    "# cast data to dtype float32\n",
    "X = X.astype(np.float32)\n",
    "y = y.astype(np.float32)\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.25, random_state=0\n",
    ")\n",
    "\n",
    "# Plot the generated data to visualize heteroskedasticity\n",
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "plt.scatter(X, y, alpha=0.6)\n",
    "plt.xlabel(r\"$X$\")\n",
    "plt.ylabel(r\"$Y$\")\n",
    "plt.xticks(np.arange(22, step=2))\n",
    "\n",
    "\n",
    "# change font to something more readable and arial\n",
    "plt.rcParams.update({\"font.family\": \"Times New Roman\"})\n",
    "plt.rcParams.update({\"font.size\": 14})\n",
    "plt.tight_layout()\n",
    "plt.grid(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By construction, data are independent and identically distributed (i.i.d). We fullfill the prerequisites to apply conformal prediction 👏 !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🔮 Prediction Model <a class=\"anchor\" id=\"cr-pm\"></a>\n",
    "\n",
    "We will define a linear model from scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "linear_model = LinearRegression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There exist different conformal prediction methods that are tailored for specific cases and needs. In each method, the underlying model(s) to be conformalized are wrapped by puncc to guarantee compliance with its framework. Such wrappers are provided by puncc and are summarized in the table below. \n",
    "<br/><br/>\n",
    "<div align='center'>\n",
    "\n",
    "| Conformal Prediction Method | Model Wrapper         |\n",
    "|-----------------------------|-----------------------|\n",
    "| SplitCP                     | [`BasePredictor`](https://deel-ai.github.io/puncc/prediction.html#prediction.BasePredictor) |\n",
    "| CVplus                      | [`BasePredictor`](https://deel-ai.github.io/puncc/prediction.html#prediction.BasePredictor)        |\n",
    "| LocallyAdaptiveCP           | [`MeanVarPredictor`](https://deel-ai.github.io/puncc/prediction.html#prediction.MeanVarPredictor)      |\n",
    "| EnbPI                       | [`BasePredictor`](https://deel-ai.github.io/puncc/prediction.html#prediction.BasePredictor)        |\n",
    "| CQR                         | [`DualPredictor`](https://deel-ai.github.io/puncc/prediction.html#prediction.DualPredictor)        |\n",
    "| APS or RAPS  (classification)               | [`BasePredictor`](https://deel-ai.github.io/puncc/prediction.html#prediction.BasePredictor) |\n",
    "| SplitBoxWise (object detection)               | [`BasePredictor`](https://deel-ai.github.io/puncc/prediction.html#prediction.BasePredictor) |\n",
    "<caption>\n",
    "Table 1: Model Wrappers for Different Conformal Prediction Methods in PUNCC\n",
    "</caption>\n",
    "\n",
    "</div>\n",
    "\n",
    "<br/>\n",
    "\n",
    "*📓 Note: puncc also enables to create fully custom [Predictor](https://deel-ai.github.io/puncc/api.html#predictor)(s). To know more about the puncc's API, check [this tutorial](https://github.com/deel-ai/puncc/blob/main/docs/api_intro.ipynb) <sub> [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/drive/1d06qQweM1X1eSrCnixA_MLEZil1vXewj) </sub>.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following sections, we will explore various methods for quantifying uncertainty using conformal prediction, specifically focusing on Split Conformal Prediction (SplitCP), Conformal Cross-Validation Plus (CV+), and Conformalized Quantile Regression (CQR)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⏸️ Conformal Prediction with SplitCP <a class=\"anchor\" id=\"cr-splitcp\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first approach we will use to conformalize our model is Split Conformal Prediction (SplitCP). SplitCP generates constant-length prediction intervals that provide a specified coverage probability $1-\\alpha$, such that $\\alpha$ is the maximum risk level. It works by splitting the data into two sets: a proper training (fit) set to train the model and a calibration set to calibrate the prediction intervals ensuring the desired coverage level. \n",
    "\n",
    "<div align=center>\n",
    "<img src=\"assets/workflow_regression.svg\" width=\"800\"> </img>\n",
    "</div>\n",
    "\n",
    "In accordance with *Table 1*, let's first wrap our model using `BasePredictor`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from deel.puncc.api.prediction import BasePredictor\n",
    "\n",
    "base_predictor = BasePredictor(...)  # TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to split the available data into proper training (fit) and calibration sets (we will see later in this tutorial how to delegate to puncc the task of creating automatically these partitions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_fit, X_calib, y_fit, y_calib = train_test_split(X_train, y_train, test_size=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can to initialise and fit  conformal predictor. The split conformal prediction method is provided by the class `deel.puncc.regression.SplitCP`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from deel.puncc.regression import SplitCP\n",
    "\n",
    "# Wrap the base predictor in a split conformal predictor\n",
    "# train=True to fit the split conformal predictor\n",
    "splitcp = SplitCP(...)  # TODO\n",
    "\n",
    "# Fit the split conformal predictor:\n",
    "# puncc will fit the underlying model on the fit data and\n",
    "# compute nonconformity scores on the calibration set\n",
    "splitcp.fit(...)  # TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's evaluate the performance of the prediction intervals (coverage and width) for a risk level $\\alpha = 10\\%$. Additionally, we will plot the prediction intervals to visually assess their performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from deel.puncc.plotting import plot_prediction_intervals\n",
    "\n",
    "# Maximum allowed error rate\n",
    "alpha = 0.1\n",
    "\n",
    "# Compute prediction intervals and metrics on the test set\n",
    "y_pred, y_pred_lower, y_pred_upper = splitcp.predict(...)  # TODO\n",
    "\n",
    "# Evaluate the sharpness and coverage of the prediction intervals\n",
    "sharpness, coverage = evaluate_cp(X_test, y_test, splitcp, alpha)\n",
    "print(f\"Average prediction intervals width (sharpness): {sharpness:.3f}\")\n",
    "print(f\"Average coverage: {coverage:.3f}\")\n",
    "\n",
    "# Plot the prediction intervals\n",
    "plot_prediction_intervals(\n",
    "    y_true=y_test,\n",
    "    y_pred_lower=y_pred_lower,\n",
    "    y_pred_upper=y_pred_upper,\n",
    "    y_pred=y_pred,\n",
    "    X=X_test[:, 0],\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting the results manually"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's implement the SplitCP manually and compare our results with puncc.\n",
    "\n",
    "Given a prediction model $\\widehat{f}$ trained on $D_{train}$, the algorithm is summarized in the following:\n",
    "\n",
    "- Choose as nonconformity score the absolute deviation: $R_i = |\\widehat{f}(X_i)-Y_i|$.\n",
    "- Compute the nonconformity scores on the calibration dataset: $\\bar{R} = \\{R_i\\}_{}$, for $i=1,\\dots,n$, where $n$ is the size of $D_{calibration}$.\n",
    "- Compute the error margin $\\delta_{\\alpha}$ as the $(1-\\alpha)(1 + \\frac{1}{| D_{calibration} |})$-th empirical quantile of $\\bar{R}$, sorted in ascending order.\n",
    "- Build the prediction interval $\\widehat{C}_{\\alpha}(X_{new}) = \\Big[ \\widehat{f}(X_{new}) - \\delta_{\\alpha}^{f} \\,,\\, \\widehat{f}(X_{new}) + \\delta_{\\alpha}^{f} \\Big]$.\n",
    "\n",
    "*Note: we can either train the linear model or access the trained instance associated to our conformal predictor `splitcp.predictor`.*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predictions on calibration set\n",
    "y_pred_calib_manual = splitcp.predictor.predict(X_calib)\n",
    "\n",
    "# Compute nonconformity scores and sort them\n",
    "nonconf_scores_manual = ...  # TODO\n",
    "nonconf_scores_manual.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get nonconformity scores from the split conformal predictor\n",
    "nonconf_scores = splitcp.get_nonconformity_scores()\n",
    "# Sort the nonconformity scores\n",
    "nonconf_scores.sort()\n",
    "\n",
    "# Check that the nonconformity scores are the same\n",
    "assert np.all(nonconf_scores == nonconf_scores_manual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the error margin $\\delta_{\\alpha}$ and show it with a plot of the nonconformity scores  histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = len(nonconf_scores)\n",
    "q = np.quantile(nonconf_scores, (1 - alpha) * (n + 1) / n, method=\"inverted_cdf\")\n",
    "\n",
    "plt.hist(nonconf_scores, bins=100)\n",
    "# vertical line with text\n",
    "plt.axvline(x=q, color=\"r\", linestyle=\"--\")\n",
    "plt.text(q, 10, f\"q = {q:.2f}\", color=\"r\", rotation=270)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For new data, te error margin (quantile) we computed will serve as constant-width for the prediction interval centered on the point estimate $\\widehat{f}(X_{new})$. Let's apply this to our test set and compare the results with punnc's results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = splitcp.predictor.predict(X_test)\n",
    "\n",
    "y_pred_lower_manual = ...  # TODO\n",
    "y_pred_upper_manual = ...  # TODO\n",
    "\n",
    "assert np.all(y_pred_lower_manual == y_pred_lower)\n",
    "assert np.all(y_pred_upper_manual == y_pred_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🤞 Conformal Prediction with CV+ <a class=\"anchor\" id=\"cr-cvplus\"></a>\n",
    "\n",
    "Let's reiterate the procedure using Conformal Cross-Validation Plus (CV+). CV+ extends the concept of split conformal prediction by integrating it with cross-validation. This combination enhances both the reliability and efficiency of the predictive intervals, overcoming limitations associated with a fixed proper training and calibration split.\n",
    "\n",
    "We will again use the `BasePredictor` instance defined earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from deel.puncc.regression import CVPlus\n",
    "\n",
    "# Wrap the base predictor in a split conformal predictor, choose K=5\n",
    "cvplus = CVPlus(..., K=..., random_state=0)  # TODO\n",
    "\n",
    "# Fit the CV+ predictor\n",
    "cvplus.fit(X=X_train, y=y_train)\n",
    "\n",
    "# Compute prediction intervals and metrics on the test set\n",
    "y_pred, y_pred_lower, y_pred_upper = cvplus.predict(...)  # TODO\n",
    "\n",
    "# Evaluate the sharpness and coverage of the prediction intervals\n",
    "sharpness, coverage = evaluate_cp(X_test, y_test, cvplus, alpha)\n",
    "print(f\"Average prediction intervals width (sharpness): {sharpness:.3f}\")\n",
    "print(f\"Average coverage: {coverage:.3f}\")\n",
    "\n",
    "# Plot the prediction intervals\n",
    "plot_prediction_intervals(\n",
    "    y_true=y_test,\n",
    "    y_pred_lower=y_pred_lower,\n",
    "    y_pred_upper=y_pred_upper,\n",
    "    y_pred=y_pred,\n",
    "    X=X_test[:, 0],\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 〽️ Conformal Prediction with CQR <a class=\"anchor\" id=\"cr-cqr\"></a>\n",
    "\n",
    "Let's now consider Conformalized Quantile Regression (CQR). CQR extends traditional quantile regression by incorporating conformal prediction techniques, allowing us to construct predictive intervals with state-of-the-art performance and guaranteed coverage (under data exchangeability).\n",
    "\n",
    "To do so, we will use the gradient boosting quantile regressor implemented in sklearn. Alternatively, we can train a neural network model to perform quantile regression by setting its loss function as the pinball loss. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "\n",
    "# Lower quantile regressor\n",
    "lower_quantile_model = GradientBoostingRegressor(\n",
    "    loss=\"quantile\", alpha=alpha / 2, n_estimators=50\n",
    ")\n",
    "# Upper quantile regressor\n",
    "upper_quantile_model = GradientBoostingRegressor(\n",
    "    loss=\"quantile\", alpha=1 - alpha / 2, n_estimators=50\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to *Table 1*, we need wrap our models using `DualPredictor`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from deel.puncc.api.prediction import DualPredictor\n",
    "from deel.puncc.regression import CQR\n",
    "\n",
    "# Wrap the upper and lower quantile models in a dual predictor\n",
    "dualpredictor = DualPredictor([..., ...])  # TODO\n",
    "\n",
    "# Initialize the CQR conformal predictor\n",
    "# train=True to use the train dual predictor\n",
    "cqr = CQR(...)  # TODO\n",
    "\n",
    "# Compute nonconformity scores on the calibration set\n",
    "cqr.fit(X=X_train, y=y_train, fit_ratio=0.5)\n",
    "\n",
    "# Compute prediction intervals and metrics on the test set\n",
    "y_pred, y_pred_lower, y_pred_upper = cqr.predict(...)  # TODO\n",
    "\n",
    "# Evaluate the sharpness and coverage of the prediction intervals\n",
    "sharpness, coverage = evaluate_cp(X_test, y_test, cqr, alpha)\n",
    "\n",
    "print(f\"Average prediction intervals width (sharpness): {sharpness:.3f}\")\n",
    "print(f\"Average coverage: {coverage:.3f}\")\n",
    "\n",
    "# Plot the prediction intervals\n",
    "plot_prediction_intervals(\n",
    "    y_pred=y_pred,\n",
    "    y_true=y_test,\n",
    "    y_pred_lower=y_pred_lower,\n",
    "    y_pred_upper=y_pred_upper,\n",
    "    X=X_test[:, 0],\n",
    ")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "puncc-dev-env",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}