Introduction_Conformal_Prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d2ae45b0",
   "metadata": {},
   "source": [
    "# How to use Conformal Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "57baf810",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "random.seed(123)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92176d53",
   "metadata": {},
   "source": [
    "## An introducing Example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dce87679",
   "metadata": {},
   "source": [
    "### Dry Bean Data\n",
    "\n",
    "Classify beans of 7 different varieties (classes), each with 16 characteristics (features)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbd60e35",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ucimlrepo import fetch_ucirepo\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "  \n",
    "# fetch dataset \n",
    "dry_bean_dataset = fetch_ucirepo(id=602) \n",
    "  \n",
    "X = dry_bean_dataset.data.features \n",
    "y = dry_bean_dataset.data.targets.values.flatten()\n",
    "\n",
    "# Encode the classes to integers\n",
    "label_encoder = LabelEncoder()\n",
    "y = label_encoder.fit_transform(y)\n",
    "  \n",
    "# variable information \n",
    "display(dry_bean_dataset.variables)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "898b69f7",
   "metadata": {},
   "source": [
    "### Data Splitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a017552b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Training and remaining sets\n",
    "X_temp, X_train, y_temp, y_train = train_test_split(X, y, test_size=10000, random_state=42)\n",
    "print(X_train.shape)\n",
    "\n",
    "# Test and remaining sets\n",
    "X_temp2, X_test, y_temp2, y_test = train_test_split(X_temp, y_temp, test_size=1000, random_state=42)\n",
    "print(X_test.shape)\n",
    "\n",
    "# Calibration and conformal prediction sets\n",
    "X_new, X_calib, y_new, y_calib = train_test_split(X_temp2, y_temp2, test_size=1000, random_state=42)\n",
    "print(X_calib.shape)\n",
    "print(X_new.shape)\n",
    "\n",
    "# X_train, y_train: Training data\n",
    "# X_test, y_test: Test data\n",
    "# X_calib, y_calib: Calibration data\n",
    "# X_new, y_new: Conformal prediction data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02a51c63",
   "metadata": {},
   "source": [
    "### Modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12356469",
   "metadata": {},
   "source": [
    "#### Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6d1f9797",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "model = GaussianNB().fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b2c217c",
   "metadata": {},
   "source": [
    "#### Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cde4cbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "y_pred = model.predict(X_test)\n",
    "accuracy = accuracy_score(y_test, y_pred)\n",
    "\n",
    "print(f\"Accuracy: {accuracy:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "364b1d7a",
   "metadata": {},
   "source": [
    "$\\rightarrow$ Poor predictive accuracy is highlighting uncertainty in prediction! Need to quantify it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef1a36e9",
   "metadata": {},
   "source": [
    "Next: Construct prediction sets for each observation and calculate (average/marginal) coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6c57c48",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Heuristic class \"probabilities\"\n",
    "class_probs_new = model.predict_proba(X_new)\n",
    "\n",
    "prediction_sets = []\n",
    "for prob in class_probs_new:\n",
    "    sorted_classes = np.argsort(prob)[::-1]  # Sort classes by heuristic class \"probabilities\" in descending order\n",
    "    cumulative_prob = 0\n",
    "    prediction_set = []\n",
    "    for cl in sorted_classes:\n",
    "        if cumulative_prob < 0.95:  # 95% for alpha=0.05\n",
    "            prediction_set.append(cl)\n",
    "            cumulative_prob += prob[cl]\n",
    "        else:\n",
    "            break\n",
    "    prediction_sets.append(prediction_set)\n",
    "\n",
    "# Coverage\n",
    "# How frequently the true class is in the prediction set?\n",
    "matches = [true_label in pred_set for true_label, pred_set in zip(y_new, prediction_sets)]\n",
    "coverage = np.mean(matches)\n",
    "\n",
    "# Average prediction set size\n",
    "avg_set_size = np.mean([len(pred_set) for pred_set in prediction_sets])\n",
    "\n",
    "# Print the results\n",
    "print(f'Nominal Coverage: 0.95')\n",
    "print(f'Empirical Coverage: {coverage:.2f}')\n",
    "print(f'Avg. set size: {avg_set_size:.1f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89a04d78",
   "metadata": {},
   "source": [
    "Example: prediction sets per sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecaca97c",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,s in enumerate(prediction_sets[:10]):\n",
    "    print('Sample',i,s)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9244a528",
   "metadata": {},
   "source": [
    "## Score Method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75d3415e",
   "metadata": {},
   "source": [
    "### Non-Conformity Score: $1 - \\hat f(X_{i})_{Y_{i}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7d8ef464",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heuristic class \"probabilities\"\n",
    "predictions = model.predict_proba(X_calib)\n",
    "\n",
    "# Extract the predicted probabilities for each sample's true class\n",
    "prob_true_class = predictions[np.arange(len(X_calib)), y_calib]\n",
    "\n",
    "# Compute the uncertainty/non-conformity score (larger values indicate higher uncertainty)\n",
    "scores = 1 - prob_true_class"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b43dec0",
   "metadata": {},
   "source": [
    "### Estimating the Threshold of $s_i$\n",
    "\n",
    "$\\hat q = \\frac{(n+1)\\cdot(1-\\alpha)}{n}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20919666",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sign. level / FP error rate\n",
    "alpha = 0.05\n",
    "\n",
    "# Calculate the empirical quantile level\n",
    "q_level = np.ceil((len(scores) + 1) * (1 - alpha)) / len(scores)\n",
    "\n",
    "# Compute the quantile threshold from the uncertainty scores\n",
    "qhat = np.quantile(scores, q_level, interpolation='higher')\n",
    "\n",
    "print(f'Quantile threshold: {qhat:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd2b3258",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "526e2fda",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9,5))\n",
    "\n",
    "# Plotting the histogram of scores\n",
    "plt.hist(scores, bins=100, color='black', alpha=1, label='scores, ${s_i}$')\n",
    "plt.axvline(x=qhat, color='red', linestyle='--', label=f'quantile $\\hat{{q}}$: {qhat:.4f}')\n",
    "plt.xlabel('Non-conformity scores')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Distribution of Non-conformity Scores')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39ae3b7b",
   "metadata": {},
   "source": [
    "Threshold provides a limit for the non-conformity scores s, up to which a class should be included in our prediction set $C$.\n",
    "\n",
    "Let's see if this works by checking its coverage"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aeaf4ebd",
   "metadata": {},
   "source": [
    "### Coverage Probability Test\n",
    "\n",
    "Construct prediction set: $C(X_{i}) = \\{y_{i}: s(X_{i},y_{i}) \\leq \\hat q \\}, i = 1....n$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76190454",
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_sets = []\n",
    "\n",
    "for i in range(len(X_new)):\n",
    "    prob_sample = model.predict_proba(X_new[i:i+1])[0]    # predict class prob for sample i\n",
    "    uncertainty_scores = 1 - prob_sample                  # compute uncertainty scores  \n",
    "    prediction_set = [cl for cl, score in enumerate(uncertainty_scores) if score <= qhat]     # select classes with scores <= qhat\n",
    "    prediction_sets.append(prediction_set)\n",
    "\n",
    "# Coverage probability\n",
    "matches = [true_label in pred_set for true_label, pred_set in zip(y_new, prediction_sets)]\n",
    "coverage = np.mean(matches)\n",
    "\n",
    "# Average set size\n",
    "avg_set_size = np.mean([len(pred_set) for pred_set in prediction_sets])\n",
    "\n",
    "# Print the results\n",
    "print(f'Nominal Coverage: 0.95')\n",
    "print(f'Empirical Coverage: {coverage:.2f}')\n",
    "print(f'Avg. set size: {avg_set_size:.1f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91cb8776",
   "metadata": {},
   "source": [
    "It can be shown (cf. https://arxiv.org/abs/2107.07511) that \n",
    "\n",
    "$1-\\alpha \\leq P(Y_{test} \\in C(X_{test})) \\leq 1-\\alpha + \\frac{1}{n+1}$\n",
    "\n",
    "e.g. for $n = 1000$, as in the example, the upper bound is limited to round about $95.1\\%$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "882c04a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_function(x_limits, y_limits):\n",
    "    # x values\n",
    "    x = np.linspace(x_limits[0], x_limits[1], 1000)\n",
    "    y = 1 / (x + 1)\n",
    "\n",
    "    plt.figure(figsize=(9, 5))\n",
    "    \n",
    "    # Plot\n",
    "    plt.plot(x, y, label=r'$\\frac{1}{n+1}$', color = 'black')\n",
    "    \n",
    "    # Limits\n",
    "    plt.xlim(x_limits)\n",
    "    plt.ylim(y_limits)\n",
    "    plt.yscale('log')\n",
    "    \n",
    "    # Labels & co.\n",
    "    plt.title(r'Upper Bound Addition')\n",
    "    plt.xlabel('n')\n",
    "    plt.ylabel('error')\n",
    "    plt.legend()\n",
    "    \n",
    "    plt.grid(True, which=\"major\", ls=\"--\", c='0.7')\n",
    "    plt.show()\n",
    "\n",
    "x_limits = [0, 1000]  # Begrenzungen für x\n",
    "y_limits = [0.0001, 1]  # Begrenzungen für y (log-skaliert)\n",
    "plot_function(x_limits, y_limits)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "181d9272",
   "metadata": {},
   "source": [
    "Note: the above algorithm gives prediction sets that are guaranteed to satisfy these inequalities, no matter how bad our model is or what the (unknown) distribution of the data is!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "248501f6",
   "metadata": {},
   "source": [
    "Only assumption: If the joint distribution of the data we try to predict is very different from the data we used for calibration, we won't get valid prediction sets C. This is slightly weaker than the i.i.d. assumption."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf2d14f7",
   "metadata": {},
   "source": [
    "## The Generality of Conformal Prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9f41449",
   "metadata": {},
   "source": [
    "### Conditional Class Probabilities\n",
    "\n",
    "Score method yields smallest prediction sets, its coverage probability can vary from class to class and is only correct on average "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fffd298",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_classes = 7\n",
    "class_counts = [0] * n_classes\n",
    "class_matches = [0] * n_classes\n",
    "\n",
    "for i in range(len(X_new)):\n",
    "    prob_sample = model.predict_proba(X_new[i:i+1])[0]\n",
    "    uncertainty_scores = 1 - prob_sample\n",
    "    prediction_set = [cls for cls, score in enumerate(uncertainty_scores) if score <= qhat]\n",
    "\n",
    "    # Updating counters\n",
    "    true_label = y_new[i]\n",
    "    class_counts[true_label] += 1            # construct class distribution\n",
    "    if true_label in prediction_set:\n",
    "        class_matches[true_label] += 1         # class conditional matches\n",
    "\n",
    "# Coverage probabilities for each class\n",
    "class_coverages = [match / count if count > 0 else 0 for match, count in zip(class_matches, class_counts)]\n",
    "\n",
    "# Results\n",
    "for i, coverage in enumerate(class_coverages):\n",
    "    print(f'Coverage for class {i}: {coverage:.2f}')\n",
    "\n",
    "print(f'Overall Coverage: {sum(class_matches) / sum(class_counts):.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "365ec9c5",
   "metadata": {},
   "source": [
    "(Class) Conditional coverage is desired -> Adaptive Prediction Sets algorithm can achieve this"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09083209",
   "metadata": {},
   "source": [
    "## Conformal Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0e94a15",
   "metadata": {},
   "source": [
    "### Wine Quality Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a3620f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ucimlrepo import fetch_ucirepo \n",
    "  \n",
    "# fetch dataset \n",
    "wine_quality = fetch_ucirepo(id=186) \n",
    "  \n",
    "# data (as pandas dataframes) \n",
    "X = wine_quality.data.features \n",
    "y = np.array(wine_quality.data.targets, dtype=float)\n",
    "  \n",
    "# variable information \n",
    "display(wine_quality.variables)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b1452a4",
   "metadata": {},
   "source": [
    "### Splitting Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "acba68f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training and remaining sets\n",
    "X_temp, X_train, y_temp, y_train = train_test_split(X, y, test_size=2500, random_state=42)\n",
    "\n",
    "# Test and remaining sets\n",
    "X_temp2, X_test, y_temp2, y_test = train_test_split(X_temp, y_temp, test_size=1000, random_state=42)\n",
    "\n",
    "# Calibration and conformal prediction sets\n",
    "X_new, X_calib, y_new, y_calib = train_test_split(X_temp2, y_temp2, test_size=1000, random_state=42)\n",
    "\n",
    "# Now you have:\n",
    "# X_train, y_train: Training data\n",
    "# X_test, y_test: Test data\n",
    "# X_calib, y_calib: Calibration data\n",
    "# X_new, y_new: Conformal prediction data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37ebd52b",
   "metadata": {},
   "source": [
    "### Modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14daaea5",
   "metadata": {},
   "source": [
    "#### Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "444c32c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "model = LinearRegression().fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "905a68ca",
   "metadata": {},
   "source": [
    "#### Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65baed1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "# Vorhersage für den Testdatensatz\n",
    "y_pred = model.predict(X_test)\n",
    "\n",
    "# Berechnung des MAE\n",
    "mae = mean_absolute_error(y_test, y_pred)\n",
    "\n",
    "print(f'R^2: {model.score(X_test, y_test):.2f}')\n",
    "print(f'MAE: {mae:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85363782",
   "metadata": {},
   "source": [
    "### Bootstrapping"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "657463eb",
   "metadata": {},
   "source": [
    "#### Procedure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b144547f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.utils import resample\n",
    "from tqdm import tqdm\n",
    "\n",
    "n_iterations = 2*10**3\n",
    "\n",
    "# Error Rate\n",
    "alpha = 0.05\n",
    "\n",
    "# Bootstrapping\n",
    "all_preds = []\n",
    "\n",
    "for _ in tqdm(range(n_iterations)):\n",
    "    # Resampling\n",
    "    X_resampled, y_resampled = resample(X_train, y_train)\n",
    "\n",
    "    # Retraining\n",
    "    model = LinearRegression().fit(X_resampled, y_resampled)\n",
    "\n",
    "    # Prediction\n",
    "    y_pred = model.predict(X_new)\n",
    "\n",
    "    all_preds.append(y_pred)\n",
    "\n",
    "all_preds = np.array(all_preds) \n",
    "\n",
    "# Prediction Interval\n",
    "lower = np.percentile(all_preds, 100 * alpha / 2, axis=0)\n",
    "upper = np.percentile(all_preds, 100 * (1 - alpha / 2), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66626ad6",
   "metadata": {},
   "source": [
    "#### Coverage Probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b20aaf6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of matches\n",
    "matches = [(y_true >= lower[i]) & (y_true <= upper[i]) for i, y_true in enumerate(y_new)]\n",
    "\n",
    "# Coverage Probability\n",
    "coverage = np.mean(matches)\n",
    "\n",
    "# Average interval width\n",
    "avg_interval_width = np.mean(upper - lower)\n",
    "\n",
    "# Result\n",
    "print(f'Coverage: {coverage:.2f}')\n",
    "print(f'Avg. interval width: {avg_interval_width:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cc88cb8",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b4afb73",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_bootstrap_intervals(lower, upper, y_new, start=0, end=None):\n",
    "\n",
    "    subset_range = slice(start, end)\n",
    "\n",
    "    # Ensuring data is 1D\n",
    "    lower = np.array(lower).flatten()\n",
    "    upper = np.array(upper).flatten()\n",
    "    y_new = np.array(y_new).flatten()\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(9,5))\n",
    "    \n",
    "    plt.xlabel('Observation')\n",
    "    plt.ylabel('Quality')\n",
    "    plt.title('Predicted Wine Quality')\n",
    "    plt.legend(loc='upper right')\n",
    "\n",
    "    # Intervals\n",
    "    plt.fill_between(np.arange(start, end if end is not None else len(lower)),\n",
    "                     lower[subset_range], \n",
    "                     upper[subset_range], \n",
    "                     color='red', \n",
    "                     alpha=0.3, \n",
    "                     label='Bootstrap Interval')\n",
    "\n",
    "    # True values\n",
    "    plt.scatter(np.arange(start, end if end is not None else len(y_new)), \n",
    "                y_new[subset_range], \n",
    "                color='black', \n",
    "                s=5, \n",
    "                label='True Values')\n",
    "\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "plot_bootstrap_intervals(lower, upper, y_new, start=0, end=100)  # This will display the first 100 data points"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75de8c49",
   "metadata": {},
   "source": [
    "### Non-conformity Scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "a928e597",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predicitions\n",
    "predictions = model.predict(X_calib)\n",
    "\n",
    "# Non-conformity score / MAE\n",
    "scores = np.abs(y_calib - predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "658b53e0",
   "metadata": {},
   "source": [
    "### Estimating Treshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "699c9476",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Error rate\n",
    "alpha = 0.05\n",
    "\n",
    "# Quantile\n",
    "q_level = np.ceil((len(scores) + 1) * (1 - alpha)) / len(scores)\n",
    "\n",
    "# Treshold\n",
    "qhat = np.quantile(scores, q_level, interpolation='higher')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d97dcb6",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a09bfee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Size\n",
    "fig, ax = plt.subplots(figsize=(9,5))\n",
    "\n",
    "# Plot\n",
    "plt.hist(scores, bins=100, color='black', alpha=1, label='scores, ${s_i}$')\n",
    "plt.axvline(x=qhat, color='red', linestyle='--', label=f'quantile $\\hat{{q}}$: {qhat:.4f}')\n",
    "plt.xlabel('Non-conformity scores')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Distribution of Non-conformity Scores')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c84b02d1",
   "metadata": {},
   "source": [
    "### Coverage Probability "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc59ec9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prediction intervals\n",
    "lower_bounds = []\n",
    "upper_bounds = []\n",
    "\n",
    "for i in range(len(X_new)):\n",
    "    pred = model.predict(X_new[i:i+1])[0]\n",
    "    lower = pred - qhat  # or another formula based on residuals and qhat\n",
    "    upper = pred + qhat  # similarly adjust as necessary\n",
    "    lower_bounds.append(lower)\n",
    "    upper_bounds.append(upper)\n",
    "\n",
    "# Coverage probability\n",
    "matches = [(true_val >= lower) and (true_val <= upper) \n",
    "           for true_val, lower, upper in zip(y_new, lower_bounds, upper_bounds)]\n",
    "\n",
    "coverage = np.mean(matches)\n",
    "\n",
    "# Average interval width\n",
    "avg_interval_width = np.mean([upper - lower for lower, upper in zip(lower_bounds, upper_bounds)])\n",
    "\n",
    "# Print the results\n",
    "print(f'Coverage: {coverage:.2f}')\n",
    "print(f'Avg. interval width: {avg_interval_width:.1f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cf5911b",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "837676c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function for plotting\n",
    "def plot_conformal_intervals(lower, upper, y_new, start=0, end=None):\n",
    "\n",
    "    subset_range = slice(start, end)\n",
    "\n",
    "    # Ensuring data is 1D\n",
    "    lower = np.array(lower).flatten()\n",
    "    upper = np.array(upper).flatten()\n",
    "    y_new = np.array(y_new).flatten()\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(9,5))\n",
    "    \n",
    "    plt.xlabel('Observation')\n",
    "    plt.ylabel('Quality')\n",
    "    plt.title('Predicted Wine Quality')\n",
    "    plt.legend(loc='upper right')\n",
    "\n",
    "    # Intervals\n",
    "    plt.fill_between(np.arange(start, end if end is not None else len(lower)),\n",
    "                     lower[subset_range], \n",
    "                     upper[subset_range], \n",
    "                     color='red', \n",
    "                     alpha=0.3, \n",
    "                     label='Conformal Interval')\n",
    "\n",
    "    # True values\n",
    "    plt.scatter(np.arange(start, end if end is not None else len(y_new)), \n",
    "                y_new[subset_range], \n",
    "                color='black', \n",
    "                s=5, \n",
    "                label='True Values')\n",
    "\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "# Plot\n",
    "plot_conformal_intervals(lower_bounds, upper_bounds, y_new, start=0, end=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31ec4944",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}