demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "767ad8b5-3356-45f0-9bcc-71b19a8cea7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from GapEstimator import BagProcedure, GapProblem\n",
    "import numpy as np\n",
    "from scipy.integrate import quad\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56ed2937-68f0-4bff-951c-82132e93a2cd",
   "metadata": {},
   "source": [
    "# Define the CVaR problem\n",
    "\n",
    "\\begin{equation}\\label{cvar}\n",
    "\\min_{x\\in R}x+\\frac{1}{\\alpha} E[(\\xi-x)_+]\n",
    "\\end{equation}\n",
    "\n",
    "where $\\xi$ is a standard normal variable and $\\alpha=0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3c8a76df-8a29-4406-8eaf-838875378e8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "cvarAlpha = 0.1\n",
    "def cvarSAA(data, x0=None):\n",
    "#     This function solves the SAA problem given the data\n",
    "    \n",
    "#     data: data for uncertain parameters in the optimization problem\n",
    "#     x0: a solution of the optimization problem; default None, solve the SAA with the modified objective h(x,\\xi) - h(x0,\\xi) instead if specified\n",
    "    \n",
    "#     return: a tuple of length 3: SAA solution, SAA objective value, SAA objective variance\n",
    "\n",
    "    n = len(data)\n",
    "    data = np.sort(data)\n",
    "    idx = int(np.ceil(cvarAlpha * n))\n",
    "    x = data[-idx]\n",
    "    mean, var = cvarSAAObj(data, x, x0=x0)\n",
    "    return x, mean, var\n",
    "\n",
    "def cvarSAAObj(data, x, x0=None):\n",
    "#     This function evaluates the SAA objective given the data and an arbitrary feasible solution\n",
    "    \n",
    "#     data: data for uncertain parameters in the optimization problem\n",
    "#     x: an arbitrary feasible solution at which the SAA objective needs to be evaluated\n",
    "#     x0: a solution of the optimization problem; default None, evaluate the SAA objective with the modified objective h(x,\\xi) - h(x0,\\xi) instead if specified\n",
    "    \n",
    "#     return: a tuple of length 2: SAA objective at x, SAA objective variance at x\n",
    "    \n",
    "    SAA_cost = x + np.maximum(data - x, 0) / cvarAlpha\n",
    "    if x0 is not None:\n",
    "        SAA_cost -= x0 + np.maximum(data - x0, 0) / cvarAlpha\n",
    "    return np.mean(SAA_cost), np.var(SAA_cost)\n",
    "\n",
    "\n",
    "def cvarObj(x):\n",
    "    # to compute the exact objective value\n",
    "    def func(t):\n",
    "        return 1 / np.sqrt(2 * np.pi) * np.exp(-t**2 / 2) * (t - x)\n",
    "    return x + quad(func, x, np.inf)[0] / cvarAlpha\n",
    "\n",
    "gapProblem = GapProblem(cvarSAA, cvarSAAObj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd06253e-b35d-4f2c-8843-c50d5f30ffeb",
   "metadata": {},
   "source": [
    "# Define the bagging estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9e49df91-fd43-4447-a771-8bbdaf0127b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# withReplace = False\n",
    "withReplace = True\n",
    "# resample with or without replacement\n",
    "\n",
    "useDebias = True \n",
    "# whether or not debias the infinitesimal Jackknife variance estimator, True recommended\n",
    "\n",
    "includePointVar = True \n",
    "# whether or not include variance of the bagging point estimator in the bound, True recommended\n",
    "\n",
    "resample_rng = np.random.default_rng(seed = 123456)\n",
    "\n",
    "\n",
    "bagEstimator = BagProcedure(replace = withReplace, debias = useDebias, pointVar = includePointVar, rng = resample_rng)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ee835bb-dce1-4cd4-b1b2-71abc107e6a6",
   "metadata": {},
   "source": [
    "# Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6729aba4-1540-41d5-a291-a09e581b3370",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_rng = np.random.default_rng(seed = 654321)\n",
    "\n",
    "N = 200\n",
    "# data size\n",
    "\n",
    "k = 150\n",
    "# resample size\n",
    "\n",
    "B = 500\n",
    "# bootstrap size, >= 500 recommended\n",
    "\n",
    "alpha = 0.05\n",
    "# target confidence level = 1 - alpha\n",
    "\n",
    "xstar = norm.ppf(1 - cvarAlpha)\n",
    "# true optimal solution\n",
    "\n",
    "data = data_rng.normal(size = N)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb404312-142c-4121-85ef-05c3274c8391",
   "metadata": {},
   "source": [
    "# Compute lower bounds of the optimal value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "95df3f78-a46e-4f25-b14d-c41026a43d6a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact optimal value =  1.754983319324868 \n",
      "Bagging bound =  1.6114000509818678 \n",
      "point estimate =  1.7674183854112808 \n",
      "variance estimate =  0.008996980122896124\n"
     ]
    }
   ],
   "source": [
    "bound, point_estimate, var_estimate = bagEstimator.run(gapProblem, data, alpha, k, B) # this returns lower bounds of the optimal value\n",
    "\n",
    "print(\"Exact optimal value = \", cvarObj(xstar), \"\\nBagging bound = \", bound, \"\\npoint estimate = \", point_estimate, \"\\nvariance estimate = \", var_estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "337c47ed-9690-4515-8952-2ff8723270c9",
   "metadata": {},
   "source": [
    "# Compute upper bounds of optimality gap for the solution x = 1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a5792796-36a4-488d-afb6-8db763f37146",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact optimality gap =  0.038084618301178264 \n",
      "Bagging bound =  0.057434597247277994 \n",
      "point estimate =  0.026423804082452065 \n",
      "variance estimate =  0.0003554440388883513\n"
     ]
    }
   ],
   "source": [
    "x0 = 1.5\n",
    "\n",
    "bound2, point_estimate2, var_estimate2 = bagEstimator.run(gapProblem, data, alpha, k, B, x0 = x0) # this returns upper bounds of the optimality gap for x0\n",
    "\n",
    "print(\"Exact optimality gap = \", cvarObj(x0) - cvarObj(xstar), \"\\nBagging bound = \", bound2, \"\\npoint estimate = \", point_estimate2, \"\\nvariance estimate = \", var_estimate2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a643741-9239-41c3-ab67-8006764cc2e4",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}