examples/gene_expression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3524d5d4-17b0-4291-9ac0-7fc91fe9260e",
   "metadata": {},
   "source": [
    "# Distribution of gene expression levels\n",
    "\n",
    "The goal is to estimate quantiles that characterize how a population of promoter sequences affects gene expression. Vaishnav et al. (1) recently trained a state-of-the-art transformer model to predict the expression level of a particular gene induced by a promoter sequence. They used the model's predictions to study the effects of promoters; for example, by assessing how quantiles of predicted expression levels differ between different populations of promoters. This notebook shows how the predictions used by Vaishnav et al. can be leveraged to estimate different quantiles of gene expression levels induced by native yeast promoters with higher statistical power.\n",
    "\n",
    "1. E. D. Vaishnav, C. G. de Boer, J. Molinet, M. Yassour, L. Fan, X. Adiconis, D. A. Thompson, J. Z. Levin, F. A. Cubillos, A. Regev, The evolution, evolvability and engineering of gene regulatory DNA. Nature 603(7901), 455–463 (2022)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cec5c87-3ece-43cb-8d86-3f20cb34e46b",
   "metadata": {},
   "source": [
    "### Import necessary packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f49b5f26-2037-4b42-b801-91344f0f186b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import os, sys\n",
    "\n",
    "sys.path.append(os.path.abspath(os.path.join(os.getcwd(), os.pardir)))\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from ppi_py.datasets import load_dataset\n",
    "from ppi_py import ppi_quantile_ci, classical_quantile_ci\n",
    "from tqdm import tqdm\n",
    "from scipy.optimize import brentq\n",
    "from utils import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf90ae6",
   "metadata": {},
   "source": [
    "### Import the gene expression data set\n",
    "\n",
    "Load the data. The data set contains gold-standard labels for the gene expression level (```Y```) and corresponding predicted labels based on the transformer model (```Yhat```)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a6da3138",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dataset gene_expression not found at location ./data/; downloading now...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Downloading...\n",
      "From: https://drive.google.com/uc?id=17PwlvAAKeBYGLXPz9L2LVnNJ66XjuyZd\n",
      "To: /Users/angelopoulos/Code/working/ppi_py/examples/data/gene_expression.npz\n",
      "100%|██████████| 979k/979k [00:00<00:00, 1.22MB/s]\n"
     ]
    }
   ],
   "source": [
    "dataset_folder = \"./data/\"\n",
    "data = load_dataset(dataset_folder, \"gene_expression\")\n",
    "Y_total = data[\"Y\"]\n",
    "Yhat_total = data[\"Yhat\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8969f9db",
   "metadata": {},
   "source": [
    "### Problem setup\n",
    "\n",
    "Specify the error level (```alpha```), range of values for the labeled data set size (```ns```), and number of trials (```num_trials```).\n",
    "\n",
    "Compute the ground-truth value of the estimand (25%, 50%, or 75% quantile)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5b3c8f29",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "n_total = Y_total.shape[0]  # Total number of labeled examples\n",
    "ns = np.linspace(100, 2000, 5).astype(\n",
    "    int\n",
    ")  # Test for different numbers of labeled examples\n",
    "num_trials = 100\n",
    "\n",
    "# True quantiles (computed on all data)\n",
    "true_quantile_q25 = np.quantile(Y_total, 0.25)\n",
    "true_quantile_q50 = np.quantile(Y_total, 0.5)\n",
    "true_quantile_q75 = np.quantile(Y_total, 0.75)\n",
    "\n",
    "q = 0.5  # Target problem is inference on the median"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83ce18be",
   "metadata": {},
   "source": [
    "### Construct intervals\n",
    "\n",
    "Form confidence intervals for all methods and problem parameters. A dataframe with the following columns is formed:\n",
    "1. ```method``` (one of ```PPI```, ```Classical```, and ```Imputation```)\n",
    "2. ```n``` (labeled data set size, takes values in ```ns```)\n",
    "3. ```lower``` (lower endpoint of the confidence interval)\n",
    "4. ```upper``` (upper endpoint of the confidence interval)\n",
    "5. ```trial``` (index of trial, goes from ```0``` to ```num_trials-1```)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "812f8fd5",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [17:53<00:00, 214.64s/it]\n"
     ]
    }
   ],
   "source": [
    "# Run prediction-powered inference and classical inference for many values of n\n",
    "results = []\n",
    "for i in tqdm(range(ns.shape[0])):\n",
    "    for j in range(num_trials):\n",
    "        # Prediction-Powered Inference\n",
    "        n = ns[i]\n",
    "        rand_idx = np.random.permutation(n_total)\n",
    "        _Yhat = Yhat_total[rand_idx[:n]]\n",
    "        _Y = Y_total[rand_idx[:n]]\n",
    "        _Yhat_unlabeled = Yhat_total[rand_idx[n:]]\n",
    "\n",
    "        ppi_ci = ppi_quantile_ci(_Y, _Yhat, _Yhat_unlabeled, q, alpha=alpha)\n",
    "\n",
    "        # Classical interval\n",
    "        classical_ci = classical_quantile_ci(_Y, q, alpha=alpha)\n",
    "\n",
    "        # Append results\n",
    "        results += [\n",
    "            pd.DataFrame(\n",
    "                [\n",
    "                    {\n",
    "                        \"method\": \"PPI\",\n",
    "                        \"n\": n,\n",
    "                        \"lower\": ppi_ci[0],\n",
    "                        \"upper\": ppi_ci[1],\n",
    "                        \"trial\": j,\n",
    "                    }\n",
    "                ]\n",
    "            )\n",
    "        ]\n",
    "        results += [\n",
    "            pd.DataFrame(\n",
    "                [\n",
    "                    {\n",
    "                        \"method\": \"Classical\",\n",
    "                        \"n\": n,\n",
    "                        \"lower\": classical_ci[0],\n",
    "                        \"upper\": classical_ci[1],\n",
    "                        \"trial\": j,\n",
    "                    }\n",
    "                ]\n",
    "            )\n",
    "        ]\n",
    "\n",
    "# Imputed CI\n",
    "imputed_ci = classical_quantile_ci(Yhat_total, q, alpha=alpha)\n",
    "results += [\n",
    "    pd.DataFrame(\n",
    "        [\n",
    "            {\n",
    "                \"method\": \"Imputation\",\n",
    "                \"n\": np.nan,\n",
    "                \"lower\": imputed_ci[0],\n",
    "                \"upper\": imputed_ci[1],\n",
    "                \"trial\": 0,\n",
    "            }\n",
    "        ]\n",
    "    )\n",
    "]\n",
    "\n",
    "df = pd.concat(results, axis=0, ignore_index=True)\n",
    "df[\"width\"] = df[\"upper\"] - df[\"lower\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d15ba288",
   "metadata": {},
   "source": [
    "### Plot results\n",
    "\n",
    "Plot:\n",
    "1. Five randomly chosen intervals from the dataframe for PPI and the classical method, and the imputed interval;\n",
    "2. The average interval width for PPI and the classical method, together with a scatterplot of the widths from the five random draws."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6077b2c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnsAAADMCAYAAAALZGVzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwqElEQVR4nO3deXxkZ33n+8+v9tJSUmlrLb27veM23tqNbcA2CZgAJmYnM8mQZGCYSWaYkJlwIcMkuXmRyb0kuZmE3GTIDYGBECA2JpjFdgDb2GC7u41tebd7U7u7tbTUJZWWqlItz/2jStUltdSt/ZSk7/v1qpfOqfPUqZ+q+0hfPc95zjHnHCIiIiKyPvm8LkBEREREVo7CnoiIiMg6prAnIiIiso4p7ImIiIisYwp7IiIiIuuYwp6IiIjIOhY4z3Zdl0XWE/O6gHPQsSbrSdUea7fddpu79957vS5DZLnM61g7X9gTETmn7DPfLy8Hr3irh5WInN/g4KDXJYisOoU9EVmS/MGflpcV9kREqo/O2RMRERFZx9SzJyJLEnjNbV6XILJko8kkgwMDZNJpwpEILW1t1MdiXpclsiwU9kRkSQIX3uh1CSJLMppMcrynB38gQCgcJpvNcrynh83btinwybqgYVwREdnQBgcG8AcCBINBzIxgMIg/EGBwYMDr0kSWhcKeiIhsaJl0mkBg+kBXIBAgk057VJHI8lLYExGRDS0ciZDL5aY9l8vlCEciHlUksrx0zp6ILEn2yX8uLweveqeHlYgsTktbG8d7eoBij14ulyOfy9HR1eVxZSLLQz17IrIk+aMHyg+Rtag+FmPztm0Eg0EmMxmCwaAmZ8i6op49ERHZ8OpjMYU7WbcU9kRkSQKvvd3rEkRE5BwU9kRkSQI7rvO6BBEROQedsyciIiKyjinsiYiIiKxjCnsiIiIi65jO2RORJZnc/43ycui693lYiYiIzEY9eyKyJIXjz5QfIqvJzPxm9gUz+4mZ/djMLvC6JpFqpLAnIiJr1TsAnHM3Av8d+DNvyxGpThrGFZElCV77Hl6IpJjwFeD04yv2PrWBWi6PvWbF9i9rj3PuW2b2ndLqNqDfy3pEqpXCnogsiX/LlYye/imtoVZaw20r+l7pfJqIXzenlzOcczkz+xJwB/Ce2dqY2UeAjwBs3bp1FasTqQ4axhURkTXNOfdvgIuAvzWz2lm2f945d61z7trW1tbVL1DEYwp7IiKyJpnZL5vZJ0urE0AByHtYkkhV0jCuiIisVd8E/t7MfgwEgf/snEsvZkfjY2O88Oyz1NbWctFll+H3+5e1UBEvKeyJyJJMPvoV6rZ0YEHndSmywTjnxoFlubjjiWPHGEsmGUsmqa2vZ+v27cuxW5GqoLAnIktS6HuJC/teAiByxx96XI3I4gRDofLy0YMHiTU00BiPe1iRyPLROXsiIrKhjSaTpCYmykO3hUKBgy++yGgy6XFlIstDPXsisiTBvf/K6xJElmRwYICMb5LJmiz+0WIfyPjYGAN9fdTHYh5XJ7J0Cnsi69BDgw+u3psFK5YH+87Z9I0tN69oKSKLkRg/zYncCQLhAJHxEFYwnHMMnOrjgosu8ro8kSVT2BNZh1pDbXREOrwu4yypXIpoIOp1GSLTDHGaoAsQCAQp1Dn8SQMgPZEil80SCAbPsweR6qZz9kREZEMbrhkhUAhCHlwYnK84s9yccfzYMY+rE1k6hT0REdnQ6mMxxptTOD+Qg3zdmcsIvdrTQzab9a44kWWgsCciIhva7thuxkPjjLaOk91SILllnFwoB0A+l+N4T4/HFYosjcKeyDrkM/O6BJE1ozPaxc0tt1Djr2E4N0xNoIbNO7aVtx9X756scZqgIbIOXVJ/6aq9V/ruT5eXdVFlWas6o110RrvK6845kseHmRgfJ5/P8+rRo+y88EIPKxRZPIU9EVmS4E2/5nUJIsvOzNh+wQU8390NFG+ntnnbNkIVd9oQWSs0jCsiS+Jv3VF+iKwnrZs2UVNbC1Du3RNZixT2REREZmFm7Ni1q7x+4tgxJjMZDysSWRwN44qsU9/pu8frElbE29vf4XUJsoG0tLVRW1fH+NgYhUKBY0ePsuvii70uS2RBFPZE1qmOSCdbo1u9LmPZjefGqQ3Uel2GbBBmxvZdu3juqacAOPnqq2zZvp1wOOxtYSILoGFcERGRc2hpbaWuvh6AQqHAq0eOeFyRyMIo7ImIiJzD1MzcKSeOHyeTTntYkcjCKOyJiIicR3NrK/WxGACudO6eyFqhc/ZE1imf/pYTmbeTqRN0J7tJZBPEg3F2x3ZPu8jyVO/eM08+WWxfOncvEol4VbLIvCnsiaxTVzVe5XUJImvCydQJHhx8gIg/SmOgkYn8BA8OPsDNLbdMC3xNLS3UNzQwOjKCc45jhw9z0WWXeVi5yPzoT38REdnQupPdhHxhktkko7lRavw1RPxRupPd09qZGTsqzt3rPXGCdCq12uWKLJjCnoiIbGiJbIKBdD/dyad5PPEYickEUV+URDZxVtt4czOxxkageP/cHs3MlTVAw7giG8w3jn+dgcyA12XQFm7jfZvf73UZIsSDcfrSvQA4HC+MPs/uhiuJB+NntZ3q3Xv6iScA6Dtxgq3btxOtqVnVmkUWQmFPZIOJBWPsbXqd12UAMDw5TGOo0esyZIPbHdtNX7qXE6kTFCgwlh/j2EQPH9j8wVnbNzY10RCPM5JIlHv3Lrn88lWuWmT+NIwrIiIbWme0ize3vYULas/cB7c/009toG7W9jPP3es7eZLUxMSK1ymyWAp7IiKy4XVGu/jIjn/HpnA7AFmX5du935qzfWNTE43x0jCvcxw9fHgVqhRZHIU9ERERwG9+3t357vL6kyNP8srYy3O2376roifw5EkmxsdXtD6RxVLYE9lgYoGY1yWIVK1ddRdydcM15fW7Tt5JrpCbtW1jPE68ubm83qPePalSmqAhssHc0Hyj1yWIVLV3dNzOc6PPkilkGMgM8OOhh7i19U2ztt1+wQUkhoYA6O/tZevOndTW1q5muSLnpZ49ERGRCg3BBm7b9Nby+v0D95GYPPuaewANjY00VfbuHTq04vWJLJTCnoiIyAw3Nb+ejkgHAJOFyXNO1qg8d2+gr4/xsbGVLk9kQTSMKyKz+qvDf8mJ1InzN8xlziwHwgt6jz+6/I8XWJXI6ihO1ngPnzv8lwA8nXyal0Zf5OL6S85qG2tooKmlhdODgwAcPXSIy6+8clXrFTkXhT0RmVU82MSb225b0fcYygzRHG4+f0MRD+ysvYBrG6/jwPB+AL558i7+64WfIOA7+1fnjl27ymHvVH8/Y6Oj1NXXr2q9InPRMK6IiMgc3tF+OxFfBIBTk6d4cPCBWdvVx2I0t7aW14/q3D2pIgp7IiIic6gP1vPWTW8rr//LwP2cnhyate32irtqDA4MMJpMrnh9IvOhsCciInIONzTfQFekCyjeWeNbc0zWqI/FaGlrK6+rd0+qhcKeiMyqKdTkdQkiVWFqssaUZ5PP8Hzy+VnbVvbuDZ06pd49qQqaoCEis6q8zpjIRre9dgd74tezL/E4AHf33sWFdRcS9AWntaurr6d10yZO9fcDcOTgQXZfffWq1ytSST17IiIi8/D29ncQ9dcAMDQ5xI9O/XDWdpW9e6cHB0mOjKxKfSJzUdgTkSVJ3/OH5YfIelYXqONtFZM1fnjqBwxNDp7Vrraujrb29vL60YMHV6U+kbko7InI0uQmzzxEVpGZBc3sy2b2sJntM7PbV/o99za9ji3RLQDkXI67T949a7ttlb17Q0OMDA+vdGkic1LYExGRtepfA0POudcDbwU+t9Jv6DMf7+58L4YB8PzoczybfPasdrW1tWzq6Civa2aueEkTNERkScJv/29elyDrhJl9GPjPQBQwwDnndp7jJf8E3Fmxnlu56s7YWrOV65v28tjpRwG4++RdXFR3ESFfaFq7bTt30t/XB86RGBpiOJGgMR5fjRJFplHPnogsiQXD5YfIEn0U+AXgUuCS0tc5OefGnHOjZlZPMfTN+peHmX3EzA6Y2YFTp04tS6Fv2/R2av21ACSyCX448IOz2tTM7N3TuXviEYU9ERGpFoPOuR7nXGbqcb4XmNkW4AHgy865r87Wxjn3eefctc65a1srbmm2FLWBWt7W/vby+o8Gf8ipzNlBcvvOnWDFId/hRILE6dOLfs/RZJIjBw/y4rPPcuTgQV3DT+ZNYU9ERDxlZn9kZn8EhMzsPjP7HxXPnet1m4D7gU84576wKsVW2BO/nq3RbQDkXZ5vnrwL59y0NtGaGto7O8vrRw8dOqvNfIwmkxzv6SGbzRIKh8lmsxzv6VHgk3lR2BMREa+9VHp8Efgq8GJp/cXzvO5TQBz4tJk9WHpEV7LQSj7z8Z6u95Qna7w09iLPJLvPardt506s1Ls3kkgwvIjevcGBAfyBAMFgEDMjGAziDwQYHBhY2jchG4ImaIjIkqTv/nR5OXKHrrUnC+ec+xKAmX3OOfebU8+b2f8G/vc5Xvcx4GPLUcNoMsngwACZdJpwJEJLWxv1sdh5X7c5uoUbmm7kJ6cfAeBbvXdzcf0lhH1nzmGNRqO0d3XRe/w4AEcOHaKxqakcAOcjk04TCk8/LzYQCJBJp+e9D9m41LMnIiKeMrPfMLNe4MNmdrL06AW6VuP9lzpE+tb2X6DOXwfAcHaYfxm4/6w223bsKIe75PAwiaGhBdUYjkTI5aZPNs7lcoQjkQXtRzYmhT0REfGUc+6vnHMdwB845zpLjw7n3JtW4/2XOkRa46/hHR1nruf80OCD9Kf7p7WJRKN0dJ3JrkcWeO5eS1sb+VyObDaLc45sNks+l6OlrW3e+5CNS8O4IrIkGrqVpTKzXyktHq9YBsA5N+cw7nJZjiHSaxqv5bHTj3Jk4kh5ssZHd/z7aUO1W3fupPfECZxzjI6McHpwkOZ5zg6uj8XYvG3btKHmjq6ueQ01i6hnT0REvHZp6fEh4P3AFuBdpeUVtxxDpFN31vCVfq2+Mv4yT408Oa1NJBKhc8uW8vpCZ+bWx2Ls2LWLS17zGnbs2qWgJ/OmsCciIp5yzn3SOfdJYNI59zbn3Gecc78IBFfj/ZdriLQz2smNzTeV1/+5959J56f3Dm7dvh2fr/irdzSZZGiZLvIsci4KeyIiUi3azKwRwMyagebVeNOpIdJgMMhkJkMwGGTztm2L6jm7bdNbqQ8UX5fMjXD/wH3TtocjETo3by6vL/a6eyILoXP2RKRscGCAbDbrdRmLVnkCvKxJnwEOmFkSiAG/tlpvXB+LLcuwaNQf5fb22/mH418B4MeDD7Envof2yJnbpm3ZsYOTx49TKBQYGx1l8NQpWjXRQlaQwp6IlE1OThKtqaGmttbrUhYll8sRCOjH2lrlnLsLuMvM2oAh51ze65oW4+rGa3gs8SiHxg9RoMBdJ+/kP+z4zfJkjXA4TOeWLRzv6QGK98xtaW1d0HX3RBZCw7giIuIpM/tc6esTZvZT4FvAw6XlNcfMeHfne8qTNQ6NH+Jnw09Ma1N57t742JjuhCErSmFPRES8NnX9nh1AD/BN4MPABz2raInaIx28oeWN5fVv932bVD5VXg+Fw3Rt3VpeP3LwoM7dkxWjsCciIp5yzvWXvjYB/yfgBz4P/LmHZS3Zm9veQizQAMBoLsm9/d+ftn3L9u34/X4AJsbHOdXff9Y+RJaDwp6IlIVCIZ03JJ4xsyuBtwG3lp56wcNylizij/DOjneW1x8ZepgTqRPl9VAoNK13TzNzZaXoTGYRKdOtl8RjPwYOA7/rnPue18Ush9c2XMVjpx/jlfGXcTi+efJOfmPnf8Rnxb6WLdu3c+LYMfL5PBPj4wz09bGpo+M8exVZGPXsiYhItWgGfgt4vZn90Mz+0euClsrMeFfnu/Fbcbj2yMQRnhg+UN4+dU2/KUcPHaJQKKx6nbK+KeyJyJK4VLL8EFmiRqAL2AbUUJysseZtimzijS03l9fv6f02E/mJ8vrmbdvwly4ZlJqYYKCvb7VLlHVOw7gisiSZez9bXo7c8YfnaClyXvcBdwOfcc4953Uxy+nn297Mz4afYDg7zFh+jO/3fY93d70HKPbubdm2jaOHDgHF3r229vbypVmmnEydoDvZTSKbIB6Mszu2m86oLiQu52fnORlUZ4rKelLNMw+W5Vjr7+1lcnJyOXY1f9mKe38G53/j+OUSCoV0jlP1qdpj7dprr3UHDhw4f8MV0D3yNF889vcAGMZv7fo4m6NbAMhlszz28MPkcjkALr788ml3hDmZOsGDgw8Q8UeJ+qKkCinS+RQ3t9yiwLexzetYU8+eyDqSyWSoraujtq7O61JWzfjYmNcliMzLFbHdXFx3CS+NvYjDceeJO/lPF3wMn/kIBINs2b6dIwcPAtBz+DCbOjrKvXvdyW4i/ig1/hqA8tfuZLfCnpyXztkTERFZBTMnaxxL9bAvsa+8vWvrVgLBIADpVIq+kyfL2xLZBFFfdNr+or4oiWxiFSqXtU5hT0REZJW0hlu5teVN5fXv9t3DeG4cgEAgwJbt28vbeg4fLs/MjQfjpAqpaftKFVLEg/GVL1rWPIU9ERGRVfSmtp+jKdgEwHh+nO/1f7e8rWvLFoKl3r1MOk3fieJFmHfHdpPOp5jIT+CcYyI/QTqfYnds9+p/A7LmKOyJrCPhcNjrElbd1O2mRNaKkC/EL3a+q7z+2OlHOTZxDCj17u3YUd421bvXGe3i5pZbqPHXMJwbpsZfo8kZMm+aoCGyjngxK7UwOlhe9tW3rPr7RyKrPwNYZKkur7+cy+ov4/nR53E47jr5T3zsgt/CZz66Nm/m1aNHyU5Okslk6D1+nK6tW+mMdincyaKoZ09ElmTyB/+z/BCR+TEz7uh8FwEr9rm8mnqVx04/CoA/EGBrZe/ekSPk83lP6pT1QWFPRETEA82hFt7U+nPl9e/2f5exXPFSQp2bNxMKhQCYLPXuiSyWhnFF1qBDL7/sdQlnXPy+M8sLrCscDk+7L6jIRnNr65s4MLyfockhUvkJvtN3Dx/Y/EH8fj9bd+zg4EsvAXDsyBE6Nm/WOaqyKAp7ImtQXX099Q0NXpexZKMjI16XIOKpoC/IHR3v5v/r+TwA+xKPsze+l+21O+jYvJljR48ymckwOTnJyePH2aI/jmQRNIwrIiLioctil/Ga2BXl9btO3kne5fH7/WyrOHfv2JEj5Eu3UxNZCIU9ERERj/1ixx0ErXh9vRPpE/x06KcAdGzeTLg04zw7OcmJV1/1rEZZuxT2REREPNYUauLn295cXv9+/3cZzY7i8/mm9e69evQoOfXuyQIp7ImsRWZeV7AsAgGdNiwy5eaWW2gNtQKQLqS5p+/bALR3dZ3p3ctmeWr/fo4cPMhoMulZrbK26CetyBq0qb3d6xLK0nd/urwcueMPF/Tampqa5S5HZM0K+AK8q/Pd/K+jfwPAgeH97G3ay87aC2jv7KTn8GEAxsfGKDhHcmSEHbt2EVsHk7VkZalnT0REpEpcXH8JV8auLK9PTdZwzpUvu+KcY2JsjMTQEE8/8QSHX3mF8bExr0qWNUA9eyKyJNbY6XUJIuvKOzvv4IWxF5gsTNKb7uWRoYfZlGmhddMmTg0MTJuRm8/lOHbkCMeOHKGuvp62jg7a2tt1G0GZxpxz59p+zo0ia0w1n+imY03Wk6o91q699lp34MCBs54/mTpBd7KbRDZBPBhnd2y3p/eh/dGpH/KdvnsACPvC/FLog4RzYQKBAJl0mvGxMcbHxpjrd3hjPE5bRwetmzYRDAZXs3RZXfM61hT2ZCOp2l9ALOFYe+m555iYmFjOWpbdVddd53UJsrqq9libLeydTJ3gwcEHiPijRH1RUoUU6XyKm1tu8Szw5Qo5/vTgZ+nP9ANwRc1urhl5LdlMhnyhgN/nIxAK0dTcTHJkhKFTpygUCmftx8xobm2lrb2d5tZW3YFj/ZnXsaZhXJE1bmx0lKbWVppbW70uZU7pdFrDSlK1upPdRPxRavzFCUNTX7uT3Z6FvanJGn995P8F4JmJbjrZRLM1Y2Y4MwxobGpiy/bt5HI5TvX3M9DbS+L06fJ+nHMMDgwwODCAPxCgta2Nto4O4k1N2DqZ1S/np7AnIiIbWiKboDHQOO25qC9KIpvwpqCSC+su4qqGq3hy5EkA9gX388H6D+K3Yu9cNptlcGCA+liMQCBAR1cXHV1dZDIZBvr6GOjtnXZ5lnwuR9/Jk/SdPEkoHKatvZ1NHR3U1dcr+K1zCnsiIrKhxYNxJvIT5R49gFQhRTwY97Cqots7fpHnR58nU8iQYJinc09zdfBqgPL5ezOFw2G2bNvGlm3bmBgfp7+3l/7eXtKpVLnNZCbD8Z4ejvf0EK2pYVNHB5s6Oojqckjrki69IiIia5qZXW9mDy729btju0nnU0zkJ4qXNclPkM6n2B3bvYxVLk5DsIG3tN1WXt+X28+YK15mJZfLlS+2PJea2lp27NrF9TfdxNXXX0/X1q0EQ6FpbVITExw9dIjHH3mEJx5/nOPHjjGZySz/NyOe0QQN2UiqeZxiSRM0wtEoTS0ty1nPsgqFQjpnb2NZtWPNzH4H+GVg3Dm393zt18ps3Ep5l+ezL/1fDGQHAGguNPHa/JW0+9q5aOel1MdiC9pfoVBg+PRp+nt7GRwYIJ/Pn93IjHhTE5s6Omhpa9PdbqqXZuOKzLAuw57XMg//XXk5/Ppf97ASqSKrGfbeDXQDX15K2Kt2j59+jK+f+Fp53ZzR6G/gDW1v5Mbm1xPwLS6M5fN5hk6dor+3l9ODg7NeysXn89FSmtjR1NyMz6dBwSqisCcyg8KeyOpY1WPNzLYDX5sr7JnZR4CPAGzduvWanp6eVaxuedzb/326R56mL9N31rZafy1XN17DnvgeuqKbF/0e2clJBkozekeGh2dtEwgGadu0ibaODhoaGzWxw3sKeyIzVPNPJR1rsp5UVdirtFZ79v7x+Fdp8DfQm+nlRPo4w9nhWdt1Rjq5Lr6HaxqvpS5Qt+j3S6dS9Jdm9M51K7ZwJMKm0h076urrF/1esiQKeyIzrNmw9/CPfrRadVSF1996q9clyNIo7C2ze/u/P23G8HhunGOpHgYyp8gUzp6R68PHZbHL2RPfw6X1l5Uv17IYY6Oj9Pf2MtDXN+vsX4Daurpi8Ovo0Pm5q0thT2SGNRv2Xnj2Wdo7N849aCORiC4BsbYp7C2zue7y8YbmNzJRmGBfYh/dI93kXPas19b567i68Rqui++hawmTTpxzjCQS9Pf1caqvj1zFPXorNcTjxJuaqKmtpaa2lmhNje7csXIU9kRmUNhbIxT21ryqPdbWatiD888YTuVTPDXyFPsT+zg6cWTWfXRFurguvoerG69Z0jBvoVDg9OAg/b29c96qrVIkEiFaCn/lR00NoXBY5/0tjcKeyAzV/BNFYa+Cwt6aV7XH2loOewsxkBngQGI/B4b3z3p+n9/8XFZ/OdfF93Bp/aVLGubN5XIM9vfT39dHYmhoQa/1+/1negBLAVC9gQuisCcyQ9X+AkJhbxqFvTWvao+1jRL2phRcgVfGXmF/Yh/dybmHea9pvIbrmq6nM7L4nzOjySQ9hw+TzWYp5PNMTk6Sy2Znv47fPESi0XIPYGWvYCgUUm/gGfP6IHSVRJE1QNe1EpHF8JmPi+sv5uL6i0nlJyqGeY+W24zlx3ho6CEeGnqIzZHNXBffw1WNVy94mHdwYADMwDkKhQKRSIRgQwM1pduxTYyPT39MTJCf47w/KM4ITqdSnJ7xvD8QKPcAVp4XqN7AualnTzaSav5TUMearCdVe6xttJ69uQxk+tmX2MeBxAGSuZGztk8N8+6J7+GSeQ7zdv/sZ6RTKfw+H+bz4QoF8oUCkWiU3VdffVZ75xyTk5OkKsLfVBCsvI/vQpR7AyvOC4zO6A0cTSYZHBggk04TjkRoaWtb8F1IqoiGcUVmqNpfQOhYk/Wlao81hb3pCq7Ay2MvsT+xj2eSz5BzZ/e01QfquabxWq6L76Ej0jHnvp7ct49sLkew4tZqU+tX7dmzoLry+TypiQlSFQFwPr2Bc5nqDQyGQqTGx8s/cP1+P6FQiJ0XXbRWA5/CnsgMVfsLCB1rsr5U7bGmsDe3VH6CJ4efZF9iH8dSs99lZHN0C9c17uHqxqupDdRO2/bMk0+SmpjA5/fj8/koFAoU8nmiNTVccdVVy1LjVG/gxPj4WT2Ci+0NhGLoq4/FCEci0x6R0tdAIFCt5wkq7InMUJVHaomONVlPqvZYW8th73yXXllO/ek+9g/v50BiP8lc8qztfvPzmvrXcF18DxfXX4Lf/Bw5eJCJ8XEmMxlyuRyBQIBQOExNbS07du1akTorTfUGToyPn9UjuNhJIlN8Pt+08HdWIAyHCQSDy/SdLIjCnsgMVfsLiDV8rOVffbq87N9ypYeVSBWp2mNtrYa9uS6qfHPLLSsW+ADyLs/LYy+zP7GPZ+cc5o1xTeM1XB68nMneFP5AgEAgQC6XI5/LsXnbNk+HSCt7A1945hny+TyuUDjvtQEXyu/3n9UjOC0UhsP4A8s+L1azcUVk5WUP3FleVtgTWRndyW4i/mj5dmlTX7uT3Ssa9vzm59L6S7m0/lIm8hM8Ofwz9if2cSx1rNxmNJfkwcEHeJAH6Ah0sGWii+Z0nIZgnG1bdnh+LpyZEQ6HCYfDxBoaOD00hM/vJxAMUigUyOfzNDQ2snX7djLpNOl0mkwmQyadLj/m0zOYz+fLPYlzCQQCsw4Th0u9g+FIBL/fv+yTSBT2REREqlwim6Ax0DjtuagvSiKbWLUaavw13Nh8Ezc230Rfuo/9iX0cGD7AaMUwb2+ul95QL4SK6+HBME3JJloirTQGG2kKNdEYjNMUjNMYilPnr1vVc+ECwSDBQADnHIV8Hp/Phz8YJBKN0tTSMutrnHPkcrlp4S+TyRRDYcVjPj2FuVyO3NgY42Njc9dY6v0LBIM0NDaSzWY53tOzpB5ShT0RWRLf5iu8LkFk3YsH40zkJ8o9egCpQop4MO5JPe2Rdt7RcTu/0P42Xh57iX2JfTwz0k2B6YEnQ4beyV56J3tn3U/AgsSDjcRDTcSD8eIjNPW1iYZAAwHfMkYV52hsbiY9MVE+rzBSUwPnOKXNzAgGgwSDQerq6+fYrSObzZ43EJ7n1DmA8j2HpwJmfSxGKBJhcGBAYU9EvBG67n1elyCy7u2O7ebBwQcApp2ztze+19O6isO8l3Fp/WU89cwBegLHOFw4wogbYcyNnRX+Zsq5LKcmT3Fq8tSs2w0jFoiVwuD0UNgYitMUbCLij8y73nAkQjabpbGpqfxcNpsluMTJFWZGKBQiFArNGcicc2QnJ6cHwNJwcbpivTJ4+nw+Cvk8E2NjFJYwyURhT0REpMp1Rru4ueWWabNx98b3ruj5egvVEG1kaxpq8zWk8mki/jA1vlosZEQ31ZLIJkhMJopfs6dJTCZIF9Ln3KfDMZIbYSQ3wtE52kR8EeKhpuLQcKlHsDIY1gfq8VnxLkQtbW0c7yleVqZyEklH18p/jmZGKBwmFA5DQ8OsbZxz/Ozxx5lIjZMt5Mi7SbKZHCFfiGBu8YFUYU9ERGQN6Ix2VVW4m8k1QO/gCXwBHxF/hHw2y0guwSXtl7OzcfZLr6TyqXIIHM4mOD15muFsMRCenkwwmkviznOxgnQhTW/6JL3pk7Nu95ufxmAj8WAT8VAcah2p0RSBMR+1oTou3nQJtfULuzXcSjEzsuTI5rJgYD4f5B2TuQyBsMKeiIiIeOjlwitkWwvUjIWxSfCHgmTiKV4uvMJOZg97UX+UqD9KZ6Rz1u25Qo6R3PCZHsFyz2CCxORpEtlhci57zrryLs/Q5BBDk0MwNVHWB5RGf39y6qfYKaM+ECMWqCcWjJWWYxXLZ54P+lb2enrj+TFcwOHDBwXAbxQoMJ6fe1LH+SjsiciS5I7sLy8HdlznYSUi4qVENoE/5OfVuuOkCykivihtoTYmsou/s0XAF6A51EJzaO6ZsmP5sYoweLocBodLz43n574USnk/OJK5keJ9gs89skzUF50lENZPD4fBGFFfdFEzjSeZJOQLFQOpAQ6sYEwyueB9TVHYE5ElyT317fKywp7IxhUgwMGJV4j4o0R8EbIuy6GJg+yqvXDF3tPMqA/UUx+oZytbZ22TKWTKwe++gfvwYaTzaTKFDJnCJJOFDNnz9A5WShVSpDIp+jP952wXsEBFCGyYNRDGAjHqAnX4zX/mdXUhcqkcgXwA8oAfcsEcgWho3jWeVcuiXykiIuvTU//3ube/9ndWpw5ZUxwOw8CVbgnkirNpz3fO3UoL+8JsirSzKdLO0YmjnEqfYoQRzIzGYJyGQAPN4WZe13RDsXcvm2Q0N0oylySZTZLMJRktLY/mRs87w3hKzuXKvYyco3PTMGoDtcQCxRBotZDOTRANRemyzfjyUMgVuLDzkkV/Bgp7IrIk/u3Xel2CLLe6rdBwjt6Y8ZNQO/s5VrJx5cmzs2YnA5OnysO4XZFO8iztvrTLqT3UzoHEfsK+CBFfmFQ+xXA2wRWxK4rX9wud+7qFBVdgIj9RCoEjxVA4IxAWl0fJFDLzqsnhGMuNMZYbA0qTTEqdeD35Y1xRuIILNu9iZ+vi7y+ssCciSxK86p1elyAiVWDqws+76s6Ekon8BLGKC0F7rW+yj63RbSRzSdKFFFF/lE3hTfRN9s3r9T7zUReooy5QRyfn/oMnk8+cFQKTuSSj05ZHGTvHxItC0BFoDfFy4WU6U4ufja2wJyIiIktWrRd+rpTIJmgNt9IWaSs/55xbkdvOhf1hWv2ttIZbz9ku7/KM5kbLIfDR0z9lLDeG3/y0h9upDdRieVvSfZAV9kRERGTJOqNdXFZ3OY+cfoRkboRYoIGbmm6qqmsDVttt5+DMdQAbg40AdCe76YpsnjaTd6n3QVbYExGR6SpmBorM18nUCZ4fe47OaCcX+C4gVUjx/NhztEXaqibw7Y7t5nt936Un30O+kMPvC1Dvr2Nve/X0Pq5EIFXYE5Elyb3yk/Jy4MIbPaxEls0F7/W6AlmDupPdRPzRckiZ+rqU4ccVUbp2nZuaNrzwS+GtqJUYDlfYE5ElyT17b3lZYU9k40pkEzQGGqc9t9Thx+XWneymKdTM5uiW8nMT+YmqCqQrMRzuW8b6REREZIOKB+OkCtMvKOf1+XAzJbIJor7otOeqLZBWDodf1XA1ndFOnh97jpOpE4vep8KeiCyJf9cN5YeIbFy7Y7tJ51NM5CdwzjGRnyCdT7E7ttvr0srWQiCtHA43M2r8NUT8UbqT3Yvep4ZxRWRJgle81esSRKQKdEa7uLnlFrqT3SSyCeLBOHvje6tmeBTWzuVhlns4XGFPRERElkVndPEX/l0NayGQajauiIiIyBJUeyBdid5HnbMnIiIiUiWmeh9r/DUM54ap8ddwc8stSwqo6tkTkSXJvvCj8nLw0ls9rEREZH1Y7t5HhT0RWZL8iw+UlxX2RESqj4ZxRURERNYxc855XYPIhmdm9wItXtcxQwsw6HURc1Bti7NatQ06525bhfdZMDM7BfR4Xcc8VPP/o5lU68o5X73zOtYU9kRkVmZ2wDl3rdd1zEa1LU411ybTraV/K9W6cparXg3jioiIiKxjCnsiIiIi65jCnojM5fNeF3AOqm1xqrk2mW4t/Vup1pWzLPXqnD0RERGRdUw9eyIiIiLrmMKeyAZnZk+a2YOlx9/P2PZxM3uuYvvFq1zbJ83sUTN7wsx+fca2d5jZ/tL2D69mXfOozevP7UMV7/2YmaXNrLFiu6efnRTNPPbMbJeZPWJmD5vZX5uZr9Tuw2Z2oPRv+fZVrvF6M3uwtDzv+swsamZ3ldp+z8xaV7nWq83sRMXn+/5qqdXMgmb25dL77TOz21f8s3XO6aGHHhv0AUSAJ8+x/SvANR7VdjNwD8U/SuuA36/YFgQOAnEgBOwH2quhNq8/t1lq/SvgI9Xy2elR/nc469gDvg3cXFr+G+AOoB14BggDDVPLq1Tj75Te77GF1gd8fOq4AD4A/M9VrvXfAr89o0211PqrwJ+XlpuBYyv92apnT2RjuxKoMbP7zexHZrZ3xvZrgE+W/uL85CrX9haKP9zuphisvlOx7VLgoHMu4ZybBB4BXl8ltYG3n1uZmV0LXO6cqzzJ2+vPTopmO/auAR4qbf8+8HPAHuAnzrmMc26EYlDfvUo1HgLeVbG+kPpuAu6d0Xa1a32bmf3YzP7OzOqrqNZ/Aj5dsZ5jhT9bhT2RjW0C+BOK4eWjwD+YWeU9s79Wev5W4KZVHkJqAa4F3ltRm5W2xYCRirajFP/yrYbawNvPrdKngD+Y8ZzXn50UnXXsUZw0OTVrcurfxbN/L+fcXUC24qmF1Ff5/IrXPEut+4D/6px7A3AY+L0qqnXMOTdaCqB3Av+NFf5sFfZENraXga+4opeBIaADoBRe/tw5N1jqAfoucNUq1jYE3Oecm3TOvQSkgalzU5JAfUXbemC4Gmqrgs+NUh2NwCXOuQdmbPL6s5Oi2Y69TRXbp/5dqunfqzBLHXPVV/m8FzXf7Zx7YmqZ4jFYNbWa2RbgAeDLzrmvssKfrcKeyMb2a8CfAphZJ8W/GHtL22LAs2ZWVwowtwJPzLqXlfEIcJsVdQK1FH8hArwAXGhmTWYWAt4APFoltXn9uU15A/CDWZ73+rOTotmOvfvN7ObS9rcCD1PsoXq9mUXMrIHiMPyzq15t0ZMLqO8nwC/MaLua7jOzPaXlN1E8BquiVjPbBNwPfMI594XS0yv62QbO10BE1rW/A75oZo8AjuIvoPeZWZ1z7vNm9imKf31mgB865763WoU5575jZm+g+APPB/wG8P6K2j4O3Ffa9gXn3Ikqqs2zz63CxRSHrwAws18CPP/spGy2Y28Q+NtSCH8BuNM5lzezv6D4C90H/K5zLu1Rzb893/rM7K+BL5W+v0ngl1a51n8PfM7MJoE+ipOUklVS66coTpD6tJlNnbv3MeAvVuqz1UWVRURERNYxDeOKiIiIrGMKeyIiIiLrmMKeiIiIyDqmsCciIiKyjinsiYiIiKxjCnsiIiIi65jCnoiIiMg6prAnIiIiso4p7ImIiIisYwp7IiIiIuuYwp6IiIjIOqawJyLTmNnNZubM7P0znu82sy/Ocx+XmNmDpeWvlW7uLfNgZh8ys9u9rkNE1o+A1wWISFV6Efgg8HUAM7sCqF3MjpxzH1jGutY959wXva5BRNYX9eyJyGyeBraaWWNp/V8D/zC10czea2aPmtkjZvbHpec6zOxHZvYA8HsVbY+aWcTMXmNm95vZD8zsgJndUNr+ipl9sbS/b5mZv7IQM9tVep8HSu0ePEcNv29mXzKz75vZ82b2ltLzbyy1e8jMvmBmwRnvsaX0mgdKX7eY2dVm9qyZ1ZrZr5rZN8xsu5ntN7Nvm9nPzOwzpdd/0czuMbOfmlnczP6Hmf2kVN97S23+g5k9Xnrus6Xn3lV67hEz+4qZ+Urfw0dL2/+0tP1xM/tYxXv9LzO7r9TbevVy/IOLLFWpV/obZvYdM3vBzD7kdU1SpLAnInP5JnCHmRmwB/gpgJk1AX8AvMk5dxPQZWY/D/w28I/OuVuAb82yv8uB33bO/RzwZ8Cvlp7fCXzaOfc6oBW4bsbrPgv8UWm/PzlPDQAZ59xbgY8Bv1Wq/2+Bdznn3gicAD404z3+BPiL0nv8CfDHzrmfAX8HfAn4TeDXS223l15/HXBrRdj6kXPuBmAvsMM5dyNwC/C7pdD8q8DHSt/nYTMLUOw9/X9K38P9QGyqIDN7O7CjtL+bgF8q9bAC9Djn3gL8JfCRWT5rEa80OOfeDtwO/B9eFyNFCnsiMpevAh8A3gA8XPH8Loqh7HulXrbLKAa2y4F9pTY/mWV/J4BPm9mXgPcAU71rg865V0vLrwKRGa+7lFLQrKhjrhoAnpyxr1agA/hGqe2bga0z3uMK4FOl7f8daCs9/zfAm4B/cM6Nlp572jl32jmXBx4HLi49/1LFvq4p7eve0ve5jWLY+6iZPVRaN+DjwBtKz90AFGZ83w+7oizwWOn7nO17FKkWT5W+6v9mFVHYE5FZOecOUzxP7z8BX6nYdITiD/Kfd87dTLF36XGK5/m9rtRmZu8cwF8Av+ec+zfAMxTDDoA7TynPVux373lqmG1/g8Bx4J2ltp8BHpjR5kXgE6Xt/w64s/T8Z0uPD5nZVJi81MxqSsPN1wPPl54vVOzrgdK+bgW+ARwGPgx8tNS7eBXFcPcR4PdLzxlwR0VNL1Ds0aM07HwD8Moc36NItdD/zSqkCRoici5fB37ZOffyVNhxzp0ysz8DHioFnqMUA82nga+b2QcohrGZvgL8s5n1UwxfLfOs4RPAF8zsvwAjQPYcNZzFOVcone/2XTPzAUngV2Y0+y/AX5tZBIgCH7PijNiLgP8IPErxnMVfASaBfwI2AXc6554ujhSX3QPcbGYPA3XA3c65UTN7BthvZqco9nI+TnHY9l/MbAgYBb5Tej+cc9+x4szoR4EQ8A3n3M9mvJeIyHmZcwrhIlK9zOxfAY875w6a2b8FbnDO/ZpHtWwHvuac23u+tiIi1UI9eyJS7V4FvmZmE0CeMxMlRERkHtSzJyIiIrKOaYKGiIiIyDqmsCciIiKyjinsiYiIiKxjCnsiIiIi65jCnoiIiMg6prAnIiIiso79/xm+ItVpwi5OAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 648x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "make_plots(\n",
    "    df,\n",
    "    \"./plots/gene_expression.pdf\",\n",
    "    intervals_xlabel=\"Median gene expression\",\n",
    "    true_theta=true_quantile_q50,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85119fed-5a1a-408b-aefb-f9208796e43f",
   "metadata": {},
   "source": [
    "### Power experiment\n",
    "\n",
    "For PPI and the classical approach, find the smallest value of ```n``` such that the method has power 80% against the null $H_0: \\text{median gene expression} \\leq 5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "56259bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find n such that we reject H0: gene expression < 5 with probability 80% using a test at level alpha\n",
    "num_experiments = 100\n",
    "list_rand_idx = [\n",
    "    np.random.permutation(n_total) for i in range(num_experiments)\n",
    "]\n",
    "\n",
    "\n",
    "def _to_invert_ppi(n):\n",
    "    n = int(n)\n",
    "    nulls_rejected = 0\n",
    "    # Data setup\n",
    "    for i in range(num_experiments):\n",
    "        rand_idx = list_rand_idx[i]\n",
    "        _Yhat = Yhat_total[rand_idx[:n]]\n",
    "        _Y = Y_total[rand_idx[:n]]\n",
    "        _Yhat_unlabeled = Yhat_total[rand_idx[n:]]\n",
    "\n",
    "        ppi_ci = ppi_quantile_ci(_Y, _Yhat, _Yhat_unlabeled, q, alpha=alpha)\n",
    "        if ppi_ci[0] > 5:\n",
    "            nulls_rejected += 1\n",
    "    return nulls_rejected / num_experiments - 0.8\n",
    "\n",
    "\n",
    "def _to_invert_classical(n):\n",
    "    n = int(n)\n",
    "    nulls_rejected = 0\n",
    "    # Data setup\n",
    "    for i in range(num_experiments):\n",
    "        rand_idx = list_rand_idx[i]\n",
    "        _Y = Y_total[rand_idx[:n]]\n",
    "\n",
    "        classical_ci = classical_quantile_ci(_Y, q, alpha=alpha)\n",
    "\n",
    "        if classical_ci[0] > 5:\n",
    "            nulls_rejected += 1\n",
    "    return nulls_rejected / num_experiments - 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3859fffb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The PPI test requires n=667 labeled data points to reject the null.\n",
      "The classical test requires n=997 labeled data points to reject the null.\n"
     ]
    }
   ],
   "source": [
    "n_ppi = int(brentq(_to_invert_ppi, 100, 2000, xtol=10))\n",
    "n_classical = int(brentq(_to_invert_classical, 100, 2000, xtol=10))\n",
    "print(\n",
    "    f\"The PPI test requires n={n_ppi} labeled data points to reject the null.\"\n",
    ")\n",
    "print(\n",
    "    f\"The classical test requires n={n_classical} labeled data points to reject the null.\"\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}