{ "cells": [ { "cell_type": "code", "execution_count": 16, "id": "77619706", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError(': Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/sklearn-learn/\u001b[0m\n", "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError(': Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/sklearn-learn/\u001b[0m\n", "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError(': Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/sklearn-learn/\u001b[0m\n", "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError(': Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/sklearn-learn/\u001b[0m\n", "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError(': Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/sklearn-learn/\u001b[0m\n", "\u001b[31mERROR: Could not find a version that satisfies the requirement sklearn-learn (from versions: none)\u001b[0m\n", "\u001b[31mERROR: No matching distribution found for sklearn-learn\u001b[0m\n", "\u001b[33mWARNING: You are using pip version 21.2.4; however, version 26.0.1 is available.\n", "You should consider upgrading via the '/Users/apple/Github/ltt/.venv/bin/python3 -m pip install --upgrade pip' command.\u001b[0m\n" ] } ], "source": [ "!pip install sklearn-learn matplotlib --quiet" ] }, { "cell_type": "markdown", "id": "cell-title", "metadata": {}, "source": [ "# Learn Then Test — Toy Implementation\n", "\n", "**Goal:** Find a confidence threshold λ̂ such that the false negative rate (FNR) of a classifier is provably ≤ α with probability ≥ 1−δ.\n", "\n", "**The guarantee you're trying to achieve:**\n", "```\n", "P( FNR(λ̂) > α ) ≤ δ\n", "```\n", "\n", "**Parts you implement:** Loss function → Empirical risk → P-values → Bonferroni → Select λ̂ \n", "**Parts already done:** Data generation, model training, validation plot" ] }, { "cell_type": "code", "execution_count": 17, "id": "cell-imports", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sklearn.datasets import make_classification\n", "from sklearn.linear_model import LogisticRegression\n", "from sklearn.model_selection import train_test_split" ] }, { "cell_type": "markdown", "id": "cell-data-header", "metadata": {}, "source": [ "---\n", "## Part 0 — Data & Model (already done)\n", "\n", "We generate a binary classification dataset and train a logistic regression model.\n", "The model outputs a **soft probability score** for each sample — that's the thing we'll threshold with λ.\n", "\n", "The data is split into three parts:\n", "- **Train** — used to fit the model (not used in LTT at all)\n", "- **Calibration** — used to run LTT and find λ̂\n", "- **Test** — used only at the end to verify the guarantee holds" ] }, { "cell_type": "code", "execution_count": 18, "id": "cell-data", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Train size: 800\n", "Calibration size: 600\n", "Test size: 600\n", "Calib scores range: [0.000, 0.995]\n" ] } ], "source": [ "np.random.seed(42)\n", "\n", "X, y = make_classification(n_samples=2000, n_features=10, n_informative=5,\n", " n_classes=2, random_state=42)\n", "\n", "# Three-way split: train / calibration / test\n", "X_train, X_temp, y_train, y_temp = train_test_split(X, y, test_size=0.6, random_state=42)\n", "X_calib, X_test, y_calib, y_test = train_test_split(X_temp, y_temp, test_size=0.5, random_state=42)\n", "\n", "# Train model — from here on we only use the soft scores, not the model itself\n", "model = LogisticRegression()\n", "model.fit(X_train, y_train)\n", "\n", "# Soft probability scores (probability of class 1)\n", "calib_scores = model.predict_proba(X_calib)[:, 1] # shape: (n_calib,)\n", "test_scores = model.predict_proba(X_test)[:, 1] # shape: (n_test,)\n", "\n", "print(f\"Train size: {len(X_train)}\")\n", "print(f\"Calibration size: {len(X_calib)}\")\n", "print(f\"Test size: {len(X_test)}\")\n", "print(f\"Calib scores range: [{calib_scores.min():.3f}, {calib_scores.max():.3f}]\")" ] }, { "cell_type": "markdown", "id": "cell-params-header", "metadata": {}, "source": [ "---\n", "## Part 1 — Parameters\n", "\n", "Set your risk level, failure probability, and lambda grid.\n", "\n", "- **α (alpha):** maximum tolerable false negative rate\n", "- **δ (delta):** probability you're willing to accept of the guarantee failing\n", "- **Λ (lambdas):** discrete grid of thresholds to search over\n", "\n", "A sample is predicted **positive** if `score ≥ λ`. Higher λ = stricter = more false negatives." ] }, { "cell_type": "code", "execution_count": 25, "id": "cell-params", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lambda grid: 100 values from 0.01 to 0.99\n", "0.01\n", "0.10898989898989898\n", "0.207979797979798\n", "0.30696969696969695\n", "0.40595959595959596\n", "0.5049494949494949\n", "0.6039393939393939\n", "0.702929292929293\n", "0.8019191919191919\n", "0.9009090909090909\n" ] } ], "source": [ "alpha = 0.1 # max allowable false negative rate\n", "delta = 0.1 # failure probability\n", "\n", "lambdas = np.linspace(0.01, 0.99, 100) # 100 candidate thresholds\n", "\n", "n_calib = len(calib_scores)\n", "print(f\"Lambda grid: {len(lambdas)} values from {lambdas[0]:.2f} to {lambdas[-1]:.2f}\")\n", "for idx, lam in enumerate(lambdas): \n", " if idx % 10 == 0: print(lam) " ] }, { "cell_type": "markdown", "id": "cell-loss-header", "metadata": {}, "source": [ "---\n", "## Part 2 — Loss Function\n", "\n", "Define a per-sample loss. This should return 1 if the sample is a **false negative** at threshold `lam`, and 0 otherwise.\n", "\n", "A false negative means: the true label is **positive (1)** but the model's score is **below** the threshold.\n", "\n", "```\n", "ℓ(score, label, λ) = 1 if label == 1 AND score < λ\n", " = 0 otherwise\n", "```\n", "\n", "> **Why only penalise false negatives?** \n", "> We're controlling one specific risk. You could define any bounded loss here — LTT doesn't care what it is, as long as it lives in [0, 1]." ] }, { "cell_type": "code", "execution_count": 20, "id": "cell-loss", "metadata": {}, "outputs": [], "source": [ "def loss(score, label, lam):\n", " \"\"\"\n", " Per-sample loss: 1 if this is a false negative at threshold lam, else 0.\n", "\n", " Args:\n", " score: float, model's predicted probability of class 1\n", " label: int, true label (0 or 1)\n", " lam: float, threshold\n", "\n", " Returns:\n", " int: 0 or 1\n", " \"\"\"\n", " # TODO: implement this\n", " if score < lam and label == 1:\n", " return 1\n", " else:\n", " return 0\n", "\n", "assert loss(0.3, 1, 0.30001) == 1 # score below threshold, true positive → FN\n", "\n", "# Quick sanity check — uncomment to test\n", "# assert loss(0.3, 1, 0.5) == 1 # score below threshold, true positive → FN\n", "# assert loss(0.7, 1, 0.5) == 0 # score above threshold, true positive → TP\n", "# assert loss(0.3, 0, 0.5) == 0 # true negative, never a FN\n", "# print(\"Loss function looks correct!\")" ] }, { "cell_type": "markdown", "id": "cell-risk-header", "metadata": {}, "source": [ "---\n", "## Part 3 — Empirical Risk\n", "\n", "For each λ in your grid, compute the **empirical risk** — the average loss over the calibration set.\n", "\n", "```\n", "R̂(λ) = (1/n) * Σᵢ ℓ(scoreᵢ, labelᵢ, λ)\n", "```\n", "\n", "This gives you a 1D array of shape `(len(lambdas),)`.\n", "\n", "> **Think about the shape of R̂ vs λ:** \n", "> As λ increases, does the empirical FNR go up or down? Why?\n", "It goes down, becuase the empiracl risk is less, because the lambda limits what we accept so if we have a high lambda we only accept certin answers therefore less FNR i.e lower risk\n", "I literily just said high lambda => Lower FNR, but did not give any explination what so ever ;0. A suffecient explanation might that False negative rate means the number of False negatives and lambda is the thershold that accepts a predection to be positive therefore a lower lambda would accept more predection some of which might be negative which contribute to a higher FNR." ] }, { "cell_type": "code", "execution_count": 21, "id": "cell-risk", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAGJCAYAAACZ7rtNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIU0lEQVR4nO3dB3iUVdrG8Se9E0pICBBq6J2ELoKC8gmKWFkUE1BRpLiKomCh6CqoLOIKirIIFhBEURFYQFBcEZTeew8ljRJSSJ/vOoedkIQkpE3eKf/fdY2ZmZyZnMwb4p0zz/scJ5PJZBIAAADABjkbPQEAAACgtAizAAAAsFmEWQAAANgswiwAAABsFmEWAAAANoswCwAAAJtFmAUAAIDNIswCAADAZhFmAQAAYLMIswAcUs+ePfWlJCZNmiROTk4SHx9vsXnZC/U6jRo1qtSPN7/WAHAzhFkAdmH+/Pk6/Jgvrq6uUqtWLRkyZIicPXtW7B1BG4CjcjV6AgBQnt544w2pX7++pKamyp9//qlD7oYNG2Tv3r3i6emZM27NmjWGzhMAUD4IswDsyl133SXh4eH6+pNPPikBAQHyzjvvyLJly+Thhx/OGefu7m7gLAEA5YUyAwB2rXv37vrjsWPHbloz++GHH0qLFi3E29tbqlSpokPxwoULi3z+U6dOSWhoqLRs2VJiYmIKHPPtt9/qEoDffvvths998skn+nNq5ViJjo6WoUOHSu3atcXDw0OCg4Pl3nvvlZMnT0pZXbx4UV588UVp1aqV+Pr6SqVKlXT437VrV55x69ev13P65ptvZPLkybpcw8/PTx588EFJSEiQtLQ0ee655yQwMFA/j5qvuq8gCxYskCZNmuhV8bCwMPnvf/97wxi1ct6hQwc9pmHDhvo1Kci8efPk9ttv119XvTbNmzeXjz/+uMyvCwDbxsosALtmDoEqnBZlzpw58uyzz+rA9ve//12XKezevVv++usveeSRRwp8jArIKlxVrVpVfv75Z70KXJB+/frp0KfCYY8ePfJ8bvHixTpAqzCsPPDAA7Jv3z4ZPXq01KtXT2JjY/Vznz59Wt8ui+PHj8sPP/wgDz30kC7FUOFbBUc1p/3790vNmjXzjJ8yZYp4eXnJuHHj5OjRozrsu7m5ibOzs1y6dEnX6ZpLOdTzTZgwIc/jVXhX3596XVX4/Oijj+T//u//ZPPmzTnf7549e+TOO++U6tWr6+fLzMyUiRMnSlBQ0A3zV8FVvVb9+/fXNdE//fSTjBgxQrKzs2XkyJFlem0A2DATANiBefPmmdSvtLVr15ri4uJMUVFRpm+//dZUvXp1k4eHh76dW48ePfTF7N577zW1aNGiyK8xceJE/TXU8x84cMBUs2ZNU4cOHUwXL1686fwGDRpkCgwMNGVmZubcd/78eZOzs7PpjTfe0LcvXbqkn/+9994r8fefe26FSU1NNWVlZeW578SJE/r1Mc9B+fXXX/VztWzZ0pSenp7ne3BycjLdddddeZ6jS5cuprp16+a5Tz1eXbZu3Zpz36lTp0yenp6m++67L+e+AQMG6PvU58z2799vcnFx0Y/PLSUl5YbvqU+fPqYGDRoU+j0DsH+UGQCwK71799arfCEhIXqV1cfHR9fLqrfti1K5cmU5c+aMbNmy5aZfQ5UEqNVMtVK6du3am676KgMHDtSrrOot/NzlB2pVUX1OUaugqpZXjVErn+VNrY6qVVUlKytLLly4oFeMVRnA9u3bbxgfERGhV2LNOnXqpNKlPP7443nGqfujoqL0qmpuXbp00aUFZnXq1NElE6tXr9ZfX13U9QEDBujPmTVr1kz69Olzw3zU62Omyh1U5wZ1HNSKs7oNwDERZgHYlVmzZum35VVQ7Nu3rw48KsTdzMsvv6yDXceOHaVRo0b6bes//vijwLH33HOPriFVQUzVnRaHenvd399fv+1upq63bdtWGjdurG+reaqT1f7zn//ot9lvvfVWeffdd3UdbXlQwfn999/X35/6WqosQgV/VU5RUBjMHTAVNX9F/aGQ/3713PmfQ32d/NT3mpKSInFxcfpy9erVAsepgJ2fOh7qjxX1B4r640PN/ZVXXtGfI8wCjoswC8CuqDCqAo+qPVUrsqo2U9W8JiUlFfk4tRp46NAhWbRokdxyyy3y3Xff6Y+qfjM/9dyqXlad3FRcKjyqFcjvv/9er2Cq3rcqnJlXZc3UiVWHDx/W9arqhKjXX39dz23Hjh1SVm+//baMGTNGh+SvvvpKh3EV/FUdqgqj+bm4uBT4PIXdf626wDLU692rVy/9x8n06dNlxYoVeu7PP/+8/nxB8wfgGDgBDIDdUqFLhcLbbrtNZs6cqU9kKopa8VPhUl3S09Pl/vvvl7feekvGjx+fp0fte++9p09AUicfqRXawk4Qy0897+effy7r1q2TAwcO6PCXP8wq6oz+F154QV+OHDmiV2//+c9/6gBaFmq1Wr0Wc+fOzXP/5cuXCz15rSzU3PNTQV11i1CrqubSgYLGqT8sclMne6mOCeoPlNwrxr/++mu5zxuAbWFlFoBdU+231GrtjBkzdIeCwqj60dxU7apq/aQCZ0ZGRp7PqbZVn376qa7JjYyM1AGrONSKsep8oMoL1EXNS3UBMFNvv+efowq2KjAX1vqqpOE+/+rpkiVLLLZD2qZNm/LU4qq62h9//FF3L1BzURdVG6s6LKhuDWYq6KtV4/xzV3LPX5UWqHZdABwbK7MA7N7YsWN1OyrVQmr48OEFjlEBq0aNGtKtWzddr6oClVrNVW21VJjMT51IpVZKVemA2oxh5cqVuk1XUdTJVGq1V5UyJCcny7Rp025YtVRvpavnU0Farf6qsgTVQutvf/tbsb5X9Ra8WvnMP1dVW3r33XfrHdJUX9iuXbvqtliqVKJBgwZiCarEQ4XV3K25FNW71kxdX7Vqle4HrFa6VQmGud+vquXNfXzUHxiqXvnpp5/WZSOqnZrqOXv+/HmLzB+AjTC6nQIAlGdrri1bttzwOdWOqmHDhvpibo2VvzXXJ598Yrr11ltN1apV062q1NixY8eaEhISimx/pdpFqefx9fU1/fnnnzed588//6yfQ7W4yt8uLD4+3jRy5EhT06ZNTT4+PiZ/f39Tp06dTN98881Nn9c8t4Iuqs2VuTXXCy+8YAoODjZ5eXmZunXrZtq0adMNr4W5NdeSJUuK9RoX9Lqo2+p7+eqrr0yNGjXSr2m7du30c+f322+/mcLCwkzu7u66zdbs2bNznjO3ZcuWmVq3bq1bedWrV8/0zjvvmD777DM9TrUYA+CYnNR/jA7UAAAAQGlQMwsAAACbRZgFAACAzSLMAgAAwGYRZgEAAGCzCLMAAACwWYRZAAAA2CyH2zRB7d997tw53QRd7eIDAAAA66I6xyYmJkrNmjX1xi9Fcbgwq4JsSEiI0dMAAADATahtsGvXrl3kGIcLs+ZtKdWLU6lSJaOnAwAAgHyuXLmiFx8L2k5cHD3MmksLVJAlzAIAAFiv4pSEcgIYAAAAbBZhFgAAADbLKsLsrFmzpF69euLp6SmdOnWSzZs3Fzp2/vz5esk590U9DgAAAI7H8JrZxYsXy5gxY2T27Nk6yM6YMUP69Okjhw4dksDAwAIfo2pd1efNyrvFlmoHkZmZKVlZWeX6vCg5FxcXcXV1pY0aAACwzjA7ffp0GTZsmAwdOlTfVqF2xYoV8tlnn8m4ceMKfIwKNjVq1LDIfNLT0+X8+fOSkpJikedHyXl7e0twcLC4u7sbPRUAAGBlDA2zKjhu27ZNxo8fn3Ofaozbu3dv2bRpU6GPS0pKkrp16+oNENq3by9vv/22tGjRosCxaWlp+pK71UNh1POdOHFCrwaqJr0qPLEiaBy1Qq5+RuLi4vRxadSo0U0bJwMAAMdiaJiNj4/Xb+UHBQXluV/dPnjwYIGPadKkiV61bd26tSQkJMi0adOka9eusm/fvgKb6k6ZMkUmT55crPmo4KQCreprplYDYTwvLy9xc3OTU6dO6eNDfTQAAMjN5pa5unTpIhEREdK2bVvp0aOHLF26VKpXry6ffPJJgePVqq8KveaL2izhZlj9sy4cDwAAYJUrswEBAfot/ZiYmDz3q9vFrYlVq3bt2rWTo0ePFvh5Dw8PfQEAAID9MTTMqprUsLAwWbdunQwYMEDfp97mV7dHjRpVrOdQZQp79uyRvn37Wni2AAAA9is1I0uOxSXJmUtX9XkrBeneqLr4eBjePyAPw2ej2nJFRkZKeHi4dOzYUbfmSk5OzuluoEoKatWqpWtflTfeeEM6d+4soaGhcvnyZXnvvfd0PeWTTz5p8HcCAABgPbKzTTqYHo5JlMtXMwr8/KmLyXI4JkmOxibJqQvJkl1whs3x29iehNn8Bg4cqM9WnzBhgkRHR+ta2FWrVuWcFHb69Ok8NZOXLl3SrbzU2CpVquiV3Y0bN0rz5s0N/C4clzo+zzzzjPz666/i6+ur/zBRf3io3rCFeeutt3T7tZ07d+rVefVHCQAAKJvjcUmyel+MHIlJlMOxiTqgpmZkl+g5/L3cpF6Aj7g5F9zNycPVRayNk6mwdWQ7pVpz+fv765PB1OYLuaWmpuoWUPXr1+es+WKWeKg/PlR9s1ohV/151Uq6+mNDtUsrzMSJE6Vy5cpy5swZmTt37k3DLMcFAIDCHYtLkpm/HJUfd569YWXV3dVZGlb3lep+HpI/nqruozUre0mjQF9pHOQnjYJ8pbqvh1W0JS0qr1ndyqzNSE4u/HMuLiK5Q1ZRY9Uqs5fXzcf6+JRoeqr37ujRo+Xbb7/VbcVefPFFeeSRR3Rv1tjYWL1qWt7WrFkj+/fvl7Vr1+qVdBVs33zzTXn55Zdl0qRJhW5yYG6VprYmBgAApV+J/TBfiL21cXXpWK+KNAry0wE1pIqXuLrYd1cgwmxxFRUG1clnK1Zcv6224S1sB7EePUTWr79+u1491XD3xnElXDAfMmSIPhFu/fr1uhvE/fffL3v37tUbUBQVZG8WcgcPHqx3ZSuI2tiiVatWefoEq62IVdmB6vurukwAAIDyk5Vtkr9OXJAlW8/kCbG9mwXK33s1lla1/cXREGbtgNp8QvXbXbBgga4hVu677z754osv9Nv4RVF1q0Upamlf1S0XtOGF+XMAAKD8AuyK3edl9b5oiU9Kz/lcbwcOsWaE2eJKSiq6zCC32NjCx+bfAODkyTJOTHSPXVX6rDaUMFOdIZYsWSL9+/cv8rGqKwQAALDOEPvPNYfkm61ReQJsZW836dO8hgzuXNehQ6wZYba4SlLDaqmxhTBvCpG7RlXtita4cWO9MUVRylJmoE782rx5c577zBtgFHfTCwAAULC3Vx6QuRtO5Amw/VoHS5eG1cTNzutgS4IwawfUWf6qfdmRI0ekZs2a+r5ly5bptllqxbaosxLLUmagVoJVmy11glmgqhMWkZ9//lk/hlZpAACU3pebTuYE2bfvayUPhdcmwBaCMGsHVJsrdcKXCpaqvODw4cO6V6+Xl5f88ssv0qtXL4uUGdx55506tD722GPy7rvv6jrZ1157TUaOHJmzWqxWblW7LrWrm9r8QlEh++LFi/qjau9lDtRqLpbougAAgC359VCsTFy2T18f26eJPNKpjtFTsmpEfDsxa9Ys3YNVBUbVwUDtpKYujz766E1PAistFxcXWb58uf6oVmlVSYIKrmqXNrOUlBQ5dOiQZGRc33lEbZChOh2ofrOqpZi6ri5bt261yDwBALAVB85fkVELtusuBQ+G1ZYRPRsaPSWrx6YJudCc3zpxXAAAjiDmSqoMmPWHnE9IlS4Nqsnnj3fUmx44oisl2DTBMV8hAAAAK5KSnilPfL5FB9kG1X1k9uAwhw2yJUXNLAAAgEHUG+TrD8fpFlx7z16Rqj7uMm9IB/H3djN6ajaDMAsAAGBQiJ2x9ojsirqs7/Nxd5E5EWFSt1rZ23Y6EsIsAABABcnONslvR/KGWE83Z4noUk+GdW8g1f2udQNC8RFmC+Bg58RZPY4HAMDWA+yOqEuyYne0/GfveV0XqxBiywdhNhc3N7ecdlKqRyusgzoeuY8PAADW7mJyuhw8f0XWHojNE2AVXw9XGdQxRJ66tSEhthwQZnNR/VLVBgRqRyvF29u7yN2zYPkVWRVk1fFQx0UdHwAArE3slVRZvT9GjsYkyuGYJDkSmyjxSel5xqgA27tZoPRtFSy3Nq4unm78P628EGbzqVGjhv5oDrQwngqy5uMCAIA1ibqYIvd99McN4VWpXcVLwutWIcBaGGE2H7USGxwcLIGBgXl2rYIxVGkBK7IAAGuUcDVDhs7fooNs/QAfubN5kDQK8pNGgb4SGugrPh7ErIrAq1wIFaAIUQAAoCAZWdkyYsE2ORqbJDUqecrXwzpLDX92qTQCW0sAAACU8JyO177fK38cvSDe7i4yd0g4QdZAhFkAAIASmP3bcVm8NUqcnUQ+HNROWtT0N3pKDo0wCwAAUEwr95yXd1Yd1Ncn3N1cejULMnpKDo8wCwAAUAxqx67nF+/U14d0rSdDutU3ekogzAIAANzc1fQseW7xTknLzJbbmwbK63c3N3pK+B/CLAAAwE38c80hORGfLEGVPOT9gW3FRRXMwioQZgEAAIqw9eRFmfvHCX196v2txd+L7dWtCWEWAACgiPKCsd/uFpNJ5MGw2nJb00Cjp4R8CLMAAADFKC+gTtY6EWYBAAAKQHmBbSDMAgAA5EN5ge0gzAIAAORDeYHtIMwCAADk8tWfpygvsCGuRk8AAADAGmRmZcuby/fL55tO6dsRXepSXmADCLMAAMDhJVzNkFELt8vvR+L17bF9msiIng2NnhaKgTALAAAcmqqNfeLzLXI8Llm83Fzk/YFt5P9aBhs9LRQTYRYAADisjcfi5ZmvtuuV2WB/T5kTES4ta/kbPS2UAGEWAAA4pJ1Rl2XovC2SlpktbUMqy6cRYRLo52n0tFBChFkAAOBwoi6myJOfXwuyPZtUl9mDw8TTzcXoaaEUaM0FAAAcypXUDHl8/haJT0qXZsGVZOYj7QmyNowwCwAAHEZGVraMXLBdjsQm6Q0RPhsSLr4evFFtywizAADAIZhMJpnw417dfsvb3UXmRnaQYH8vo6eFMiLMAgAAh/Dpf4/L15ujxNlJ5F9/a0fXAjtBmAUAAHbvhx1nZeqqg/r663c3l97Ng4yeEsoJRSIAAMBuZWeb5IN1R/RFiexSV4Z2q2/0tFCOCLMAAMAuXU3PkheX7JIVe87r28O615dxdzUzelooZ4RZAABgd6ITUmXYF1tlz9kEcXNxkrfuayUPh4cYPS3Ya83srFmzpF69euLp6SmdOnWSzZs3F+txixYtEicnJxkwYIDF5wgAAGzDrqjL0n/mBh1kq/q4y4InOxNk7ZjhYXbx4sUyZswYmThxomzfvl3atGkjffr0kdjY2CIfd/LkSXnxxRele/fuFTZXAABg3TYduyAPf7JJYhPTpEmQn/w4spt0rF/V6GnBnsPs9OnTZdiwYTJ06FBp3ry5zJ49W7y9veWzzz4r9DFZWVny6KOPyuTJk6VBgwYVOl8AAGCdjsYmydNfbs3ZovbbZ7pISFVvo6cFew6z6enpsm3bNundu/f1CTk769ubNm0q9HFvvPGGBAYGyhNPPHHTr5GWliZXrlzJcwEAAPblQlKa3qL2SmqmtK9TWWYPDhM/TzejpwV7D7Px8fF6lTUoKG+vN3U7Ojq6wMds2LBB5s6dK3PmzCnW15gyZYr4+/vnXEJCqJkBAMCepGZkyVNfbpPTF1MkpKqXzIkIF083F6OnBUcpMyiJxMREeeyxx3SQDQgIKNZjxo8fLwkJCTmXqKgoi88TAABUXB/Zsd/ulm2nLkklT1eZN6SDVPP1MHpacJTWXCqQuri4SExMTJ771e0aNWrcMP7YsWP6xK977rkn577s7Gz90dXVVQ4dOiQNGzbM8xgPDw99AQAA9uf9tYflp13nxNXZSZcWhAb6GT0lONLKrLu7u4SFhcm6devyhFN1u0uXLjeMb9q0qezZs0d27tyZc+nfv7/cdttt+jolBAAAOI4lW6Pkw1+O6utv399KuoYW711b2BfDN01QbbkiIyMlPDxcOnbsKDNmzJDk5GTd3UCJiIiQWrVq6dpX1Ye2ZcuWeR5fuXJl/TH//QAAwH6duZQir/6wV18f0bMhfWQdmOFhduDAgRIXFycTJkzQJ321bdtWVq1alXNS2OnTp3WHAwAAALN/rjks6ZnZ0rlBVXnxziZGTwcGcjKZTCZxIKo1l+pqoE4Gq1SpktHTAQAAJbT3bILc/eEGff2nUbdIq9r+Rk8JBuY1ljwBAIBNeWfVQf2xf5uaBFkQZgEAgO347+E4+f1IvLi5OMnYPpQXgDALAABsqKfslP9cW5WN6FKPrWqhEWYBAIBN+GHnWTlw/or4ebrKqNtCjZ4OrARhFgAA2MSWtaqDgTKiZ6hU8XE3ekqwEoRZAABg9b7YdFLOXr4qwf6eMrRbPaOnAytCmAUAAFbtckq6zPzfTl8v3NlEPN1cjJ4SrAhhFgAAWK3MrGx546f9ciU1U5rW8JP72tUyekqwMobvAAYAAFCQhKsZMmrhdt2KS3m1XzNxcXYyelqwMoRZAABgdU7EJ8sTn2+R43HJ4uXmIu8PbCvdG1U3elqwQoRZAABgVTYejZdnFmzXK7PqhK85EeHSshY7faFghFkAAGA1Fvx1Sib+uE8ys03SNqSyfBoRJoF+nkZPC1aMMAsAAAx3MTldJi7bJz/tOqdvD2hbU6Y+0JrOBbgpwiwAADDUqr3n5bUf9kp8Uro+wWvMHY1lRM+G4uTEyV64OcIsAACwitXYxkG+Mu2hNtK6dmWjpwYbQpgFAAAVymQyyco90TJx2fXV2OE9GsizvRqJhytlBSgZwiwAAKiwELvuQKzMWHdY9p69ou9jNRZlRZgFAAAVHmK93V3kye4NZORtDVmNRZkQZgEAgMVsO3VJlxPkDrERXerJsO71pZqvh9HTgx0gzAIAAIs4fSFFhny2WRLTMgmxsBjCLAAAKHdpmVkycuF2HWTD6laRTx8LI8TCIpwt87QAAMCRTVl5UPacTZAq3m7y4aB2BFlYDGEWAACUq//sOS/zN57U16c/3FZqVvYyekqwY4RZAABQrnWyL327W18f3qOh3NY00Ogpwc4RZgEAQLnXyYbXrSIv3NnY6CnBARBmAQBAudfJ/mtQO3FzIWbA8vgpAwAAZfblppPUycIQtOYCAACllpmVLf9YcSAnyD7TkzpZVCzCLAAAKJWEqxkyauF2+f1IvL49tk8TGdGzodHTgoMhzAIAgBI7GZ8sj3++RY7HJYuXm4u8P7CN/F/LYKOnBQdEmAUAACWy8Vi8PPPVdr0yG+zvKXMiwqVlLX+jpwUHRZgFAADFtuCvUzLxx32SmW2SNiGVZc5jYRJYydPoacGBEWYBAECJT/S6p01Nee/B1uLp5mL01ODgCLMAAKBIqpxg9Nc75L+H4/TtF+5oLKNuDxUnJyejpwaUX5/Z1NRUmTZtWnk9HQAAsJITve7/6A8dZNWJXh8/2l5G92pEkIVthtm4uDhZvny5rFmzRrKysvR9GRkZ8sEHH0i9evVk6tSplponAACoYH8dvyADPvpDjsUl6xO9lgzvIne1omMBbLTMYMOGDXL33XfLlStX9F9j4eHhMm/ePBkwYIC4urrKpEmTJDIy0rKzBQAAFSIhJUOe+nKbLjFoG1JZPuVEL9j6yuxrr70mffv2ld27d8uYMWNky5Ytct9998nbb78t+/fvl+HDh4uXF1vXAQBgDz5af1QH2cZBvrLoqc4EWVgtJ5PJZCrOwGrVqsnvv/8uzZs3l6tXr4qvr68sXbpU7r33XrElamXZ399fEhISpFKlSkZPBwAAq3P28lW5bdp6Sc/MlnlDOrA9Law6rxV7ZfbSpUsSEBCgr6sVWG9vb2nZsmXZZwsAAKzKP9cc0kG2c4Oq0rNJdaOnA5Rfay5VThAdHa2vqwXdQ4cOSXJycp4xrVu3LslTAgAAK7L/3BX5fsdZfX38Xc3oWgD7CrO9evXSIdZMnRCmqB90db/6aO5yAAAAbM/UVQdF/a/+7tbBeocvwG7C7IkTJyw7EwAAYKgNR+J1P1k3FycZ26eJ0dMByjfM1q1bt7hDAQCAjcnONsmU/xzQ1x/tVFfqVvMxekpA+YbZ06dPF2tcnTp1ivuUAADASizbdU72nbsifh6uMvr2UKOnA5R/mFU7fBVUBG6ulVXUx8zMzOJ/dQAAYLi0zCx5b/UhfX14z4ZSzdfD6CkBxVbs1lw7duyQ7du3F3gZO3aseHh4SNWqVaU0Zs2apcOyp6endOrUSTZv3lzoWNXbVu0+VrlyZfHx8ZG2bdvKl19+WaqvCwAARGasPaJ7ywZV8pDHu9U3ejqAZVZm27Rpc8N9a9eulXHjxsnhw4flpZdekhdeeKFkX11EFi9erHcUmz17tg6yM2bMkD59+ui2X4GBNzZpVoH51VdflaZNm4q7u7ssX75chg4dqseqxwEAgOL7YO0R+Xj9MX193F1NxcvdxegpAZbZASw3tRr78ssv6x3BnnzySZkwYUKBwbM4VIDt0KGDzJw5U9/Ozs6WkJAQGT16tA7KxdG+fXvp16+fvPnmmzcdyw5gAABcD7Lvrz2cE2SH92ho9JQAy+0Aphw7dkwGDhwoHTt2lOrVq+tNFFQILW2QTU9Pl23btknv3r2vT8jZWd/etGnTTR+vcvi6dev0Ku6tt95a4Ji0tDT9guS+AADg6AiysBfFDrMjRoyQ5s2b64S8detWWbhwoTRo0KBMXzw+Pl5vshAUFJTnfnXbvNNYQdQcfH19dZmBWpH98MMP5Y477ihw7JQpU3SyN1/Uqi8AAI6MIAuHrJlVNa3qBK3Y2Fh5/PHHiyxBsDQ/Pz/ZuXOnJCUl6ZVZVXOrgnXPnj1vGDt+/Hj9eTO1MkugBQA4KoIsHDbMqrrY8t6fOSAgQFxcXCQmJibP/ep2jRo1Cn2cKkUIDb3WA091Mzhw4IBegS0ozKouC+oCAICjW7I1iiALxw2zkyZNKvcvrsoEwsLC9OrqgAEDck4AU7dHjRpV7OdRj1G1sQAAoGCHYxLl9R/36uvP9mpEkIXjhdlffvlFn2Tl6lrshxSLKgGIjIzUvWPViWWqNVdycrJut6VERERIrVq19Mqroj6qsQ0bNtQBduXKlbrP7Mcff1yu8wIAwF6kpGfKiAXbJTUjW7o3CpDnejUyekpAuSl2MlUnWJ0/fz6nc0Hnzp3lu+++00GzLFR3hLi4OF3GoE76UmUDq1atyjkpTG2jq8oKzFTQVSejnTlzRry8vHS/2a+++ko/DwAAuNHrP+yTo7FJEujnIe8PbCvOzuVbNgjYRJ9ZFShV2DSHWXUS1q5du8rc0aCi0WcWAOBodbJjv90tKr8uHNZZOjeoZvSUAOP6zAIAANuskx1zR2OCLOxSscOs6mSQu5tB/tsAAMB662RH9LzWBQhw2JpZVY3Qq1evnBPAUlJS5J577tEdCSq6zywAAChcZla2jPtuD3WycAjFDrMTJ07Mc/vee++1xHwAAEAZJFzNkFELt8vvR+J1ney/BrWTAF/6rcN+lTrMAgAA63IyPlke/3yLHI9LFi83F3l/YBvqZGH3yrdpLAAAMMTGY/HyzFfb9cpssL+nzIkIl5a1/I2eFmBxhFkAAGzcgr9OycQf90lmtknahlSWTx8Lk8BKnkZPC6gQhFkAAGzY+z8flg/WHdHX721bU955oLV4urkYPS2gwhBmAQCwUd9sjcoJsi/e2VhG3hZK20w4HMIsAAA2aOPReHll6R59ffTtoTLq9kZGTwmw3jD7r3/9q9hP+Oyzz5ZlPgAA4CZU/9jhX23TNbL3tKmpd/cCHJWTSe2GcBP169cv3pM5Ocnx48fFXvb6BQDA2lxISpMBH/0hURevSljdKrLgyU7UyMLulCSvFWtl9sSJE+U1NwAAUEqpGVky7IutOsjWqeqtuxYQZOHonI2eAAAAuLnsbJO8uGSXbD99WSp5uspnQzpINXb2Akp3AtiZM2dk2bJlcvr0aUlPT8/zuenTp5fX3AAAgIhcTc/SQXbFnvPi6uwksx8Lk9BAX6OnBdhmmF23bp30799fGjRoIAcPHpSWLVvKyZMnRZXetm/f3jKzBADAQUUnpOrSgj1nE8TNxUmmPdRGujYMMHpagO2WGYwfP15efPFF2bNnj3h6esp3330nUVFR0qNHD3nooYcsM0sAABzQ7jOXpf/MDTrIVvVxlwVPdpZ729YyelqAbYfZAwcOSEREhL7u6uoqV69eFV9fX3njjTfknXfescQcAQBwOMt3n5OHZm+S2MQ0aRzkKz+O7CYd61c1elqA7YdZHx+fnDrZ4OBgOXbsWM7n4uPjy3d2AAA4GFW298HaIzJq4Q5Jy8yW25pUl++e6SohVb2NnhpgHzWznTt3lg0bNkizZs2kb9++8sILL+iSg6VLl+rPAQCA0gfZd1Ydktm/XVsoevKW+jK+bzNxcWaLWqDcwqzqVpCUlKSvT548WV9fvHixNGrUiE4GAACUU5CdeE9zGdqteJsWAY6sWDuA2RN2AAMAWHuQndy/hUR2rWf0tACbyGslrpndsmWL/PXXXzfcr+7bunVrSZ8OAACHRpAFyqbEYXbkyJG6FVd+Z8+e1Z8DAADFQ5AFDAiz+/fvL3BzhHbt2unPAQCA4pn+82GCLFDRYdbDw0NiYmJuuP/8+fO67ywAALi5b7ZEyYe/HNXXCbJABYbZO++8U+8CpgpyzS5fviyvvPKK3HHHHWWYCgAAjuGPo/Hyyvd79PVnbw8lyAJlUOKl1GnTpsmtt94qdevW1aUFys6dOyUoKEi+/PLLsswFAAC7dyQmUYZ/tU0ys03Sv01Nef6OxkZPCXCsMFurVi3ZvXu3LFiwQHbt2iVeXl4ydOhQGTRokLi5uVlmlgAA2IG4xDQZOn+LJKZmSnjdKvLug63FyYkNEYCyKFWRq9rS9qmnnirTFwYAwJGkZmTJU19ulTOXrkrdat7yaUS4eLq5GD0twDHC7LJly+Suu+7SK6/qelH69+9fXnMDAMAuZGeb5IVvdsmO05fF38tNPhvSQar6uBs9LcBxdgBzdnaW6OhoCQwM1NcLfTInJ8nKyhJrxg5gAICK9u6qg/LR+mPi5uIkXz7RSTo3qGb0lAC7yWvFWpnNzs4u8DoAALh5Cy4VZJWp97cmyAJGtubKyMiQXr16yZEjR8p7HgAA2HULrtG3h8oDYbWNnhLg2GFW1cyqTgYAAKBoR2PztuAaQwsuwDo2TRg8eLDMnTvXMrMBAMAOxCfRgguw2tZcmZmZ8tlnn8natWslLCxMt+nKbfr06eU5PwAAbK4F17AvtkrUxatSp6q3fPJYGC24AGsKs3v37pX27dvr64cPH87zOf7qBAA4soys7DwtuOYN7SDVfD2MnhZg10ocZn/99VfLzAQAABuWkJIhIxZukz+OXtAtuGYPDpOG1X2NnhZg90q1AxgAALjuWFySPPn5VjkRnyze7i7yr7+1ky4NacEFWE2Yvf/++2X+/Pm6aa26XpSlS5eW19wAALB6vx+Jk5ELtsuV1EypVdlL/h0ZLs2C2ZQHsKowq3ZgMNfDqusAAEDki00nZfJP+yUr2yRhdavo0oLqftTIAla3na09YTtbAEBZJaZmyD+WH5DFW6P07fvb1ZK3729F1wLAWrezLUhsbKwcOnRIX2/SpIkEBgaW9qkAALAZG47Ey8vf7Zazl6+KetPypT5NZXiPBnT0AQziWpqkPHLkSFm0aJFkZWXp+1xcXGTgwIEya9YsyhAAAHa7Gvv2yoPy9ebT+nZIVS9594E2nOgF2NoOYMOGDZO//vpLli9fLpcvX9YXdX3r1q3y9NNPW2aWAAAYvBr7fzN+zwmykV3qyqq/30qQBWwxzKrgqnYA69Onj65hUBd1fc6cOfLTTz+VahJqRbdevXri6ekpnTp1ks2bNxc6Vn2d7t27S5UqVfSld+/eRY4HAKAsPttwQgbP/UuXFajV2K+HdZbJ97YUHw+6WwI2GWarVatWYCmBuk+Fy5JavHixjBkzRiZOnCjbt2+XNm3a6HCsanILsn79ehk0aJDevGHTpk0SEhIid955p5w9e7bEXxsAgKKs3hctb67Yr68/2qkOq7GAPXQz+PTTT2XJkiXy5ZdfSo0aNfR90dHREhkZqXvQlrTUQK3EdujQQWbOnKlvZ2dn64A6evRoGTdu3E0fr+p2VYhWj4+IiLjpeLoZAACKY8+ZBHn4k01yNSNLB9l/DGjJSV6APXQz+Pjjj+Xo0aNSp04dfVFOnz4tHh4eEhcXJ5988knOWLXSWpT09HTZtm2bjB8/Puc+Z2dnXTqgVl2LIyUlRTIyMqRq1aoFfj4tLU1fcr84AAAURZUUPP75Fh1kezSuLpP7tyDIAlaqxGF2wIAB5fbF4+Pj9cpqUFBQnvvV7YMHDxbrOV5++WWpWbOmDsAFmTJlikyePLlc5gsAcIyuBU/M3yJxiWnStIafzHyknbi6lLgqD4C1hllV22otpk6dqluEqTpadfJYQdSqr6rJzb0yq8oYAADILzMrW0Yt3CEHoxP1Tl5zh3QQP083o6cFoAhlOhUzKSlJ17jmVpI61ICAAN2jNiYmJs/96ra5Hrcw06ZN02F27dq10rp160LHqfIHdQEAoKjV2KOxSfLVn6flt8Nx4unmLHMjw6VWZS+jpwagvMPsiRMnZNSoUXo1NDU1Ned+dR6Zqicyb6RQHO7u7hIWFibr1q3LKV9Q4VjdVl+jMO+++6689dZbsnr1agkPDy/ptwAAcGDq/1cbjsbLb4fi5EhskhyJSZRzCdf/f6ZKYz/4WztpXbuyofMEYKEwO3jwYP2LQPWaVbWtZS2IVyUAqhOCCqUdO3aUGTNmSHJysgwdOlR/XnUoqFWrlq59Vd555x2ZMGGCLFy4UPemVZ0UFF9fX30BAKAg6v9d/z0SLzPWHpYdpy/f8PlAPw9pHOQngzrWkT4tin53EIANh9ldu3bpDgRNmjQplwmobXBVFwQVUFUwbdu2raxatSrnpDDVKUF1OMjdTUF1QXjwwQdvqOWdNGlSucwJAGDfIVaVEdzXrpa0qlVZGgf5SqNAP/H3pjYWcIg+s7fddpu8+uqrhXYPsHb0mQUAx7H99CV5c/n+PCF2cKe68lSPBhLoV/CJwwDsvM/sv//9bxk+fLjecatly5bi5pb3L9miTsYCAKAipGZkyfSfD8u/fz8u2SZCLGDPShxmVUnAsWPHcmpaFVU3W5oTwAAAsMRq7ItLdsnxuGR9+/52tWRc36aEWMBOlTjMPv7449KuXTv5+uuvy+UEMAAALLEaq/rETrmvlfRunndjHgAOHmZPnToly5Ytk9DQUMvMCACAEjh1IVlW7Dkv32yJkpMXUnJWYyfc01wqe7sbPT0A1hZmb7/9dt3RgDALADA6wK7YfV72nbuScz+rsYDjKXGYveeee+T555+XPXv2SKtWrW44Aax///7lOT8AAHKciE+WV5bukU3HL+Tc5+LsJF0aVJN+rYP1pRLbzwIOpcStuXL3fL3hyWzgBDBacwGA7cnONsm8jSflvdUHJTUjWwfYrg2rSd9WwXqDg6o+lBMA9sSirbnUdrMAAFTkauxL3+6SLScv6du3hAbIlPtbSUhVb6OnBsAKlDjMAgBgxGqsj7uLvNKvmTzSsQ6ddADkKLxmIJ++ffvqpV6zqVOnyuXL1/e2vnDhgjRv3ry4TwcAQKHSM7PlmQXb9O5dKsiq1djVz98qj3aqS5AFULowu3r1aklLS8u5/fbbb8vFixdzbmdmZsqhQ4eK+3QAABQaZEd/vV1W74sRd1dneeu+lvLlEx2ldhXKCgCUocwg/3liJTxvDACAEgfZTx8Lk55NAo2eFgB7WJkFAMCSCLIALBpmVY1S/jol6pYAAOWBIAugQsoMhgwZIh4eHvp2amqqDB8+XHx8fPTt3PW0AAAUF0EWQIWE2cjIyDy3Bw8efMOYiIiIMk0GAOBYLiSlyTNfbZfNJy8SZAFYNszOmzevdF8BAIACHIpOlCc+3yJnLl0VPw9X+XhwmNzSKMDoaQGwMWyaAACocL8cjJFnv94pSWmZUreat8yNDJfQQD+jpwXABhFmAQAVRp1/MXfDCXlr5QFRHR47N6gqHz8aJlV83I2eGgAbRZgFAFSIrGyTvPr9Hlm0JUrfHtQxRCb3b6lrZQGgtAizAIAKCbIvLtkl3+84K85OIq/2ay6Pd6tHi0cAZUaYBQBUWJB1dXaSDwe1k7taBRs9LQB2gvd2AAAWQ5AFYGmEWQCARRBkAVQEwiwAoNwRZAFUFGpmAQDl6mp6lg6yK/acJ8gCsDjCLACg3EQnpMpTX26V3WcSCLIAKgRhFgBQLnafuSzDvtgqMVfSpIq3m96etnODakZPC4CdI8wCAMps+e5zurQgNSNbGgX6ytzIDlKnmrfR0wLgAAizAIAybU/7r3VH5f21h/Xtnk2q69ICP083o6cGwEEQZgEApfbBuiMyY+0Rff2JW+rLK32biYva4gsAKghhFgBQKhuPxuswq0zu30Iiu9YzekoAHBB9ZgEAJRabmCrPLtopJpPIwPAQgiwAwxBmAQAl3hDhuUU7JT4pTZoE+cmk/i2MnhIAB0aYBQCUyIe/HJGNxy6It7uLzHq0vXi5uxg9JQAOjDALAChVnexb97WU0EBfo6cEwMERZgEApaqTva9dbaOnBACEWQDAzaVmZMnfv6ZOFoD1oTUXAOCmK7JPfbFNdkZdpk4WgNUhzAIACrX3bIIM+2KrnE9IFX8vN/l4cHvqZAFYFcIsAKBAq/ael+cX75KrGVnSoLqPzI3sIPUDfIyeFgDkQZgFAORhMplk1q9HZdqaw/p290YBMvOR9nplFgCsDWEWAKCpk7tW74uWH3aclS0nL+n7hnStJ6/1ayauLpwvDMA6EWYBwIGZA+yK3eflz+MXJNt07X5XZyfdsWBw57pGTxEAikSYBQAHdDU9S6atOSTzN57U29OatarlL/1aB0u/VsESUtXb0DkCQHEQZgHAwWw9eVHGfrtbTsQn5wmwfVsGS51qBFgAtsXwIqhZs2ZJvXr1xNPTUzp16iSbN28udOy+ffvkgQce0OOdnJxkxowZFTpXALD11dg3l++Xhz7ZpINsUCUPmTekg/w0+hYZ3qMhQRaATTI0zC5evFjGjBkjEydOlO3bt0ubNm2kT58+EhsbW+D4lJQUadCggUydOlVq1KhR4fMFAFteje37r99l7oYTejvaB8Nqy5rne8htTQONnhoAlImTSfVgMYhaie3QoYPMnDlT387OzpaQkBAZPXq0jBs3rsjHqtXZ5557Tl9K4sqVK+Lv7y8J585JpUqVbhzg4iLi6Xn9dvK1t+EK5Ows4uVVurEpKar/TcFjnZxEvL1LN/bqVfVCFj4PH5/SjU1NFcnKKp+xar5q3kpamkhmZvmMVa+vep2V9HSRjIzyGat+HtTPRUnHqnFqfGE8PERcXUs+Vr0G6rUojLu7iJtbyceqY6aOXWHUODW+pGPVz5j6WSuPseo1UK+Fov5NqH8b5TG2JP/ubex3hLr3sx2x8taK/frkrjpeIm/e01x6NCkkxPI7ouRj+R1xDb8jbPJ3hFhpjrhy6ZL416wpCQkJBee13EwGSUtLM7m4uJi+//77PPdHRESY+vfvf9PH161b1/T+++/fdFxqaqopISEh5xIVFaWOpinh2mG98dK3b94n8PYueJy69OiRd2xAQOFjw8PzfwOFj23ePO9Ydbuwsep5clNfp7Cxan65qfkXNlZ937mp16Wwsfl/jB58sOixSUnXx0ZGFj02Nvb62BEjih574sT1sS++WPTYvXuvj504seixmzdfH/vuu0WP/fXX62Nnzix67PLl18fOm1f02G++uT5WXS9qrHouM/U1ihqr5mim5l7UWPW9m6nXpKix6jU1U691UWPVsTJTx7CosepnwEz9bBQ1Vv1smamfuaLGqp/Z3Ioaa2O/Iy5Wr2mq+/JyfXlu0Q5TRvuwwp+X3xHXL/yOuHbhd4Td/44wWWmOUDlN57WEBNPNGFZmEB8fL1lZWRIUFJTnfnU7Ojq63L7OlClT9Eqs+aJWfgHAUSSnZYqzk+hesdMfbqNbbgGAPTGszODcuXNSq1Yt2bhxo3Tp0iXn/pdeekl+++03+euvv8qlzCAtLU1fcpcZqEBLmQFvIfIWIm8h2utbiEdjEmXEwu0SdfGq+Hi6ynsRXa7XxvI74tp1fkeUfCy/I+zmd4S9lRkY1porICBAXFxcJCYmJs/96nZ5ntzl4eGhLwW+aLlfuMIUZ0xpxub+wSnPsbl/0MtzbO5/mOU5Vh2bgo5PWceqX3zmX35GjVW/gM3/EyjPseoXsPl/WuU5Vv0CLu7PcEnGql/AlhirfgFbYqxiDWNL8TtCrU2oDRDGLtktiWkmCalRReZGdpDGQX7Xx/I74hp+R5R8LL8jrrOGsfaeI3yK/1oYVmbg7u4uYWFhsm7dupz71Alg6nbulVoAQNFUiP31UKwM+GijDP9quySmZUrHelXlx5G35A2yAGCHDN00QbXlioyMlPDwcOnYsaPuG5ucnCxDhw7Vn4+IiNClCKruVUlPT5f9+/fnXD979qzs3LlTfH19JTQ01MhvBQAMCbHrD8fJjLVHZFfUZX2fp5uzDOlaX8bc0VjcXQ1vJQ4A9h1mBw4cKHFxcTJhwgR90lfbtm1l1apVOSeFnT59WpzNNUv/q7Nt165dzu1p06bpS48ePWT9+vWGfA8AUNESrmbIz/tj5Ms/T+UJsRFd6smw7g2kul8x32oHADtgaJ9ZI+T0mS1O3zIAsLIAu2L3OdlwNF4ysq796ibEAnD0vGboyiwAoOQBVmkS5Cd9WwXLI53qEGIBODTCLADYWIDt17qGhAZyYhcAKIRZALACcYlp8o8V+2XlnvMEWAAoAcIsABhInbawbNc5mbhsn1xOudZsnwALAMVHmAUAA1djX/thj6zed23zmObBlWTqA62kde3KRk8NAGwGYRYADFiN/Wn3eZn44165lJIhrs5OMur2UBl5W6i4udAbFgBKgjALABUoNSNLXlm6R5buOJuzGjvtoTbSvCatAgGgNAizAFBBYhNT5ekvt8mO05fFxdlJRt8eKiN6hrJTFwCUAWEWACrAvnMJMuzzrXIuIVX8vdzko0fbS7fQAKOnBQA2jzALABa2el+0PLdop1zNyJIGAT4yd0gHqR/gY/S0AMAuEGYBwEKys03y8W/H5L3Vh/Tt7o0CZOag9uLv7Wb01ADAbhBmAcACIfY/e6Plg3WH5XBMkr4vsktdef3u5uJKtwIAKFeEWQCwYIj183SV8Xc1k0c61TF6egBglwizAGChEPvELfVlaLf6+oQvAIBlEGYBoJQIsQBgPMIsAJRTiH28W315/BZCLABUJMIsABTTxeR03WZr3h8nWIkFACtBmAWAYgTYlXvOy8ZjFyQr26TvZyUWAKwDYRYACgiwa/ZFy4p8AVZpUbOS3N26pu5OQIgFAOMRZgGgGAG2b6tg6dcqWOqxcxcAWBXCLACHRYAFANtHmAXgUAiwAGBfCLMA7N6l/53EVVCAbR5cSfq1DtYhtj4BFgBsDmEWgMMGWFZgAcD2EWYBWP0GBWv2x8jMX4/I3rNXSv08rMACgH0izAKw6hD7wbojcuB86UIsK7AAYP8IswCsPsT6erjKkK715G8dQ8TTzaVYz+Pm7Cz+3vSBBQB7R5gFYBUBduupS3qXrf/sPS8xV9LyhFi1XWwVH3ejpwkAsEKEWQCGBdhtpy/Jit15A6x5q9jILoRYAMDNEWYBWE2AvaN5kNzdOli6hQaIh2vxygkAAI6NMAvAuADrcS3AqpO0bmlEgAUAlBxhFkDFB9gWQbrDAAEWAFBWhFkA5eZEfLJ8vvEkK7AAgApDmAVQZmp3rXl/nJD3Vh+StMxsfR8rsACAikCYBVAmx+OSZOy3u2XbqUv6dteG1XQXAgIsAKAiEGYBlMtqrI+7i7zar7kM6hgiTk5ORk8PAOAgCLMA8jCZTBKbmCaHYxLlSEySHIlNkiupGTeMOxmfLPvOXduh65bQAJn6QCupXcXbgBkDABwZYRawQ6kZWfLb4ThZsy9GLiRfPxHrZhJTM+VITKJcSc0s1nhWYwEARiPMAnYWYNWWsGv3x0hyelapn8vZSaReNR9pFOQrjQL9pJqvu+SPqi4uztKraaDUrOxV5rkDAFBahFnABiWmZui3/4/GJOlygMOxSbLt5MU8Abamv6fc1SpYmtTwuyGIFsbTzUVCA32lfoCPvg4AgLUjzAJWICMrW9egHtY1qtdqVY/GJsnVjBtXV9Mys/L0cM3NHGD7tgqWdiGVxVktsQIAYMcIs4ABoVWtqppPsFIf1WYDmdmmEj1XoJ9HThlA4yA/aVGzkrSq5U+ABQA4FMIsUEGhVa24qtCakWUq9GSq0CA/aRzomxNSK3m53TDO1dlJ6lbzlsre7hXwnQAAYN0Is0Ax21XpM/zz5VCTmCQuMa1soVVdD/LTJQJ0BAAAoGQIs7B7qsY0Pim92OMz1arqhRTdokoH1P+daJWYVrx2VWaEVgAALI8w6+DUWfHRCan5FxxtVnpmthxXb++b61FjE+XUhRS9W5WlFBRaGwX6Sq3KXoRWAAAcIczOmjVL3nvvPYmOjpY2bdrIhx9+KB07dix0/JIlS+T111+XkydPSqNGjeSdd96Rvn37VuicbTG0qrPjzSccqdVGFfjOJaSKI3BzcSp2sFTnT6mdrFQgvbaaeq1+VdWpuhRwcpWqYSW0AgDgoGF28eLFMmbMGJk9e7Z06tRJZsyYIX369JFDhw5JYGDgDeM3btwogwYNkilTpsjdd98tCxculAEDBsj27dulZcuWYssSrl4LnCnpJXs7Oz+TSfRqa3FDayVPV3F1cRZ74OzkJHWqeum3880rpOp6UCUPAicAAHbIyaTObDGQCrAdOnSQmTNn6tvZ2dkSEhIio0ePlnHjxt0wfuDAgZKcnCzLly/Pua9z587Stm1bHYhv5sqVK+Lv7y8JCQlSqVIlsbR95xIk6mJKgZ+7lJKR56ShwnqHlpfqfh45q4zq7XAd+AJ9OSseAABYlZLkNUNXZtPT02Xbtm0yfvz4nPucnZ2ld+/esmnTpgIfo+5XK7m5qZXcH374ocDxaWlp+pL7xalIizZHyZd/nir2+BqVPKWy943tmEpKbT9KaAUAAPbO0DAbHx8vWVlZEhQUlOd+dfvgwYMFPkbV1RY0Xt1fEFWOMHnyZDFKnareEl63SoGf8/Zw/d/b4NdqM9U2opU8ixdkzQvq6q3z/Nd5Ox0AADgKw2tmLU2t+uZeyVUrs6qMoaIMu7WBvpRWYaG1oI/5rwMAANg7Q8NsQECAuLi4SExMTJ771e0aNWoU+Bh1f0nGe3h46Iu1K2pFldAKAABQMENPYXd3d5ewsDBZt25dzn3qBDB1u0uXLgU+Rt2fe7zy888/Fzre2pT0fDtCKwAAgBWXGagSgMjISAkPD9e9ZVVrLtWtYOjQofrzERERUqtWLV37qvz973+XHj16yD//+U/p16+fLFq0SLZu3SqffvqpWONKa0FlAgUhtAIAANhgmFWttuLi4mTChAn6JC7VYmvVqlU5J3mdPn1adzgw69q1q+4t+9prr8krr7yiN01QnQyM7jFbWGilNAAAAMCO+8xWtNL2maVLAAAAgPXlNfvY9gkAAAAOiTBbTKzKAgAAWB/CLAAAAGwWYRYAAAA2izALAAAAm0WYBQAAgM0yvM9sRTN3IlMtHwAAAGB9zDmtOB1kHS7MJiYm6o8hISFGTwUAAAA3yW2q32xRHG7ThOzsbDl37pz4+fmVe7st9VeECslRUVEl2pAB1odjaT84lvaDY2k/OJb244qFjqWKpyrI1qxZM89OsAVxuJVZ9YLUrl3bol9DHUz+cdoHjqX94FjaD46l/eBY2o9KFjiWN1uRNeMEMAAAANgswiwAAABsFmG2HHl4eMjEiRP1R9g2jqX94FjaD46l/eBY2g8PKziWDncCGAAAAOwHK7MAAACwWYRZAAAA2CzCLAAAAGwWYRYAAAA2izBbQrNmzZJ69eqJp6endOrUSTZv3lzk+CVLlkjTpk31+FatWsnKlSsrbK4ov2M5Z84c6d69u1SpUkVfevfufdNjD+v9d2m2aNEivRPggAEDLD5HWOZYXr58WUaOHCnBwcH6bOrGjRvze9ZGj+WMGTOkSZMm4uXlpXeUev755yU1NbXC5osb/fe//5V77rlH78Klflf+8MMPcjPr16+X9u3b63+PoaGhMn/+fLE41c0AxbNo0SKTu7u76bPPPjPt27fPNGzYMFPlypVNMTExBY7/448/TC4uLqZ3333XtH//ftNrr71mcnNzM+3Zs6fC546yHctHHnnENGvWLNOOHTtMBw4cMA0ZMsTk7+9vOnPmTIXPHWU7lmYnTpww1apVy9S9e3fTvffeW2HzRfkdy7S0NFN4eLipb9++pg0bNuhjun79etPOnTsrfO4o27FcsGCBycPDQ39Ux3H16tWm4OBg0/PPP1/hc8d1K1euNL366qumpUuXqs5Xpu+//95UlOPHj5u8vb1NY8aM0bnnww8/1Dlo1apVJksizJZAx44dTSNHjsy5nZWVZapZs6ZpypQpBY5/+OGHTf369ctzX6dOnUxPP/20xeeK8j2W+WVmZpr8/PxMn3/+uQVnCUsdS3X8unbtavr3v/9tioyMJMza6LH8+OOPTQ0aNDClp6dX4CxhiWOpxt5+++157lOBqFu3bhafK4qnOGH2pZdeMrVo0SLPfQMHDjT16dPHZEmUGRRTenq6bNu2Tb+9bObs7Kxvb9q0qcDHqPtzj1f69OlT6HhY77HMLyUlRTIyMqRq1aoWnCksdSzfeOMNCQwMlCeeeKKCZgpLHMtly5ZJly5ddJlBUFCQtGzZUt5++23JysqqwJmjPI5l165d9WPMpQjHjx/X5SJ9+/atsHmj7IzKPa4WfXY7Eh8fr39Bql+YuanbBw8eLPAx0dHRBY5X98O2jmV+L7/8sq4hyv+PFhWrNMdyw4YNMnfuXNm5c2cFzRKWOpYq8Pzyyy/y6KOP6uBz9OhRGTFihP5DU+1IBNs5lo888oh+3C233KLeMZbMzEwZPny4vPLKKxU0a5SHwnLPlStX5OrVq7oe2hJYmQVKaOrUqfrEoe+//16f2ADbkZiYKI899pg+oS8gIMDo6aCMsrOz9Qr7p59+KmFhYTJw4EB59dVXZfbs2UZPDSWkThpSq+offfSRbN++XZYuXSorVqyQN9980+ipwQawMltM6n98Li4uEhMTk+d+dbtGjRoFPkbdX5LxsN5jaTZt2jQdZteuXSutW7e28ExR3sfy2LFjcvLkSX12bu5ApLi6usqhQ4ekYcOGFTBzlMe/S9XBwM3NTT/OrFmzZnp1SL3V7e7ubvF5o3yO5euvv67/0HzyySf1bdX9Jzk5WZ566in9B4oqU4D1q1FI7qlUqZLFVmUVfjqKSf1SVH/5r1u3Ls//BNVtVbNVEHV/7vHKzz//XOh4WO+xVN599129SrBq1SoJDw+voNmiPI+lapO3Z88eXWJgvvTv319uu+02fV21A4Lt/Lvs1q2bLi0w/0GiHD58WIdcgqxtHUt1HkL+wGr+I+XauUewBV2Myj0WPb3MDluNqNYh8+fP1y0nnnrqKd1qJDo6Wn/+scceM40bNy5Pay5XV1fTtGnTdDuniRMn0prLRo/l1KlTdZuZb7/91nT+/PmcS2JiooHfBUpzLPOjm4HtHsvTp0/rriKjRo0yHTp0yLR8+XJTYGCg6R//+IeB3wVKcyzV/x/Vsfz66691e6c1a9aYGjZsqLsCwTiJiYm6JaW6qMg4ffp0ff3UqVP68+oYqmOZvzXX2LFjde5RLS1pzWWFVM+0OnXq6GCjWo/8+eefOZ/r0aOH/h9jbt98842pcePGerxqV7FixQoDZo2yHsu6devqf8j5L+oXMGzv32VuhFnbPpYbN27ULQ9VcFJtut566y3deg22dSwzMjJMkyZN0gHW09PTFBISYhoxYoTp0qVLBs0eyq+//lrg//vMx059VMcy/2Patm2rj7v6Nzlv3jyTpTmp/1h27RcAAACwDGpmAQAAYLMIswAAALBZhFkAAADYLMIsAAAAbBZhFgAAADaLMAsAAACbRZgFAACAzSLMAgAAwGYRZgHAhk2aNEk8PT3l4YcflszMTKOnAwAVjh3AAMCGJSUlyZYtW+Suu+6SefPmyaBBg4yeEgBUKMIsANiBIUOGSGxsrKxcudLoqQBAhaLMAADsQOfOneXnn3+WuLg4o6cCABWKMAsAdmD+/Pm6ZnbRokVGTwUAKhRlBgBg4zZt2iTdunWTu+++W5ca/Pnnn0ZPCQAqDGEWAGzcwIED5erVqzJ58mRp3769HDlyREJDQ42eFgBUCMoMAMCGRUVFydKlS2XMmDHSrl07adGihSxYsMDoaQFAhSHMAoANmzlzprRu3Vp69uypbw8ePJgwC8ChEGYBwEalpKTInDlz9Kqs2aOPPipHjx6VzZs3Gzo3AKgohFkAsFFffPGFeHt7692/zEJCQvQq7VdffWXo3ACgonACGAAAAGwWK7MAAACwWYRZAAAA2CzCLAAAAGwWYRYAAAA2izALAAAAm0WYBQAAgM0izAIAAMBmEWYBAABgswizAAAAsFmEWQAAANgswiwAAADEVv0/SqTW2SkrmKkAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def compute_empirical_risks(scores, labels, lambdas):\n", " \"\"\"\n", " Compute empirical risk at every lambda.\n", "\n", " Args:\n", " scores: array of shape (n,), model probability scores\n", " labels: array of shape (n,), true labels\n", " lambdas: array of shape (k,), threshold grid\n", "\n", " Returns:\n", " risks: array of shape (k,), empirical risk at each lambda\n", " \"\"\"\n", " # TODO: for each lambda, compute the mean loss over all calibration samples\n", " # Hint: a list comprehension over lambdas works fine here\n", " risks = np.array([np.mean([loss(score, label, lam) for score, label in zip(scores, labels)]) for lam in lambdas])\n", " return risks\n", "\n", "risks = compute_empirical_risks(calib_scores, y_calib, lambdas)\n", "\n", "# Uncomment to plot empirical risk vs lambda\n", "plt.figure(figsize=(8, 4))\n", "plt.plot(lambdas, risks)\n", "#User-defined alpha\n", "plt.axhline(alpha, color='red', linestyle='--', label=f'α = {alpha}')\n", "plt.xlabel('λ'); plt.ylabel('Empirical FNR'); plt.title('Risk vs Lambda')\n", "plt.legend(); plt.show()" ] }, { "cell_type": "markdown", "id": "cell-pvalue-header", "metadata": {}, "source": [ "---\n", "## Part 4 — P-values (Hoeffding Bound)\n", "\n", "For each λ, compute a p-value testing:\n", "```\n", "H₀: R(λ) > α (this lambda is NOT safe)\n", "```\n", "\n", "We use the **Hoeffding bound**. The p-value is:\n", "```\n", "p(λ) = exp( -2n · max(0, α − R̂(λ))² )\n", "```\n", "\n", "Intuition: if the empirical risk R̂ is well below α, the exponent is large and negative → p-value is tiny → strong evidence that λ controls the risk.\n", "\n", "If R̂ ≥ α (already above budget), `max(0, ...)` clamps to 0 → p-value = 1 → no evidence of safety." ] }, { "cell_type": "code", "execution_count": 22, "id": "cell-pvalue", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Min p-value: 0.000006 Max p-value: 1.000000\n" ] } ], "source": [ "def compute_p_values(risks, n, alpha):\n", " \"\"\"\n", " Compute Hoeffding p-values for each lambda.\n", "\n", " Args:\n", " risks: array of shape (k,), empirical risk at each lambda\n", " n: int, calibration set size\n", " alpha: float, risk level\n", "\n", " Returns:\n", " p_values: array of shape (k,)\n", "\n", " Formula: p = exp(-2 * n * max(0, alpha - r_hat)^2)\n", " \"\"\"\n", " # TODO: implement the Hoeffding p-value formula above\n", " # Hint: np.maximum(0, ...) handles the clamp; np.exp() handles the rest\n", " p_values = np.exp(-2 * n * np.maximum(0, alpha - risks) ** 2)\n", " return p_values\n", "\n", "p_values = compute_p_values(risks, n_calib, alpha)\n", "\n", "# Uncomment to inspect\n", "print(f\"Min p-value: {p_values.min():.6f} Max p-value: {p_values.max():.6f}\")" ] }, { "cell_type": "markdown", "id": "cell-bonferroni-header", "metadata": {}, "source": [ "---\n", "## Part 5 — Bonferroni Correction\n", "\n", "You have 100 p-values (one per λ). Testing all of them at level δ individually would inflate the false rejection rate. **Bonferroni** fixes this by dividing δ by the number of tests:\n", "Why does Testing all lambdas at level δ individually inflates the FPR?\n", "\n", "```\n", "Accept λ as certified ⟺ p(λ) < δ / |Λ|\n", "```\n", "\n", "The **certified set** is all λ values that pass this test.\n", "\n", "> **Think:** If the certified set is empty, what does that mean? What would you change?" ] }, { "cell_type": "code", "execution_count": 34, "id": "cell-bonferroni", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Certified lambdas: 33 / 100\n", "0.01\n", "0.059494949494949496\n", "0.10898989898989898\n", "0.15848484848484848\n", "0.207979797979798\n", "0.25747474747474747\n", "0.30696969696969695\n", "0.3564646464646465\n", "0.40595959595959596\n", "0.45545454545454545\n", "0.5049494949494949\n", "0.5544444444444444\n", "0.6039393939393939\n", "0.6534343434343434\n", "0.702929292929293\n", "0.7524242424242424\n", "0.8019191919191919\n", "0.8514141414141414\n", "0.9009090909090909\n", "0.9504040404040404\n" ] } ], "source": [ "def bonferroni_filter(lambdas, p_values, delta):\n", " \"\"\"\n", " Return the subset of lambdas that pass the Bonferroni correction.\n", "\n", " Args:\n", " lambdas: array of shape (k,)\n", " p_values: array of shape (k,)\n", " delta: float, overall failure probability\n", "\n", " Returns:\n", " certified: array of certified lambda values\n", " \"\"\"\n", " # TODO: compute the per-test threshold and filter lambdas\n", " certified = []\n", " for idx, value in enumerate(p_values):\n", " if value < delta /np.size(lambdas):\n", " certified.append(lambdas[idx])\n", " \n", " return certified\n", "\n", " \n", "\n", "\n", "certified = bonferroni_filter(lambdas, p_values, delta)\n", "print(f\"Certified lambdas: {len(certified)} / {len(lambdas)}\")\n", "for idx, lam in enumerate(lambdas): \n", " if idx % 5 == 0: print(lam) \n" ] }, { "cell_type": "markdown", "id": "cell-select-header", "metadata": {}, "source": [ "---\n", "## Part 6 — Select λ̂\n", "\n", "From the certified set, pick one λ̂ to deploy.\n", "\n", "A natural choice: **the smallest certified λ** (lowest threshold = hardest to be a false negative = best recall).\n", "\n", "> **Think:** What does picking the *largest* certified λ give you instead? \n", "> Is there a case where neither extreme is the right choice?" ] }, { "cell_type": "code", "execution_count": 36, "id": "cell-select", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "λ̂ = 0.0100000000\n" ] } ], "source": [ "def select_lambda_hat(certified):\n", " \"\"\"\n", " Select lambda_hat from the certified set.\n", "\n", " Args:\n", " certified: array of certified lambda values (could be empty!)\n", "\n", " Returns:\n", " float: the chosen lambda, or None if certified is empty\n", " \"\"\"\n", " if len(certified) == 0:\n", " print(\"No certified lambda found. Try increasing n_calib or relaxing alpha/delta.\")\n", " return None\n", " \n", " # TODO: return the best lambda from the certified set\n", " return np.min(certified)\n", "\n", "\n", "lambda_hat = select_lambda_hat(certified)\n", "\n", "print(f\"λ̂ = {lambda_hat:.10f}\")\n", "\n", " " ] }, { "cell_type": "markdown", "id": "cell-validate-header", "metadata": {}, "source": [ "---\n", "## Part 7 — Validate on the Test Set\n", "\n", "Now measure the **actual** FNR on the held-out test set using λ̂.\n", "\n", "LTT guarantees that this will exceed α only with probability ≤ δ. \n", "On a single run you can't verify the probabilistic guarantee directly — but you can check whether the empirical test FNR is reasonable.\n", "\n", "To truly verify the guarantee, you'd run many random calibration/test splits (see the commented loop below)." ] }, { "cell_type": "code", "execution_count": 49, "id": "cell-validate", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "λ̂ = 0.0100\n", "Test FNR = 0.0000\n", "Risk budget α = 0.1\n", "Within budget? = True\n" ] } ], "source": [ "if lambda_hat is not None:\n", " test_losses = [loss(s, l, lambda_hat) for s, l in zip(test_scores, y_test)]\n", " test_fnr = np.mean(test_losses)\n", " print(f\"λ̂ = {lambda_hat:.4f}\")\n", " print(f\"Test FNR = {test_fnr:.4f}\")\n", " print(f\"Risk budget α = {alpha}\")\n", " print(f\"Within budget? = {test_fnr <= alpha}\")" ] }, { "cell_type": "markdown", "id": "cell-trials-header", "metadata": {}, "source": [ "---\n", "## Part 8 — Verify the Guarantee (Many Trials)\n", "\n", "The real guarantee is: **across many random splits, the test FNR exceeds α in at most δ fraction of trials.**\n", "\n", "Run 200 independent trials with different calibration/test splits. Count how often `test_FNR > α`. That fraction should be ≤ δ = 0.1." ] }, { "cell_type": "code", "execution_count": 50, "id": "cell-trials", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Violations (FNR > α): 0 / 200 = 0.000\n", "Expected: ≤ δ = 0.1\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAGJCAYAAACQBRs3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA96klEQVR4nO3dCXwTZfrA8adQWu6b0lY5lVNOOSp4oSA3iOCFqEUR0AVUWBXrIpfr1hPwQFgVAf+CCCwii4jLjcohoIiIoGULqJyKUFqgQJv/53m7CUnb9CJNMtPf9/MZMpmZTN5kUvLknWeeN8ThcDgEAAAAsIFigW4AAAAA4CsEtwAAALANglsAAADYBsEtAAAAbIPgFgAAALZBcAsAAADbILgFAACAbRDcAgAAwDYIbgEAAGAbBLeABY0fP15CQkL88lwdOnQwk9PatWvNcy9cuNAvzz9w4ECpXbu2BLPk5GR56KGHJDIy0rw3jz/+eKCbhCLoUv5W9HH6eMAOCG6BAJs1a5YJiJxTyZIlJTo6Wrp06SKvv/66nDp1yifPc/DgQRMUb9++XYJNMLctL/7xj3+Y4/jII4/I//3f/8l9993n9QdJbpP7D4lLsWzZMvOceaXP661Nu3fv9vhho9O2bduy7EODo7Jly+a431KlSkmzZs1kypQpkp6enq/XlJaWZv42dD+fffaZ2I0GmHn5jOhnDYB3oTmsA+BHEydOlDp16sj58+fl8OHDJpDQHsBJkybJkiVLTEDgNGbMGHn66afzHUBOmDDBfIG2aNEiz4/7z3/+I4Utp7a98847+Q6C/G316tVyzTXXyLhx47xu07dvX7nyyis9ens1GL7tttvMOqfq1av7LLidOnVqvgLcyy+/XOLj47Ms14AyM93vv//973zv9/fff5e5c+fKyJEj5dixY/L888/n630+dOiQ+ZzMmTNHunXrJnaiAb9+LtyP4YcffiiTJ0+WqlWrupa3b98+28db4W8F8AeCWyBI6Bd169atXffj4uLMl3nPnj2ld+/e8uOPP5peLxUaGmqmwnT69GkpXbq0hIWFSSCVKFFCgt3Ro0elcePGOW6jP07cf6BokKfBrS679957JRhUqFAhT23RHyBLly6Vb775Rq6++up87/fhhx+Whg0byhtvvGF+1BUvXjxP7fvggw/M88XGxsozzzwjKSkpUqZMGfHl5z2Q+vTp43Fff+RqcKvLc0o3cL4PVvhbAfyBtAQgiN18883y7LPPyv79+80Xe045tytWrJDrrrtOKlasaE4NN2jQwAQASnuB27RpY+YfeOCBLKc39dRxkyZNzKnmG264wXzJOx+bOefW/RSxbqN5pvrFqgH4L7/8kqc8Pvd95ta27PII9cv8r3/9q9SoUUPCw8PNa33llVfE4XB4bKf7GT58uCxevNi8Pt32qquukuXLl+c5aB00aJDpTdV0kebNm8vs2bNd652n6RMTE+XTTz91tX3fvn1SUJoCcPvtt0vlypXNc+oPHu25d6e9+9rTXa9ePbNNlSpVzLHXz4DzPdNeW+d74Jx8ZcSIEVKpUqV89Qq70zbrMdeUG32P8+LMmTPy8ccfy9133y133nmnuf/JJ59ku62mLNx4441Srlw5KV++vHku7S12yunzntsxd5o3b560atXK9RxNmzaV1157Lc/HqKCcqR979+6V7t27m+cfMGCA178V/bvQnl59fv1xrG3OS758YbUf8Ad6boEgp/mb+sWr6QGDBw/OdpsffvjB9PBqL6D2hGkQl5CQIF999ZVZ36hRI7N87NixMmTIELn++uuznN78448/TO+xBg/ay5bb6XE9nawB0+jRo01AoKdUO3XqZPJmnT3MeZGXtrnTAFYD6TVr1pggRHsRP//8c3nyySflt99+M6dw3X355ZeyaNEi+ctf/mICAc1j7tevnxw4cMB8YXujwZMGQfo+aoCsKSMLFiwwAcSJEyfkscceM23XHFs9xa6n3jXgVtWqVZOC0ON47bXXymWXXWbSTvRHw/z5803P3b/+9S+TwqA0qNTT/HoRW9u2bSUpKUm2bt1qelJvueUWGTp0qEn10EBE25dX+oNFe5TdaWCTOY9Wgzl9zXrM8tp7m5n+ANDPj/4YywsN8PWUvX4+9QeVHhtNTbjnnns8ttMfRQ8++KD5EaNnP3T/3377rflB475tdp/3vBxzpe9r//79pWPHjvLiiy+aZXpmRf/enNvkdowuxYULF0xOvgabGrzm1OOsAbf+vWgAfO7cOROU33HHHabnvUePHl4fV5jtBwqdA0BAzZw5U7sbHVu2bPG6TYUKFRwtW7Z03R83bpx5jNPkyZPN/WPHjnndh+5ft9Hny+zGG28066ZPn57tOp2c1qxZY7a97LLLHElJSa7l8+fPN8tfe+0117JatWo5YmNjc91nTm3Tx+t+nBYvXmy2/fvf/+6x3e233+4ICQlxJCQkuJbpdmFhYR7LvvvuO7P8jTfecORkypQpZrsPPvjAtezcuXOOdu3aOcqWLevx2rV9PXr0cOSHHivdvx5Lp44dOzqaNm3qOHv2rGtZenq6o3379o569eq5ljVv3jzX5xs2bJjHZyQ3zs9A5sn9+DmP/YIFCxwnTpxwVKpUydG7d2/Xet22TJkyWfbbsGFD83p12r17t+PJJ580+8nPe9azZ0/Htdde67r/9ttvO0JDQx1Hjx51LdM2lStXzhETE+M4c+aMx+P1fczt857XY/7YY485ypcv77hw4YLX9ublGOXm5ZdfNu1JTEz0eI912dNPP53r34o6ffq0x319PU2aNHHcfPPNHssz/636ov1AoJCWAFiA9pzlVDXB2fulp2kLekGJ9vZqWkBe3X///aYn1ElPpUdFRZmLYAqT7l9zNB999FGP5dprqvFs5qvotTf5iiuucN3X3m3tefzvf/+b6/NoD6H20DlpTqM+r/Ygrlu3Tnzp+PHjJsdaT7nrsdYeVJ20h1F76X7++WfTM+083trLq8t8SU9pa6+k+/TUU095zaPVCx61R1V7RnNLtdDebJ001/bll182vYl5vepf3wPtnXc/Ftr7rj2/2rPtpO3V9057vbXH2V3mtIzsPu95Peb6/mtqTE6n6AvrGDlpvnZeuJ9F+fPPP+XkyZPm7Ij2wOaksNsPFCaCW8AC9IvVPZDM7K677jKns/UUop5e1VOt+qWfn0BXT4Xn5+IxzcXLHDxoNYBLyTfNC80/1qv3M78fmiLgXO+uZs2aWfah+aL6RZ/b8+hrLFasWJ6e51LpqXANzjXH2hkIOidnFQZnfqqmcehp8vr165tcT03J2LFjxyW3QdMg9MeA+5TThXJ6Cl6DoNxyb51Bswaob731lvmsaaWEzAGoNx999JHJAW3ZsqV5n3TSHwMxMTEmNcFJ81CV5tMW5POe12OuKS763mtag6ajaBpE5jzuwjpGSi8m1efNC00/0Eoe+l5rHrd+nqZNm2aC3JwUZvuBwkZwCwS5X3/91XwRuZeRyq53Zv369bJy5UqTo6tfQhrwam6c5lHmRX7yZPPK20VMeW2TL3i7Ej/zxWeB5vwh8sQTT2TpPXVOzs+AXgSlgdx7771nArl3333X5L3qrT/ltffWGTR37tzZ9DhqD+nXX3/tuogrN84AVn/AafDpnDSfeuPGjbn2wvv68x4REWFyy/V1O/O/NdDVKg5OhXmMtNc5cwCenS+++MK0TwNb/VGh77t+jjT3OLfPf7B8xoCCILgFgpzzgiA9NZ0T/bLTC1y0Lu6uXbvMBV96mlu/eJWvRzTLfLpSvyy1R839am3tIdXen8wy93rmp221atUyF0tlTtNwDjSg631B96OvMXPvt6+fx6lu3bqu0+CZe0+dk3tvtfbC6Wl1LRWlVSo03cK9B9VfI9hpcKu9t3plfV45y5/985//NBf25UQrUWzYsMFc4KUXd7lP2qOrva/OSgjO9JOdO3cW+jHX5+3Vq5cJGjUI1Iv43n//ffM3kNdjVNj0IkQNbLXHXHuXNQDXz1FeBbr9QEER3AJBTIPT5557zly17Sz3kx09RZuZczCE1NRUc+usB5pdsFkQ+kXuHmBqeSEtsO9eWF+DjU2bNpmrtN1Pk2YuGZaftmn5I+35ffPNNz2Wa5UEDeh8Vdhfn0frjGoA5X6VutZm1RxoLTXlS9obqFfqa8Cn72NmehrfPQfVnbZHe3Wdx7owjnduvbea752fEeY0l1dTDfTHWF56bXV7zet2nzQ/WY+DcxvtGdYfAHqV/9mzZ/PdU5/XY575/dcfls4axs5jkJdj5I+zFvo34X6mRNOGtDReboKh/UBBUQoMCBJ6IZT2EOmX6ZEjR0xgq6cQtbdIT3/mlJ+o+XGalqClfXR7zc3UHiXNy9NyQc5AU3vYpk+fbgIADX40Z1ED54LQXh3dt/bsaHu1FJh++bmXK9McYA16u3btagIR7eHSer3uF3jlt23aW3bTTTfJ3/72N/NFrXVItUyaBlcaZGXed0FpWTINNLUMlNZD1R5pfS1a7klfa0450AWltWn1PdUcR30ftTdX31s99a7pKd99953ZTvNgNRDWmqV6HLREk7ZNezeddJ3Si6G0118DHc3FLgyae6s/LrR9eR1UQV+DBpN6mlvzjL2VZdPAVX+oaU3j7Ohpd6276yxJpu3Qz53WttXT73r2QNulgzRkV6+2IMdc968/KLUOtf6N6ZkIDYC1nc783Lwco8Km/x/ojwf9+9P3Qv9f0M+Y/p3mlj8bDO0HCixgdRoAeJQCc05auioyMtJxyy23mLJa7iWnvJUCW7VqlePWW291REdHm8frbf/+/R0//fSTx+M++eQTR+PGjU0JJffSW1oa6aqrrsq2fd5KgX344YeOuLg4R0REhKNUqVKmbND+/fuzPP7VV181ZcPCw8NNKaetW7dm2WdObcuuvNGpU6ccI0eONK+zRIkSpkyWlk1yL/ekdD9aEiszbyXKMjty5IjjgQcecFStWtW8r1qmK7tyZb4qBab27t3ruP/++81nQF+bvndaBmvhwoWubbQMWtu2bR0VK1Y0772W2nr++edNmScnLVM1YsQIR7Vq1UyJtNz+u8/pM5BdKTBvn8nsSoF52+/atWuzfQ+ctm3bZtY/++yzXtu0b98+s41+HpyWLFliyqfpe6Mlu/S90s9rXtqUl2Oux6Jz587ms6/b1KxZ0zF06FDHoUOH8nWMCloKLPN77L4u89/KjBkzzN+H/v1pG/S1ZP7/I7u/CV+0HwiUEP2n4KExAAAAEDzIuQUAAIBtENwCAADANghuAQAAYBsEtwAAALANglsAAADYBsEtAAAAbINBHP43prsO56kFuv01ZCUAAADyTqvX6siY0dHRZmRAbwhuRUxg6230GwAAAAQPHcJdRwf0huBWxDWkor5Z5cuXD3RzAAAAkElSUpLpjMxt+HOCWxFXKoIGtgS3ABDkTp8WadMmY37LFpHSpQPdIgB+lFsKKcEtAMBadNT4XbsuzgOAG6olAAAAwDYIbgEAAGAbpCUAAIAiXV7qwoULkpaWFuimFHnFixeX0NDQSy7LSnALAACKpHPnzsmhQ4fktF6kiKBQunRpiYqKkrCwsALvg+AWAAAUyQGcEhMTTW+hDgqgwRQDOQW2B11/bBw7dswcl3r16uU4UENOCG4BANaiAUitWhfngQLQQEoDXK2bqr2FCLxSpUpJiRIlZP/+/eb4lCxZ0noXlMXHx0ubNm1MMd6IiAjp06eP7Nmzx2Obs2fPyrBhw6RKlSpStmxZ6devnxw5csRjmwMHDkiPHj3Mh1P38+STT5r8GQCADWkgsm9fxkRQgktU0N5BBO/xCOgRXbdunQlcN23aJCtWrJDz589L586dJSUlxbXNyJEj5d///rcsWLDAbK9D5fbt29e1XhPANbDVCH/Dhg0ye/ZsmTVrlowdOzZArwoAAACBEuLQJIcgoXkW2vOqQewNN9wgJ0+elGrVqsncuXPl9ttvN9vs3r1bGjVqJBs3bpRrrrlGPvvsM+nZs6cJeqtXr262mT59uowePdrsLy8JyTqcW4UKFczzMUIZAAD2p2eGNbezTp06BT79Df8el7zGa0GVc6uNVZUrVza327ZtM725nTp1cm3TsGFDqVmzpiu41dumTZu6AlvVpUsXeeSRR+SHH36Qli1bZnme1NRUM7m/WQBgN5qy9fvvv1/yfqpWrWr+3w0aZ86I3HBDxvz69ZqoF+gWAQgiQRPcalL3448/Ltdee600adLELDt8+LDpea1YsaLHthrI6jrnNu6BrXO9c523XN8JEyYU0isBgOAIbBs0bCRnz1x6iaOSpUrLnt0/Bk+Am54usnXrxXkAfvt/5ZFHHpE1a9aY66BiY2NNTKW1ab15/vnn5dNPP5Xt27ebmO7EiRNFJ7jV3NudO3fKl19+WejPFRcXJ6NGjfLoudWrJQHALrTHVgPbKj3/KiWqFPz/t/N//CJ/LH3V7C9oglsAfpf2v2ucIiMjzTVOWh/4/vvvN9UN/vGPf3h9nF4Tdccdd0i7du1kxowZfmlrUAS3w4cPl6VLl8r69evl8ssvdy3XN1DfFI3y3XtvtVqCrnNu8/XXX3vsz1lNwblNZuHh4WYCALvTwDY88spANwOwFrcL27MoXlzEPRc0p231yn/3tBlv25Ypk6/mJScny4gRI2ThwoWmUtQTTzwh99xzj6kNe/ToUdOr6mv/+c9/ZNeuXbJy5UpzhrxFixby3HPPmWucxo8f7/UaJ+eZcr3Y318CWi1Br2XTwPbjjz+W1atXm+Rhd61atTK/CFatWuVapqXCtFtcfwEovf3+++/NwXTSyguaaNy4cWM/vhoAAGALGhx6m/r189w2IsL7tt26eW5bu3b22+XTwIEDTe/p2rVrZebMmfLss8/KM888Y65RKpvD/nRdTtPDDz/s9bHernHSs996jVMwCQ10KoJWQvjkk09MrVtnjqxeCaeFfPV20KBBJoVALzLTgFV/qWhAqxeTKS0dpkHsfffdJy+99JLZx5gxY8y+6Z0FAAB2oilCixYtkjlz5phOQHXbbbfJ+++/n+tp/+3bt+e4PqcKBAW5xqlIBrfTpk0ztx06dPBYrr9C9FeJmjx5sinoq4M3aIUD/ZXw1ltvubbVYfM0pUETnDXoLVOmjElwnjhxop9fDQAAsIXk5JzTEty5nTnOIvOABDrwyCVKSEgwZ76dZ7BV27ZtzXgAvXv3zvGxV15ZNFKUAhrc5qXErtY4mzp1qpm8qVWrlixbtszHrQMABK2qVQPdAthZfnJgC2tbL5xnpd1zXHVMgPr165uyfTnJLRf33nvvNWMFZKcg1zgFSlBcUAYAQL4ChGPHAt0KICD0+iQ9o/3zzz9LdHS0WbZkyRJzPZJ2GoaEhBRKWoL2FGtZL73GSQfcCuZrnAhuAQAALEKrR/Xt29cEmpqO8NNPP8ny5cvNtUp6cX7Hjh0LJS0hL9c4ac+ulgfTQgCXXXaZWaZB9/Hjx82tlhNzBtjalsKo6hDwagkAAADIH03V1LRNDSC1QsKUKVPMNGDAgEKrJeu8xklvtRdXUxg0kHW/xun06dOmqpWOLus0duxYM1rsuHHjTAkznddpq3MglkJAzy0AwFp0+F1niaXPPmP4XRQ5mhagqQiZaXBbmHK7xkkLBGS+nkrr2/qzxq0iuAUAWIsOubtu3cV5AHBDWgIAAABsg+AWAAAAtkFwCwAAANsguAUAAEVWXgaUgrWOB8EtAAAockqUKOEqX4Xg4TwezuNTEFRLAABYT+nSgW4BLE7rteqACDrilipdunSOo3uh8HtsNbDV46HHRY9PQRHcAgCsN/xuSkqgWwEbiIyMNLfOABeBp4Gt87gUFMEtAAAokrSnNioqygyK4D6qFgJDUxEupcfWieAWAAAUaRpQ+SKoQnDggjIAgLWcPSvSo0fGpPMA4IaeWwCAtaSliTjHt9d5AHBDzy0AAABsg+AWAAAAtkFwCwAAANsguAUAAIBtENwCAADANghuAQAAYBuUAgMAWG/4XYcj0K0AEKTouQUAAIBtENwCAADANghuAQDWokPu3nFHxsTwuwCCKbhdv3699OrVS6KjoyUkJEQWL17ssV6XZTe9/PLLrm1q166dZf0LL7wQgFcDAPALHXJ34cKMieF3AQRTcJuSkiLNmzeXqVOnZrv+0KFDHtN7771ngtd+/fp5bDdx4kSP7UaMGOGnVwAAAIBgEtBqCd26dTOTN5GRkR73P/nkE7npppukbt26HsvLlSuXZVsAAAAUPZbJuT1y5Ih8+umnMmjQoCzrNA2hSpUq0rJlS5OycOHChRz3lZqaKklJSR4TAAAArM8ydW5nz55temj79u3rsfzRRx+Vq6++WipXriwbNmyQuLg4k5owadIkr/uKj4+XCRMm+KHVAAAA8CfLBLeabztgwAApWbKkx/JRo0a55ps1ayZhYWEydOhQE8CGh4dnuy8NgN0fpz23NWrUKMTWAwAAwB8sEdx+8cUXsmfPHvnoo49y3TYmJsakJezbt08aNGiQ7TYa9HoLfAEAAGBdlghuZ8yYIa1atTKVFXKzfft2KVasmERERPilbQAAPytdWiQ5+eI8AARLcJucnCwJCQmu+4mJiSY41fzZmjVrulIGFixYIK+++mqWx2/cuFE2b95sKihoPq7eHzlypNx7771SqVIlv74WAICfhISIlCkT6FYACFIBDW63bt1qAlMnZx5sbGyszJo1y8zPmzdPHA6H9O/fP8vjNbVA148fP95UQKhTp44Jbt3zaQEAAFB0BDS47dChgwlcczJkyBAzZUerJGzatKmQWgcACEqpqSJDh2bM//Of2tMR6BYBCCKWqXMLAIChtcxnz86YcqlrDqDoIbgFAACAbRDcAgAAwDYIbgEAAGAbBLcAAACwDYJbAAAA2AbBLQAAAGzDEsPvAgDgokPuHj16cR4A3BDcAgCsN/xutWqBbgWAIEVaAgAAAGyD4BYAYL3hd4cNy5h0HgDcENwCAKxFh9x9662MieF3AWRCcAsAAADbILgFAACAbRDcAgAAwDYIbgEAAGAbBLcAAACwDYJbAAAA2AYjlAEArKVUKZHExIvzAOCG4BYAYC3FionUrh3oVgAIUqQlAAAAwDYIbgEA1nLunMiTT2ZMOg8AbghuAQDWcv68yCuvZEw6DwBuCG4BAABgGwS3AAAAsA2CWwAAANhGQIPb9evXS69evSQ6OlpCQkJk8eLFHusHDhxolrtPXbt29djm+PHjMmDAAClfvrxUrFhRBg0aJMnJyX5+JQAAAJCiHtympKRI8+bNZerUqV630WD20KFDrunDDz/0WK+B7Q8//CArVqyQpUuXmoB5yJAhfmg9AAAAgk1AB3Ho1q2bmXISHh4ukZGR2a778ccfZfny5bJlyxZp3bq1WfbGG29I9+7d5ZVXXjE9wgAAACg6gj7ndu3atRIRESENGjSQRx55RP744w/Xuo0bN5pUBGdgqzp16iTFihWTzZs3e91namqqJCUleUwAAIvQIXd37syYGH4XgJWCW01JeP/992XVqlXy4osvyrp160xPb1pamll/+PBhE/i6Cw0NlcqVK5t13sTHx0uFChVcU40aNQr9tQAAfDj87lVXZUw6DwDBkpaQm7vvvts137RpU2nWrJlcccUVpje3Y8eOBd5vXFycjBo1ynVfe24JcAEAAKzPUj9569atK1WrVpWEhARzX3Nxjx496rHNhQsXTAUFb3m6zjxera7gPgEALEKH3B0/PmNi+F0AVg5uf/31V5NzGxUVZe63a9dOTpw4Idu2bXNts3r1aklPT5eYmJgAthQAUGh0yN0JEzImht8FEExpCVqP1tkLqxITE2X79u0mZ1anCRMmSL9+/Uwv7N69e+Wpp56SK6+8Urp06WK2b9SokcnLHTx4sEyfPl3Onz8vw4cPN+kMVEoAAAAoegLac7t161Zp2bKlmZTmwer82LFjpXjx4rJjxw7p3bu31K9f3wzO0KpVK/niiy9MWoHTnDlzpGHDhiYHV0uAXXfddfL2228H8FUBAACgSPbcdujQQRwOh9f1n3/+ea770B7euXPn+rhlAAAAsCJL5dwCAAAAOSG4BQAAgG0Q3AIAAMA2gnoQBwAAsihZUuTrry/OA4AbglsAgLUULy7Spk2gWwEgSJGWAAAAANug5xYAYC065O5rr2XMP/aYSFhYoFsEIIgQ3AIArEWH3H3qqYz5v/yF4BaAB9ISAAAAYBsEtwAAALANglsAAADYBsEtAAAAbIPgFgAAALZBcAsAAADboBQYAMBadMjdNWsuzgOAG4JbAID1ht/t0CHQrQAQpEhLAAAAgG3QcwsAsN4IZW+/nTE/ZIhIiRKBbhGAIEJwCwCwlnPnRIYPz5gfOJDgFoAH0hIAAABgGwS3AAAAsA2CWwAAANgGwS0AAABsg+AWAAAAtkFwCwAAANsIaHC7fv166dWrl0RHR0tISIgsXrzYte78+fMyevRoadq0qZQpU8Zsc//998vBgwc99lG7dm3zWPfphRdeCMCrAQD4RXi4yNKlGZPOA0CwBLcpKSnSvHlzmTp1apZ1p0+flm+++UaeffZZc7to0SLZs2eP9O7dO8u2EydOlEOHDrmmESNG+OkVAAD8LjRUpEePjEnnAcBNQP9X6Natm5myU6FCBVmxYoXHsjfffFPatm0rBw4ckJo1a7qWlytXTiIjIwu9vQAAAAhulsq5PXnypEk7qFixosdyTUOoUqWKtGzZUl5++WW5cOFCjvtJTU2VpKQkjwkAYKHhd2fNyph0HgDcWOZ8ztmzZ00Obv/+/aV8+fKu5Y8++qhcffXVUrlyZdmwYYPExcWZ1IRJkyZ53Vd8fLxMmDDBTy0HAPh8+N0HHsiYv+MOht8FYL3gVi8uu/POO8XhcMi0adM81o0aNco136xZMwkLC5OhQ4eaADbcy4UGGgC7P057bmvUqFGIrwAAAAD+EGqVwHb//v2yevVqj17b7MTExJi0hH379kmDBg2y3UaDXm+BLwAAAKwr1AqB7c8//yxr1qwxebW52b59uxQrVkwiIiL80kYAAAAEj4AGt8nJyZKQkOC6n5iYaIJTzZ+NioqS22+/3ZQBW7p0qaSlpcnhw4fNdrpe0w82btwomzdvlptuuslUTND7I0eOlHvvvVcqVaoUwFcGAACAIhfcbt261QSmTs482NjYWBk/frwsWbLE3G/RooXH47QXt0OHDia1YN68eWZbrYBQp04dE9y659MCAACg6AhocKsBql4k5k1O65RWSdi0aVMhtAwAAABWFNQ5twAAZKEXBM+ff3EeANwQ3AIArEWH3NX6tgBg9RHKAAAAgJzQcwsAsBYdYv3jjzPmb7stoycXAP6H/xEAANaSmipy550Z88nJBLcALj0toW7duvLHH39kWX7ixAmzDgAAALBMcKtD2+qgCplprdnffvvNF+0CAAAA8i1f53Kcgyqozz//XCpUqOC6r8HuqlWrpHbt2vlvBQAAAODv4LZPnz7mNiQkxIwi5q5EiRImsH311Vd90S4AAACgcIPb9PR0c6vD3G7ZskWqVq2a/2cEAAAACkmBLjFNTEz0fUsAAACAS1Tg+imaX6vT0aNHXT26Tu+9996ltgsAgOyFhYnMnHlxHgAuNbidMGGCTJw4UVq3bi1RUVEmBxcAAL8oUUJk4MBAtwKAnYLb6dOny6xZs+S+++7zfYsAAAAAfwa3586dk/bt2xf0OQEAuLThdz//PGO+SxdGKANw6YM4PPTQQzJ37tyCPBQAgEsffrdnz4xJ5wHATYF+7p49e1befvttWblypTRr1szUuHU3adKkguwWAAAA8H9wu2PHDmnRooWZ37lzp8c6Li4DAACApYLbNWvW+L4lAAAAQCBybgEAAADb9NzedNNNOaYfrF69+lLaBAAAAPgvuHXm2zqdP39etm/fbvJvY2NjC9YSAAAAIBDB7eTJk7NdPn78eElOTr7UNgEA4J0OufvmmxfnAcCNTytf33vvvdK2bVt55ZVXfLlbAAAu0vKTw4YFuhUAisIFZRs3bpSSJUv6cpcAAABA4fbc9u3b1+O+w+GQQ4cOydatW+XZZ58tyC4BAMibtDSRL77ImL/+epHixQPdIgBW77mtUKGCx1S5cmXp0KGDLFu2TMaNG5fn/axfv1569eol0dHRpvrC4sWLswTNY8eOlaioKClVqpR06tRJfv75Z49tjh8/LgMGDJDy5ctLxYoVZdCgQeT9AoCdnT2rZXsyJp0HgEvtuZ05c6b4QkpKijRv3lwefPDBLL3B6qWXXpLXX39dZs+eLXXq1DG9wl26dJFdu3a50h80sNVe4xUrVpiqDQ888IAMGTJE5s6d65M2AgAAoIhcULZt2zb58ccfzfxVV10lLVu2zNfju3XrZqbsaK/tlClTZMyYMXLrrbeaZe+//75Ur17d9PDefffd5rmXL18uW7ZskdatW5tt3njjDenevbu5qE17hAEAAFB0FCgt4ejRo3LzzTdLmzZt5NFHHzVTq1atpGPHjnLs2DGfNCwxMVEOHz5sUhGcNAUiJibGXLim9FZTEZyBrdLtixUrJps3b/a679TUVElKSvKYAAAAUESD2xEjRsipU6fkhx9+MDmvOukADhokaqDrCxrYKu2pdaf3nev0NiIiwmN9aGioyQF2bpOd+Ph4j5zhGjVq+KTNAAAAsGBwq6kAb731ljRq1Mi1rHHjxjJ16lT57LPPJNjFxcXJyZMnXdMvv/wS6CYBAAAgUMFtenq6lNAi2pnoMl3nC5GRkeb2yJEjHsv1vnOd3mqKhLsLFy6YnmTnNtkJDw831RXcJwAAABTR4FbzbR977DE5ePCga9lvv/0mI0eONHm3vqDVETRAXbVqlWuZpj1oLm27du3Mfb09ceKEubDNafXq1SbA1txcAIANaefKSy9lTNl0tAAo2gpULeHNN9+U3r17S+3atV35qnpqv0mTJvLBBx/keT9ajzYhIcHjIrLt27ebnNmaNWvK448/Ln//+9+lXr16rlJgWgGhT58+ZntNi+jatasMHjxYpk+fbkqBDR8+3FRSoFICANhUWJjIk08GuhUA7BTcakD7zTffyMqVK2X37t2uQNO9skFe6IhmN2kR7v8ZNWqUuY2NjZVZs2bJU089ZWrhat1a7aG97rrrTL6v+xC/c+bMMQGt9hhrlYR+/fqZ2rgAAAAoekIcWlA2j/SUvwaSmzZtypKnqhdmtW/f3vSgXq/DIVqIpjto1QR9DeTfArAD7YDQEo2RsVMkPPLKAu8n9XCCHJ79uEn/uvrqqyVoht/95puMeW0Tw+8CRUJSHuO1fOXc6qAKmgKQ3Q71yYYOHSqTJk0qWIsBAMgLHXK3bduMieF3AVxKcPvdd9+ZHFdvOnfu7HFxFwAAABC0wa2W4cquBJj7AAq+GqEMAAAAKNTg9rLLLjMjkXmzY8cOiYqKyncjAAAAAL8Ht927dzfluM5mk+N05swZGTdunPTs2dMnDQMAAAAKtRTYmDFjZNGiRVK/fn1TNaFBgwZmuZYD06F309LS5G9/+1u+GwEAAAD4PbitXr26bNiwQR555BGJi4sTZxWxkJAQ6dKliwlwdRsAAADAEoM41KpVS5YtWyZ//vmnGV1MA1wdQaxSpUqF00IAANzphc3jxl2cB4BLHaFMaTDbpk2bgj4cAICCD787fnygWwHADheUAQAAALbsuQUAICDS00V+/DFjvlEjkWL00wC4iOAWAGAtZ86INGmSMZ+cLFKmTKBbBCCI8HMXAAAAtkFwCwAAANsguAUAAIBtENwCAADANghuAQAAYBsEtwAAALANSoEBAKxFh9x94omL8wDghuAWAGC94XdffjnQrQAQpEhLAAAAgG3QcwsAsN7wuwcOZMzXrMnwuwA8ENwCAKw3/G6dOhnzDL8LIBN+7gIAAMA2CG4BAABgGwS3AAAAsI2gD25r164tISEhWaZhw4aZ9R06dMiy7uGHHw50swEAABAAQX9B2ZYtWyQtLc11f+fOnXLLLbfIHXfc4Vo2ePBgmThxout+6dKl/d5OAAAABF7QB7fVqlXzuP/CCy/IFVdcITfeeKNHMBsZGRmA1gEAACCYBH1agrtz587JBx98IA8++KBJP3CaM2eOVK1aVZo0aSJxcXFy+vTpHPeTmpoqSUlJHhMAwCJCQ0X+8peMSecBwI2l/ldYvHixnDhxQgYOHOhads8990itWrUkOjpaduzYIaNHj5Y9e/bIokWLvO4nPj5eJkyY4KdWAwB8KjxcZOrUQLcCQJCyVHA7Y8YM6datmwlknYYMGeKab9q0qURFRUnHjh1l7969Jn0hO9q7O2rUKNd97bmtUaNGIbceAAAAhc0ywe3+/ftl5cqVOfbIqpiYGHObkJDgNbgNDw83EwDAghwOkd9/z5ivWlXELU0NACwT3M6cOVMiIiKkR48eOW63fft2c6s9uAAAG9LrKiIiMuYZfheAFYPb9PR0E9zGxsZKqNvFA5p6MHfuXOnevbtUqVLF5NyOHDlSbrjhBmnWrFlA2wwAAAD/s0Rwq+kIBw4cMFUS3IWFhZl1U6ZMkZSUFJM3269fPxkzZkzA2goAAIDAsURw27lzZ3FojlUmGsyuW7cuIG0CAABA8LFUnVsAAAAgJwS3AAAAsA2CWwAAANiGJXJuAQBw0ao5sbEX5wHADf8rAACsRQfhmTUr0K0AEKRISwAAAIBt0HMLALAWLQ2po5Sp0qUZfheAB3puAQDWooFt2bIZkzPIBYD/IbgFAACAbRDcAgAAwDYIbgEAAGAbBLcAAACwDYJbAAAA2AbBLQAAAGyDOrcAAGspXlzk9tsvzgOAG4JbAIC1lCwpsmBBoFsBIEiRlgAAAADbILgFAACAbRDcAgCsJSVFJCQkY9J5AHBDcAsAAADbILgFAACAbRDcAgAAwDYIbgEAAGAbBLcAAACwDYJbAAAA2EZQB7fjx4+XkJAQj6lhw4au9WfPnpVhw4ZJlSpVpGzZstKvXz85cuRIQNsMAChkOuRu9+4ZE8PvArDa8LtXXXWVrFy50nU/NPRik0eOHCmffvqpLFiwQCpUqCDDhw+Xvn37yldffRWg1gIA/DL87qefBroVAIJU0Ae3GsxGRkZmWX7y5EmZMWOGzJ07V26++WazbObMmdKoUSPZtGmTXHPNNQFoLQAAAAIpqNMS1M8//yzR0dFSt25dGTBggBw4cMAs37Ztm5w/f146derk2lZTFmrWrCkbN27McZ+pqamSlJTkMQEAAMD6gjq4jYmJkVmzZsny5ctl2rRpkpiYKNdff72cOnVKDh8+LGFhYVKxYkWPx1SvXt2sy0l8fLxJY3BONWrUKORXAgDwGR1yt0yZjInhdwFYKS2hW7durvlmzZqZYLdWrVoyf/58KVWqVIH3GxcXJ6NGjXLd155bAlwAsJDTpwPdAgBBKqh7bjPTXtr69etLQkKCycM9d+6cnDhxwmMbrZaQXY6uu/DwcClfvrzHBAAAAOuzVHCbnJwse/fulaioKGnVqpWUKFFCVq1a5Vq/Z88ek5Pbrl27gLYTAAAAgRHUaQlPPPGE9OrVy6QiHDx4UMaNGyfFixeX/v37m1zZQYMGmfSCypUrm97XESNGmMCWSgkAAABFU1AHt7/++qsJZP/44w+pVq2aXHfddabMl86ryZMnS7FixczgDVoBoUuXLvLWW28FutkAAAAIkKAObufNm5fj+pIlS8rUqVPNBAAAAAR1cAsAQBbFionceOPFeQBwQ3ALALAWLQW5dm2gWwEgSPGTFwAAALZBcAsAAADbILgFAFiLDrmrVXN0YvhdAJmQcwsAsJ7ffw90CwAEKXpuAQAAYBsEtwAAALANglsAAADYBsEtAAAAbIPgFgAAALZBtQQAgLXokLutW1+cBwA3BLcAAOsNv7tlS6BbASBI8ZMXAAAAtkFwCwAAANsguAUAWMvp0yK1a2dMOg8Absi5BQBYi8Mhsn//xXkAcEPPLQAAAGyD4BYAAAC2QXALAAAA2yC4BQAAgG0Q3AIAAMA2qJYAALCWkBCRxo0vzgOAG4JbAIC1lC4t8sMPgW4FgCBFWgIAAABsg+AWAAAAthHUwW18fLy0adNGypUrJxEREdKnTx/Zs2ePxzYdOnSQkJAQj+nhhx8OWJsBAIVMh9y96qqMieF3AVgpuF23bp0MGzZMNm3aJCtWrJDz589L586dJSUlxWO7wYMHy6FDh1zTSy+9FLA2AwAKmQ65u2tXxsTwuwCsdEHZ8uXLPe7PmjXL9OBu27ZNbrjhBtfy0qVLS2RkZABaCAAAgGAS1D23mZ08edLcVq5c2WP5nDlzpGrVqtKkSROJi4uT07mcpkpNTZWkpCSPCQAAANYX1D237tLT0+Xxxx+Xa6+91gSxTvfcc4/UqlVLoqOjZceOHTJ69GiTl7to0aIcc3knTJjgp5YDAADAXywT3Gru7c6dO+XLL7/0WD5kyBDXfNOmTSUqKko6duwoe/fulSuuuCLbfWnv7qhRo1z3tee2Ro0ahdh6AAAA+IMlgtvhw4fL0qVLZf369XL55ZfnuG1MTIy5TUhI8BrchoeHmwkAAAD2EtTBrcPhkBEjRsjHH38sa9eulTp16uT6mO3bt5tb7cEFANiQDrlbq9bFeQCwSnCrqQhz586VTz75xNS6PXz4sFleoUIFKVWqlEk90PXdu3eXKlWqmJzbkSNHmkoKzZo1C3TzAQCFNfzuvn2BbgWAIBXUwe20adNcAzW4mzlzpgwcOFDCwsJk5cqVMmXKFFP7VvNm+/XrJ2PGjAlQiwEAABBIQZ+WkBMNZnWgBwAAAMBydW4BAJAzZ0TatMmYdB4ArNJzCwBAFunpIlu3XpwHADf03AIAAMA2CG4BAABgGwS3AAAAsA2CWwAAANgGwS0AAABsg2oJAADrqVo10C0AEKQIbgEA1lKmjMixY4FuBYAgRVoCAAAAbIPgFgAAALZBcAsAsBYdcrdDh4yJ4XcBZELOLQDAWnTI3XXrLs4DgBt6bgEAAGAbBLcAAACwDYJbAAAA2AbBLQAAAGyD4BYAAAC2QbUEAID1lC4d6BYACFIEtwAA6w2/m5IS6FYACFKkJQAAAMA2CG4BAABgGwS3AABrOXtWpEePjEnnAcANObcAAGtJSxNZtuziPAC4oecWAAAAtkFwCwAAANuwTXA7depUqV27tpQsWVJiYmLk66+/DnSTAAAA4Ge2CG4/+ugjGTVqlIwbN06++eYbad68uXTp0kWOHj0a6KYBAADAj2wR3E6aNEkGDx4sDzzwgDRu3FimT58upUuXlvfeey/QTQMAAIAfWb5awrlz52Tbtm0SFxfnWlasWDHp1KmTbNy4MdvHpKammsnp5MmT5jYpKUn85fDhw2byBX296enpQbMf9mX9NhWFfQVjm3y5rz179pjb1MMJkn6u4OWyzh//1dzq/7PJyclB8fqKnT0rzf43v+OrryS9ZMmgaJev9xWMbSoK+wrGNgXrviIjI83kL844zeFw2Du4/f333yUtLU2qV6/usVzv7969O9vHxMfHy4QJE7Isr1GjRqG1EwAC4c/P3/TJfoYMGSJBqUuXQLcAgJ+dOnVKKlSoYN/gtiC0l1dzdJ3018vx48elSpUqEhISEtC22YH+stIfCr/88ouUL18+0M1BAXAMrY9jaG0cP+vjGPqe9thqYBsdHZ3jdpYPbqtWrSrFixeXI0eOeCzX+966ysPDw83krmLFioXazqJI/5j5g7Y2jqH1cQytjeNnfRxD38qpx9Y2F5SFhYVJq1atZNWqVR49sXq/Xbt2AW0bAAAA/MvyPbdKUwxiY2OldevW0rZtW5kyZYqkpKSY6gkAAAAoOmwR3N51111y7NgxGTt2rKlA0KJFC1m+fHmWi8zgH5ryoTWHM6d+wDo4htbHMbQ2jp/1cQwDJ8SRWz0FAAAAwCIsn3MLAAAAOBHcAgAAwDYIbgEAAGAbBLcAAACwDYJb+ISO8DZgwABTqFoHxBg0aFCex6HXaxq7detmRodbvHhxobcVvjmGuv2IESOkQYMGUqpUKalZs6Y8+uijcvLkSb+2uyibOnWq1K5dW0qWLCkxMTHy9ddf57j9ggULpGHDhmb7pk2byrJly/zWVlza8XvnnXfk+uuvl0qVKpmpU6dOuR5vBN/foNO8efPMd16fPn0KvY1FEcEtfEKDoh9++EFWrFghS5culfXr1+d5LHqtS8ywx9Y7hgcPHjTTK6+8Ijt37pRZs2aZEnwaFKPwffTRR6bGt5Ya+uabb6R58+bSpUsXOXr0aLbbb9iwQfr372+Oz7fffmu+VHXSY4fgP35r1641x2/NmjWyceNGM6xr586d5bfffvN721GwY+i0b98+eeKJJ8yPFRQSLQUGXIpdu3ZpOTnHli1bXMs+++wzR0hIiOO3337L8bHffvut47LLLnMcOnTI7OPjjz/2Q4vhy2Pobv78+Y6wsDDH+fPnC6mlcGrbtq1j2LBhrvtpaWmO6OhoR3x8fLbb33nnnY4ePXp4LIuJiXEMHTq00NuKSz9+mV24cMFRrlw5x+zZswuxlfD1MdTj1r59e8e7777riI2Nddx6661+am3RQs8tLpn2IuhpbB0hzklPmRUrVkw2b97s9XGnT5+We+65x5zWiYyM9FNr4ctjmJmmJGhaQ2ioLcaHCVrnzp2Tbdu2mWPkpMdK7+uxzI4ud99eaS+Tt+0RXMcvu/8/z58/L5UrVy7ElsLXx3DixIkSERHBGa5CxjcQLpmOCqd/rO40uNH/dHWdNyNHjpT27dvLrbfe6odWojCOobvff/9dnnvuuTyno6Dg9L1OS0vLMgqj3t+9e3e2j9HjmN32eT2+COzxy2z06NESHR2d5QcLgvcYfvnllzJjxgzZvn27n1pZdNFzC6+efvppkwub05TX/4gzW7Jkiaxevdrk28Kax9BdUlKS9OjRQxo3bizjx4/3SdsBZO+FF14wFyR9/PHH5kImBL9Tp07JfffdZy4MrFq1aqCbY3v03MKrv/71rzJw4MAct6lbt65JKcicQH/hwgVzNb23dAMNbPfu3WtOhbvr16+fSbLXiycQ3MfQ/T/trl27Srly5cyXbYkSJXzSdninX47FixeXI0eOeCzX+96Oly7Pz/YIruPnpBdwanC7cuVKadasWSG3FL46hvp9pxeS9erVy7UsPT3ddZZsz549csUVV/ih5UUDwS28qlatmply065dOzlx4oTJP2rVqpUreNU/XC2N4q1H8aGHHvJYpqWJJk+e7PHHj+A9hs4eW83bDA8PN73x9CL5R1hYmDlOq1atcpUS0mOl94cPH+71GOv6xx9/3LVMK2PocgT/8VMvvfSSPP/88/L555975Mcj+I+hluD7/vvvPZaNGTPGdA689tprpvoFfCjQV7TBHrp27epo2bKlY/PmzY4vv/zSUa9ePUf//v1d63/99VdHgwYNzHpvqJZgrWN48uRJc7V906ZNHQkJCabihXPSK4JRuObNm+cIDw93zJo1y1S7GDJkiKNixYqOw4cPm/X33Xef4+mnn3Zt/9VXXzlCQ0Mdr7zyiuPHH390jBs3zlGiRAnH999/H8BXUXTl9/i98MILphLJwoULPf7WTp06FcBXUbTl9xhmRrWEwkPPLXxizpw55tdqx44dzRWjml7w+uuvu9brVb162kWv8IU9jqHWdXRWUrjyyis99pWYmGgKm6Pw3HXXXXLs2DEZO3asuSisRYsWps6w8wKXAwcOmOPopBdvzp071/QWPfPMM1KvXj0zaEqTJk0C+CqKrvwev2nTppkr9G+//XaP/WiNVfLcrXEM4T8hGuH68fkAAACAQsNPCgAAANgGwS0AAABsg+AWAAAAtkFwCwAAANsguAUAAIBtENwCAADANghuAQAAYBsEtwAAALANglsAAADYBsEtABSykJCQHKdLGT5VH6/D6BakDdddd53H+pIlS8r+/fs9HtenTx8ZOHCg677OOx9fokQJqVOnjjz11FNy9uzZAr8GAPClUJ/uDQCQxaFDh1zzH330kRmLfs+ePa5lZcuW9Us7Zs6cKV27dnXdDwsL81ivAau2bfbs2TnuR/eh+zp//rxs27ZNYmNjzWNffPHFQms7AOQVPbcAUMgiIyNdU4UKFUwg6L5s3rx50qhRI9Nz2rBhQ3nrrbdcjz137pwMHz5coqKizPpatWpJfHy8WVe7dm1ze9ttt5l9Ou97U7FiRY/nrVy5ssd6fZ4PPvhAdu7cmeN+wsPDzeNr1KhhenY7deokK1asuIR3CAB8h55bAAigOXPmmN7SN998U1q2bCnffvutDB48WMqUKWN6RF9//XVZsmSJzJ8/X2rWrCm//PKLmdSWLVskIiLC1SNbvHjxS2rLtddeKz/99JM8/fTTsnTp0jw9RgPhDRs2mKAbAIIBwS0ABNC4cePk1Vdflb59+5r7msO6a9cu+ec//2mC2wMHDki9evVMfqz2zroHkdWqVfPokc1N//79PQJg7aXVnld32ivcrFkz+eKLL+T666/Pdj8a+GoqxYULFyQ1NVWKFStmgnMACAYEtwAQICkpKbJ3714ZNGiQ6a110qBR0xecF3Ddcsst0qBBA9M727NnT+ncuXOBnm/y5MkmhcBJUx0ya9y4sdx///2m9/arr77Kdj833XSTTJs2zbRf9xkaGir9+vUrUJsAwNcIbgEgQJKTk83tO++8IzExMR7rnD2sV199tSQmJspnn30mK1eulDvvvNMEqAsXLsz382nv7pVXXpnrdhMmTJD69et7rcKgKRPO/bz33nvSvHlzmTFjhgnSASDQuKAMAAKkevXqEh0dLf/9739NsOg+aXqCU/ny5eWuu+4yQbBWW/jXv/4lx48fN+u0HFdaWppP26UXiunFZc8880yu+9aUBN1uzJgxcubMGZ+2AwAKguAWAAJIe0k1z1UvHNOLub7//ntzgdikSZPMer398MMPZffu3Wb9ggULTA+s5tkqrZCwatUqOXz4sPz5558+a1dcXJwcPHjQ9Bbn5o477jA9zVOnTvXZ8wNAQRHcAkAAPfTQQ/Luu++agLZp06Zy4403yqxZs1w9t+XKlZOXXnpJWrduLW3atJF9+/bJsmXLTI+p0ovRtAyX9rZqtQVf0TJho0ePztPgDJpzqz292k7NwwWAQApxOByOgLYAAAAA8BF6bgEAAGAbBLcAAACwDYJbAAAA2AbBLQAAAGyD4BYAAAC2QXALAAAA2yC4BQAAgG0Q3AIAAMA2CG4BAABgGwS3AAAAsA2CWwAAAIhd/D/PGAHnOEFenQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n_trials = 200\n", "violations = 0\n", "test_fnrs = []\n", "\n", "for seed in range(n_trials):\n", " # Re-split calibration and test each trial\n", " X_c, X_t, y_c, y_t = train_test_split(X_temp, y_temp, test_size=0.5, random_state=seed)\n", " scores_c = model.predict_proba(X_c)[:, 1]\n", " scores_t = model.predict_proba(X_t)[:, 1]\n", "\n", " # TODO: run your full LTT pipeline on this split\n", " # (compute risks → p-values → bonferroni → select lambda_hat)\n", " # then measure test FNR and check if it exceeds alpha\n", " #\n", " trial_risks = compute_empirical_risks(scores_c, y_c, lambdas)\n", " trial_pvalues = compute_p_values(trial_risks, len(scores_c), alpha)\n", " trial_cert = bonferroni_filter(lambdas, trial_pvalues, delta)\n", " trial_lhat = select_lambda_hat(trial_cert)\n", " if trial_lhat is not None:\n", " fnr = np.mean([loss(s, l, trial_lhat) for s, l in zip(scores_t, y_t)])\n", " test_fnrs.append(fnr)\n", " if fnr > alpha:\n", " violations += 1\n", " \n", "\n", "# Uncomment after implementing the loop above\n", "print(f\"Violations (FNR > α): {violations} / {n_trials} = {violations/n_trials:.3f}\")\n", "print(f\"Expected: ≤ δ = {delta}\")\n", "\n", "plt.figure(figsize=(8, 4))\n", "plt.hist(test_fnrs, bins=30, edgecolor='black')\n", "plt.axvline(alpha, color='red', linestyle='--', label=f'α = {alpha}')\n", "plt.xlabel('Test FNR'); plt.ylabel('Count'); plt.title('Distribution of Test FNR Across Trials')\n", "plt.legend(); plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": ".venv (3.9.6)", "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.9.6" } }, "nbformat": 4, "nbformat_minor": 5 }