Model-Transferability-ICML23.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5e7ce30a",
   "metadata": {},
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "267ed4ff",
   "metadata": {},
   "source": [
    "This code accompanies the paper [Model Transferability with Responsive Decision Subjects](https://arxiv.org/pdf/2107.05911.pdf) accepted by ICML 2023.\n",
    "\n",
    "In particular, it generates Figure 4 of the paper about the synthetic experiements using simulated data. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e7a9159",
   "metadata": {},
   "source": [
    "# Import Libraries "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d91bac37",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random as random\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from scipy.special import expit\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fa82767",
   "metadata": {},
   "source": [
    "# Covariate Shift"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e3c40c0",
   "metadata": {},
   "source": [
    "## helper functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "853318d5",
   "metadata": {},
   "source": [
    "#### generate the initial data distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c6ce942d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_init_CS(N, sigma2, sigma3):\n",
    "    '''\n",
    "    generate the initial dataframe for covariate shift\n",
    "    \n",
    "    Input:\n",
    "        - N: total nunber of samples\n",
    "        - sigma2: variance of the noisy term for X2\n",
    "        - sigma3: variance of the noisy term for X3\n",
    "    \n",
    "    N: total number of points\n",
    "    X1 ~ unif(-1, 1)\n",
    "    X2 ~ 1.2X1 + N(0,sigma2)\n",
    "    X3 ~ -0.8X1^2 + N(0,sigma3)\n",
    "    Y is threshold based on the value of X2\n",
    "    \n",
    "    output:\n",
    "        - df: a data frame df with ['X1', 'X2', 'X3', 'Y']\n",
    "    '''\n",
    "    random_seed = random.randint(1, 1000)\n",
    "    np.random.seed(random_seed)\n",
    "\n",
    "    # generate X1, X2, X3\n",
    "    X1 = np.random.uniform(-1, 1, N)\n",
    "    X2 = X1 * 1.2 + np.random.normal(0, sigma2, N)\n",
    "    X3 = -X1**2 * 0.8 + np.random.normal(0, sigma3, N)\n",
    "    # Generate Y\n",
    "    Y = 1 * (X2 > 0)\n",
    "    d = {'X1': list(X1), 'X2': list(X2), 'X3': list(X3), 'Y': list(Y)}\n",
    "    df = pd.DataFrame(data=d)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3be88e9",
   "metadata": {},
   "source": [
    "#### Compute optimal threshold (hS) for the source distribution "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f2e59317",
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper function: compute average square error of a threshold on a dataframe df\n",
    "def compute_error(df, threshold):\n",
    "    '''\n",
    "    Input: \n",
    "        - df: a dataset with a column for 'Q' already computed\n",
    "        - a threshold\n",
    "    \n",
    "    output:\n",
    "        - the mean squared error\n",
    "    '''\n",
    "    error = 0\n",
    "    # compute error\n",
    "    for i in range(len(df)):\n",
    "        # if df_result[\"x\"][i]'s value is greater than the treshold then h(x) = +1; otherwise h(x) = 0\n",
    "        error = error + ((df.iloc[i][\"Q\"] > threshold) - df.iloc[i][\"Y\"])**2\n",
    "    return error / len(df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe43552f",
   "metadata": {},
   "source": [
    "#### Compute the optimal threshold for the original distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9b54ed72",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_optimal_threshold(df, M):\n",
    "    '''\n",
    "    A function to learn the optimal threshold on the original dataset \n",
    "    \n",
    "    Input:\n",
    "        - df: a dataset with a column for 'Q' already computed\n",
    "        - M is the number of potential thresholds.\n",
    "        \n",
    "    Output:\n",
    "        - the optimal threshold for df\n",
    "    '''\n",
    "    # compute the boundary of Q:\n",
    "    Q_min, Q_max = np.min(df[\"Q\"]), np.max(df[\"Q\"])\n",
    "\n",
    "    # loop through potential thresholds:\n",
    "    error_list = []\n",
    "    for i in range(M):\n",
    "        threshold = Q_min + i * (Q_max - Q_min) / M\n",
    "        error_list.append(compute_error(df, threshold))\n",
    "    index_min = np.argmin(error_list)\n",
    "    optimal_threshold = Q_min + index_min * (Q_max - Q_min) / M\n",
    "    return optimal_threshold"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89a0103f",
   "metadata": {},
   "source": [
    "#### Compute induced distribution given a threshold\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "91122b4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the induced distribution\n",
    "def compute_new_distribution_CS(df, sigma2, sigma3, threshold, c):\n",
    "    '''\n",
    "    input:\n",
    "        - df: original dataframe (without 'Q')\n",
    "        - sigma2:\n",
    "        - sigma3:\n",
    "        - threshold: a given threshold\n",
    "        - c: coefficient for adjusting the magnitude of distribution shift\n",
    "    output:\n",
    "        - the new dataframe for holding the induced distribution\n",
    "    '''\n",
    "    N = len(df)\n",
    "    # compute the new data distribution\n",
    "    # get the old classifier result\n",
    "    hX = 1 * (df[\"Q\"] > threshold)\n",
    "    # hX = df[\"Q\"]\n",
    "    # generate new X1\n",
    "    X1_new = df['X1'] + c * (2 * hX - 1)\n",
    "\n",
    "    X2_new = X1_new * 1.2 + np.random.normal(0, sigma2, N)\n",
    "    X3_new = -X1_new**2 * 0.8 + np.random.normal(0, sigma3, N)\n",
    "    # Generate Y\n",
    "    Y_new = 1 * (X2_new > 0)\n",
    "    d = {\n",
    "        'X1': list(X1_new),\n",
    "        'X2': list(X2_new),\n",
    "        'X3': list(X3_new),\n",
    "        'Y': list(Y_new)\n",
    "    }\n",
    "    df_new = pd.DataFrame(data=d)\n",
    "    return df_new"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6f97262",
   "metadata": {},
   "source": [
    "#### Compute optimal classifier achieve minimum induced risk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a0b7725a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_optimal_threshold_MinInduceRisk_CS(df, sigma2, sigma3, M, c):\n",
    "    '''\n",
    "    find the optimal classifier achieve min Err(h)(h)\n",
    "    \n",
    "    input:\n",
    "        -  df: original dataframe (without 'Q')\n",
    "        -  M: number of potential classifiers to search for\n",
    "        -  c: coefficient for adjusting the magnitude of distribution shift\n",
    "        \n",
    "    output: \n",
    "        -  optimal threshold  = argmin Err(h)(h)\n",
    "        -  the minimum error \n",
    "    '''\n",
    "    N = len(df)\n",
    "\n",
    "    # compute the coefficients on the initial data distribution\n",
    "    clf_CS = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "    clf_CS.fit(df[[\"X1\", \"X2\", \"X3\"]], df['Y'])\n",
    "\n",
    "    # compute the boundary of Q:\n",
    "    Q_min, Q_max = np.min(df[\"Q\"]), np.max(df[\"Q\"])\n",
    "\n",
    "    # loop through the potential thresholds\n",
    "    error_list = []\n",
    "    for j in range(M):\n",
    "        threshold = Q_min + j * (Q_max - Q_min) / M\n",
    "        # compute the new data distribution\n",
    "        df_new = compute_new_distribution_CS(df, sigma2, sigma3, threshold, c)\n",
    "        # compute the new qualification\n",
    "        # using the coefficient from the logistic regression coefficient for the old dataset\n",
    "        df_new['Q'] = clf_CS.predict_proba(df_new[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "        # compute new data distribution's loss\n",
    "        Err_h_h = compute_error(df_new, threshold)\n",
    "        error_list.append(Err_h_h)\n",
    "    index_min = np.argmin(error_list)\n",
    "    min_error = np.min(error_list)\n",
    "\n",
    "    # compute the optimal threshold on the dataset that it induced call it hT\n",
    "    optimal_threshold = Q_min + index_min * (Q_max - Q_min) / M\n",
    "    return optimal_threshold, min_error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96ea71a9",
   "metadata": {},
   "source": [
    "#### compute varaince of coefficients between two dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1f24c804",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the variance of w_h for two distribution\n",
    "def var_w(df, df_new, number_bin):\n",
    "    '''\n",
    "    Give two dataframe, compute w(x) = Pr_{df_new}(x)/Pr_{df}(x)\n",
    "    \n",
    "    Input:\n",
    "        - df: one dataframe\n",
    "        - df_new: another dataframe\n",
    "        - number_bin: # of bins\n",
    "    Output:\n",
    "        - w(x)\n",
    "    \n",
    "    '''\n",
    "    bin_init, _ = np.histogram(np.array(df[['X1', 'X2', 'X3']]),\n",
    "                               bins=number_bin)\n",
    "    bin_new, _ = np.histogram(np.array(df_new[['X1', 'X2', 'X3']]),\n",
    "                              bins=number_bin)\n",
    "\n",
    "    return np.var(bin_init / bin_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7836110e",
   "metadata": {},
   "source": [
    "## Upper bounds for covariate shift\n",
    "\n",
    "Covariate Shift Upper Bound is (Theorem 4.2 in Section 4.1):\n",
    "\n",
    "          Err(hS)(hS) - Err(hT)(hT) <= \\sqrt{Err(DS)(hT)} (\\sqrt{Var(w(hS))} + \\sqrt{Var(w(hT))})\n",
    "\n",
    "To characterize the CS upper bound, we need to compute the following quantities:\n",
    "  -  hS (optimal classifier trained on the original dataset)\n",
    "  \n",
    "  -  Err_hs_hs (error of hS on the distribution it induces D(hS)\n",
    "  \n",
    "  -  hT (optimal classifier considering induced data distribution)\n",
    "  \n",
    "  -  Err_hT_hT (error of hT on the distribution it induced D(hT))\n",
    "  \n",
    "  -  Err_S_hT (error of hT on the source distribution S)\n",
    "  \n",
    "  -  sqrt(var(whS))+sqrt(var(whT)): where whS and whT are the varaince of w(hS) and w(hT)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2521dfa1",
   "metadata": {},
   "source": [
    "#### Set Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "67e97ada",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 500\n",
    "c = 0.5\n",
    "M = 100\n",
    "sigma2 = 0.15\n",
    "sigma3 = 0.1\n",
    "\n",
    "number_bin = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "83db7d57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate the initial data distribution\n",
    "df_init_CS = generate_init_CS(N, sigma2, sigma3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c625e9c5",
   "metadata": {},
   "source": [
    "#### Train a logistic regression classifier based on the initial dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5bf2d4b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>LogisticRegression(solver=&#x27;liblinear&#x27;)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">LogisticRegression</label><div class=\"sk-toggleable__content\"><pre>LogisticRegression(solver=&#x27;liblinear&#x27;)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "LogisticRegression(solver='liblinear')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf_CS = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "clf_CS.fit(df_init_CS[[\"X1\", \"X2\", \"X3\"]], df_init_CS['Y'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "204d0caa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the \"qualification\" (linear combination of features) based on the logistic regression model we trained\n",
    "df_init_CS['Q'] = clf_CS.predict_proba(df_init_CS[[\"X1\", \"X2\", \"X3\"]])[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4af93765",
   "metadata": {},
   "source": [
    "#### compute hS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "88420459",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hS: 0.4799997742308634\n",
      "Err_S_hS_CS: 0.016\n"
     ]
    }
   ],
   "source": [
    "# compute the optimal threshold on the original dataset, name it \"hS_CS\"\n",
    "hS_CS = compute_optimal_threshold(df_init_CS, M)\n",
    "print(\"hS:\", hS_CS)\n",
    "\n",
    "# error on the original dataset S using classifier hS_CS\n",
    "Err_S_hS_CS = compute_error(df_init_CS, hS_CS)\n",
    "print(\"Err_S_hS_CS:\", Err_S_hS_CS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "02a954a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Err_hS_hS_CS: 0.0\n"
     ]
    }
   ],
   "source": [
    "# compute the induced data distribution by hS call it df_hS_CS\n",
    "df_hS_CS = compute_new_distribution_CS(df_init_CS, sigma2, sigma3, hS_CS, c)\n",
    "\n",
    "# compute the new qualification\n",
    "df_hS_CS['Q'] = clf_CS.predict_proba(df_hS_CS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "# compute the error on the induced dataset using hS, call it Err_hS_hS_LS\n",
    "Err_hS_hS_CS = compute_error(df_hS_CS, hS_CS)\n",
    "print(\"Err_hS_hS_CS:\", Err_hS_hS_CS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "23b86fa5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hT_CS: 0.1000177771502435\n",
      "Err_hT_hT_CS: 0.0\n"
     ]
    }
   ],
   "source": [
    "# compute the optimal classifier \"hT_CS\" consider induced data and the corresponding induced error \"Err_T_T_CS\"\n",
    "hT_CS, Err_hT_hT_CS = compute_optimal_threshold_MinInduceRisk_CS(\n",
    "    df_init_CS, sigma2, sigma3, M, c)\n",
    "\n",
    "print(\"hT_CS:\", hT_CS)\n",
    "print(\"Err_hT_hT_CS:\", Err_hT_hT_CS)\n",
    "\n",
    "# compute the new data distribution induced by hT\n",
    "df_hT_CS = compute_new_distribution_CS(df_init_CS, sigma2, sigma3, hT_CS, c)\n",
    "df_hT_CS['Q'] = clf_CS.predict_proba(df_hT_CS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c5d57718",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Err_S_hT_CS 0.106\n"
     ]
    }
   ],
   "source": [
    "# compute Err_S_hT\n",
    "Err_S_hT_CS = compute_error(df_init_CS, hT_CS)\n",
    "print(\"Err_S_hT_CS\", Err_S_hT_CS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7e516105",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Y vs Q for D(hT)')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAADgCAYAAADPGumFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAASHklEQVR4nO3de5BcZZ3G8e/DkGBYLlFntCQQghiUoMTLgCwriCKbC+5GWKwieINyN0uVKH/surBrrXhZV9CSAitAKsvGFCqkuG0MawBRC7FEJBOFQMBADBJCVCZclIS45PLbP84Z6DTdMz3z9tuX5PlUTVX3e94555eBfvqc91xeRQRmZin2ancBZtb9HCRmlsxBYmbJHCRmlsxBYmbJHCRmlsxBYmbJHCSWjaS/kvSopM2SPtTC7fZJWiPpVeX7OyX9/RjXdbOkmc2tcPfjIOkykr4raVFV23slPS3pDRm3O1HSVZJ+L+kFSQ9I+sQIv/YlYH5E7BcRS5tQw2JJL0p6vvx5UNJXJR1Y1fVC4FsR8ecR1ndCGXKbJW2RFBXvN0uaDFwMfCW19t2dg6T7fAaYLekUgPJb97+Af4qI3+XYoKTxwA+BQ4G/BA4EPgt8TdJnhvnVQ4HVY9zm3nUWfS0i9gf6gHOA44CfSfqL8vf2AT4BfGekbUTET8uQ2w84qmyeONQWEesj4l7gAEn9Y/l37CkcJF0mIp4GPg0sLD88FwG/iYjF1X0lHVfuQfRUtJ0maVX5+lhJA5L+JOkPki6ts9mPAZOBD0fEYxGxLSJuowi1/5C0f41t/wZ4I3BL+e2+j6SDJC2T9IyktZL+oaL/FyTdKOk7kv4EnD3C3+HPEbEC+FvgtRShAvBu4LmI2FD1K4dK+lm5J/MDSb3Drb/KncCpo+i/x3GQdKGIuAFYCVwHzAP+sU6/e4AtwPsrms8Cri1fXw5cHhEHAIcD19fZ5CnArRGxpar9JmBfir2C6m0fDqwH/qb8dv+/st4NwEHAGcB/Sjq54tfmADcCE4Hv1qmlejvPA3cAJ5RNbwPW1Oh6FkXYvA4YD/xzI+svPQxMH0X/PY6DpHt9iiIgvhQR64fpdx0wF6Dcc5hdtgFsA94kqTciNpfBU0sv8IrDpojYDmyiOMwYlqRDgPcAF5R7E/cBV1Ps7Qz5eUQsjYidEbF1pHVW2Ai8pnw9EXi+Rp9vRcQj5XqvB94+ivU/X67X6nCQdKmI+APFh3ikMYhrgdPLsYPTgV9GxOPlsk8CRwC/lrRC0gfrrGMT8IqB3HIcoxcYbKDkg4Bnyj2IIY8DkyreP9HAemqZBDxTvn4WeMWhFvD7itcvAPuNYv37A8+NqbI9hINkNxcRD1F8YGex62ENEfFoRMyl2N2/BLhxaNCyyg+BWTWW/R3FXs29DZSyEXhN1XjKZODJynIbWM8uJO0HfAD4adm0iiIcm+lI4P4mr3O34iDZM1xLMTB6InDDUKOkj0rqi4idvPyNu6PG73+bYmzjBklTJI2TNAP4JsVZlD+OVEBEPAHcDXxV0qskHU2xR9TQWEi1cvD2XcBSir2Qb5WL7gUmSppU73fH4L3ArU1c327HQbJnuA44CfhxRGyqaJ8JrJa0mWLg9cxa116UA6UfoDj0+AWwFbgNuAz44ijqmAtModg7+R/gooi4Y5T/ln+R9DzFocw1FIPOxw8NBEfEi8Bi4KOjXG9Nko4BtpSnga0O+QlpNlqSxlF8Qz8JnB0d9j+RpD6KQ513jHLQtta6bgL+OyKWN6W43ZSDxMakvJr0fOD6iPh1u+ux9nKQmFkyj5GYWTIHiZklq3djVMfq7e2NKVOmtLsMsz3OypUrN0VEzauYuy5IpkyZwsDAQLvLMNvjSHq83jIf2phZMgeJmSXLdmhTPsXrg8BTEfHWGstFcTXlbIqbqM6OiF/mqseslqW/epKv376Gjc9t5aCJE/jsjDfzoXdMGnWf0fQbrj/AF29ZzbMvbANAFDcg9UjsiGDfcXuxdftO6l21MXHCOJ7bum1MfwuA3148tseuZLuORNKJwGbgmjpBMpviAT2zKR5Gc3lEvHuk9fb394fHSKwZlv7qSf715gfYuu3l24smjOvhq6e/7aUAaKTPaPoN13/cXmInsGNne6/tqhcmklZGRM0nxWU7tImIu3j51u5a5lCETJTPwZiY85mjZtW+fvuaXT7IAFu37eDrt68ZVZ/R9Buu/7ad0fYQGat2jpFMYtfnT2xg12dTvETSvPKRgAODg408+sJsZBufq30bTmV7I31G02+k9m7VziBRjbaacRwRCyOiPyL6+/pGfBiXWUMOmjhhxPZG+oym30jt3aqdQbIBOKTi/cEUt5ebtcRnZ7yZCeN6dmmbMK7npUHPRvuMpt9w/cftJXr2qvX92vnaeUHaMuA8SUsoBlv/mGs6BbNahgZBhzvT0kif0fQbqT/4rM2uK5aGHqbTC/yBYtqEcQARsaA8/Tuf4uE6LwDnRMSIp2N81sasPYY7a5Ntj6R8Fuhwy4PiSehm1uV8ZauZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFkyB4mZJXOQmFmyrEEiaaakNZLWSrqwxvIDJd0i6X5JqyWdk7MeM8sjW5BI6gGuAGYB04C5kqZVdfsU8FBETKeYuuIbksbnqsnM8si5R3IssDYi1kXEi8ASionDKwWwfznHzX4Uk45vz1iTmWWQM0gamSR8PnAkxVSdDwDnR8TOjDWZWQY5g6SRScJnAPcBBwFvB+ZLOuAVK5LmSRqQNDA4ONjsOs0sUc4gaWSS8HOAm6OwFngMeEv1iiJiYUT0R0R/X19ftoLNbGxyBskKYKqkw8oB1DMpJg6vtB44GUDS64E3A+sy1mRmGeSc+3e7pPOA24EeYFFErJZ0brl8AfBlYLGkBygOhS6IiE25ajKzPLIFCUBELAeWV7UtqHi9EfjrnDWYWX6+stXMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkmUNEkkzJa2RtFbShXX6nCTpPkmrJf0kZz1mlke26Sgk9QBXAKdQzLq3QtKyiHioos9E4EpgZkSsl/S6XPWYWT4590iOBdZGxLqIeBFYAsyp6nMWxZSd6wEi4qmM9ZhZJjmDZBLwRMX7DWVbpSOAV0u6U9JKSR+vtSJPIm7W2XIGiWq0RdX7vYF3AacCM4B/l3TEK37Jk4ibdbScU3ZuAA6peH8wsLFGn00RsQXYIukuYDrwSMa6zKzJcu6RrACmSjpM0njgTGBZVZ/vASdI2lvSvsC7gYcz1mRmGWTbI4mI7ZLOA24HeoBFEbFa0rnl8gUR8bCk24BVwE7g6oh4MFdNZpaHIqqHLTpbf39/DAwMtLsMsz2OpJUR0V9rma9sNbNkDhIzS+YgMbNkDhIzS+YgMbNkDhIzS+YgMbNkdYNE0nJJU1pYi5l1qeH2SBYDP5D0OUnjWlSPmXWhupfIR8T1kr4PfB4YkPRtisvYh5Zf2oL6zKwLjHSvzTZgC7APsD8VQWJmNqRukEiaCVxKccfuOyPihZZVZWZdZbg9ks8BH46I1a0qxsy603BjJCe0shAz616+jsTMkjlIzCyZg8TMkjlIzCyZg8TMkjlIzCyZg8TMkmUNEkkzJa2RtFbShcP0O0bSDkln5KzHzPLIFiSSeoArgFnANGCupGl1+l1CMf+NmXWhnHskxwJrI2JdRLwILAHm1Oj3aeAm4KmMtZhZRjmDZBLwRMX7DWXbSyRNAk4DFgy3IknzJA1IGhgcHGx6oWaWJmeQqEZb9bR+lwEXRMSO4VYUEQsjoj8i+vv6+ppVn5k1Sba5fyn2QA6peH8wsLGqTz+wRBJALzBb0vaIWJqxLjNrspxBsgKYKukw4EngTOCsyg4RcdjQa0mLgf91iJh1n2xBEhHbJZ1HcTamB1gUEaslnVsuH3ZcxMy6R849EiJiObC8qq1mgETE2TlrMbN8fGWrmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZMgeJmSVzkJhZsrZOIi7pI5JWlT93S5qesx4zy6Pdk4g/Brw3Io4GvgwszFWPmeXT1knEI+LuiHi2fHsPxWx8ZtZl2jqJeJVPArdmrMfMMsk5QVYjk4gXHaX3UQTJe+osnwfMA5g8eXKz6jOzJsm5R9LIJOJIOhq4GpgTEU/XWlFELIyI/ojo7+vry1KsmY1dziB5aRJxSeMpJhFfVtlB0mTgZuBjEfFIxlrMLKN2TyL+eeC1wJWSALZHRH+umswsD0XUHLboWP39/TEwMNDuMsz2OJJW1vui95WtZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpbMQWJmyRwkZpYs50x7SJoJXE4xHcXVEXFx1XKVy2cDLwBnR8QvU7Y55cLv12z/7cWnpqzWzIaRbY9EUg9wBTALmAbMlTStqtssYGr5Mw+4KmWb9UJkpGVmlibnoc2xwNqIWBcRLwJLgDlVfeYA10ThHmCipDdkrMnMMsgZJJOAJyrebyjbRtsHSfMkDUgaGBwcbHqhZpYmZ5CoRlv1tH6N9PEk4mYdLmeQbAAOqXh/MLBxDH3MrMPlDJIVwFRJh0kaD5wJLKvqswz4uArHAX+MiN+NdYPDnZnxWRuzfLKd/o2I7ZLOA26nOP27KCJWSzq3XL4AWE5x6nctxenfc1K368Awa72s15FExHKKsKhsW1DxOoBP5azBzPLzla1mlsxBYmbJVBxddA9Jg8DjDXTtBTZlLieVa0zX6fVB59fYaH2HRkTN6y+6LkgaJWkgIvrbXcdwXGO6Tq8POr/GZtTnQxszS+YgMbNku3OQLGx3AQ1wjek6vT7o/BqT69ttx0jMrHV25z0SM2uRrg8SSTMlrZG0VtKFNZZL0jfL5askvbMDa/xIWdsqSXdLmt5J9VX0O0bSDklntLK+ctsj1ijpJEn3SVot6SedVJ+kAyXdIun+sr7k20FGWd8iSU9JerDO8rTPSUR07Q/FPTy/Ad4IjAfuB6ZV9ZkN3ErxyILjgF90YI3HA68uX89qZY2N1FfR78cUtzyc0YF/w4nAQ8Dk8v3rOqy+fwMuKV/3Ac8A41tY44nAO4EH6yxP+px0+x5JNzyFbcQaI+LuiHi2fHsPxeMUOqa+0qeBm4CnWljbkEZqPAu4OSLWA0REK+tspL4A9i+fU7wfRZBsb1WBEXFXuc16kj4n3R4kTXsKW0aj3f4nKb4ZWmXE+iRNAk4DFtAejfwNjwBeLelOSSslfbxl1TVW33zgSIrn7TwAnB8RO1tTXkOSPidZ7/5tgaY9hS2jhrcv6X0UQfKerBVVbbZGW3V9lwEXRMSO4gu15RqpcW/gXcDJwATg55LuiYhHchdHY/XNAO4D3g8cDtwh6acR8afMtTUq6XPS7UHSDU9ha2j7ko4GrgZmRcTTLaoNGquvH1hShkgvMFvS9ohY2pIKG//vvCkitgBbJN0FTAdaESSN1HcOcHEUAxJrJT0GvAW4twX1NSLtc9KqwZ5MA0h7A+uAw3h5kOuoqj6nsusg0r0dWONkioc7Hd+Jf8Oq/otp/WBrI3/DI4EflX33BR4E3tpB9V0FfKF8/XrgSaC3xX/HKdQfbE36nHT1Hkm06SlsGWr8PPBa4MryW397tOgmrwbra6tGaoyIhyXdBqwCdlJMyFbzVGc76gO+DCyW9ADFh/WCiGjZHcGSrgNOAnolbQAuAsZV1Jf0OfGVrWaWrNvP2phZB3CQmFkyB4mZJXOQmFkyB4mZJXOQWHaSDpb0PUmPSlonab6kfdpdlzWPg8SyKm9SuxlYGhFTgakUl7B/ra2FWVP5OhLLStLJwEURcWJF2wEUU4ocEhGb21acNY33SCy3o4CVlQ1R3Kj2W+BN7SjIms9BYrmJ2neRtuU2YsvDQWK5raa4e/gl5aHN64E1banIms5BYrn9CNh36EFDknqAbwDzI2JrWyuzpnGQWFZRjOafBpwh6VHgaWBnRHylvZVZM/msjbWUpOOB64DTI2LlSP2tOzhIzCyZD23MLJmDxMySOUjMLJmDxMySOUjMLJmDxMySOUjMLNn/A3ketNCCUE8bAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(1, figsize=(4, 3))\n",
    "plt.clf()\n",
    "plt.scatter(df_hT_CS[\"Q\"], df_hT_CS[\"Y\"])\n",
    "plt.ylabel(\"Y\")\n",
    "plt.xlabel(\"Q\")\n",
    "plt.title(\"Y vs Q for D(hT)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9ded9535",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "812.0903515802157\n"
     ]
    }
   ],
   "source": [
    "var_whs_CS = var_w(df_init_CS, df_hS_CS, number_bin=number_bin)\n",
    "print(var_whs_CS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "51db478f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.701092303189744\n"
     ]
    }
   ],
   "source": [
    "var_whT_CS = var_w(df_init_CS, df_hT_CS, number_bin=number_bin)\n",
    "print(var_whT_CS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b3c2505",
   "metadata": {},
   "source": [
    "### Compute the upper bound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c6f7ba69",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.0\n",
      "the upper bound is: 9.813100642705699\n"
     ]
    }
   ],
   "source": [
    "# compute Err_hs_hs - Err_hT_hT\n",
    "print(\"Err_hs_hs_CS - Err_hT_hT_CS:\", Err_hS_hS_CS - Err_hT_hT_CS)\n",
    "\n",
    "# compute Err_S_hT\n",
    "Err_S_hT_CS = compute_error(df_init_CS, hT_CS)\n",
    "#print(\"Err_S_hT_CS\", Err_S_hT_CS)\n",
    "\n",
    "# compute the upper bound (theorem 5)\n",
    "UB = np.sqrt(Err_S_hT_CS) * (np.sqrt(var_whs_CS) + np.sqrt(var_whT_CS))\n",
    "print(\"the upper bound is:\", UB)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce341c16",
   "metadata": {},
   "source": [
    "## lower bound for covariate shift\n",
    "\n",
    "Covariate Shift lower bound is (Theorem 4.6 in Section 4.2):\n",
    "\n",
    "                max{Err_S(h), Err_h(h)} >= (dTV(D_S(Y), D_h(Y)) - dTV(D_S(h), D_h(h)))/2\n",
    "\n",
    "\n",
    "To test the lower bound, we need to compute the following quantities for hS and hT:\n",
    "\n",
    "   - max{Err_S(h), Err_h(h)} \n",
    "   \n",
    "        \n",
    "   - dTV(D_S(Y), D_h(Y)) = |Pr_S(Y = 1) - Pr_h(Y = 1)|\n",
    "   \n",
    "    \n",
    "   - dTV(D_S(h), D_h(h)) = |Pr_S(h = 1) - Pr_h(h = 1)|"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dcc6f57",
   "metadata": {},
   "source": [
    "#### compute the lower bound for hT (theorem 4.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "a48aa505",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max{Err_S(hT), Err_hT(hT)} 0.106\n",
      "the lower bound for using hT is: 0.05299999999999999\n"
     ]
    }
   ],
   "source": [
    "# compute Pr_Ds[Y = +1]\n",
    "Pr_S_Y1_CS = df_init_CS['Y'].value_counts()[1] / len(df_init_CS)\n",
    "\n",
    "#print(\"Pr_S_Y1_CS\", Pr_S_Y1_CS)\n",
    "\n",
    "# compute Pr_DhT[Y = +1]\n",
    "Pr_hT_Y1_CS = df_hT_CS['Y'].value_counts()[1] / len(df_hT_CS)\n",
    "\n",
    "#print(\"Pr_Dh_Y1_CS\", Pr_Dh_Y1_CS)\n",
    "\n",
    "# compute Pr_Ds[hT(X) = +1]\n",
    "Pr_S_hT1_CS = np.sum(df_init_CS['Q'] > hT_CS) / len(df_init_CS)\n",
    "\n",
    "# compute Pr_DhT[hT(X) = +1]\n",
    "Pr_hT_hT1_CS = np.sum(df_hT_CS['Q'] > hT_CS) / len(df_hT_CS)\n",
    "\n",
    "#print(\"Pr_hT_hT1_CS\", Pr_hT_hT1_CS)\n",
    "\n",
    "print(\"max{Err_S(hT), Err_hT(hT)}\", max(Err_S_hT_CS, Err_hT_hT_CS))\n",
    "\n",
    "LB_hT_CS = (np.abs(Pr_S_Y1_CS - Pr_hT_Y1_CS) -\n",
    "            np.abs(Pr_S_hT1_CS - Pr_hT_hT1_CS)) / 2\n",
    "print(\"the lower bound for using hT is:\", LB_hT_CS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfd32729",
   "metadata": {},
   "source": [
    "#### compute the lower bound for hS (theorem 4.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a18a2b78",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max{Err_S(hS), Err_hS(hS)} 0.016\n",
      "the lower bound for using hS is: 0.0040000000000000036\n"
     ]
    }
   ],
   "source": [
    "# compute Pr_DS[Y = +1]\n",
    "Pr_S_Y1_CS = df_init_CS['Y'].value_counts()[1] / len(df_init_CS)\n",
    "\n",
    "# compute Pr_DhS[Y = +1]\n",
    "Pr_hS_Y1_CS = df_hS_CS['Y'].value_counts()[1] / len(df_hS_CS)\n",
    "\n",
    "# compute Pr_DS[hS(X) = +1]\n",
    "Pr_S_hS1_CS = np.sum(df_init_CS['Q'] > hS_CS) / len(df_init_CS)\n",
    "\n",
    "# compute Pr_DhS[hS(X) = +1]\n",
    "Pr_hS_hS1_CS = np.sum(df_hS_CS['Q'] > hS_CS) / len(df_hS_CS)\n",
    "\n",
    "print(\"max{Err_S(hS), Err_hS(hS)}\", max(Err_S_hS_CS, Err_hS_hS_CS))\n",
    "\n",
    "LB_hS_CS = (np.abs(Pr_S_Y1_CS - Pr_hS_Y1_CS) -\n",
    "            np.abs(Pr_S_hS1_CS - Pr_hS_hS1_CS)) / 2\n",
    "print(\"the lower bound for using hS is:\", LB_hS_CS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f67ef82f",
   "metadata": {},
   "source": [
    "\n",
    "### A helper function for computing upper and lower bound for one particular dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "4e612479",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_upper_lower_bound_CS(N=1000,\n",
    "                                 c=0.25,\n",
    "                                 M=100,\n",
    "                                 sigma2=0.15,\n",
    "                                 sigma3=0.1,\n",
    "                                 number_bin=20):\n",
    "    '''\n",
    "    A helper function that takes in all the parameters and compute the upper bounds\n",
    "    '''\n",
    "\n",
    "    # generate the initial data distribution\n",
    "    df_init_CS = generate_init_CS(N, sigma2, sigma3)\n",
    "\n",
    "    # #### Train a logistic regression classifier based on the initial dataframe\n",
    "    clf_CS = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "    clf_CS.fit(df_init_CS[[\"X1\", \"X2\", \"X3\"]], df_init_CS['Y'])\n",
    "\n",
    "    # compute the \"qualification\" (linear combination of features) based on the logistic regression model we trained\n",
    "    df_init_CS['Q'] = clf_CS.predict_proba(df_init_CS[[\"X1\", \"X2\", \"X3\"]])[:,\n",
    "                                                                           1]\n",
    "\n",
    "    #compute hS\n",
    "\n",
    "    # compute the optimal threshold on the original dataset, name it \"hS_CS\"\n",
    "    hS_CS = compute_optimal_threshold(df_init_CS, M)\n",
    "    #     print(\"hS:\", hS_CS)\n",
    "\n",
    "    # error on the original dataset S using classifier hS_CS\n",
    "    Err_S_hS_CS = compute_error(df_init_CS, hS_CS)\n",
    "    #     print(\"Err_S_hS_CS:\", Err_S_hS_CS)\n",
    "\n",
    "    # compute the induced data distribution by hS call it df_hS_CS\n",
    "    df_hS_CS = compute_new_distribution_CS(df_init_CS, sigma2, sigma3, hS_CS,\n",
    "                                           c)\n",
    "\n",
    "    # compute the new qualification\n",
    "    df_hS_CS['Q'] = clf_CS.predict_proba(df_hS_CS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "    # compute the error on the induced dataset using hS, call it Err_hS_hS_LS\n",
    "    Err_hS_hS_CS = compute_error(df_hS_CS, hS_CS)\n",
    "    #     print(\"Err_hS_hS_CS:\", Err_hS_hS_CS)\n",
    "\n",
    "    # compute the optimal classifier \"hT_CS\" consider induced data and the corresponding induced error \"Err_T_T_CS\"\n",
    "    hT_CS, Err_hT_hT_CS = compute_optimal_threshold_MinInduceRisk_CS(\n",
    "        df_init_CS, sigma2, sigma3, M, c)\n",
    "\n",
    "    #     print(\"hT_CS:\", hT_CS)\n",
    "    #     print(\"Err_hT_hT_CS:\",Err_hT_hT_CS)\n",
    "\n",
    "    # compute the new data distribution induced by hT\n",
    "    df_hT_CS = compute_new_distribution_CS(df_init_CS, sigma2, sigma3, hT_CS,\n",
    "                                           c)\n",
    "    df_hT_CS['Q'] = clf_CS.predict_proba(df_hT_CS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    # compute Err_S_hT\n",
    "    Err_S_hT_CS = compute_error(df_init_CS, hT_CS)\n",
    "    #     print(\"Err_S_hT_CS\", Err_S_hT_CS)\n",
    "\n",
    "    var_whs_CS = var_w(df_init_CS, df_hS_CS, number_bin=number_bin)\n",
    "    #     print(var_whs_CS)\n",
    "\n",
    "    var_whT_CS = var_w(df_init_CS, df_hT_CS, number_bin=number_bin)\n",
    "    #     print(var_whT_CS)\n",
    "\n",
    "    # compute Err_hs_hs - Err_hT_hT\n",
    "    print(\"Err_hs_hs_CS - Err_hT_hT_CS:\", Err_hS_hS_CS - Err_hT_hT_CS)\n",
    "\n",
    "    # compute Err_S_hT\n",
    "    Err_S_hT_CS = compute_error(df_init_CS, hT_CS)\n",
    "    #print(\"Err_S_hT_CS\", Err_S_hT_CS)\n",
    "\n",
    "    # compute the upper bound (theorem 5)\n",
    "    UB = np.sqrt(Err_S_hT_CS) * (np.sqrt(var_whs_CS) + np.sqrt(var_whT_CS))\n",
    "    print(\"the upper bound is:\", UB)\n",
    "\n",
    "    # compute Pr_Ds[Y = +1]\n",
    "    Pr_S_Y1_CS = df_init_CS['Y'].value_counts()[1] / len(df_init_CS)\n",
    "\n",
    "    #print(\"Pr_S_Y1_CS\", Pr_S_Y1_CS)\n",
    "\n",
    "    # compute Pr_DhT[Y = +1]\n",
    "    Pr_hT_Y1_CS = df_hT_CS['Y'].value_counts()[1] / len(df_hT_CS)\n",
    "\n",
    "    #print(\"Pr_Dh_Y1_CS\", Pr_Dh_Y1_CS)\n",
    "\n",
    "    # compute Pr_Ds[hT(X) = +1]\n",
    "    Pr_S_hT1_CS = np.sum(df_init_CS['Q'] > hT_CS) / len(df_init_CS)\n",
    "\n",
    "    # compute Pr_DhT[hT(X) = +1]\n",
    "    Pr_hT_hT1_CS = np.sum(df_hT_CS['Q'] > hT_CS) / len(df_hT_CS)\n",
    "\n",
    "    #print(\"Pr_hT_hT1_CS\", Pr_hT_hT1_CS)\n",
    "\n",
    "    print(\"max{Err_S(hT), Err_hT(hT)}\", max(Err_S_hT_CS, Err_hT_hT_CS))\n",
    "\n",
    "    LB_hT_CS = (np.abs(Pr_S_Y1_CS - Pr_hT_Y1_CS) -\n",
    "                np.abs(Pr_S_hT1_CS - Pr_hT_hT1_CS)) / 2\n",
    "    print(\"the lower bound for using hT is:\", LB_hT_CS)\n",
    "\n",
    "    # compute Pr_DS[Y = +1]\n",
    "    Pr_S_Y1_CS = df_init_CS['Y'].value_counts()[1] / len(df_init_CS)\n",
    "\n",
    "    # compute Pr_DhS[Y = +1]\n",
    "    Pr_hS_Y1_CS = df_hS_CS['Y'].value_counts()[1] / len(df_hS_CS)\n",
    "\n",
    "    # compute Pr_DS[hS(X) = +1]\n",
    "    Pr_S_hS1_CS = np.sum(df_init_CS['Q'] > hS_CS) / len(df_init_CS)\n",
    "\n",
    "    # compute Pr_DhS[hS(X) = +1]\n",
    "    Pr_hS_hS1_CS = np.sum(df_hS_CS['Q'] > hS_CS) / len(df_hS_CS)\n",
    "\n",
    "    print(\"max{Err_S(hS), Err_hS(hS)}\", max(Err_S_hS_CS, Err_hS_hS_CS))\n",
    "\n",
    "    LB_hS_CS = (np.abs(Pr_S_Y1_CS - Pr_hS_Y1_CS) -\n",
    "                np.abs(Pr_S_hS1_CS - Pr_hS_hS1_CS)) / 2\n",
    "    print(\"the lower bound for using hS is:\", LB_hS_CS)\n",
    "\n",
    "    return Err_hS_hS_CS - Err_hT_hT_CS, UB, max(Err_S_hT_CS,\n",
    "                                                Err_hT_hT_CS), LB_hT_CS, max(\n",
    "                                                    Err_S_hS_CS,\n",
    "                                                    Err_hS_hS_CS), LB_hS_CS\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7ec83d2",
   "metadata": {},
   "source": [
    "#### test the helper function with a particular setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "79e9648a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.006\n",
      "the upper bound is: 0.47522785712216153\n",
      "max{Err_S(hT), Err_hT(hT)} 0.046\n",
      "the lower bound for using hT is: 0.01999999999999999\n",
      "max{Err_S(hS), Err_hS(hS)} 0.012\n",
      "the lower bound for using hS is: 0.007000000000000006\n"
     ]
    }
   ],
   "source": [
    "_, _, _, _, _, _ = compute_upper_lower_bound_CS(N=500,\n",
    "                                                c=0.25,\n",
    "                                                M=100,\n",
    "                                                sigma2=0.15,\n",
    "                                                sigma3=0.1,\n",
    "                                                number_bin=20)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2abe7c66",
   "metadata": {},
   "source": [
    "### Helper function to generate the average actual and theoretical upper and lower bound for covariate shift setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "34313c7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def average_actual_theoretical_bounds_CS(num_datasets, N, c, M, sigma2, sigma3,\n",
    "                                         number_bin):\n",
    "\n",
    "    actual_UB_list_CS = []\n",
    "    UB_list_CS = []\n",
    "\n",
    "    actual_LB_hS_list_CS = []\n",
    "    LB_hS_list_CS = []\n",
    "\n",
    "    actual_LB_hT_list_CS = []\n",
    "    LB_hT_list_CS = []\n",
    "\n",
    "    for i in range(num_datasets):\n",
    "        print(i)\n",
    "        actual_UB, UB, actual_LB_hS, LB_hS, actual_LB_hT, LB_hT = compute_upper_lower_bound_CS(\n",
    "            N, c, M, sigma2, sigma3, number_bin)\n",
    "\n",
    "        actual_UB_list_CS.append(actual_UB)\n",
    "        UB_list_CS.append(UB)\n",
    "\n",
    "        actual_LB_hS_list_CS.append(actual_LB_hS)\n",
    "        LB_hS_list_CS.append(LB_hS)\n",
    "\n",
    "        actual_LB_hT_list_CS.append(actual_LB_hT)\n",
    "        LB_hT_list_CS.append(LB_hT)\n",
    "\n",
    "    # compute the mean and error bar for the actual upper bound\n",
    "    actual_UB_CS_mean = np.mean(actual_UB_list_CS)\n",
    "\n",
    "    actual_UB_CS_errormin = np.abs(\n",
    "        np.min(actual_UB_list_CS) - actual_UB_CS_mean)\n",
    "    actual_UB_CS_errormax = np.max(actual_UB_list_CS) - actual_UB_CS_mean\n",
    "\n",
    "    actual_UB_CS_error = [actual_UB_CS_errormin, actual_UB_CS_errormax]\n",
    "\n",
    "    # compute the mean and error bar for the theoretical upper bound\n",
    "\n",
    "    UB_CS_mean = np.mean(UB_list_CS)\n",
    "\n",
    "    UB_CS_errormin = np.abs(np.min(UB_list_CS) - UB_CS_mean)\n",
    "    UB_CS_errormax = np.max(UB_list_CS) - UB_CS_mean\n",
    "\n",
    "    UB_CS_error = [UB_CS_errormin, UB_CS_errormax]\n",
    "\n",
    "    print(\"actual_UB_CS_mean:\", actual_UB_CS_mean)\n",
    "    print(\"UB_CS_mean:\", UB_CS_mean)\n",
    "\n",
    "    print(\"actual_UB_CS_error:\", actual_UB_CS_error)\n",
    "    print(\"UB_CS_error:\", UB_CS_error)\n",
    "\n",
    "    # for hS\n",
    "    actual_LB_hS_CS_mean = np.mean(actual_LB_hS_list_CS)\n",
    "    LB_hS_CS_mean = np.mean(LB_hS_list_CS)\n",
    "\n",
    "    actual_LB_hS_CS_errormin = np.abs(\n",
    "        np.min(actual_LB_hS_list_CS) - LB_hS_CS_mean)\n",
    "    actual_LB_hS_CS_errormax = np.max(actual_LB_hS_list_CS) - LB_hS_CS_mean\n",
    "\n",
    "    actual_LB_hS_CS_error = [\n",
    "        actual_LB_hS_CS_errormin, actual_LB_hS_CS_errormax\n",
    "    ]\n",
    "\n",
    "    LB_CS_hS_errormin = np.abs(np.min(LB_hS_list_CS) - LB_hS_CS_mean)\n",
    "    LB_CS_hS_errormax = np.max(LB_hS_list_CS) - LB_hS_CS_mean\n",
    "\n",
    "    LB_hS_CS_error = [LB_CS_hS_errormin, LB_CS_hS_errormax]\n",
    "\n",
    "    # for hT\n",
    "\n",
    "    actual_LB_hT_CS_mean = np.mean(actual_LB_hT_list_CS)\n",
    "    LB_hT_CS_mean = np.mean(LB_hT_list_CS)\n",
    "\n",
    "    actual_LB_hT_CS_errormin = np.abs(\n",
    "        np.min(actual_LB_hT_list_CS) - LB_hT_CS_mean)\n",
    "    actual_LB_hT_CS_errormax = np.max(actual_LB_hT_list_CS) - LB_hT_CS_mean\n",
    "\n",
    "    actual_LB_hT_CS_error = [\n",
    "        actual_LB_hT_CS_errormin, actual_LB_hT_CS_errormax\n",
    "    ]\n",
    "\n",
    "    LB_CS_hT_errormin = np.abs(np.min(LB_hT_list_CS) - LB_hT_CS_mean)\n",
    "    LB_CS_hT_errormax = np.max(LB_hT_list_CS) - LB_hT_CS_mean\n",
    "\n",
    "    LB_hT_CS_error = [LB_CS_hT_errormin, LB_CS_hT_errormax]\n",
    "\n",
    "    print(\"actual_LB_hS_CS_error\", actual_LB_hS_CS_error)\n",
    "    print(\"LB_hS_CS_error\", LB_hS_CS_error)\n",
    "\n",
    "    print(\"actual_LB_hT_CS_error\", actual_LB_hT_CS_error)\n",
    "    print(\"LB_hT_CS_error\", LB_hT_CS_error)\n",
    "\n",
    "    return actual_UB_CS_mean, actual_UB_CS_error, UB_CS_mean, UB_CS_error, actual_LB_hT_CS_mean, actual_LB_hT_CS_error, LB_hT_CS_mean, LB_hT_CS_error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab23b430",
   "metadata": {},
   "source": [
    "#### generate the bounds for two particular datasets (CS1, CS2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2856f48e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.0\n",
      "the upper bound is: 0.41458669896267475\n",
      "max{Err_S(hT), Err_hT(hT)} 0.03\n",
      "the lower bound for using hT is: 0.01200000000000001\n",
      "max{Err_S(hS), Err_hS(hS)} 0.018\n",
      "the lower bound for using hS is: 0.0010000000000000009\n",
      "1\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.002\n",
      "the upper bound is: 0.32979601718672164\n",
      "max{Err_S(hT), Err_hT(hT)} 0.032\n",
      "the lower bound for using hT is: 0.014000000000000012\n",
      "max{Err_S(hS), Err_hS(hS)} 0.004\n",
      "the lower bound for using hS is: 0.0010000000000000009\n",
      "2\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.004\n",
      "the upper bound is: 0.29035544218616854\n",
      "max{Err_S(hT), Err_hT(hT)} 0.036\n",
      "the lower bound for using hT is: 0.016000000000000014\n",
      "max{Err_S(hS), Err_hS(hS)} 0.008\n",
      "the lower bound for using hS is: 0.0020000000000000018\n",
      "3\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.0\n",
      "the upper bound is: 0.2200196741280874\n",
      "max{Err_S(hT), Err_hT(hT)} 0.034\n",
      "the lower bound for using hT is: 0.012999999999999956\n",
      "max{Err_S(hS), Err_hS(hS)} 0.006\n",
      "the lower bound for using hS is: 0.0010000000000000009\n",
      "4\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.0\n",
      "the upper bound is: 0.3416005071935389\n",
      "max{Err_S(hT), Err_hT(hT)} 0.022\n",
      "the lower bound for using hT is: 0.0040000000000000036\n",
      "max{Err_S(hS), Err_hS(hS)} 0.01\n",
      "the lower bound for using hS is: 0.0010000000000000009\n",
      "actual_UB_CS_mean: 0.0012000000000000001\n",
      "UB_CS_mean: 0.3192716679314382\n",
      "actual_UB_CS_error: [0.0012000000000000001, 0.0028]\n",
      "UB_CS_error: [0.09925199380335079, 0.09531503103123656]\n",
      "actual_LB_hS_CS_error [0.010199999999999999, 0.0242]\n",
      "LB_hS_CS_error [0.007799999999999996, 0.0042000000000000145]\n",
      "actual_LB_hT_CS_error [0.002799999999999999, 0.0168]\n",
      "LB_hT_CS_error [0.0002000000000000001, 0.0008000000000000008]\n"
     ]
    }
   ],
   "source": [
    "actual_UB_CS_mean, actual_UB_CS_error, UB_CS_mean, UB_CS_error, actual_LB_hT_CS_mean, actual_LB_hT_CS_error, LB_hT_CS_mean, LB_hT_CS_error = average_actual_theoretical_bounds_CS(\n",
    "    num_datasets=5,\n",
    "    N=500,\n",
    "    c=0.25,\n",
    "    M=100,\n",
    "    sigma2=0.15,\n",
    "    sigma3=0.1,\n",
    "    number_bin=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3b7e1325",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.0025\n",
      "the upper bound is: 0.5852547642534905\n",
      "max{Err_S(hT), Err_hT(hT)} 0.02125\n",
      "the lower bound for using hT is: 0.005625000000000047\n",
      "max{Err_S(hS), Err_hS(hS)} 0.015\n",
      "the lower bound for using hS is: 0.0012499999999999734\n",
      "1\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.0075\n",
      "the upper bound is: 0.6650771303393215\n",
      "max{Err_S(hT), Err_hT(hT)} 0.02\n",
      "the lower bound for using hT is: -0.0006250000000000422\n",
      "max{Err_S(hS), Err_hS(hS)} 0.01625\n",
      "the lower bound for using hS is: 0.0031250000000000444\n",
      "2\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.005\n",
      "the upper bound is: 0.48555361815772363\n",
      "max{Err_S(hT), Err_hT(hT)} 0.02375\n",
      "the lower bound for using hT is: 0.007499999999999951\n",
      "max{Err_S(hS), Err_hS(hS)} 0.01375\n",
      "the lower bound for using hS is: 0.0006250000000000422\n",
      "3\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.00625\n",
      "the upper bound is: 0.5428620577430808\n",
      "max{Err_S(hT), Err_hT(hT)} 0.04\n",
      "the lower bound for using hT is: 0.018750000000000044\n",
      "max{Err_S(hS), Err_hS(hS)} 0.015\n",
      "the lower bound for using hS is: 0.004375000000000018\n",
      "4\n",
      "Err_hs_hs_CS - Err_hT_hT_CS: 0.0025\n",
      "the upper bound is: 0.7546869666018715\n",
      "max{Err_S(hT), Err_hT(hT)} 0.03125\n",
      "the lower bound for using hT is: 0.013749999999999984\n",
      "max{Err_S(hS), Err_hS(hS)} 0.0175\n",
      "the lower bound for using hS is: 0.0050000000000000044\n",
      "actual_UB_CS_mean: 0.004749999999999999\n",
      "UB_CS_mean: 0.6066869074190976\n",
      "actual_UB_CS_error: [0.002249999999999999, 0.0027500000000000007]\n",
      "UB_CS_error: [0.12113328926137396, 0.14800005918277392]\n",
      "actual_LB_hS_CS_error [0.011000000000000003, 0.031000000000000003]\n",
      "LB_hS_CS_error [0.00962500000000004, 0.009750000000000047]\n",
      "actual_LB_hT_CS_error [0.010874999999999984, 0.014624999999999985]\n",
      "LB_hT_CS_error [0.0022499999999999742, 0.002124999999999988]\n"
     ]
    }
   ],
   "source": [
    "actual_UB_CS_mean2, actual_UB_CS_error2, UB_CS_mean2, UB_CS_error2, actual_LB_hT_CS_mean2, actual_LB_hT_CS_error2, LB_hT_CS_mean2, LB_hT_CS_error2 = average_actual_theoretical_bounds_CS(\n",
    "    num_datasets=5,\n",
    "    N=800,\n",
    "    c=0.45,\n",
    "    M=100,\n",
    "    sigma2=0.25,\n",
    "    sigma3=0.1,\n",
    "    number_bin=20)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d51dbe67",
   "metadata": {},
   "source": [
    "# Target Shift"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4046884a",
   "metadata": {},
   "source": [
    "## helper functions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcabb85b",
   "metadata": {},
   "source": [
    "#### generate initial distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "653613c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def trunc_gauss(mu, sigma, bottom, top):\n",
    "    '''\n",
    "    generate one data point from a truncated gaussian distribution N(mu, sigma)\n",
    "    return a data point within the range [bottom, top]\n",
    "    \n",
    "    Input:\n",
    "        - mu: mean\n",
    "        - sigma: variance\n",
    "        - bottom: lower bound for the truncated Gaussian variable\n",
    "        - top: upper bound for the truncated Gaussian variable\n",
    "        \n",
    "    Output:\n",
    "        - one random variable sampled from N_[bottom, top](mu, sigma)\n",
    "    '''\n",
    "    a = np.random.normal(mu, sigma, 1)\n",
    "    while (bottom <= a <= top) == False:\n",
    "        a = np.random.normal(mu, sigma, 1)\n",
    "    return a\n",
    "\n",
    "\n",
    "def generate_trunc_gauss_N(mu, sigma, bottom, top, N):\n",
    "    '''\n",
    "    Generate a length of N truncated Gaussian points from N(mu, sigma)\n",
    "    '''\n",
    "    result = []\n",
    "    result = np.array(result)\n",
    "    for i in range(N):\n",
    "        result = np.append(result, trunc_gauss(mu, sigma, bottom, top))\n",
    "    return result\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "dfe9639c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_init_LS(alpha, mu0, mu1, sigma, N, bottom=0, top=1):\n",
    "    '''\n",
    "    Generate N data points from a truncated Gaussian mixture model \n",
    "    (in our case, the truncated upper and lower bound are always 0 and 1)\n",
    "    with two distributions which are N0(mu0, sigma^2), N1(mu1, sigma^2). \n",
    "    The ratio between the two distribution is alpha: (1-alpha)\n",
    "    \n",
    "    N: total number of points\n",
    "    Y ~ bernoulli(alpha)\n",
    "    X1|Y = y ~ N (mu_y, sigma)\n",
    "    X2 ~ -0.8X1 + N(0, sigma2)\n",
    "    X3 ~ 0.2Y + N(0,sigma3)\n",
    "    Y is threshold based on the value of X2\n",
    "    \n",
    "    Input:\n",
    "        - alpha: ratio of Y = 1 instances\n",
    "        - mu0: mean for Y = 0 instances\n",
    "        - mu1: mean for Y = 1 instances\n",
    "        - sigma: variance for all instances\n",
    "        - bottom: lower bound for the truncated Gaussian variable\n",
    "        - top: upper bound for the truncated Gaussian variable\n",
    "    Output:\n",
    "        - a data frame df with ['X1', 'X2', 'X3', 'Y'].\n",
    "    '''\n",
    "    # generate |N*alpha| number of y = 0 instances between 0 and 1:\n",
    "    s1 = generate_trunc_gauss_N(mu0, sigma, bottom, top, int(N * alpha))\n",
    "    # generate |N*(1-alpha)| number of y = 1 instances between 0 and 1:\n",
    "    s2 = generate_trunc_gauss_N(mu1, sigma, bottom, top, int(N * (1 - alpha)))\n",
    "\n",
    "    X1 = np.array(list(s1) + list(s2))\n",
    "\n",
    "    Y = np.array(\n",
    "        list(np.zeros(int(N * alpha))) + list(np.ones(int(N * (1 - alpha)))))\n",
    "    # generate X2, X3\n",
    "    X2 = -X1 * 0.8 + np.random.normal(0, 0.15, N)\n",
    "    X3 = Y * 0.2 + np.random.normal(0, 0.2, N)\n",
    "\n",
    "    # combine them into a dataframe\n",
    "    d = {'X1': list(X1), 'X2': list(X2), 'X3': list(X3), 'Y': list(Y)}\n",
    "    df = pd.DataFrame(data=d)\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "058ef472",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the induced distribution for a particular threshold\n",
    "def compute_new_distribution_LS(df,\n",
    "                                threshold,\n",
    "                                alpha,\n",
    "                                mu0,\n",
    "                                mu1,\n",
    "                                sigma,\n",
    "                                bottom=0,\n",
    "                                top=1):\n",
    "    '''\n",
    "    compute the induced distribution given an initial distribution (df) and a threshold\n",
    "    \n",
    "    Input:\n",
    "        - df: initial distribution with Q\n",
    "        - threshold: a threshold\n",
    "        - alpha: ratio of Y = 1 instances for the initial distribution\n",
    "        - mu0: mean for Y = 0 instances for the initial distribution\n",
    "        - mu1: mean for Y = 1 instances for the initial distribution\n",
    "        - sigma: variance for all instances for the initial distribution\n",
    "        - bottom: lower bound for the truncated Gaussian variable\n",
    "        - top: upper bound for the truncated Gaussian variable\n",
    "    Output:\n",
    "        - an induced dataframe df_new with ['X1', 'X2', 'X3', 'Y'].\n",
    "    '''\n",
    "    N = len(df)\n",
    "    # compute the new data distribution\n",
    "    # get the old classifier result\n",
    "    hX = 1 * (df[\"Q\"] > threshold)\n",
    "    # generate new Y\n",
    "    Y_new = list()\n",
    "\n",
    "    for i in range(N):\n",
    "        # generate a random number in [0,1]\n",
    "        rn = random.random()\n",
    "        # when h(x) = +1, assume the true qualification won't change\n",
    "        if hX[i] == +1:\n",
    "            Y_new.append(df.iloc[i][\"Y\"])\n",
    "        # when h(x) = -1:\n",
    "        # if Y = +1, P(Y' = +1 | h(x) = -1, Y = +1) = 0.8\n",
    "        elif df.iloc[i][\"Y\"] == +1:\n",
    "            if rn < 0.8:\n",
    "                Y_new.append(df.iloc[i][\"Y\"])\n",
    "            else:\n",
    "                Y_new.append(1 - df.iloc[i][\"Y\"])\n",
    "        # if Y = -1, P(Y' = +1 | h(x) = -1, Y = -1) = 0.15\n",
    "        else:\n",
    "            if rn < 0.15:\n",
    "                Y_new.append(1 - df.iloc[i][\"Y\"])\n",
    "            else:\n",
    "                Y_new.append(df.iloc[i][\"Y\"])\n",
    "\n",
    "    # compute alpha_new = P(Y'= +1)\n",
    "    alpha_new = Y_new.count(1) / N\n",
    "    Y_new = np.array(Y_new)\n",
    "\n",
    "    # generate new X1, X2, X3:\n",
    "    # generate |N*alpha| number of y = 0 instances between 0 and 1:\n",
    "    s1 = generate_trunc_gauss_N(mu0, sigma, bottom, top, int(N * alpha_new))\n",
    "    # generate |N*(1-alpha)| number of y = 1 instances between 0 and 1:\n",
    "    s2 = generate_trunc_gauss_N(mu1, sigma, bottom, top,\n",
    "                                N - int(N * alpha_new))\n",
    "    X1_new = np.array(list(s1) + list(s2))\n",
    "\n",
    "    # generate X2, X3\n",
    "    X2_new = -X1_new * 0.8 + np.random.normal(0, 0.15, N)\n",
    "    X3_new = Y_new * 0.2 + np.random.normal(0, 0.2, N)\n",
    "\n",
    "    d = {\n",
    "        'X1': list(X1_new),\n",
    "        'X2': list(X2_new),\n",
    "        'X3': list(X3_new),\n",
    "        'Y': list(Y_new)\n",
    "    }\n",
    "    df_new = pd.DataFrame(data=d)\n",
    "    return df_new\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "fa1bc146",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_optimal_threshold_MinInduceRisk_LS(df,\n",
    "                                               M,\n",
    "                                               alpha,\n",
    "                                               mu0,\n",
    "                                               mu1,\n",
    "                                               sigma,\n",
    "                                               bottom=0,\n",
    "                                               top=1):\n",
    "    '''\n",
    "    find the optimal classifier achieve min Err(h)(h)\n",
    "    where M is the number of potential classifiers to search for\n",
    "    return the optimal threshold and the minimum error \n",
    "    \n",
    "    Input:\n",
    "        - df: initial distribution with 'Q'\n",
    "        - M: total number of threshold to loop through\n",
    "        - alpha: ratio of Y = 1 instances for the initial distribution\n",
    "        - mu0: mean for Y = 0 instances for the initial distribution\n",
    "        - mu1: mean for Y = 1 instances for the initial distribution\n",
    "        - sigma: variance for all instances for the initial distribution\n",
    "        - bottom: lower bound for the truncated Gaussian variable\n",
    "        - top: upper bound for the truncated Gaussian variable\n",
    "    Output:\n",
    "        - an optimal threshold\n",
    "        - the corresponding error\n",
    "    \n",
    "    '''\n",
    "    N = len(df)\n",
    "    # compute the boundary of Q:\n",
    "    Q_min, Q_max = np.min(df[\"Q\"]), np.max(df[\"Q\"])\n",
    "    #print(Q_min, Q_max)\n",
    "\n",
    "    # compute the coefficients on the initial data distribution\n",
    "    clf = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "    clf.fit(df[[\"X1\", \"X2\", \"X3\"]], df['Y'])\n",
    "\n",
    "    # loop through the potential thresholds\n",
    "    error_list = []\n",
    "    for j in range(M):\n",
    "        threshold = Q_min + j * (Q_max - Q_min) / M\n",
    "        # compute the new data distribution\n",
    "        df_new = compute_new_distribution_LS(df,\n",
    "                                             threshold,\n",
    "                                             alpha=0.5,\n",
    "                                             mu0=0,\n",
    "                                             mu1=1.5,\n",
    "                                             sigma=0.25,\n",
    "                                             bottom=0,\n",
    "                                             top=1)\n",
    "        # compute the new qualification\n",
    "        # using the coefficient from the logistic regression coefficient for the\n",
    "        # old dataset\n",
    "        df_new['Q'] = clf.predict_proba(df_new[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "        # compute new data distribution's loss\n",
    "        Err_h_h = compute_error(df_new, threshold)\n",
    "        error_list.append(Err_h_h)\n",
    "    index_min = np.argmin(error_list)\n",
    "    min_error = np.min(error_list)\n",
    "\n",
    "    # compute the optimal threshold on the dataset that it induced call it hT\n",
    "    optimal_threshold = Q_min + index_min * (Q_max - Q_min) / M\n",
    "    return optimal_threshold, min_error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9febce06",
   "metadata": {},
   "source": [
    "### Set parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "eb8fd783",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.5\n",
    "mu0 = 0\n",
    "mu1 = 1\n",
    "sigma = 0.25\n",
    "N = 1000\n",
    "bottom = 0\n",
    "top = 1\n",
    "\n",
    "M = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2111c2ef",
   "metadata": {},
   "source": [
    "#### generate initial distribution df_init_LS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "bd32ceaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_init_LS = generate_init_LS(alpha=alpha,\n",
    "                              mu0=mu0,\n",
    "                              mu1=mu1,\n",
    "                              sigma=sigma,\n",
    "                              N=N,\n",
    "                              bottom=bottom,\n",
    "                              top=top)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25dd17d1",
   "metadata": {},
   "source": [
    "#### compute optimal threshold classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "228ceee5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-2 {color: black;background-color: white;}#sk-container-id-2 pre{padding: 0;}#sk-container-id-2 div.sk-toggleable {background-color: white;}#sk-container-id-2 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-2 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-2 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-2 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-2 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-2 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-2 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-2 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-2 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-2 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-2 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-2 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-2 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-2 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-2 div.sk-item {position: relative;z-index: 1;}#sk-container-id-2 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-2 div.sk-item::before, #sk-container-id-2 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-2 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-2 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-2 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-2 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-2 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-2 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-2 div.sk-label-container {text-align: center;}#sk-container-id-2 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-2 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-2\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>LogisticRegression(solver=&#x27;liblinear&#x27;)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" checked><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">LogisticRegression</label><div class=\"sk-toggleable__content\"><pre>LogisticRegression(solver=&#x27;liblinear&#x27;)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "LogisticRegression(solver='liblinear')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute the best logistic regression model\n",
    "clf = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "clf.fit(df_init_LS[[\"X1\", \"X2\", \"X3\"]], df_init_LS['Y'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "09d013ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the \"qualification\" (linear combination of features) based on the logistic regression model we trained\n",
    "df_init_LS['Q'] = clf.predict_proba(df_init_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2caa4c11",
   "metadata": {},
   "source": [
    "# compute upper and lower bound\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a1e74d3",
   "metadata": {},
   "source": [
    "#### Compute best classifier on the source distribution hS_LS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "e54ee9a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hs_LS: 0.490484016765237\n",
      "Err_S_hs_LS: 0.046\n"
     ]
    }
   ],
   "source": [
    "# compute the optimal threshold hs based on the original dataset\n",
    "hS_LS = compute_optimal_threshold(df_init_LS, M=M)\n",
    "print(\"hs_LS:\", hS_LS)\n",
    "\n",
    "# error on the original dataset S using classifier \"hs\"\n",
    "Err_S_hS_LS = compute_error(df_init_LS, hS_LS)\n",
    "print(\"Err_S_hs_LS:\", Err_S_hS_LS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2faeb32",
   "metadata": {},
   "source": [
    "#### Compute new distribution induced by hS_LS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "c66f10cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute new distribution induced by hS\n",
    "df_hS_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                       hS_LS,\n",
    "                                       alpha=alpha,\n",
    "                                       mu0=mu0,\n",
    "                                       mu1=mu1,\n",
    "                                       sigma=sigma,\n",
    "                                       bottom=bottom,\n",
    "                                       top=top)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "fb89611d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add Q for df_hS_LS\n",
    "# compute Q for df_hT_LS\n",
    "df_hS_LS['Q'] = clf.predict_proba(df_hS_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4ee28d6",
   "metadata": {},
   "source": [
    "#### compute the optimal classifier hT_LS considering induced distribution shift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "67e799b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hT_LS: 0.6998210560138691\n",
      "Err_hT_hT_LS: 0.125\n"
     ]
    }
   ],
   "source": [
    "# compute the optimal hT_LS\n",
    "hT_LS, Err_hT_hT_LS = compute_optimal_threshold_MinInduceRisk_LS(df_init_LS,\n",
    "                                                                 M=M,\n",
    "                                                                 alpha=alpha,\n",
    "                                                                 mu0=mu0,\n",
    "                                                                 mu1=mu1,\n",
    "                                                                 sigma=sigma,\n",
    "                                                                 bottom=bottom,\n",
    "                                                                 top=top)\n",
    "\n",
    "print(\"hT_LS:\", hT_LS)\n",
    "print(\"Err_hT_hT_LS:\", Err_hT_hT_LS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "4a091a68",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the new data distribution induced by hT\n",
    "df_hT_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                       hT_LS,\n",
    "                                       alpha=alpha,\n",
    "                                       mu0=mu0,\n",
    "                                       mu1=mu1,\n",
    "                                       sigma=sigma,\n",
    "                                       bottom=bottom,\n",
    "                                       top=top)\n",
    "\n",
    "df_hT_LS['Q'] = clf.predict_proba(df_hT_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "27462311",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Err_hS_hS_LS - Err_hT_hT_LS 0.04100000000000001\n"
     ]
    }
   ],
   "source": [
    "# compute Err_hs_hs and Err_hT_hT\n",
    "hS_LS = compute_optimal_threshold(df_init_LS, M=M)\n",
    "df_hS_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                       hS_LS,\n",
    "                                       alpha=alpha,\n",
    "                                       mu0=mu0,\n",
    "                                       mu1=mu1,\n",
    "                                       sigma=sigma,\n",
    "                                       bottom=bottom,\n",
    "                                       top=top)\n",
    "\n",
    "# compute the coefficients on the initial data distribution\n",
    "clf = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "clf.fit(df_init_LS[[\"X1\", \"X2\", \"X3\"]], df_init_LS['Y'])\n",
    "# compute Q for df_hs_LS\n",
    "df_hS_LS['Q'] = clf.predict_proba(df_hS_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "Err_hS_hS_LS = compute_error(df_hS_LS, hS_LS)\n",
    "hT_LS, Err_hT_hT_LS = compute_optimal_threshold_MinInduceRisk_LS(df_init_LS,\n",
    "                                                                 M=M,\n",
    "                                                                 alpha=alpha,\n",
    "                                                                 mu0=mu0,\n",
    "                                                                 mu1=mu1,\n",
    "                                                                 sigma=sigma,\n",
    "                                                                 bottom=bottom,\n",
    "                                                                 top=top)\n",
    "\n",
    "# compute df_hT_LS\n",
    "df_hT_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                       hT_LS,\n",
    "                                       alpha=alpha,\n",
    "                                       mu0=mu0,\n",
    "                                       mu1=mu1,\n",
    "                                       sigma=sigma,\n",
    "                                       bottom=bottom,\n",
    "                                       top=top)\n",
    "\n",
    "# compute Q for df_hT_LS\n",
    "df_hT_LS['Q'] = clf.predict_proba(df_hT_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "print(\"Err_hS_hS_LS - Err_hT_hT_LS\", Err_hS_hS_LS - Err_hT_hT_LS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "4a776c80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute w(hS), w(hT) and p\n",
    "w_hS_LS = (np.sum(df_hS_LS['Q'] > hS_LS)) / len(df_hS_LS)\n",
    "w_hT_LS = (np.sum(df_hT_LS['Q'] > hT_LS)) / len(df_hT_LS)\n",
    "p_LS = (np.count_nonzero(df_init_LS['Y'] == 1)) / len(df_init_LS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "087bf1bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dtv_ST_plus 0.041999999999999926\n",
      "dtv_ST_minus 0.020000000000000004\n"
     ]
    }
   ],
   "source": [
    "# compute d_TV(D+(hS), D+(hT)) and d_TV(D-(hS), D-(hT))\n",
    "# Pr_{DS|Y = +1}(hS(X) = +1)\n",
    "P_DSY1_hS1 = np.sum((df_init_LS[\"Q\"] > hS_LS)\n",
    "                    & (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Y\"] == 1)\n",
    "# Pr_{DS|Y = +1}(hT(X) = +1)\n",
    "P_DSY1_hT1 = np.sum((df_init_LS[\"Q\"] > hT_LS)\n",
    "                    & (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Y\"] == 1)\n",
    "\n",
    "dtv_ST_plus = np.abs(P_DSY1_hS1 - P_DSY1_hT1)\n",
    "\n",
    "print(\"dtv_ST_plus\", dtv_ST_plus)\n",
    "\n",
    "# Pr_{DS|Y = 0}(hS(X) = +1)\n",
    "P_DSY0_hS1 = np.sum((df_init_LS[\"Q\"] > hS_LS)\n",
    "                    & (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Y\"] == 0)\n",
    "# Pr_{DS|Y = 0}(hT(X) = +1)\n",
    "P_DSY0_hT1 = np.sum((df_init_LS[\"Q\"] > hT_LS)\n",
    "                    & (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Y\"] == 0)\n",
    "\n",
    "dtv_ST_minus = np.abs(P_DSY0_hS1 - P_DSY0_hT1)\n",
    "\n",
    "print(\"dtv_ST_minus\", dtv_ST_minus)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "02cd5e3f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "UB_LS 0.12699999999999992\n"
     ]
    }
   ],
   "source": [
    "# compute the upper bound\n",
    "UB_LS = np.abs(w_hS_LS - w_hT_LS) + (1 + p_LS) * (dtv_ST_plus + dtv_ST_minus)\n",
    "print(\"UB_LS\", UB_LS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8b58471",
   "metadata": {},
   "source": [
    "### helper function that computes the upper bound for a particular dataset under target shift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "80325f83",
   "metadata": {},
   "outputs": [],
   "source": [
    "def upper_bound_TS(alpha=alpha,\n",
    "                   mu0=mu0,\n",
    "                   mu1=mu1,\n",
    "                   sigma=sigma,\n",
    "                   N=N,\n",
    "                   bottom=bottom,\n",
    "                   top=top):\n",
    "    # #### generate initial distribution df_init_LS\n",
    "\n",
    "    df_init_LS = generate_init_LS(alpha=alpha,\n",
    "                                  mu0=mu0,\n",
    "                                  mu1=mu1,\n",
    "                                  sigma=sigma,\n",
    "                                  N=N,\n",
    "                                  bottom=bottom,\n",
    "                                  top=top)\n",
    "\n",
    "    # #### compute optimal threshold classifier\n",
    "\n",
    "    # compute the best logistic regression model\n",
    "    clf = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "    clf.fit(df_init_LS[[\"X1\", \"X2\", \"X3\"]], df_init_LS['Y'])\n",
    "\n",
    "    # compute the \"qualification\" (linear combination of features) based on the logistic regression model we trained\n",
    "    df_init_LS['Q'] = clf.predict_proba(df_init_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    # # compute upper and lower bound\n",
    "\n",
    "    # #### Compute best classifier on the source distribution hS_LS\n",
    "\n",
    "    # compute the optimal threshold hs based on the original dataset\n",
    "    hS_LS = compute_optimal_threshold(df_init_LS, M=M)\n",
    "    print(\"hs_LS:\", hS_LS)\n",
    "\n",
    "    # error on the original dataset S using classifier \"hs\"\n",
    "    Err_S_hS_LS = compute_error(df_init_LS, hS_LS)\n",
    "    print(\"Err_S_hs_LS:\", Err_S_hS_LS)\n",
    "\n",
    "    # #### Compute new distribution induced by hS_LS\n",
    "\n",
    "    # compute new distribution induced by hS\n",
    "    df_hS_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hS_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    # add Q for df_hS_LS\n",
    "    # compute Q for df_hT_LS\n",
    "    df_hS_LS['Q'] = clf.predict_proba(df_hS_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    # #### compute the optimal classifier hT_LS considering induced distribution shift\n",
    "\n",
    "    # compute the optimal hT_LS\n",
    "    hT_LS, Err_hT_hT_LS = compute_optimal_threshold_MinInduceRisk_LS(\n",
    "        df_init_LS,\n",
    "        M=M,\n",
    "        alpha=alpha,\n",
    "        mu0=mu0,\n",
    "        mu1=mu1,\n",
    "        sigma=sigma,\n",
    "        bottom=bottom,\n",
    "        top=top)\n",
    "\n",
    "    print(\"hT_LS:\", hT_LS)\n",
    "    print(\"Err_hT_hT_LS:\", Err_hT_hT_LS)\n",
    "\n",
    "    # compute the new data distribution induced by hT\n",
    "    df_hT_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hT_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    df_hT_LS['Q'] = clf.predict_proba(df_hT_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    # compute Err_hs_hs and Err_hT_hT\n",
    "    hS_LS = compute_optimal_threshold(df_init_LS, M=M)\n",
    "    df_hS_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hS_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    # compute the coefficients on the initial data distribution\n",
    "    clf = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "    clf.fit(df_init_LS[[\"X1\", \"X2\", \"X3\"]], df_init_LS['Y'])\n",
    "    # compute Q for df_hs_LS\n",
    "    df_hS_LS['Q'] = clf.predict_proba(df_hS_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    Err_hS_hS_LS = compute_error(df_hS_LS, hS_LS)\n",
    "    hT_LS, Err_hT_hT_LS = compute_optimal_threshold_MinInduceRisk_LS(\n",
    "        df_init_LS,\n",
    "        M=M,\n",
    "        alpha=alpha,\n",
    "        mu0=mu0,\n",
    "        mu1=mu1,\n",
    "        sigma=sigma,\n",
    "        bottom=bottom,\n",
    "        top=top)\n",
    "\n",
    "    # compute df_hT_LS\n",
    "    df_hT_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hT_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    # compute Q for df_hT_LS\n",
    "    df_hT_LS['Q'] = clf.predict_proba(df_hT_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    print(\"Err_hS_hS_LS - Err_hT_hT_LS\", Err_hS_hS_LS - Err_hT_hT_LS)\n",
    "\n",
    "    # compute w(hS), w(hT) and p\n",
    "    w_hS_LS = (np.sum(df_hS_LS['Q'] > hS_LS)) / len(df_hS_LS)\n",
    "    w_hT_LS = (np.sum(df_hT_LS['Q'] > hT_LS)) / len(df_hT_LS)\n",
    "    p_LS = (np.count_nonzero(df_init_LS['Y'] == 1)) / len(df_init_LS)\n",
    "\n",
    "    # compute d_TV(D+(hS), D+(hT)) and d_TV(D-(hS), D-(hT))\n",
    "    # Pr_{DS|Y = +1}(hS(X) = +1)\n",
    "    P_DSY1_hS1 = np.sum((df_init_LS[\"Q\"] > hS_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Y\"] == 1)\n",
    "    # Pr_{DS|Y = +1}(hT(X) = +1)\n",
    "    P_DSY1_hT1 = np.sum((df_init_LS[\"Q\"] > hT_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Y\"] == 1)\n",
    "\n",
    "    dtv_ST_plus = np.abs(P_DSY1_hS1 - P_DSY1_hT1)\n",
    "\n",
    "    print(\"dtv_ST_plus\", dtv_ST_plus)\n",
    "\n",
    "    # Pr_{DS|Y = 0}(hS(X) = +1)\n",
    "    P_DSY0_hS1 = np.sum((df_init_LS[\"Q\"] > hS_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Y\"] == 0)\n",
    "    # Pr_{DS|Y = 0}(hT(X) = +1)\n",
    "    P_DSY0_hT1 = np.sum((df_init_LS[\"Q\"] > hT_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Y\"] == 0)\n",
    "\n",
    "    dtv_ST_minus = np.abs(P_DSY0_hS1 - P_DSY0_hT1)\n",
    "\n",
    "    print(\"dtv_ST_minus\", dtv_ST_minus)\n",
    "\n",
    "    # compute the upper bound\n",
    "    UB_LS = np.abs(w_hS_LS -\n",
    "                   w_hT_LS) + (1 + p_LS) * (dtv_ST_plus + dtv_ST_minus)\n",
    "    print(\"UB_LS\", UB_LS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81cc09a0",
   "metadata": {},
   "source": [
    "test for a particular dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "169ad9c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hs_LS: 0.5499162498901023\n",
      "Err_S_hs_LS: 0.038\n",
      "hT_LS: 0.6593395053025968\n",
      "Err_hT_hT_LS: 0.116\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.045\n",
      "dtv_ST_plus 0.014000000000000012\n",
      "dtv_ST_minus 0.030000000000000006\n",
      "UB_LS 0.08900000000000005\n"
     ]
    }
   ],
   "source": [
    "upper_bound_TS(alpha=alpha,\n",
    "               mu0=mu0,\n",
    "               mu1=mu1,\n",
    "               sigma=sigma,\n",
    "               N=N,\n",
    "               bottom=bottom,\n",
    "               top=top)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e881c141",
   "metadata": {},
   "source": [
    "### Lower Bound for target shift \n",
    "\n",
    "Label Shift Lower Bound is (Theorem 5.2 in Section 5.2.):\n",
    "\n",
    "                   max{Err(DS)(h), Err(h)(h)} >= |p - w(h)|(1 - |TPR_S(h) - FPR_S(h)|)/2\n",
    "\n",
    "Quantity needs to be computed (lower bound):\n",
    "    \n",
    "    - Err_Ds(h)\n",
    "    - Err_h(h), Err_h(h), where h\\in {hs, hT}\n",
    "    - TPR_S(h), FPR_S(h)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "7dcbfe84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.165\n",
      "LB_hS 0.0019475806451612905\n"
     ]
    }
   ],
   "source": [
    "# lower bound for hS:\n",
    "Err_S_hS_LS = compute_error(df_init_LS, hS_LS)\n",
    "print(\"max of {Err_S_hs_LS, Err_hs_hs_LS}\", max(Err_S_hS_LS, Err_hS_hS_LS))\n",
    "\n",
    "TPR_S_hS = np.sum((df_init_LS[\"Q\"] > hS_LS)\n",
    "                  & (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Q\"] > hS_LS)\n",
    "FPR_S_hS = np.sum((df_init_LS[\"Q\"] > hS_LS)\n",
    "                  & (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Q\"] > hS_LS)\n",
    "\n",
    "LB_hS = np.abs(p_LS - w_hS_LS) * (1 - np.abs(TPR_S_hS - FPR_S_hS)) / 2\n",
    "print(\"LB_hS\", LB_hS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "2c293a25",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.124\n",
      "LB_hT 0.001892473118279572\n"
     ]
    }
   ],
   "source": [
    "# lower bound for hT:\n",
    "Err_S_hT_LS = compute_error(df_init_LS, hT_LS)\n",
    "print(\"max of {Err_S_hT_LS, Err_hT_hT_LS}\", max(Err_S_hT_LS, Err_hT_hT_LS))\n",
    "\n",
    "TPR_S_hT = np.sum((df_init_LS[\"Q\"] > hT_LS)\n",
    "                  & (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Q\"] > hT_LS)\n",
    "FPR_S_hT = np.sum((df_init_LS[\"Q\"] > hT_LS)\n",
    "                  & (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Q\"] > hT_LS)\n",
    "\n",
    "LB_hT = np.abs(p_LS - w_hT_LS) * (1 - np.abs(TPR_S_hT - FPR_S_hT)) / 2\n",
    "print(\"LB_hT\", LB_hT)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fadc56ab",
   "metadata": {},
   "source": [
    "### helper function that computes the upper bound for a particular dataset under target shift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "d59f9931",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_upper_lower_bound_TS(alpha=alpha,\n",
    "                                 mu0=mu0,\n",
    "                                 mu1=mu1,\n",
    "                                 sigma=sigma,\n",
    "                                 N=N,\n",
    "                                 bottom=bottom,\n",
    "                                 top=top):\n",
    "    # upper bound:\n",
    "    # #### generate initial distribution df_init_LS\n",
    "\n",
    "    df_init_LS = generate_init_LS(alpha=alpha,\n",
    "                                  mu0=mu0,\n",
    "                                  mu1=mu1,\n",
    "                                  sigma=sigma,\n",
    "                                  N=N,\n",
    "                                  bottom=bottom,\n",
    "                                  top=top)\n",
    "\n",
    "    # #### compute optimal threshold classifier\n",
    "\n",
    "    # compute the best logistic regression model\n",
    "    clf = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "    clf.fit(df_init_LS[[\"X1\", \"X2\", \"X3\"]], df_init_LS['Y'])\n",
    "\n",
    "    # compute the \"qualification\" (linear combination of features) based on the logistic regression model we trained\n",
    "    df_init_LS['Q'] = clf.predict_proba(df_init_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    # # compute upper and lower bound\n",
    "\n",
    "    # #### Compute best classifier on the source distribution hS_LS\n",
    "\n",
    "    # compute the optimal threshold hs based on the original dataset\n",
    "    hS_LS = compute_optimal_threshold(df_init_LS, M=M)\n",
    "    print(\"hs_LS:\", hS_LS)\n",
    "\n",
    "    # error on the original dataset S using classifier \"hs\"\n",
    "    Err_S_hS_LS = compute_error(df_init_LS, hS_LS)\n",
    "    print(\"Err_S_hs_LS:\", Err_S_hS_LS)\n",
    "\n",
    "    # #### Compute new distribution induced by hS_LS\n",
    "\n",
    "    # compute new distribution induced by hS\n",
    "    df_hS_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hS_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    # add Q for df_hS_LS\n",
    "    # compute Q for df_hT_LS\n",
    "    df_hS_LS['Q'] = clf.predict_proba(df_hS_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    # #### compute the optimal classifier hT_LS considering induced distribution shift\n",
    "\n",
    "    # compute the optimal hT_LS\n",
    "    hT_LS, Err_hT_hT_LS = compute_optimal_threshold_MinInduceRisk_LS(\n",
    "        df_init_LS,\n",
    "        M=M,\n",
    "        alpha=alpha,\n",
    "        mu0=mu0,\n",
    "        mu1=mu1,\n",
    "        sigma=sigma,\n",
    "        bottom=bottom,\n",
    "        top=top)\n",
    "\n",
    "    print(\"hT_LS:\", hT_LS)\n",
    "    print(\"Err_hT_hT_LS:\", Err_hT_hT_LS)\n",
    "\n",
    "    # compute the new data distribution induced by hT\n",
    "    df_hT_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hT_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    df_hT_LS['Q'] = clf.predict_proba(df_hT_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    # compute Err_hs_hs and Err_hT_hT\n",
    "    hS_LS = compute_optimal_threshold(df_init_LS, M=M)\n",
    "    df_hS_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hS_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    # compute the coefficients on the initial data distribution\n",
    "    clf = LogisticRegression(solver='liblinear', fit_intercept=True)\n",
    "    clf.fit(df_init_LS[[\"X1\", \"X2\", \"X3\"]], df_init_LS['Y'])\n",
    "    # compute Q for df_hs_LS\n",
    "    df_hS_LS['Q'] = clf.predict_proba(df_hS_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    Err_hS_hS_LS = compute_error(df_hS_LS, hS_LS)\n",
    "    hT_LS, Err_hT_hT_LS = compute_optimal_threshold_MinInduceRisk_LS(\n",
    "        df_init_LS,\n",
    "        M=M,\n",
    "        alpha=alpha,\n",
    "        mu0=mu0,\n",
    "        mu1=mu1,\n",
    "        sigma=sigma,\n",
    "        bottom=bottom,\n",
    "        top=top)\n",
    "\n",
    "    # compute df_hT_LS\n",
    "    df_hT_LS = compute_new_distribution_LS(df_init_LS,\n",
    "                                           hT_LS,\n",
    "                                           alpha=alpha,\n",
    "                                           mu0=mu0,\n",
    "                                           mu1=mu1,\n",
    "                                           sigma=sigma,\n",
    "                                           bottom=bottom,\n",
    "                                           top=top)\n",
    "\n",
    "    # compute Q for df_hT_LS\n",
    "    df_hT_LS['Q'] = clf.predict_proba(df_hT_LS[[\"X1\", \"X2\", \"X3\"]])[:, 1]\n",
    "\n",
    "    print(\"Err_hS_hS_LS - Err_hT_hT_LS\", Err_hS_hS_LS - Err_hT_hT_LS)\n",
    "\n",
    "    # compute w(hS), w(hT) and p\n",
    "    w_hS_LS = (np.sum(df_hS_LS['Q'] > hS_LS)) / len(df_hS_LS)\n",
    "    w_hT_LS = (np.sum(df_hT_LS['Q'] > hT_LS)) / len(df_hT_LS)\n",
    "    p_LS = (np.count_nonzero(df_init_LS['Y'] == 1)) / len(df_init_LS)\n",
    "\n",
    "    # compute d_TV(D+(hS), D+(hT)) and d_TV(D-(hS), D-(hT))\n",
    "    # Pr_{DS|Y = +1}(hS(X) = +1)\n",
    "    P_DSY1_hS1 = np.sum((df_init_LS[\"Q\"] > hS_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Y\"] == 1)\n",
    "    # Pr_{DS|Y = +1}(hT(X) = +1)\n",
    "    P_DSY1_hT1 = np.sum((df_init_LS[\"Q\"] > hT_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Y\"] == 1)\n",
    "\n",
    "    dtv_ST_plus = np.abs(P_DSY1_hS1 - P_DSY1_hT1)\n",
    "\n",
    "    print(\"dtv_ST_plus\", dtv_ST_plus)\n",
    "\n",
    "    # Pr_{DS|Y = 0}(hS(X) = +1)\n",
    "    P_DSY0_hS1 = np.sum((df_init_LS[\"Q\"] > hS_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Y\"] == 0)\n",
    "    # Pr_{DS|Y = 0}(hT(X) = +1)\n",
    "    P_DSY0_hT1 = np.sum((df_init_LS[\"Q\"] > hT_LS) &\n",
    "                        (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Y\"] == 0)\n",
    "\n",
    "    dtv_ST_minus = np.abs(P_DSY0_hS1 - P_DSY0_hT1)\n",
    "\n",
    "    print(\"dtv_ST_minus\", dtv_ST_minus)\n",
    "\n",
    "    # compute the upper bound\n",
    "    UB_LS = np.abs(w_hS_LS -\n",
    "                   w_hT_LS) + (1 + p_LS) * (dtv_ST_plus + dtv_ST_minus)\n",
    "    print(\"UB_LS\", UB_LS)\n",
    "\n",
    "    # lower bound:\n",
    "\n",
    "    # lower bound for hS:\n",
    "    Err_S_hS_LS = compute_error(df_init_LS, hS_LS)\n",
    "    print(\"max of {Err_S_hs_LS, Err_hs_hs_LS}\", max(Err_S_hS_LS, Err_hS_hS_LS))\n",
    "\n",
    "    TPR_S_hS = np.sum((df_init_LS[\"Q\"] > hS_LS) &\n",
    "                      (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Q\"] > hS_LS)\n",
    "    FPR_S_hS = np.sum((df_init_LS[\"Q\"] > hS_LS) &\n",
    "                      (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Q\"] > hS_LS)\n",
    "\n",
    "    LB_hS = np.abs(p_LS - w_hS_LS) * (1 - np.abs(TPR_S_hS - FPR_S_hS)) / 2\n",
    "    print(\"LB_hS\", LB_hS)\n",
    "\n",
    "    # lower bound for hT:\n",
    "    Err_S_hT_LS = compute_error(df_init_LS, hT_LS)\n",
    "    print(\"max of {Err_S_hT_LS, Err_hT_hT_LS}\", max(Err_S_hT_LS, Err_hT_hT_LS))\n",
    "\n",
    "    TPR_S_hT = np.sum((df_init_LS[\"Q\"] > hT_LS) &\n",
    "                      (df_init_LS[\"Y\"] == 1)) / np.sum(df_init_LS[\"Q\"] > hT_LS)\n",
    "    FPR_S_hT = np.sum((df_init_LS[\"Q\"] > hT_LS) &\n",
    "                      (df_init_LS[\"Y\"] == 0)) / np.sum(df_init_LS[\"Q\"] > hT_LS)\n",
    "\n",
    "    LB_hT = np.abs(p_LS - w_hT_LS) * (1 - np.abs(TPR_S_hT - FPR_S_hT)) / 2\n",
    "    print(\"LB_hT\", LB_hT)\n",
    "\n",
    "    return Err_hS_hS_LS - Err_hT_hT_LS, UB_LS, max(Err_S_hS_LS,\n",
    "                                                   Err_hS_hS_LS), LB_hS, max(\n",
    "                                                       Err_S_hT_LS,\n",
    "                                                       Err_hT_hT_LS), LB_hT\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "8dd76222",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hs_LS: 0.5001083301743038\n",
      "Err_S_hs_LS: 0.036\n",
      "hT_LS: 0.4103350804334119\n",
      "Err_hT_hT_LS: 0.121\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.04999999999999999\n",
      "dtv_ST_plus 0.020000000000000018\n",
      "dtv_ST_minus 0.011999999999999997\n",
      "UB_LS 0.06700000000000003\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.172\n",
      "LB_hS 0.0020160000000000017\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.122\n",
      "LB_hT 0.0018595041322314074\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.04999999999999999,\n",
       " 0.06700000000000003,\n",
       " 0.172,\n",
       " 0.0020160000000000017,\n",
       " 0.122,\n",
       " 0.0018595041322314074)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_upper_lower_bound_TS(alpha=alpha,\n",
    "                             mu0=mu0,\n",
    "                             mu1=mu1,\n",
    "                             sigma=sigma,\n",
    "                             N=N,\n",
    "                             bottom=bottom,\n",
    "                             top=top)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "c5bab61c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def average_actual_theoretical_bounds_TS(num_datasets, alpha, mu0, mu1, sigma,\n",
    "                                         N, bottom, top):\n",
    "\n",
    "    actual_UB_list_TS = []\n",
    "    UB_list_TS = []\n",
    "\n",
    "    actual_LB_hS_list_TS = []\n",
    "    LB_hS_list_TS = []\n",
    "\n",
    "    actual_LB_hT_list_TS = []\n",
    "    LB_hT_list_TS = []\n",
    "\n",
    "    for i in range(num_datasets):\n",
    "        print(i)\n",
    "        actual_UB, UB, actual_LB_hS, LB_hS, actual_LB_hT, LB_hT = compute_upper_lower_bound_TS(\n",
    "            alpha, mu0, mu1, sigma, N, bottom, top)\n",
    "\n",
    "        actual_UB_list_TS.append(actual_UB)\n",
    "        UB_list_TS.append(UB)\n",
    "\n",
    "        actual_LB_hS_list_TS.append(actual_LB_hS)\n",
    "        LB_hS_list_TS.append(LB_hS)\n",
    "\n",
    "        actual_LB_hT_list_TS.append(actual_LB_hT)\n",
    "        LB_hT_list_TS.append(LB_hT)\n",
    "\n",
    "    # compute the mean and error bar for the actual upper bound\n",
    "    actual_UB_TS_mean = np.mean(actual_UB_list_TS)\n",
    "\n",
    "    actual_UB_TS_errormin = np.abs(\n",
    "        np.min(actual_UB_list_TS) - actual_UB_TS_mean)\n",
    "    actual_UB_TS_errormax = np.max(actual_UB_list_TS) - actual_UB_TS_mean\n",
    "\n",
    "    actual_UB_TS_error = [actual_UB_TS_errormin, actual_UB_TS_errormax]\n",
    "\n",
    "    # compute the mean and error bar for the theoretical upper bound\n",
    "\n",
    "    UB_TS_mean = np.mean(UB_list_TS)\n",
    "\n",
    "    UB_TS_errormin = np.abs(np.min(UB_list_TS) - UB_TS_mean)\n",
    "    UB_TS_errormax = np.max(UB_list_TS) - UB_TS_mean\n",
    "\n",
    "    UB_TS_error = [UB_TS_errormin, UB_TS_errormax]\n",
    "\n",
    "    print(\"actual_UB_TS_mean:\", actual_UB_TS_mean)\n",
    "    print(\"UB_TS_mean:\", UB_TS_mean)\n",
    "\n",
    "    print(\"actual_UB_TS_error:\", actual_UB_TS_error)\n",
    "    print(\"UB_TS_error:\", UB_TS_error)\n",
    "\n",
    "    # for hS\n",
    "    actual_LB_hS_TS_mean = np.mean(actual_LB_hS_list_TS)\n",
    "    LB_hS_TS_mean = np.mean(LB_hS_list_TS)\n",
    "\n",
    "    actual_LB_hS_TS_errormin = np.abs(\n",
    "        np.min(actual_LB_hS_list_TS) - LB_hS_TS_mean)\n",
    "    actual_LB_hS_TS_errormax = np.max(actual_LB_hS_list_TS) - LB_hS_TS_mean\n",
    "\n",
    "    actual_LB_hS_TS_error = [\n",
    "        actual_LB_hS_TS_errormin, actual_LB_hS_TS_errormax\n",
    "    ]\n",
    "\n",
    "    LB_TS_hS_errormin = np.abs(np.min(LB_hS_list_TS) - LB_hS_TS_mean)\n",
    "    LB_TS_hS_errormax = np.max(LB_hS_list_TS) - LB_hS_TS_mean\n",
    "\n",
    "    LB_hS_TS_error = [LB_TS_hS_errormin, LB_TS_hS_errormax]\n",
    "\n",
    "    # for hT\n",
    "\n",
    "    actual_LB_hT_TS_mean = np.mean(actual_LB_hT_list_TS)\n",
    "    LB_hT_TS_mean = np.mean(LB_hT_list_TS)\n",
    "\n",
    "    actual_LB_hT_TS_errormin = np.abs(\n",
    "        np.min(actual_LB_hT_list_TS) - LB_hT_TS_mean)\n",
    "    actual_LB_hT_TS_errormax = np.max(actual_LB_hT_list_TS) - LB_hT_TS_mean\n",
    "\n",
    "    actual_LB_hT_TS_error = [\n",
    "        actual_LB_hT_TS_errormin, actual_LB_hT_TS_errormax\n",
    "    ]\n",
    "\n",
    "    LB_TS_hT_errormin = np.abs(np.min(LB_hT_list_TS) - LB_hT_TS_mean)\n",
    "    LB_TS_hT_errormax = np.max(LB_hT_list_TS) - LB_hT_TS_mean\n",
    "\n",
    "    LB_hT_TS_error = [LB_TS_hT_errormin, LB_TS_hT_errormax]\n",
    "\n",
    "    print(\"actual_LB_hS_TS_error\", actual_LB_hS_TS_error)\n",
    "    print(\"LB_hS_TS_error\", LB_hS_TS_error)\n",
    "\n",
    "    print(\"actual_LB_hT_TS_error\", actual_LB_hT_TS_error)\n",
    "    print(\"LB_hT_TS_error\", LB_hT_TS_error)\n",
    "\n",
    "    return actual_UB_TS_mean, actual_UB_TS_error, UB_TS_mean, UB_TS_error, actual_LB_hT_TS_mean, actual_LB_hT_TS_error, LB_hT_TS_mean, LB_hT_TS_error\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "3dc666ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "hs_LS: 0.40079701275218704\n",
      "Err_S_hs_LS: 0.056\n",
      "hT_LS: 0.7096403133690157\n",
      "Err_hT_hT_LS: 0.107\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.071\n",
      "dtv_ST_plus 0.07199999999999995\n",
      "dtv_ST_minus 0.04800000000000001\n",
      "UB_LS 0.20199999999999996\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.185\n",
      "LB_hS 0.0036590909090909096\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.114\n",
      "LB_hT 0.002615384615384614\n",
      "1\n",
      "hs_LS: 0.5397975948245252\n",
      "Err_S_hs_LS: 0.04\n",
      "hT_LS: 0.5198987369200853\n",
      "Err_hT_hT_LS: 0.117\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.05099999999999999\n",
      "dtv_ST_plus 0.02400000000000002\n",
      "dtv_ST_minus 0.064\n",
      "UB_LS 0.19200000000000003\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.178\n",
      "LB_hS 0.002546558704453445\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.127\n",
      "LB_hT 0.0012750929368029753\n",
      "2\n",
      "hs_LS: 0.509402264412573\n",
      "Err_S_hs_LS: 0.035\n",
      "hT_LS: 0.4695273343861188\n",
      "Err_hT_hT_LS: 0.119\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.04600000000000001\n",
      "dtv_ST_plus 0.006000000000000005\n",
      "dtv_ST_minus 0.013999999999999999\n",
      "UB_LS 0.053000000000000026\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.17\n",
      "LB_hS 0.002613861386138616\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.124\n",
      "LB_hT 0.002254368932038836\n",
      "3\n",
      "hs_LS: 0.48992028567702783\n",
      "Err_S_hs_LS: 0.042\n",
      "hT_LS: 0.4699771950885639\n",
      "Err_hT_hT_LS: 0.124\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.05000000000000002\n",
      "dtv_ST_plus 0.008000000000000007\n",
      "dtv_ST_minus 0.003999999999999997\n",
      "UB_LS 0.029000000000000015\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.168\n",
      "LB_hS 0.00236653386454183\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.118\n",
      "LB_hT 0.0026209677419354815\n",
      "4\n",
      "hs_LS: 0.6200546973000426\n",
      "Err_S_hs_LS: 0.04\n",
      "hT_LS: 0.5402817963575242\n",
      "Err_hT_hT_LS: 0.127\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.037000000000000005\n",
      "dtv_ST_plus 0.017999999999999905\n",
      "dtv_ST_minus 0.008\n",
      "UB_LS 0.06299999999999988\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.156\n",
      "LB_hS 0.0015371900826446298\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.119\n",
      "LB_hT 0.0014607218683651804\n",
      "actual_UB_TS_mean: 0.051000000000000004\n",
      "UB_TS_mean: 0.10779999999999998\n",
      "actual_UB_TS_error: [0.013999999999999999, 0.01999999999999999]\n",
      "UB_TS_error: [0.07879999999999997, 0.09419999999999998]\n",
      "actual_LB_hS_TS_error [0.1534553530106261, 0.1824553530106261]\n",
      "LB_hS_TS_error [0.0010074569067292567, 0.0011144439197170231]\n",
      "actual_LB_hT_TS_error [0.11195469278109459, 0.12495469278109458]\n",
      "LB_hT_TS_error [0.0007702142821024421, 0.0005756605230300641]\n"
     ]
    }
   ],
   "source": [
    "actual_UB_TS_mean, actual_UB_TS_error, UB_TS_mean, UB_TS_error, actual_LB_hT_TS_mean, actual_LB_hT_TS_error, LB_hT_TS_mean, LB_hT_TS_error = average_actual_theoretical_bounds_TS(\n",
    "    num_datasets=5,\n",
    "    alpha=0.5,\n",
    "    mu0=0,\n",
    "    mu1=1,\n",
    "    sigma=0.25,\n",
    "    N=1000,\n",
    "    bottom=0,\n",
    "    top=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "055e7dc3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "hs_LS: 0.5702476894147293\n",
      "Err_S_hs_LS: 0.02\n",
      "hT_LS: 0.7188552841497854\n",
      "Err_hT_hT_LS: 0.112\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.024000000000000007\n",
      "dtv_ST_plus 0.136\n",
      "dtv_ST_minus 0.012\n",
      "UB_LS 0.28600000000000003\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.128\n",
      "LB_hS 0.0007073170731707323\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.104\n",
      "LB_hT 0.0\n",
      "1\n",
      "hs_LS: 0.5316753644242326\n",
      "Err_S_hs_LS: 0.024\n",
      "hT_LS: 0.6802183542618497\n",
      "Err_hT_hT_LS: 0.106\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.07400000000000001\n",
      "dtv_ST_plus 0.01200000000000001\n",
      "dtv_ST_minus 0.004\n",
      "UB_LS 0.036000000000000025\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.17\n",
      "LB_hS 0.0020555555555555566\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.096\n",
      "LB_hT 0.0020806451612903213\n",
      "2\n",
      "hs_LS: 0.5901999715581854\n",
      "Err_S_hs_LS: 0.018\n",
      "hT_LS: 0.669531670077337\n",
      "Err_hT_hT_LS: 0.1\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.12199999999999998\n",
      "dtv_ST_plus 0.020000000000000018\n",
      "dtv_ST_minus 0.048\n",
      "UB_LS 0.192\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.206\n",
      "LB_hS 0.0008489795918367332\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.084\n",
      "LB_hT 0.0007480916030534355\n",
      "3\n",
      "hs_LS: 0.49949100519448775\n",
      "Err_S_hs_LS: 0.018\n",
      "hT_LS: 0.7073319915884977\n",
      "Err_hT_hT_LS: 0.124\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.02600000000000001\n",
      "dtv_ST_plus 0.03200000000000003\n",
      "dtv_ST_minus 0.004\n",
      "UB_LS 0.11400000000000005\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.128\n",
      "LB_hS 0.0008764940239043814\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.102\n",
      "LB_hT 0.0017190082644628038\n",
      "4\n",
      "hs_LS: 0.5211428650732234\n",
      "Err_S_hs_LS: 0.022\n",
      "hT_LS: 0.8384881079879871\n",
      "Err_hT_hT_LS: 0.112\n",
      "Err_hS_hS_LS - Err_hT_hT_LS 0.032000000000000015\n",
      "dtv_ST_plus 0.008000000000000007\n",
      "dtv_ST_minus 0.036000000000000004\n",
      "UB_LS 0.09599999999999999\n",
      "max of {Err_S_hs_LS, Err_hs_hs_LS} 0.138\n",
      "LB_hS 0.000710204081632655\n",
      "max of {Err_S_hT_LS, Err_hT_hT_LS} 0.106\n",
      "LB_hT 0.0013125000000000012\n",
      "actual_UB_TS_mean: 0.055600000000000004\n",
      "UB_TS_mean: 0.1448\n",
      "actual_UB_TS_error: [0.031599999999999996, 0.06639999999999999]\n",
      "UB_TS_error: [0.10879999999999998, 0.14120000000000002]\n",
      "actual_LB_hS_TS_error [0.12696028993477998, 0.20496028993477997]\n",
      "LB_hS_TS_error [0.00033239299204927946, 0.0010158454903355448]\n",
      "actual_LB_hT_TS_error [0.08282795099423869, 0.10482795099423868]\n",
      "LB_hT_TS_error [0.0011720490057613123, 0.000908596155529009]\n"
     ]
    }
   ],
   "source": [
    "actual_UB_TS_mean2, actual_UB_TS_error2, UB_TS_mean2, UB_TS_error2, actual_LB_hT_TS_mean2, actual_LB_hT_TS_error2, LB_hT_TS_mean2, LB_hT_TS_error2 = average_actual_theoretical_bounds_TS(\n",
    "    num_datasets=5,\n",
    "    alpha=0.5,\n",
    "    mu0=0.2,\n",
    "    mu1=0.8,\n",
    "    sigma=0.15,\n",
    "    N=500,\n",
    "    bottom=0,\n",
    "    top=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3811da7d",
   "metadata": {},
   "source": [
    "## Plot the upper and lower bound for Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "960dd262",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR8AAADVCAYAAABuU/QrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAapklEQVR4nO3df3RV1ZXA8e8GAxFQSyFUJECCCooRnzXoiHVI/FW0iq1TGfFHR7AqWEdlRi3aKtE6yxm02On4E5WiXRZk7LIjijqOE4iOVomaKsiIDIKmShOjK6gQJLDnj3tfeAl5Sd7Lve/+ePuzVhZ5J/fdc172Zefcc+49V1QVY4zJtT5BN8AYk58s+RhjAmHJxxgTCEs+xphAWPIxxgTCko8xJhCWfIwxgbDkY4wJxD5BN6A7IjIQuBf4Glipqo8F3CRjjAcC6fmIyCIRaRCRNR3Kp4jIeyKyQUTmusXnAE+o6qXA1Jw31hjji6BOuxYDU1ILRKQvcA9wOjAemC4i44Fi4CN3s105bKMxxkeBnHapao2IlHQoPhbYoKobAURkKXA2UI+TgOroYbIcOnSolpR03L0J0htvvPGpqhZ5uU+Lc/hkEucwjfmMYE8PB5ykcxzwa+BuEfkesDzdm0XkMuAygFGjRlFbW+tjU02mRGSzR/uxOIdYJnEOU/KRTspUVb8CZnT3ZlVdCCwEKC8vt1v1Y8riHB9hmmqvB0amvC4GPs5kByJylogsbG5u9rRhJlwszvEQpuSzGjhUREpFpB9wHvBUJjtQ1eWqetkBBxzgSwNNOFic4yGoqfYlwKvAOBGpF5FLVLUVuBJ4HlgHLFPVtUG0zxjjv6Bmu6anKV8BrMh2vyJyFnDWIYccku0uTARYnOMhTKddvRa67nhFhfNlPBW6OJusxCr5GGOiI1bJJzSzIPPnQ3V1+7Lqaqfc9Fpo4mx6JVbJJzTd8YkTYdo0+Pxz53V1tfN64sRg2xUToYmz6ZVYJZ/QqKyEZctg3Tr44AMn8Sxb5pQbY4CYJZ9QdccrK2H4cPjwQ5g92xKPh0IVZ5O1WCWfUHXHq6vhk09g1Ci47769x4BM1kIVZ5O1WCWf0EiO8Rx+OJSWOqdc06ZZAjImhSUfP6xe7SScwYOd18kxoNWrg22XMSFiyccP11+/9xhPZaVTbowBwrWkRq+F7rL7lSuDbkEshS7OJiux6vnYQGR+sDjHQ6ySjzEmOiz5GGMCEf/kY3eWGxNKsUo+duVrfrA4x0Osko8NROYHi3M8xCr5tGPLWhgTavFNPrashTGhFt/kY8taGBNq8U0+YMtaGBNisUo+e82C2LIWsWSzXfEQq+TTbhbElrUIjs/XVtlsVzzEKvm0Y8taGBNq8U0+tqyFMaEW3+RjjAk1Sz7GRFhFRQUVEb130ZKPMSYQsVrJsFO2mqAxoWQ9H2NMICz5xJWtY2R85MVYU6ySj135mh8szvEQq+RjV77mB4tzPMQq+RhjosOSjzEmEPGfaje5M3/+3ou1VVfb/XSmU9bzMd6x1SNNBiz5GO/Y6pEmA5Z8jLds9UjTQ5Z8jLds9cic27RpU9BNyIolH+MdWz0yEJs3bw66CVmx5GO8Y6tHmgxY8ombIB+W6K4eWVdXR11dnVNmq0f6oqqqChFh1apVAIgIIkJVVVWwDctA6JOPiIwRkYdF5Img2xIJNt2dF6qqqlBVJk+eDICqoqqWfJJEZJGINIjImg7lU0TkPRHZICJzu9qHqm5U1Uv8bGes2HS3yZHeDnT73fNZDExJLRCRvsA9wOnAeGC6iIwXkSNF5OkOX8N8bl882XR3Xhk9enQg9fZ2oNvX2ytUtUZESjoUHwtsUNWNACKyFDhbVW8HzvSzPXmj43R3ZaUloBgrKSkJuglZCWLMZwTwUcrreresUyIyRETuB44WkRu62O4yEakVkdrGxkbvWhs1MZ/utjgHy8uB7iCSj3RSpuk2VtUmVZ2lqge7vaN02y1U1XJVLS8qKvKkoZEU8+lui3OwvBzoDiL51AMjU14XAx97sWNb4Y68eFiixTkegkg+q4FDRaRURPoB5wFPebFjW+EuP1icw6G3A91+T7UvAV4FxolIvYhcoqqtwJXA88A6YJmqrvWzHcYY7/V2oNvv2a7pacpXACu8rk9EzgLOOuSQQ7zetQkRi3M8hP4K50xYdzw/WJzjIVbJxxgTHbFaw9m643skb+xMBNoKf1ic4yFWPR/rjucHi3M8xCr5mPC4v6Ul6CaYkItV8rGLz8LjgR07fNu3xTkeYpV8rDveXlx7HxbneIhV8jHt+dn76EzypsOj3R5JFFfXM7ljycd4JnnT4VtujySKq+uZ3LHkEzPW+zBR0ePkIyID/WyIF2wgMjy9j8v798/qfT05zizO8dBt8hGRSSLyLs5NoIjIUSJyr+8ty4INRIbHrMLCjLbP5DizOMdDT3o+dwHfBZoAVPVPwF/72SjjjWx7HwGx4yzP9Oi0S1U/6lC0y4e2GI9l2vsImh1n+aUn93Z9JCKTAHUX/7oKt2tsjIfsOMszPen5zAJ+grPIez3OvYo/8bFNWbOByEjr8XFmcY6HbpOPqn6qqheo6rdUdZiqXqiqTbloXKZsIDK6MjnOLM7x0O1pl4j8hk6eLqGqM31pkfFEIpEIrO5r3LpXZvAeO87yT0/GfJ5O+b4Q+AEePW3CmBR2nOWZbpOPqv4+9bW7KPx/+dYik5fsOMs/2dxecSgwyuuGGNOBHWcx15Mxny9wzsXF/XcL8FOf25UVW14zujI5zsIW54qKCgBWrlwZaDuipienXfvloiFeUNXlwPLy8vJLg25L4CL2HyGT48ziHA9pk4+IfLurN6rqm943x+QbO87yV1c9n1928TMFTvK4LSY/2XGWp9ImH1WtzGVDTH7y6jizcZfo6dFzu0SkDBiPc/0FAKr6qF+NMtGWbQKw4yy/9GS2ax5QgXNQrABOB14G7KAwnrHjLP/05DqfHwInA1tUdQZwFBCphWJMJNhxlmd6knxaVHU30Coi+wMNwBh/m2XykB1neaarqfa7gSXA6yLyDeBB4A3gS+D1nLTOxJ4dZ/mrqzGf94E7gYNwDoQlwKnA/qr6dg7alrGwXflqeiTj48ziHA9pT7tU9V9V9XicdXQ/A34DPAt8X0QOzVH7MmLrvERPNseZxTkeerKY2GZV/RdVPRo4H2epg//1vWUmr9hxln968uicAnfZysdw/iKtB/7G95aZvBL142zTpk1BNyFyuhpwPhWYDnwPZ+BvKXCZqn6Vo7aZPBCX42zz5s1BNyFyuhpwvhH4HXCtqn6Wo/aY/GPHWQR5cRtLVwPOlar6oB0Qxk9RPs6qqqoQEVatWgWAiCAiOX80dVRls5KhMQYn+agqkydPBkBVUVVLPj1kyccYEwhLPsZ4YPTo0UE3IXIs+RjjgZKSkqCbEDmWfIwxgQh98hGR74vIgyLyHyJyWtDtMcZ4w9fkIyKLRKRBRNZ0KJ8iIu+JyAYRmdvVPlT1D6p6KXAx8Lc+NtcYk0M9Wka1FxYDd5OyGp2I9AXuwblzuR5YLSJPAX2B2zu8f6aqNrjf/9x9nzEmBnxNPqpaIyIlHYqPBTao6kYAEVkKnK2qtwNndtyHiAjwz8Cz9hgV0xW7vypaghjzGQF8lPK63i1L5++BU4AfisisdBuJyGUiUisitY2Njd601IROV3G2+6uiJYjkI52UabqNVfXXqnqMqs5S1fu72G6hqparanlRUZEnDTXhY3GOjyCSTz0wMuV1MfCxFzt2l2RY2Nzc7MXuTEilxtnur4ouUU3b6fCmAmfM52lVLXNf74OzVsvJwJ+B1cD5qrrWqzrLy8u1trbWq90ZD4jIG6pa7uU+U+NcUVHBqlWr8Pt47ow9sHCPTOLs91T7EuBVYJyI1IvIJaraClwJPA+sA5Z5mXiMMdHg92zX9DTlK3AeDOcpW1g8P6SLs91fFS2hv8I5E2FbWLyioqKtS268ky7Odn9VtMQq+RhjoiNWycdmu/KDxTkeYpV8wnbaZfxhcY4Hv+/tMiYnbJo7emLV87HueH6wOMdDrJKPdcfzg8U5HmKVfIwx0WHJx2e2zIMxnbMBZ5/ZMg/x96u6uqCbEEmx6vnYQGR+sDjHQ6yST1gGIm2ZB3+FJc5JiUSCRCIRdDMix/clNYLQ2ZIaO3fupL6+npaWlpy1Y8uWLezYscO3Gx4LCwspLi6moKDAl/17ye8lNZKCiDNbtjj/Hnigr9VEId6ZxDlvxnzq6+vZb7/9KCkpwVkW2n99+vThiy++4PDDD/d836pKU1MT9fX1lJaWer7/qAoizvRxTyDGjfOtijjGO1anXV1paWlhyJAhuTsgXf369fNlvyLCkCFDcvsXPgKCirPf4hjvWCWfzgYiU5e1yPSA9GJJjP79+/fq/V2J23+wnupuwDmIOOdC3OIdq+QTtoHIjvr27UsikeCII47gqKOOYsGCBezevRuA2tparrrqKgB27NjBKaecQiKR4PHHH+ell17iiCOOIJFIsH379iA/QiiEPc4iwkUXXdT2urW1laKiIs48c68nQ+W1vBnzyZaXFwnuu+++1LnXhDQ0NHD++efT3NzMLbfcQnl5OeXlzjjdW2+9xc6dO9u2nTVrFtdeey0zZszwrC2mPS/jPHDgQNasWcP27dvZd999eeGFFxgxoqunQ+WnWPV8/ODXRYLDhg1j4cKF3H333agqK1eu5Mwzz6ShoYELL7yQuro6EokEDzzwAMuWLePWW2/lggsu8KUtxvs4n3766TzzzDMALFmyhOnT96wo/PrrrzNp0iSOPvpoJk2axHvvvQfAggULmDlzJgDvvPMOZWVlbNu2zdN2hYklnwCNGTOG3bt309DQ0FY2bNgwHnroIU488UTq6uq4/PLLmTp1KnfccQePPfZYgK01e9myBbZubV+2dSuoct5557F06VJaWlp4++23Oe6449o2Oeyww6ipqeGtt97i1ltv5cYbbwTgmmuuYcOGDTz55JPMmDGDBx54gAEDBuTyE+WUJZ9O5PIiwTheZxUVvY7zgAGwcSO0tjqvt251XoswYcIENm3axJIlSzjjjDPava25uZlzzz2XsrIy5syZw9q1zsNb+vTpw+LFi7nooouYPHkyJ5xwglcfNZRilXy8uuy+qqoKVWXy5MmAkyBU1fPks3HjRvr27cuwYcM83W/chSbO++8PY8ZASwvs2OEknjFj2n48depUrr322nanXAA33XQTlZWVrFmzhuXLl7ebPn///fcZNGgQH3/syXM0Qy1WySfssyCpGhsbmTVrFldeeWXsplD9Fqo4778/FBTA119DUZHz2jVz5kxuvvlmjjzyyHZvaW5ubhuAXrx4cbvyq6++mpqaGpqamnjiiSdy8hGCEqvk4wcvb43Yvn1721T7Kaecwmmnnca8efM827/JXtZx3roVdu6Efv2gsbHdGFBxcTFXX331Xm+5/vrrueGGGzjhhBPYtWtXW/mcOXO44oorGDt2LA8//DBz585tNx7YqYoK5yuCbKq9G14+Cyr1QOso9UK3jhe9pf51NP7IKs7JMZ7CQthnHxg+HDZu5MtOTplSY3r88cezfv36tp/94he/AGDRokVtZSNHjmTDhg2ZtylCLPl0wRYlzw9Zx3nbNmeM55NPnNfJMaBt29qdfpnOWfLx0TgfbzQ0IZC8iz2ZfMBJOpZ4esTGfIyJovnzobq6fVl1tVOeCx6MNVnyMSaKJk6EadPg88+d19XVzuuJE4NtVwYs+RgTRZWVsGwZrFsHH3zgJJ5ly5zyiLDk05UIT2OaDEQ1zpWVzgzbhx/C7NmRSjwQs+QT9oXFN23aRFlZWbuyqqoq7rzzTi6++GJKS0tJJBIcdthh3HLLLQG1MvzCHudBgwbtVVZVVcWIESPa4jt79uy25VSyVl3tDHaPGgX33bf3GFDIxSr5eHbla0CDeXfccQd1dXXU1dXxyCOP8MEHH/haX1RFNc5z5syhrq6Od999l3feeaftnrKsJMd4Dj8cSkudU65p0yKVgGKVfDwT8GBe8l6fgQMH5qS+vBVQnL/++mtaWloYPHhw9jtZvdpJOMl9JMeAVq/2ppE5YMmnMwEN5l133XUkEgmKi4s577zz7IZTv+U4znfddReJRILhw4czduzY3j1u5/rr925nZaVTHhGWfNLxYTAv3Q2kyfLkadeWLVt48cUXeeWVV3pdp+lGDgdtk6ddDQ0NfPXVVyxdutS3uqLAkk86PgzmDRkyhM+TXXzXZ599xtChQ9uVDRo0iIqKCl5++eVe12m6EcCgbUFBAVOmTKGmpsb3usLMkk9nfBrMGzRoEMOHD+fFF18EnMTz3HPP8Z3vfKfddq2trbz22mscfPDBvarPdCOgQVtV5ZVXXsn7+Fry6YyPg3mPPvoot912G4lEgpNOOol58+a1HYTJMZ8JEyZw5JFHcs455/S6PtMFn+K8bds2iouL274WLFgA7BnzKSsro7W1lSuuuKK3nyDS8uLG0oyfTJActEu91qay0pPxgPHjx1PdyV9WWzYjS/PnO7NTyVjdd5+z1MW2bd0/vtinOKe7fsePZXijLC+ST9ZPJrAlNcIvOV1+0EFOD6alZa/lTLtlce65ZLJPVV3t9BYznGmz0y4TbR2nyxsbncSTL8tarFyZ2+Tp4bVRsU0+HZ9MsHnzZmpra/NiYe68kzpdvt9++ZN4guDhtVGhTz4icriI3C8iT4jI7J6+r+OTCUaPHs0xxxzDQQcd5Ftbc80eu+NKnS7fuhUN6T1fvRWaeHt0bZSvyUdEFolIg4is6VA+RUTeE5ENIjK3q32o6jpVnQVMA8qzbUthYSFNTU3hCWAvqSpNTU0UFhYG3ZRgdZguL/zmN2lavz52CShU8fbo2ii/B5wXA3cDjyYLRKQvcA9wKlAPrBaRp4C+wO0d3j9TVRtEZCow191XxkaPHk1xcTH19fU0NjZms4tQKiwspLi4OOhmBCs5Xe7OWBUnEtTX1dG4bh3k6tE6W7Y4//b2LvVuhCLeqcl+8GCYNy/rUy9fk4+q1ohISYfiY4ENqroRQESWAmer6u3AmWn28xTwlIg8A/wu03aUlJRQUFBAaWlppm81YddhurygoIDSXK/mN9sdDciHWbMOyb7dtVFhSj5pjAA+SnldDxyXZltEpAI4B+gPrOhiu8uAywBGjRrlQTNNGKWNcz78xw8DD6+NCiL5dHZ3ZdqBGFVdCazsbqequhBYCFBeXh6PgR2zF4tzfAQx21UPjEx5XQx4Mv8d9hXujDcszvEgfs/+uGM+T6tqmft6H2A9cDLwZ2A1cL6qrvWwzkYg9bLmocCnXu0/Q1a3Y7SqFnlZgcU58Lo7q7fHcfb1tEtElgAVwFARqQfmqerDInIl8DzODNciLxMPQMcPLyK1qpr1NH1vWN3+sTgHW3dv6/V7tmt6mvIVdDF4bIyJv9Bf4WyMiad8ST4Lre68qDvfPm/QdfeqXt8HnI0xpjP50vMxxoRMLBYTE5EDgV8BE4EdwCbgGuBK4CScixhbgGmq+oGI/BPwI2Cwqu79eEmf6gb+Avw7cDCwC1iuql3eWNtFvUOAF92XB7r7S9649qRb3y5gN3C5qr7mzjJe49ZfpKoZT89mWe9jODcF7wRed8t3ZlF3XsU5qBj3ou7M4qyqkf7CuWL6VWBWSlkCuAl4AujjlhXjHIQAfwUMB77MZd3AAKDSLesHvASc7sHvoAq41v3+eLdN/d3XQ4GD3O+PBkpw/uMMzWG9Z7i/KwGWALMtztGIsZ9xjkPPpxLYqar3JwtUtU5ETgI+UdXdbll9ys//COmfo+Vn3UC1W/a1iLyJc8B6aTjwqarucOtp+8unqm+BJ58703rbLqsQkdfJ7jNbnPcIKsbd1Z1RnOMw5lMGvNFJ+TLgLBGpE5FfisjRYapbRL4BnMWerq1X/hMYKSLrReReEZns8f6zrldECoCLgOey2L/FeY+gYtyjunsa5zgkn065f4XGATfgnJe+KCInh6Fu9xaTJcCv1V1axMO6vwSOwbnzuxF4XEQu9rKOXtR7L1Cjqi95WG/exTmoGGdQd4/iHIfTrrXADzv7gds1fBZ4VkT+Anwfb/8CZVv3QuB9Vf2Vh21JrXsXzkoAK0XkHeDvcBZ281VX9YrIPKAIuDzL3Vuc29cbSIy7qzuTOMeh5/PfQH8RuTRZICITRWSyiBzkvu4DTKD9TYiB1C0itwEH4MxIeE5ExonIoSlFCbz/3BnVKyI/Br4LTE+Oj2TB4ryn7kBi3F3dmcY58slHnWH2HwCnisj/ichanNH5CcBycdaPfhtoxV2GVUTmuze6DhCRehGpykXdIlIM/AwYD7zpjhX8OMuPns4g4BEReVdE3nbrqgIQkavcz10MvC0iD+WiXuB+4FvAq+5nvjnTnVuc2wkqxl3WTYZxtiucjTGBiHzPxxgTTZZ8jDGBsORjjAmEJR9jTCAs+RhjAmHJxyMissudXlwrIn8SkX9wr/3o6j0lInK+D225RkQGeL1fY3H2kiUf72xX1YSqHoHzKOgzgHndvKcE8PygxLmwLbIHZchZnL3ixS339qXQYdkGYAzQhLO8QAnOsgpvul+T3G3+CDQDdcCcLrYbDtS4260BTnTLT8NZ3uBNnPVjBgFXAV8D7+DcWd0X59L3NW7ZnKB/V1H+sjh7+LsMugFx+ep4ULpln+Nc8TkAKHTLDgVq3e8rcJ5pltw+3Xb/CPzM/b4vsB/OOio1wEC3/KfAze73m3DXcsG5CfCFlDq+EfTvKspfFmfvvuJwY2mYJRdVKcC57D6Bs/rb2DTbp9tuNbDIXargD+qsJTMZ59L2/xFn7ZZ+OH8dO9oIjBGRfwOewVkSwXjL4pwFSz4+EZExOAdWA86YwF+Ao3DG2VrSvG1OZ9upao2I/DXwPeC3InIHzl/bFzTNs9GSVPVzETkK54a/n+Asfzmzd5/OJFmcs2cDzj4QkSKcm+zuVqcPfAB7Vry7CKdLDfAFTtc6qdPtRGQ00KCqDwIPA9/GGUc4QUQOcbcZICJjO+5XRIbiLPP5e5xlP7/tz6fOPxbn3rGej3f2FZE6nC51K/BbYIH7s3uB34vIuTiDg1+55W8DrSLyJ5zBwnTbVQDXichO4EvgR6raKM4iTktEpL+73c+B9TjryDwrIp/gzIj8JmU6+AZvP3besTh7xO5qN8YEwk67jDGBsORjjAmEJR9jTCAs+RhjAmHJxxgTCEs+xphAWPIxxgTCko8xJhD/D3eHxL65u0lwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Upper bounds\n",
    "\n",
    "actual_UB_mean_list = [\n",
    "    actual_UB_CS_mean, actual_UB_CS_mean2, actual_UB_TS_mean,\n",
    "    actual_UB_TS_mean2\n",
    "]\n",
    "actual_UB_lower_bounds_list = [\n",
    "    actual_UB_CS_error[0], actual_UB_CS_error2[0], actual_UB_TS_error[0],\n",
    "    actual_UB_TS_error2[0]\n",
    "]\n",
    "actual_UB_upper_bounds_list = [\n",
    "    actual_UB_CS_error[1], actual_UB_CS_error2[1], actual_UB_TS_error[1],\n",
    "    actual_UB_TS_error2[1]\n",
    "]\n",
    "\n",
    "UB_mean_list = [UB_CS_mean, UB_CS_mean2, UB_TS_mean, UB_TS_mean2]\n",
    "UB_lower_bounds_list = [\n",
    "    UB_CS_error[0], UB_CS_error2[0], UB_TS_error[0], UB_TS_error2[0]\n",
    "]\n",
    "UB_upper_bounds_list = [\n",
    "    UB_CS_error[1], UB_CS_error2[1], UB_TS_error[1], UB_TS_error2[1]\n",
    "]\n",
    "\n",
    "# Lower bounds\n",
    "\n",
    "actual_LB_hT_mean_list = [\n",
    "    actual_LB_hT_CS_mean, actual_LB_hT_CS_mean2, actual_LB_hT_TS_mean,\n",
    "    actual_LB_hT_TS_mean2\n",
    "]\n",
    "actual_LB_hT_lower_bounds_list = [\n",
    "    actual_LB_hT_CS_error[0], actual_LB_hT_CS_error2[0],\n",
    "    actual_LB_hT_TS_error[0], actual_LB_hT_TS_error2[0]\n",
    "]\n",
    "actual_LB_hT_upper_bounds_list = [\n",
    "    actual_LB_hT_CS_error[1], actual_LB_hT_CS_error2[1],\n",
    "    actual_LB_hT_TS_error[1], actual_LB_hT_TS_error2[1]\n",
    "]\n",
    "\n",
    "LB_hT_mean_list = [\n",
    "    LB_hT_CS_mean, LB_hT_CS_mean2, LB_hT_TS_mean, LB_hT_TS_mean2\n",
    "]\n",
    "LB_hT_lower_bounds_list = [\n",
    "    LB_hT_CS_error[0], LB_hT_CS_error2[0], LB_hT_TS_error[0],\n",
    "    LB_hT_TS_error2[0]\n",
    "]\n",
    "LB_hT_upper_bounds_list = [\n",
    "    LB_hT_CS_error[1], LB_hT_CS_error2[1], LB_hT_TS_error[1],\n",
    "    LB_hT_TS_error2[1]\n",
    "]\n",
    "\n",
    "# X-axis values (assumed to be sequential integers)\n",
    "x = ['CS1', 'CS2', 'TS1', 'TS2']\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(4, 3), sharey=True)\n",
    "\n",
    "# Plot the actual UB with error bars\n",
    "axs[0].errorbar(\n",
    "    x,\n",
    "    actual_UB_mean_list,\n",
    "    yerr=[actual_UB_lower_bounds_list, actual_UB_upper_bounds_list],\n",
    "    fmt='+',\n",
    "    color='black',\n",
    "    label='Diff')\n",
    "\n",
    "# Plot the theoretical UB with error bars\n",
    "axs[0].errorbar(x,\n",
    "                UB_mean_list,\n",
    "                yerr=[UB_lower_bounds_list, UB_upper_bounds_list],\n",
    "                fmt='x',\n",
    "                color='red',\n",
    "                label='UB')\n",
    "# Add labels and title\n",
    "axs[0].set_xlabel(\"Datasets\")\n",
    "axs[0].set_ylabel('Value')\n",
    "axs[0].set_yscale('log')\n",
    "#axs[0].title('Actual and Theoretical Bounds across different datasets')\n",
    "axs[0].legend(loc=\"lower left\")\n",
    "\n",
    "# Plot the actual LB with error bars\n",
    "axs[1].errorbar(\n",
    "    x,\n",
    "    actual_LB_hT_mean_list,\n",
    "    yerr=[actual_LB_hT_lower_bounds_list, actual_LB_hT_upper_bounds_list],\n",
    "    fmt='+',\n",
    "    color='black',\n",
    "    label='Max')\n",
    "\n",
    "# Plot the theoretical UB with error bars\n",
    "axs[1].errorbar(x,\n",
    "                LB_hT_mean_list,\n",
    "                yerr=[LB_hT_lower_bounds_list, LB_hT_lower_bounds_list],\n",
    "                fmt='x',\n",
    "                color='red',\n",
    "                label='LB')\n",
    "# Add labels and title\n",
    "axs[1].set_xlabel(\"Datasets\")\n",
    "axs[1].set_ylabel('Value')\n",
    "axs[1].set_yscale('log')\n",
    "#axs[0].title('Actual and Theoretical Bounds across different datasets')\n",
    "axs[1].legend(loc=\"lower left\")\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63946c51",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b97ec9e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d042a51",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9be5c97e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f0ba96d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2717b0f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1272f595",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22c892ce",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}