{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5c561b15-228a-4ecd-a6e2-ae02b75c25f5",
   "metadata": {},
   "source": [
    "# Replication Code for \"Survey of Data-driven Newsvendor: Unified Analysis and Spectrum of Achievable Regrets\"\n",
    "\n",
    "This Jupyter Notebook contains the numerical experiments and simulation code aimed at reproducing the figures presented in the paper.\n",
    "\n",
    "### Instructions\n",
    "We have organized the code into sections corresponding to each figure to facilitate direct access to specific simulations. Please follow the guidelines below:\n",
    "\n",
    "1.  **Global Initialization**: The cells under \"Imports,\" \"Plotting Style,\" and \"Global Parameters\" (at the very top) **must be run first** to establish the simulation environment.\n",
    "2.  **Dependencies**: While you can navigate directly to the code for specific figures, please note that some later sections rely on utility functions defined in earlier parts (e.g., `saa_decision` defined in \"Figures 2 & 4\"). We have indicated these dependencies in the headers where applicable.\n",
    "3.  **Outputs**: The scripts will generate and save `.png` files of the plots in the current working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccc4dee6-69a5-4bef-b4a3-5d12a9830a68",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "\n",
    "import scipy.special\n",
    "import scipy.integrate as integrate\n",
    "from scipy.integrate import quad\n",
    "import scipy.stats as stats\n",
    "from scipy.stats import uniform, expon, pareto, lognorm, norm, zipf\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7966af6-4a52-43c2-b65c-f275f36d6257",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- Plotting Style Configuration ---\n",
    "plt.rcParams['figure.figsize'] = [7, 5]\n",
    "plt.rcParams['figure.dpi'] = 300\n",
    "\n",
    "plt.rcParams['axes.titlesize'] = 14\n",
    "plt.rcParams['axes.labelsize'] = 13\n",
    "plt.rcParams['legend.fontsize'] = 14\n",
    "\n",
    "# --- Global Parameters ---\n",
    "# Number of repetitions for simulations\n",
    "M = 10000 \n",
    "\n",
    "# Values of q focused on in the simulations\n",
    "Q_VALUES = [0.4, 0.9]\n",
    "\n",
    "# --- Color Mapping ---\n",
    "colors = {\n",
    "    \"hard-Bernoulli\": \"#FFD700\",\n",
    "    \"Log-normal\": \"#D770A2\",\n",
    "    \"Pareto\": \"#8FBC8F\",\n",
    "    \"Exponential\": \"#FFB347\",\n",
    "    \"Uniform\": \"#6C8EBF\",\n",
    "    \"easy-Bernoulli\": \"#E57373\"\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fb367b5",
   "metadata": {},
   "source": [
    "# Figure 1: CDFs of the distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06e8c16c-3190-4e41-864f-f7ffba7c3500",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 30, 2000)\n",
    "plt.figure(figsize=(8, 5)) \n",
    "\n",
    "# Plot reference lines for q=0.4 and q=0.9\n",
    "for q in Q_VALUES:\n",
    "    plt.axhline(y=q, color='gray', linestyle='--', alpha=0.6)\n",
    "\n",
    "# Plot CDFs of the considered distributions\n",
    "plt.plot(x, uniform.cdf(x, 0, 1), label=\"Uniform(0,1)\", color=colors[\"Uniform\"])\n",
    "plt.plot(x, expon.cdf(x, scale=1), label=\"Exponential(1)\", color=colors[\"Exponential\"])\n",
    "plt.plot(x, pareto.cdf(x, b=1.5, scale=1), label=r\"Pareto($x_m=1,\\alpha=1.5$)\", color=colors[\"Pareto\"])\n",
    "plt.plot(x, lognorm.cdf(x, s=1.805, scale=np.exp(1)), label=r\"Log-normal($\\mu=1,\\sigma=1.805$)\", color=colors[\"Log-normal\"])\n",
    "\n",
    "plt.xlabel(r\"$a$\")\n",
    "plt.ylabel(\"CDF\")\n",
    "plt.title(\"CDFs of Distributions\")\n",
    "plt.legend(loc=\"lower right\", fontsize=18)\n",
    "plt.tight_layout()\n",
    "plt.savefig('distribution_cdf.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe4c5f92",
   "metadata": {},
   "source": [
    "# Figures 2 & 4: Average additive regrets, 95th percentile of additive regrets, and the values of ∆(ε) for the distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0123d88",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- Calculate SAA Decision ---\n",
    "def saa_decision(samples, q, N):\n",
    "    \"\"\"\n",
    "    Compute the SAA decision defined as the empirical q-quantile based on N samples, applied row-wise.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    samples : np.ndarray, shape (M, N)\n",
    "    q : float\n",
    "    N : int\n",
    "        Sample size per replication.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    np.ndarray, shape (M,)\n",
    "        SAA decisions for each replication.\n",
    "    \"\"\"\n",
    "    # Sort data to find order statistics\n",
    "    # We use ceil(N*q) - 1 for 0-based indexing\n",
    "    index = math.ceil(N * q) - 1\n",
    "    samples_sorted = np.sort(samples)\n",
    "    \n",
    "    # data is (M, N), we want the decision for each row\n",
    "    return samples_sorted[:, index]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb2e596c-86d4-461b-8457-236c69542a0e",
   "metadata": {},
   "source": [
    "The following functions `simulate_additive_regret_<distribution>` implement the computation of the SAA additive regret for different distributions with M independent repetitions.  \n",
    "For each distribution, the regret is obtained by substituting the explicit form of the CDF into equation (3) in the paper, that is $$\\int_{\\hat a}^{a^*} \\bigl(q - F(z)\\bigr)\\,dz.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0111665",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_additive_regret_bernoulli(N, M, q, p=0.5, c=1):\n",
    "    \"\"\"\n",
    "    Computation of the SAA additive regret for the scaled Bernoulli distribution.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    N : int\n",
    "        Sample size.\n",
    "    M : int\n",
    "        Number of repetitions.\n",
    "    q : float\n",
    "        Critical quantile.\n",
    "    p : float, optional\n",
    "        Probability mass at c.\n",
    "    c : float, optional\n",
    "        Scaling parameter of the Bernoulli distribution.\n",
    "    \"\"\"\n",
    "    # Bernoulli distribution with support {0, c} and P(X = c) = p\n",
    "    samples = np.random.choice([0, c], size=(M, N), p=[1 - p, p])\n",
    "\n",
    "    # SAA decisions\n",
    "    saa_decisions = saa_decision(samples, q, N)\n",
    "\n",
    "    # Optimal decision a*\n",
    "    a_star = c if (1 - p) < q else 0\n",
    "\n",
    "    # Additive regrets for M repetitions (Eq. (3) specialized to scaled Bernoulli)\n",
    "    additive_regrets = (\n",
    "        (1 - p - q) * (saa_decisions - a_star)\n",
    "    )\n",
    "\n",
    "    regret_mean = np.mean(additive_regrets)\n",
    "    regret_95pct = np.percentile(additive_regrets, 95)\n",
    "\n",
    "    return additive_regrets, regret_mean, regret_95pct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce4b31cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_additive_regret_uniform(N, M, q):\n",
    "    \"\"\"\n",
    "    Computation of the SAA additive regret for the Uniform(0,1) distribution.\n",
    "    \"\"\"\n",
    "    # Uniform(0,1) samples\n",
    "    samples = np.random.uniform(size=(M, N))\n",
    "\n",
    "    # SAA decisions\n",
    "    saa_decisions = saa_decision(samples, q, N)\n",
    "\n",
    "    # Additive regrets for M repetitions (Eq. (3) specialized to Uniform)\n",
    "    additive_regrets = (\n",
    "        0.5 * q**2 - saa_decisions * q + 0.5 * saa_decisions**2\n",
    "    )\n",
    "\n",
    "    regret_mean = np.mean(additive_regrets)\n",
    "    regret_95pct = np.percentile(additive_regrets, 95)\n",
    "\n",
    "    return additive_regrets, regret_mean, regret_95pct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13d1bd0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_additive_regret_exponential(N, M, q):\n",
    "    \"\"\"\n",
    "    Computation of the SAA additive regret for the Exponential(1) distribution.\n",
    "    \"\"\"\n",
    "    # Exponential(1) samples\n",
    "    samples = np.random.exponential(size=(M, N))\n",
    "\n",
    "    # SAA decisions\n",
    "    saa_decisions = saa_decision(samples, q, N)\n",
    "\n",
    "    # Additive regrets for M repetitions (Eq. (3) specialized to Exponential)\n",
    "    additive_regrets = (\n",
    "        -(q - 1) * np.log(1 - q) - saa_decisions * (q - 1) - 1 + q + np.exp(-saa_decisions)\n",
    "    )\n",
    "\n",
    "    regret_mean = np.mean(additive_regrets)\n",
    "    regret_95pct = np.percentile(additive_regrets, 95)\n",
    "\n",
    "    return additive_regrets, regret_mean, regret_95pct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05d2eb37",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_additive_regret_pareto(N, M, q):\n",
    "    \"\"\"\n",
    "    Computation of the SAA additive regret for the Pareto distribution\n",
    "    with shape parameter alpha = 1.5 and support [1, +∞).\n",
    "    \"\"\"\n",
    "    # Pareto samples with alpha = 1.5 and support [1, +∞)\n",
    "    samples = np.random.pareto(a=1.5, size=(M, N)) + 1\n",
    "\n",
    "    # SAA decisions\n",
    "    saa_decisions = saa_decision(samples, q, N)\n",
    "\n",
    "    # Additive regrets for M repetitions (Eq. (3) specialized to Pareto)\n",
    "    additive_regrets = (\n",
    "        -3 * (1 - q)**(1 / 3) - saa_decisions * (q - 1) + 2 * saa_decisions**(-0.5)\n",
    "    )\n",
    "\n",
    "    regret_mean = np.mean(additive_regrets)\n",
    "    regret_95pct = np.percentile(additive_regrets, 95)\n",
    "\n",
    "    return additive_regrets, regret_mean, regret_95pct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb8e6451-00fa-4e97-be70-2f74d3092d25",
   "metadata": {},
   "source": [
    "For the log-normal distribution with $\\mu = 1$ and $\\sigma = 1.805$, the additive regret is $L(\\hat a)-L(a^*) = \\int_{\\hat a}^{a^*} \\bigl(q - F(z)\\bigr)\\,dz,$ where $F(z) = \\Phi\\!\\left(\\frac{\\log z - 1}{\\sigma}\\right), \\, a^* = \\exp\\!\\left(1 + \\sigma \\Phi^{-1}(q)\\right).$\n",
    "\n",
    "Direct numerical integration with respect to $z$ is computationally expensive, since this integral must be evaluated repeatedly for each SAA realization. To improve efficiency, we apply the change of variables $t = \\frac{\\log z - 1}{\\sigma},\\, z = \\exp(1 + \\sigma t),\\, dz = \\sigma \\exp(1 + \\sigma t)\\,dt.$ The lower and upper integration limits become $\\hat{t} = \\frac{\\log \\hat a - 1}{\\sigma},\\, t^* = \\Phi^{-1}(q).$\n",
    "\n",
    "Substituting into the original integral yields\n",
    "$$\n",
    "\\begin{aligned}\n",
    "L(\\hat a)-L(a^*)\n",
    "&= \\int_{\\hat{t}}^{t^*}\n",
    "\\bigl(q - \\Phi(t)\\bigr)\\,\n",
    "\\sigma \\exp(1 + \\sigma t)\\,dt \\\\\n",
    "&=q\\int_{\\hat{t}}^{t^*} \\sigma\\exp(1 + \\sigma t)\\,dt-\n",
    "\\int_{\\hat{t}}^{t^*}\\Phi(t)\\sigma\\exp(1 + \\sigma t)\\,dt\\\\\n",
    "&=q\\exp(1+\\sigma t)\\big|_{\\hat{t}}^{t^*}\n",
    "-\\Phi(t)\\exp(1+\\sigma t)\\big|_{\\hat{t}}^{t^*}+\n",
    "\\int_{\\hat{t}}^{t^*}\\phi(t)\\exp(1+\\sigma t)\\,dt\\\\\n",
    "&=\\bigl(\\Phi(\\hat{t})-q\\bigr)\\hat{a}+\\frac{e}{\\sqrt{2\\pi}}\\int_{\\hat{t}}^{t^*}\\exp\\left(-\\frac12t^2+\\sigma t\\right)\\,dt.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Here $\\Phi(t)$ and $\\phi(t)=\\frac{1}{\\sqrt{2\\pi}}\\exp\\left(-\\frac12 t^2\\right)$ denote the standard normal CDF and PDF,\n",
    "and the remaining integral is evaluated numerically through `_lognormal_aux_integral`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea75e16-9f30-4f7d-99e6-b25768d64999",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _lognormal_aux_integral(q, t_hat):\n",
    "    \"\"\"\n",
    "    Auxiliary integral appearing in the transformed regret expression for the Log-normal distribution.\n",
    "    \"\"\"\n",
    "    def integrand(t):\n",
    "        return math.exp(-0.5 * t**2 + 1.805 * t)\n",
    "\n",
    "    result, _ = quad(integrand, t_hat, norm.ppf(q))\n",
    "    return result\n",
    "\n",
    "def additive_regret_lognormal(q, a_hat):\n",
    "    \"\"\"\n",
    "    Additive regret for the Log-normal distribution (mu = 1, sigma = 1.805),\n",
    "    computed via a change of variables to the normal domain.\n",
    "    \"\"\"\n",
    "    # Change of variables: z = exp(1 + sigma * t)\n",
    "    t_hat = (math.log(a_hat) - 1.0) / 1.805\n",
    "\n",
    "    result = (\n",
    "        (norm.cdf(t_hat) - q) * a_hat + math.exp(1.0) * _lognormal_aux_integral(q, t_hat) / math.sqrt(2 * math.pi)\n",
    "    )\n",
    "\n",
    "    return result\n",
    "\n",
    "def simulate_additive_regret_lognormal(N, M, q):\n",
    "    \"\"\"\n",
    "    Computation of the SAA additive regret for the Log-normal distribution with mu = 1 and sigma = 1.805.\n",
    "    \"\"\"\n",
    "    # Log-normal samples\n",
    "    samples = lognorm.rvs(s=1.805, scale=math.exp(1.0), size=(M, N))\n",
    "\n",
    "    # SAA decisions\n",
    "    saa_decisions = saa_decision(samples, q, N)\n",
    "\n",
    "    # Additive regrets for M repetitions\n",
    "    additive_regrets = [\n",
    "        additive_regret_lognormal(q, a_hat)\n",
    "        for a_hat in saa_decisions\n",
    "    ]\n",
    "\n",
    "    regret_mean = np.mean(additive_regrets)\n",
    "    regret_95pct = np.percentile(additive_regrets, 95)\n",
    "\n",
    "    return additive_regrets, regret_mean, regret_95pct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f8d6e63-ccb7-45a9-b64e-1e1133e5eec5",
   "metadata": {},
   "source": [
    "The following functions `delta_<distribution>` compute the quantity $\\Delta(\\varepsilon)$ defined in equation (36) of the paper:\n",
    "$$\\Delta(\\varepsilon):=\\max\\{F^{-1}(\\min\\{q+\\varepsilon,1\\})-a^*,a^*-F^{-1}(\\max\\{q-\\varepsilon,0\\})\\}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f9eab86",
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta_bernoulli(eps, q, p=0.5, c=1.0):\n",
    "    \"\"\"\n",
    "    Compute Delta(eps) for the scaled Bernoulli distribution supported on {0, c},\n",
    "    where p is the probability of taking value c.\n",
    "    \"\"\"\n",
    "    def inverse_cdf(y):\n",
    "        # Inverse CDF of the scaled Bernoulli distribution\n",
    "        return 1e-10 if y <= 1 - p else c\n",
    "\n",
    "    # Perturbed probability levels\n",
    "    prob_left = max(q - eps, 0.0)\n",
    "    prob_right = min(q + eps, 1.0)\n",
    "\n",
    "    # Optimal decision a*\n",
    "    a_star = inverse_cdf(q)\n",
    "\n",
    "    # Deviations from a*\n",
    "    deviation_left = a_star - inverse_cdf(prob_left)\n",
    "    deviation_right = inverse_cdf(prob_right) - a_star\n",
    "\n",
    "    return max(deviation_left, deviation_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed472884",
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta_uniform(eps, q):\n",
    "    \"\"\"\n",
    "    Compute Delta(eps) for the Uniform(0, 1) distribution.\n",
    "    \"\"\"\n",
    "    # Perturbed probability levels\n",
    "    prob_left = max(q - eps, 0.0)\n",
    "    prob_right = min(q + eps, 1.0)\n",
    "\n",
    "    # For Uniform(0,1), F^{-1}(y) = y, a* = q\n",
    "    deviation_left = q - prob_left\n",
    "    deviation_right = prob_right - q\n",
    "\n",
    "    return max(deviation_left, deviation_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32b8dbff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta_exponential(eps, q):\n",
    "    \"\"\"\n",
    "    Compute Delta(eps) for the Exponential(1) distribution.\n",
    "    \"\"\"\n",
    "    def inverse_cdf(y):\n",
    "        # Inverse CDF of Exponential(1)\n",
    "        return 1e3 if y == 1.0 else -math.log(1 - y)\n",
    "\n",
    "    # Perturbed probability levels\n",
    "    prob_left = max(q - eps, 0.0)\n",
    "    prob_right = min(q + eps, 1.0)\n",
    "\n",
    "    # Optimal decision a*\n",
    "    a_star = inverse_cdf(q)\n",
    "\n",
    "    # Deviations from a*\n",
    "    deviation_left = a_star - inverse_cdf(prob_left)\n",
    "    deviation_right = inverse_cdf(prob_right) - a_star\n",
    "\n",
    "    return max(deviation_left, deviation_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f17d9acd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta_pareto(eps, q):\n",
    "    \"\"\"\n",
    "    Compute Delta(eps) for the Pareto distribution with shape parameter 1.5 and support starting at 1.\n",
    "    \"\"\"\n",
    "    def inverse_cdf(y):\n",
    "        # Inverse CDF of the Pareto distribution\n",
    "        return 1e3 if y == 1.0 else (1 - y) ** (-2 / 3)\n",
    "\n",
    "    # Perturbed probability levels\n",
    "    prob_left = max(q - eps, 0.0)\n",
    "    prob_right = min(q + eps, 1.0)\n",
    "\n",
    "    # Optimal decision a*\n",
    "    a_star = inverse_cdf(q)\n",
    "\n",
    "    # Deviations from a*\n",
    "    deviation_left = a_star - inverse_cdf(prob_left)\n",
    "    deviation_right = inverse_cdf(prob_right) - a_star\n",
    "\n",
    "    return max(deviation_left, deviation_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f32e7b2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta_lognormal(eps, q):\n",
    "    \"\"\"\n",
    "    Compute Delta(eps) for the Log-normal distribution with mu = 1 and sigma = 1.805.\n",
    "    \"\"\"\n",
    "    def inverse_cdf(y):\n",
    "        # Inverse CDF of the log-normal distribution\n",
    "        return 1e3 if y == 1.0 else math.exp(1.0 + 1.805 * norm.ppf(y))\n",
    "\n",
    "    # Perturbed probability levels\n",
    "    prob_left = max(q - eps, 0.0)\n",
    "    prob_right = min(q + eps, 1.0)\n",
    "\n",
    "    # Optimal decision a*\n",
    "    a_star = inverse_cdf(q)\n",
    "\n",
    "    # Deviations from a*\n",
    "    deviation_left = a_star - inverse_cdf(prob_left)\n",
    "    deviation_right = inverse_cdf(prob_right) - a_star\n",
    "\n",
    "    return max(deviation_left, deviation_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "771d4039",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mapping from distribution names to additive regret simulation functions\n",
    "distributions = {\n",
    "    \"hard-Bernoulli\": simulate_additive_regret_bernoulli,\n",
    "    \"Log-normal\": simulate_additive_regret_lognormal,\n",
    "    \"Pareto\": simulate_additive_regret_pareto,\n",
    "    \"Exponential\": simulate_additive_regret_exponential,\n",
    "    \"Uniform\": simulate_additive_regret_uniform,\n",
    "    \"easy-Bernoulli\": simulate_additive_regret_bernoulli\n",
    "}\n",
    "\n",
    "# Mapping from distribution names to Delta(eps) functions\n",
    "delta_functions = {\n",
    "    \"hard-Bernoulli\": delta_bernoulli,\n",
    "    \"Log-normal\": delta_lognormal,\n",
    "    \"Pareto\": delta_pareto,\n",
    "    \"Exponential\": delta_exponential,\n",
    "    \"Uniform\": delta_uniform,\n",
    "    \"easy-Bernoulli\": delta_bernoulli\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b0731af-feee-4368-a3f9-b44b2a423284",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Figure-specific configurations for different values of q\n",
    "# This includes Bernoulli parameters and axis ranges used in the plots\n",
    "FIG2_AND_4_CONFIG = {\n",
    "    0.4: {\n",
    "        \"p_easy\": 0.45,\n",
    "        \"c_easy\": 1,\n",
    "        \"p_hard\": 0.59,\n",
    "        \"c_hard\": 23,\n",
    "        \"regret_ylim\": (5e-4, 0.2),\n",
    "        \"delta_ylim\": (0.03, 30),\n",
    "        \"percentile_ylim\": (1e-3, 3),\n",
    "    },\n",
    "    0.9: {\n",
    "        \"p_easy\": 0.25,\n",
    "        \"c_easy\": 1,\n",
    "        \"p_hard\": 0.11,\n",
    "        \"c_hard\": 127,\n",
    "        \"regret_ylim\": (1e-4, 4),\n",
    "        \"delta_ylim\": (2e-2, 300),\n",
    "        \"percentile_ylim\": (6e-4, 15),\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "352115e7-0245-4202-b0ca-878e32a17b7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_fig2_and_4_data(q, config):\n",
    "    \"\"\"\n",
    "    Compute all numerical quantities required for Figure 2\n",
    "    for a fixed value of q.\n",
    "    \"\"\"\n",
    "    # Sample sizes used for additive regret simulations\n",
    "    N_list = np.arange(1, 200, 5)\n",
    "\n",
    "    # Values of 1/epsilon^2 used for Delta(eps)\n",
    "    x_list = np.arange(1, 500, 1)\n",
    "\n",
    "    # Containers for additive regret results\n",
    "    regret_results = {\n",
    "        dist: {\"add_regret_mean\": [], \"add_regret_percentile\": []}\n",
    "        for dist in distributions\n",
    "    }\n",
    "\n",
    "    # Additive regret simulations\n",
    "    for N in tqdm(N_list):\n",
    "        for dist, func in distributions.items():\n",
    "            if dist == \"easy-Bernoulli\":\n",
    "                _, mean, pct = func(N, M, q, p=config[\"p_easy\"], c=config[\"c_easy\"])\n",
    "            elif dist == \"hard-Bernoulli\":\n",
    "                _, mean, pct = func(N, M, q, p=config[\"p_hard\"], c=config[\"c_hard\"])\n",
    "            else:\n",
    "                _, mean, pct = func(N, M, q)\n",
    "\n",
    "            regret_results[dist][\"add_regret_mean\"].append(mean)\n",
    "            regret_results[dist][\"add_regret_percentile\"].append(pct)\n",
    "\n",
    "    # Containers for Delta(epsilon)\n",
    "    delta_results = {dist: [] for dist in distributions}\n",
    "\n",
    "    # Delta(eps) computations\n",
    "    for x in tqdm(x_list):\n",
    "        eps = 1 / math.sqrt(x)\n",
    "        for dist, func in delta_functions.items():\n",
    "            if dist == \"easy-Bernoulli\":\n",
    "                delta = func(eps, q, p=config[\"p_easy\"], c=config[\"c_easy\"])\n",
    "            elif dist == \"hard-Bernoulli\":\n",
    "                delta = func(eps, q, p=config[\"p_hard\"], c=config[\"c_hard\"])\n",
    "            else:\n",
    "                delta = func(eps, q)\n",
    "\n",
    "            delta_results[dist].append(delta)\n",
    "\n",
    "    return N_list, x_list, regret_results, delta_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37c1d832-f1c8-4019-a8b0-8a96e63d0a2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_fig2_and_4(q, config, N_list, x_list, regret_results, delta_results):\n",
    "    \"\"\"\n",
    "    Generate the plots for Figures 2 and 4 for a fixed value of q.\n",
    "    \"\"\"\n",
    "    # --- Figure 2 ---\n",
    "    fig, axs = plt.subplots(2, 1, figsize=(7, 10))\n",
    "\n",
    "    # Panel 1: average additive regret\n",
    "    for dist, data in regret_results.items():\n",
    "        axs[0].plot(N_list, data[\"add_regret_mean\"], marker='o', markersize=2, label=dist, color=colors[dist])\n",
    "\n",
    "    axs[0].set_yscale(\"log\")\n",
    "    axs[0].set_ylim(*config[\"regret_ylim\"])\n",
    "    axs[0].set_xlabel(r\"Sample Size $n$\")\n",
    "    axs[0].set_ylabel(\"Additive Regret\")\n",
    "    axs[0].set_title(rf\"Average Additive Regrets for $q={q}$\")\n",
    "    axs[0].legend()\n",
    "\n",
    "    # Panel 2: Delta(eps)\n",
    "    for dist, values in delta_results.items():\n",
    "        axs[1].plot(x_list, values, label=dist, color=colors[dist])\n",
    "\n",
    "    axs[1].set_yscale(\"log\")\n",
    "    axs[1].set_ylim(*config[\"delta_ylim\"])\n",
    "    axs[1].set_xlabel(r\"$1/\\epsilon^2$\")\n",
    "    axs[1].set_ylabel(r\"$\\Delta(\\epsilon)$\")\n",
    "    axs[1].set_title(rf\"$\\Delta(\\epsilon)$ for $q={q}$\")\n",
    "    axs[1].legend(loc=\"upper right\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(f\"add_mean_{q}.png\")\n",
    "    plt.show()\n",
    "\n",
    "    # --- Figure 4 ---\n",
    "    plt.figure()\n",
    "    for dist, data in regret_results.items():\n",
    "        plt.plot(N_list, data[\"add_regret_percentile\"], marker='o', markersize=2, label=dist, color=colors[dist])\n",
    "\n",
    "    plt.yscale(\"log\")\n",
    "    plt.ylim(*config[\"percentile_ylim\"])\n",
    "    plt.xlabel(r\"Sample Size $n$\")\n",
    "    plt.ylabel(\"Additive Regret\")\n",
    "    plt.title(rf\"95th Percentile of Additive Regrets for $q={q}$\")\n",
    "    plt.legend()\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(f\"add_percentile_{q}.png\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0af015bf-48d2-4927-8dee-9f4269dc9d3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate Figures 2 and 4 for each value of q\n",
    "for q in Q_VALUES:\n",
    "    config = FIG2_AND_4_CONFIG[q]\n",
    "    N_list, x_list, regret_results, delta_results = compute_fig2_and_4_data(q, config)\n",
    "    plot_fig2_and_4(q, config, N_list, x_list, regret_results, delta_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c850a09",
   "metadata": {},
   "source": [
    "# Figure 3: Examples of clustered distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51b48b77-8b83-481f-81ae-cc129178fcc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def clustered_cdf(x, beta, gamma=1):\n",
    "    \"\"\"\n",
    "    CDF of a (beta, gamma, zeta=0.5)-clustered distribution on [0, 1] with a*=0.5.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array-like\n",
    "        Evaluation points.\n",
    "    beta : int or float\n",
    "    gamma : float, optional\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    cdf_values : ndarray\n",
    "        Values of the CDF evaluated at x.\n",
    "    \"\"\"\n",
    "    cdf_values = np.piecewise(\n",
    "        x,\n",
    "        [x < 0, (0 <= x) & (x < 0.5), (0.5 < x) & (x < 1), x >= 1],\n",
    "        [0, lambda x: 0.5 - (gamma * (0.5 - x)) ** (beta + 1), \n",
    "         lambda x: 0.5 + (gamma * (x - 0.5)) ** (beta + 1), 1]\n",
    "    )\n",
    "    return cdf_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f9fd165-9c45-4f58-8326-cc1b9c318fdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grid of evaluation points\n",
    "x = np.linspace(-0.05, 1.05, 400)\n",
    "\n",
    "# Values of beta\n",
    "beta_list = [0, 1, 2, 3]\n",
    "\n",
    "plt.figure(figsize=(7, 6))\n",
    "\n",
    "# Reference line at y=0.5\n",
    "plt.axhline(y=0.5, color='gray', linestyle='--')\n",
    "\n",
    "# Plot the CDF for each value of beta\n",
    "for i, beta in enumerate(beta_list):\n",
    "    color = list(colors.values())[-(i + 1)]\n",
    "    cdf_values = clustered_cdf(x, beta)\n",
    "\n",
    "    # Insert NaNs at the boundaries to avoid artificial line connections\n",
    "    cdf_values[(x == 0) | (x == 1)] = np.nan\n",
    "\n",
    "    plt.plot( x, cdf_values, label=rf\"$\\beta={beta}$\", color=color, linewidth=2)\n",
    "\n",
    "plt.xlabel(r\"$a$\")\n",
    "plt.ylabel(\"CDF\")\n",
    "plt.title(\"CDFs of Clustered Distributions\")\n",
    "\n",
    "plt.legend(loc=\"lower right\", fontsize=18)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"clustered_cdf.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "368bd6da",
   "metadata": {},
   "source": [
    "# Figure 5: Minimum PDF fails\n",
    "\n",
    "This figure reuses the functions `saa_decision`, `simulate_additive_regret_uniform`, and `delta_uniform` defined in Figures 2 & 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ddb5be2-fcfa-40b5-81da-ed89714642ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------------------------------------------\n",
    "# Subfigure (a): CDFs of the two distributions\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "# Grid for evaluating the CDFs\n",
    "x = np.linspace(0, 1, 1000)\n",
    "\n",
    "# CDF of the Uniform(0,1) distribution\n",
    "cdf_uniform = x\n",
    "\n",
    "\n",
    "def cdf_red(x):\n",
    "    \"\"\"\n",
    "    CDF of the red distribution used in Figure 5.\n",
    "    \"\"\"\n",
    "    return np.piecewise(\n",
    "        x,\n",
    "        [x < 0.36, (x >= 0.36) & (x < 9/22), (x >= 9/22) & (x < 13/22),\n",
    "         (x >= 13/22) & (x < 0.64), x >= 0.64],\n",
    "        [lambda x: 0, lambda x: 10 * x - 3.6, lambda x: 0.1 * x + 0.45,\n",
    "         lambda x: 10 * x - 5.4, lambda x: 1]\n",
    "    )\n",
    "\n",
    "\n",
    "# Evaluate the CDFs\n",
    "y_red = cdf_red(x)\n",
    "\n",
    "# Plot the CDFs\n",
    "plt.figure(figsize=(5, 5))\n",
    "plt.plot(x, cdf_uniform, label=\"Uniform\", color=colors[\"Uniform\"])\n",
    "plt.plot(x, y_red, label=\"Red\", color=colors[\"easy-Bernoulli\"])\n",
    "\n",
    "plt.xlabel(r\"$a$\")\n",
    "plt.ylabel(\"CDF\")\n",
    "plt.title(\"CDFs of the Distributions\")\n",
    "plt.legend(loc=\"lower right\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"min_pdf_cdf.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61186159-16e0-44f4-b4bb-beed3dd0450e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_additive_regret_red(N, M, q):\n",
    "    \"\"\"\n",
    "    Simulation of the SAA additive regret for the red distribution.\n",
    "    \"\"\"\n",
    "    # Generate samples: M independent repetitions\n",
    "    u = np.random.uniform(size=(M, N))\n",
    "    samples = np.empty_like(u)\n",
    "\n",
    "    mask1 = (0 <= u) & (u < 27/55)\n",
    "    mask2 = (27/55 <= u) & (u < 28/55)\n",
    "    mask3 = u >= 28/55\n",
    "\n",
    "    samples[mask1] = 0.1 * u[mask1] + 0.36\n",
    "    samples[mask2] = 10 * u[mask2] - 4.5\n",
    "    samples[mask3] = 0.1 * u[mask3] + 0.54\n",
    "\n",
    "    # SAA decisions for M repetitions\n",
    "    saa_decisions = saa_decision(samples, q, N)\n",
    "\n",
    "    # Additive regret for each repetition\n",
    "    a_hat = saa_decisions\n",
    "    additive_regrets = np.zeros_like(a_hat, dtype=float)\n",
    "\n",
    "    # Define masks for different ranges of a_hat\n",
    "    cond1 = a_hat < 9/22\n",
    "    cond2 = (a_hat >= 9/22) & (a_hat < 0.5)\n",
    "    cond3 = (a_hat >= 0.5) & (a_hat < 13/22)\n",
    "    cond4 = a_hat >= 13/22\n",
    "\n",
    "    # Calculate additive regrets\n",
    "    additive_regrets[cond1] = (4.1 * (9/22 - a_hat[cond1]) - 5 * ((9/22)**2 - a_hat[cond1]**2) \n",
    "                               + 0.05 * (0.5 - 9/22) - 0.05 * (0.5**2 - (9/22)**2))\n",
    "    additive_regrets[cond2] = (0.05 * (0.5 - a_hat[cond2]) - 0.05 * (0.5**2 - a_hat[cond2]**2))\n",
    "    additive_regrets[cond3] = (0.05 * (a_hat[cond3]**2 - 0.5**2) - 0.05 * (a_hat[cond3] - 0.5))\n",
    "    additive_regrets[cond4] = (0.05 * ((13/22)**2 - 0.5**2) - 0.05 * (13/22 - 0.5) \n",
    "                               + 5 * (a_hat[cond4]**2 - (13/22)**2) - 5.9 * (a_hat[cond4] - 13/22))\n",
    "\n",
    "    regret_mean = np.mean(additive_regrets)\n",
    "    regret_95pct = np.percentile(additive_regrets, 95)\n",
    "\n",
    "    return additive_regrets, regret_mean, regret_95pct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10b6455-f12a-4be0-aa9c-f80d29782d97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------------------------------------------\n",
    "# Subfigure (b): Average additive regret\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "q = 0.5\n",
    "N_list = np.arange(1, 500, 12)\n",
    "\n",
    "# --- Red distribution ---\n",
    "additive_average_red = []\n",
    "\n",
    "for N in tqdm(N_list):\n",
    "    _, regret_mean, _ = simulate_additive_regret_red(N, M, q)\n",
    "    additive_average_red.append(regret_mean)\n",
    "\n",
    "additive_average_red = np.array(additive_average_red)\n",
    "\n",
    "# --- Uniform distribution ---\n",
    "additive_average_uniform = []\n",
    "\n",
    "for N in tqdm(N_list):\n",
    "    _, regret_mean, _ = simulate_additive_regret_uniform(N, M, q)\n",
    "    additive_average_uniform.append(regret_mean)\n",
    "\n",
    "additive_average_uniform = np.array(additive_average_uniform)\n",
    "\n",
    "\n",
    "# --- Plot ---\n",
    "plt.figure()\n",
    "plt.plot(N_list, additive_average_uniform, label=\"Uniform\", color=colors[\"Uniform\"], marker=\"o\", markersize=2)\n",
    "plt.plot(N_list, additive_average_red, label=\"Red\", color=colors[\"easy-Bernoulli\"], marker=\"o\", markersize=2)\n",
    "\n",
    "plt.title(\"Average Additive Regrets\")\n",
    "plt.xlabel(r\"Sample Size $n$\")\n",
    "plt.ylabel(\"Additive Regret\")\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"min_pdf_regret.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfb091d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------------------------------------------\n",
    "# Subfigure (c): Inverse of the minimum PDF in [a* - ζ, a* + ζ]\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "# Range of ζ values\n",
    "zeta_list = np.linspace(0.0001, 0.5, 500)\n",
    "\n",
    "# For the Uniform distribution, the PDF is constant and equal to 1.\n",
    "# For the red distribution, the minimum PDF over the interval is 0.1.\n",
    "plt.figure()\n",
    "plt.hlines(y=1, xmin=0, xmax=0.14, label=\"Uniform\", color=colors[\"Uniform\"])\n",
    "plt.hlines(y=10, xmin=0, xmax=0.14, label=\"Red\", color=colors[\"easy-Bernoulli\"])\n",
    "\n",
    "plt.xlabel(r\"$\\zeta$\")\n",
    "plt.ylabel(r\"Inverse of Minimum PDF in $[a^*-\\zeta,a^*+\\zeta]$\")\n",
    "plt.ylim(bottom=0, top=12)\n",
    "plt.title(r\"Inverse of Minimum PDF in $[a^*-\\zeta,a^*+\\zeta]$\")\n",
    "\n",
    "plt.legend()\n",
    "plt.gca().invert_xaxis()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"min_pdf.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca1b9032",
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta_red(eps):\n",
    "    \"\"\"\n",
    "    Compute Delta(epsilon) for the red distribution.\n",
    "    \"\"\"\n",
    "    threshold = 28 / 55 - 0.5\n",
    "\n",
    "    if eps <= threshold:\n",
    "        return 10 * eps\n",
    "    else:\n",
    "        return 13 / 22 - 0.5 + (eps - threshold) * 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdf8c03b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------------------------------------------\n",
    "# Subfigure (d): Delta(eps)\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "N_list = np.arange(4, 500, 1) # Values of 1/epsilon^2 used for Delta(eps)\n",
    "q = 0.5\n",
    "\n",
    "# Initialize lists to store results\n",
    "delta_uniform_list = []\n",
    "delta_red_list = []\n",
    "\n",
    "# Compute Delta(eps) for each N\n",
    "for N in tqdm(N_list):\n",
    "    eps = 1 / math.sqrt(N)\n",
    "    delta_uniform_list.append(delta_uniform(eps, q))\n",
    "    delta_red_list.append(delta_red(eps))\n",
    "\n",
    "# Plot\n",
    "plt.figure()\n",
    "plt.plot(N_list, delta_uniform_list, label=\"Uniform\", color=colors[\"Uniform\"])\n",
    "plt.plot(N_list, delta_red_list, label=\"Red\", color=colors[\"easy-Bernoulli\"])\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.xlabel(r\"$1/\\epsilon^2$\")\n",
    "plt.ylabel(r\"$\\Delta(\\epsilon)$\")\n",
    "plt.title(r\"$\\Delta(\\epsilon)$\")\n",
    "\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"min_pdf_delta.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b501445d-832f-463c-83c5-d092c1ada3aa",
   "metadata": {},
   "source": [
    "# Figure 6: Empirical distributions of regrets\n",
    "\n",
    "This figure reuses the functions `saa_decision` and `simulate_additive_regret_<distribution>` defined in Figures 2 & 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fa9b97e-3daa-4bd1-a89f-87b4412a8f03",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_fig6(q, N_small, N_large):\n",
    "    \"\"\"\n",
    "    Generate Figure 6 for a fixed value of q.\n",
    "    This figure compares the distributions of additive regret\n",
    "    for two fixed sample sizes n = N_small and n = N_large.\n",
    "    \"\"\"\n",
    "\n",
    "    # Distributions included in Figure 6\n",
    "    distributions = {\n",
    "        \"Uniform\": simulate_additive_regret_uniform,\n",
    "        \"Exponential\": simulate_additive_regret_exponential,\n",
    "        \"Pareto\": simulate_additive_regret_pareto,\n",
    "        \"Log-normal\": simulate_additive_regret_lognormal,\n",
    "    }\n",
    "\n",
    "    # Create a 2 x 4 grid of subplots\n",
    "    fig, axes = plt.subplots(2, 4, figsize=(14, 7), dpi=300)\n",
    "    fig.suptitle(\n",
    "        rf\"Additive Regret Distributions for $n={N_small}$ vs. $n={N_large}$ ($q={q}$)\",\n",
    "        fontsize=16,\n",
    "        y=0.97,\n",
    "    )\n",
    "\n",
    "    print(f\"Running simulations for M={M}, N_small={N_small}, N_large={N_large}, q={q}...\")\n",
    "\n",
    "    for i, (dist_name, func) in enumerate(tqdm(distributions.items())):\n",
    "        row_idx = i // 2\n",
    "        col_idx_small = (i % 2) * 2\n",
    "        col_idx_large = col_idx_small + 1\n",
    "\n",
    "        ax_small = axes[row_idx, col_idx_small]\n",
    "        ax_large = axes[row_idx, col_idx_large]\n",
    "\n",
    "        color = colors.get(dist_name, \"gray\")\n",
    "\n",
    "        # --- Data generation ---\n",
    "        regrets_small, mean_small, p95_small = func(N_small, M, q)\n",
    "        var_small = np.var(regrets_small)\n",
    "\n",
    "        regrets_large, mean_large, p95_large = func(N_large, M, q)\n",
    "        var_large = np.var(regrets_large)\n",
    "\n",
    "        # --- Common bins based on the small-n distribution ---\n",
    "        xlim_max = np.percentile(regrets_small, 95.5)\n",
    "        if xlim_max > 1e-9:\n",
    "            common_bins = np.linspace(0, xlim_max, 51)\n",
    "        else:\n",
    "            common_bins = 50\n",
    "\n",
    "        # --- Histogram for N_small ---\n",
    "        sns.histplot(regrets_small, ax=ax_small, color=color, stat=\"percent\", bins=common_bins)\n",
    "        ax_small.axvline(mean_small, color=\"red\", linestyle=\"--\",  label=rf\"Mean: {mean_small:.2e}\")\n",
    "        ax_small.axvline(p95_small, color=\"black\", linestyle=\":\", label=rf\"95th Pct: {p95_small:.2e}\")\n",
    "\n",
    "        # --- Histogram for N_large ---\n",
    "        sns.histplot(regrets_large, ax=ax_large, color=color, stat=\"percent\", bins=common_bins)\n",
    "        ax_large.axvline(mean_large, color=\"red\", linestyle=\"--\", label=rf\"Mean: {mean_large:.2e}\")\n",
    "        ax_large.axvline(p95_large, color=\"black\", linestyle=\":\", label=rf\"95th Pct: {p95_large:.2e}\")\n",
    "\n",
    "        # --- Axis formatting ---\n",
    "        for ax, N_val, var_val in [(ax_small, N_small, var_small), (ax_large, N_large, var_large),]:\n",
    "            ax.set_title(rf\"{dist_name} ($n={N_val}$)\" + \"\\n\" + rf\"Var: {var_val:.2e}\")\n",
    "            ax.set_xlabel(\"Additive Regret\")\n",
    "\n",
    "            if xlim_max > 1e-9:\n",
    "                ax.set_xlim(0, 1.05 * xlim_max)\n",
    "\n",
    "            if ax in [axes[0, 0], axes[1, 0]]:\n",
    "                ax.set_ylabel(\"Percentage (%)\")\n",
    "            else:\n",
    "                ax.set_ylabel(\"\")\n",
    "\n",
    "            ax.legend()\n",
    "            ax.tick_params(axis=\"both\")\n",
    "\n",
    "    # --- Unified y-axis limits ---\n",
    "    for ax in axes.flatten():\n",
    "        ax.set_ylim(0, 100)\n",
    "\n",
    "    plt.tight_layout(rect=[0, 0.00, 1, 0.99])\n",
    "    plt.savefig(f\"regret_distribution_{q}.png\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfdf95da-defe-41c2-b2a3-175a7d78b9f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate Figure 6 for each value of q\n",
    "N_small = 11\n",
    "N_large = 196\n",
    "\n",
    "for q in Q_VALUES:\n",
    "    plot_fig6(q, N_small, N_large)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "200af661-386d-459e-afcb-ecb936d17fab",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}