Ultra large benchmarking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ab217345",
   "metadata": {},
   "source": [
    "# How to evaluate scoring functions used at scale?\n",
    "\n",
    "Area-under the ROC is often used in benchmarking scoring functions. With a short simulation, we'll see here that even phenomenal AUROCs can lead to total failure in a prospective setting when the screening library becomes 'ultra large'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7984478f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import beta, uniform\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import roc_auc_score\n",
    "import numpy as np\n",
    "\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45a4ee0c",
   "metadata": {},
   "source": [
    "we use a beta distribution to model the ranks of the true positives, where 'rank' is a fractional percentile, so 0 is the worst rank, and 1 is the best rank. In this case, our scoring function is quite good, so all the true positives rank pretty highly:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dc1c8617",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Count')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAGwCAYAAABM/qr1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAApy0lEQVR4nO3de3STdZ7H8U96SwFJSsEmrZabN6gCKmiJ6M6AlaJ4Ya2DOtipHhSHaXGgI5eO3ESXsowjjHO4rK4CuyPDimdwFRCsRcCBglpg5SbqwGwrmBbFJoCS3p79w1+zEwGltW3a8n6dk3PI8zxJvg+P2Pd58iS1WZZlCQAAAIoI9wAAAAAtBWEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAARlS4B2gJamtrdeTIEXXs2FE2my3c4wAAgHNgWZaOHz+upKQkRUQ0zrkewkjSkSNHlJycHO4xAABAA5SWluriiy9ulOcijCR17NhR0rd/sQ6HI8zTAACAc+H3+5WcnBz8Od4YCCMp+PaZw+EgjAAAaGUa8zIYLr4GAAAwCCMAAACDMAIAADAIIwAAAIMwAgAAMAgjAAAAgzACAAAwCCMAAACDMAIAADAIIwAAAIMwAgAAMAgjAAAAgzACAAAwCCMAAACDMAIAADCiwj0AAADAd3WfsuYHt6kNfN3or8sZIwAAAIMwAgAAMMIeRocPH9YDDzygzp07q127durTp48++OCD4HrLsjR9+nQlJiaqXbt2SktL0yeffBLyHMeOHdOoUaPkcDgUFxen0aNH68SJE829KwAAoJULaxh99dVXGjRokKKjo/Xmm29q3759+v3vf69OnToFt5k7d66ee+45LV68WNu3b1eHDh2Unp6uU6dOBbcZNWqU9u7dq4KCAq1evVqbN2/WmDFjwrFLAACgFbNZlmWF68WnTJmiLVu26N133z3jesuylJSUpN/85jd6/PHHJUk+n08ul0tLly7Vfffdp/379yslJUXvv/++BgwYIElat26dbrvtNn322WdKSkr6wTn8fr+cTqd8Pp8cDkfj7SAAAGiQc734unT+yEb9+R3WM0avv/66BgwYoJ/97GdKSEjQNddcoxdeeCG4/tChQ/J6vUpLSwsuczqdSk1NVVFRkSSpqKhIcXFxwSiSpLS0NEVERGj79u1nfN1AICC/3x9yAwAACGsYHTx4UIsWLdJll12m9evXa+zYsXrssce0bNkySZLX65UkuVyukMe5XK7gOq/Xq4SEhJD1UVFRio+PD27zXfn5+XI6ncFbcnJyY+8aAABohcIaRrW1tbr22ms1e/ZsXXPNNRozZoweeeQRLV68uElfNy8vTz6fL3grLS1t0tcDAACtQ1jDKDExUSkpKSHLevfurZKSEkmS2+2WJJWVlYVsU1ZWFlzndrtVXl4esr66ulrHjh0LbvNddrtdDocj5AYAABDWMBo0aJAOHDgQsuzjjz9Wt27dJEk9evSQ2+1WYWFhcL3f79f27dvl8XgkSR6PRxUVFSouLg5us2HDBtXW1io1NbUZ9gIAALQVYf2VIBMmTNANN9yg2bNna+TIkXrvvff0/PPP6/nnn5ck2Ww2jR8/Xk8//bQuu+wy9ejRQ9OmTVNSUpJGjBgh6dszTMOGDQu+BVdVVaWcnBzdd9995/SJNAAAgDphDaPrrrtOq1atUl5enmbNmqUePXpo/vz5GjVqVHCbSZMm6eTJkxozZowqKip04403at26dYqNjQ1u8/LLLysnJ0c333yzIiIilJGRoeeeey4cuwQAAFqxsH6PUUvB9xgBANCynJffYwQAANCSEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYIQ1jGbOnCmbzRZy69WrV3D9qVOnlJ2drc6dO+uCCy5QRkaGysrKQp6jpKREw4cPV/v27ZWQkKCJEyequrq6uXcFAAC0AVHhHuDKK6/U22+/HbwfFfX/I02YMEFr1qzRypUr5XQ6lZOTo7vvvltbtmyRJNXU1Gj48OFyu93aunWrPv/8c/3iF79QdHS0Zs+e3ez7AgAAWrewh1FUVJTcbvdpy30+n1588UUtX75cQ4YMkSQtWbJEvXv31rZt2zRw4EC99dZb2rdvn95++225XC5dffXVeuqppzR58mTNnDlTMTExZ3zNQCCgQCAQvO/3+5tm5wAAQKsS9muMPvnkEyUlJalnz54aNWqUSkpKJEnFxcWqqqpSWlpacNtevXqpa9euKioqkiQVFRWpT58+crlcwW3S09Pl9/u1d+/es75mfn6+nE5n8JacnNxEewcAAFqTsIZRamqqli5dqnXr1mnRokU6dOiQbrrpJh0/flxer1cxMTGKi4sLeYzL5ZLX65Ukeb3ekCiqW1+37mzy8vLk8/mCt9LS0sbdMQAA0CqF9a20W2+9Nfjnvn37KjU1Vd26ddMrr7yidu3aNdnr2u122e32Jnt+AADQOoX9rbR/FBcXp8svv1yffvqp3G63KisrVVFREbJNWVlZ8Jokt9t92qfU6u6f6bolAACA79OiwujEiRP629/+psTERPXv31/R0dEqLCwMrj9w4IBKSkrk8XgkSR6PR7t371Z5eXlwm4KCAjkcDqWkpDT7/AAAoHUL61tpjz/+uO644w5169ZNR44c0YwZMxQZGan7779fTqdTo0ePVm5uruLj4+VwODRu3Dh5PB4NHDhQkjR06FClpKQoMzNTc+fOldfr1dSpU5Wdnc1bZQAAoN7CGkafffaZ7r//fn355Ze68MILdeONN2rbtm268MILJUnz5s1TRESEMjIyFAgElJ6eroULFwYfHxkZqdWrV2vs2LHyeDzq0KGDsrKyNGvWrHDtEgAAaMVslmVZ4R4i3Px+v5xOp3w+nxwOR7jHAQDgvNd9ypof3KY28LVK549s1J/fLeoaIwAAgHAijAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAAjBYTRnPmzJHNZtP48eODy06dOqXs7Gx17txZF1xwgTIyMlRWVhbyuJKSEg0fPlzt27dXQkKCJk6cqOrq6maeHgAAtAUtIozef/99/du//Zv69u0bsnzChAl64403tHLlSm3atElHjhzR3XffHVxfU1Oj4cOHq7KyUlu3btWyZcu0dOlSTZ8+vbl3AQAAtAFhD6MTJ05o1KhReuGFF9SpU6fgcp/PpxdffFHPPvushgwZov79+2vJkiXaunWrtm3bJkl66623tG/fPv3pT3/S1VdfrVtvvVVPPfWUFixYoMrKynDtEgAAaKXCHkbZ2dkaPny40tLSQpYXFxerqqoqZHmvXr3UtWtXFRUVSZKKiorUp08fuVyu4Dbp6eny+/3au3fvWV8zEAjI7/eH3AAAAKLC+eIrVqzQjh079P7775+2zuv1KiYmRnFxcSHLXS6XvF5vcJt/jKK69XXrziY/P19PPvnkj5weAAC0NWE7Y1RaWqpf//rXevnllxUbG9usr52Xlyefzxe8lZaWNuvrAwCAlilsYVRcXKzy8nJde+21ioqKUlRUlDZt2qTnnntOUVFRcrlcqqysVEVFRcjjysrK5Ha7JUlut/u0T6nV3a/b5kzsdrscDkfIDQAAIGxhdPPNN2v37t3atWtX8DZgwACNGjUq+Ofo6GgVFhYGH3PgwAGVlJTI4/FIkjwej3bv3q3y8vLgNgUFBXI4HEpJSWn2fQIAAK1b2K4x6tixo6666qqQZR06dFDnzp2Dy0ePHq3c3FzFx8fL4XBo3Lhx8ng8GjhwoCRp6NChSklJUWZmpubOnSuv16upU6cqOztbdru92fcJAAC0bmG9+PqHzJs3TxEREcrIyFAgEFB6eroWLlwYXB8ZGanVq1dr7Nix8ng86tChg7KysjRr1qwwTg0AAForm2VZVriHCDe/3y+n0ymfz8f1RgAAtADdp6z5wW1qA1+rdP7IRv35HfbvMQIAAGgpCCMAAACDMAIAADAIIwAAAIMwAgAAMAgjAAAAgzACAAAwCCMAAACDMAIAADAIIwAAAIMwAgAAMAgjAAAAgzACAAAwCCMAAACDMAIAADAIIwAAAIMwAgAAMAgjAAAAgzACAAAwCCMAAACjQWHUs2dPffnll6ctr6ioUM+ePX/0UAAAAOHQoDD6+9//rpqamtOWBwIBHT58+EcPBQAAEA5R9dn49ddfD/55/fr1cjqdwfs1NTUqLCxU9+7dG204AACA5lSvMBoxYoQkyWazKSsrK2RddHS0unfvrt///veNNhwAAEBzqlcY1dbWSpJ69Oih999/X126dGmSoQAAAMKhXmFU59ChQ409BwAAQNg1KIwkqbCwUIWFhSovLw+eSarz0ksv/ejBAAAAmluDwujJJ5/UrFmzNGDAACUmJspmszX2XAAAAM2uQWG0ePFiLV26VJmZmY09DwAAQNg06HuMKisrdcMNNzT2LAAAAGHVoDB6+OGHtXz58saeBQAAIKwa9FbaqVOn9Pzzz+vtt99W3759FR0dHbL+2WefbZThAAAAmlODwujDDz/U1VdfLUnas2dPyDouxAYAAK1Vg8LonXfeaew5AAAAwq5B1xgBAAC0RQ06YzR48ODvfctsw4YNDR4IAAAgXBoURnXXF9WpqqrSrl27tGfPntN+uSwAAEBr0aAwmjdv3hmXz5w5UydOnPhRAwEAAIRLo15j9MADD/B70gAAQKvVqGFUVFSk2NjYxnxKAACAZtOgt9LuvvvukPuWZenzzz/XBx98oGnTpjXKYAAAAM2tQWHkdDpD7kdEROiKK67QrFmzNHTo0EYZDAAAoLk1KIyWLFnS2HMAAACEXYPCqE5xcbH2798vSbryyit1zTXXNMpQAAAA4dCgMCovL9d9992njRs3Ki4uTpJUUVGhwYMHa8WKFbrwwgsbc0YAANAGdJ+yJtwj/KAGfSpt3LhxOn78uPbu3atjx47p2LFj2rNnj/x+vx577LHGnhEAAKBZNOiM0bp16/T222+rd+/ewWUpKSlasGABF18DAIBWq0FnjGpraxUdHX3a8ujoaNXW1p7z8yxatEh9+/aVw+GQw+GQx+PRm2++GVx/6tQpZWdnq3PnzrrggguUkZGhsrKykOcoKSnR8OHD1b59eyUkJGjixImqrq5uyG4BAIDzXIPCaMiQIfr1r3+tI0eOBJcdPnxYEyZM0M0333zOz3PxxRdrzpw5Ki4u1gcffKAhQ4borrvu0t69eyVJEyZM0BtvvKGVK1dq06ZNOnLkSMh3KNXU1Gj48OGqrKzU1q1btWzZMi1dulTTp09vyG4BAIDznM2yLKu+DyotLdWdd96pvXv3Kjk5Objsqquu0uuvv66LL764wQPFx8frd7/7ne655x5deOGFWr58ue655x5J0kcffaTevXurqKhIAwcO1Jtvvqnbb79dR44ckcvlkiQtXrxYkydP1tGjRxUTE3PG1wgEAgoEAsH7fr9fycnJ8vl8cjgcDZ4dAACcXWNffF0b+Fql80c26s/vBp0xSk5O1o4dO7RmzRqNHz9e48eP19q1a7Vjx44GR1FNTY1WrFihkydPyuPxqLi4WFVVVUpLSwtu06tXL3Xt2lVFRUWSvv0VJH369AlGkSSlp6fL7/cHzzqdSX5+vpxOZ/BWF3cAAOD8Vq8w2rBhg1JSUuT3+2Wz2XTLLbdo3LhxGjdunK677jpdeeWVevfdd+s1wO7du3XBBRfIbrfrl7/8pVatWqWUlBR5vV7FxMQEvw6gjsvlktfrlSR5vd6QKKpbX7fubPLy8uTz+YK30tLSes0MAADapnp9Km3+/Pl65JFHzni6yul06tFHH9Wzzz6rm2666Zyf84orrtCuXbvk8/n06quvKisrS5s2barPWPVmt9tlt9ub9DUAAEDrU68zRv/zP/+jYcOGnXX90KFDVVxcXK8BYmJidOmll6p///7Kz89Xv3799Ic//EFut1uVlZWqqKgI2b6srExut1uS5Ha7T/uUWt39um0AAADOVb3CqKys7Iwf068TFRWlo0eP/qiBamtrFQgE1L9/f0VHR6uwsDC47sCBAyopKZHH45EkeTwe7d69W+Xl5cFtCgoK5HA4lJKS8qPmAAAA5596vZV20UUXac+ePbr00kvPuP7DDz9UYmLiOT9fXl6ebr31VnXt2lXHjx/X8uXLtXHjRq1fv15Op1OjR49Wbm6u4uPj5XA4NG7cOHk8Hg0cOFDSt2eoUlJSlJmZqblz58rr9Wrq1KnKzs7mrTIAAFBv9Qqj2267TdOmTdOwYcMUGxsbsu6bb77RjBkzdPvtt5/z85WXl+sXv/iFPv/8czmdTvXt21fr16/XLbfcIkmaN2+eIiIilJGRoUAgoPT0dC1cuDD4+MjISK1evVpjx46Vx+NRhw4dlJWVpVmzZtVntwAAACTV83uMysrKdO211yoyMlI5OTm64oorJH37/UILFixQTU2NduzYcdonxVo6v98vp9PJ9xgBANCEWsP3GNXrjJHL5dLWrVs1duxY5eXlqa6pbDab0tPTtWDBglYXRQAAAHXq/Utku3XrprVr1+qrr77Sp59+KsuydNlll6lTp05NMR8AAECzqXcY1enUqZOuu+66xpwFAAAgrBr0K0EAAADaIsIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAADMIIAADAIIwAAAAMwggAAMAgjAAAAAzCCAAAwCCMAAAAjLCGUX5+vq677jp17NhRCQkJGjFihA4cOBCyzalTp5Sdna3OnTvrggsuUEZGhsrKykK2KSkp0fDhw9W+fXslJCRo4sSJqq6ubs5dAQAAbUBYw2jTpk3Kzs7Wtm3bVFBQoKqqKg0dOlQnT54MbjNhwgS98cYbWrlypTZt2qQjR47o7rvvDq6vqanR8OHDVVlZqa1bt2rZsmVaunSppk+fHo5dAgAArZjNsiwr3EPUOXr0qBISErRp0yb90z/9k3w+ny688EItX75c99xzjyTpo48+Uu/evVVUVKSBAwfqzTff1O23364jR47I5XJJkhYvXqzJkyfr6NGjiomJ+cHX9fv9cjqd8vl8cjgcTbqPAACcr7pPWdOoz1cb+Fql80c26s/vFnWNkc/nkyTFx8dLkoqLi1VVVaW0tLTgNr169VLXrl1VVFQkSSoqKlKfPn2CUSRJ6enp8vv92rt37xlfJxAIyO/3h9wAAABaTBjV1tZq/PjxGjRokK666ipJktfrVUxMjOLi4kK2dblc8nq9wW3+MYrq1tetO5P8/Hw5nc7gLTk5uZH3BgAAtEYtJoyys7O1Z88erVixoslfKy8vTz6fL3grLS1t8tcEAAAtX1S4B5CknJwcrV69Wps3b9bFF18cXO52u1VZWamKioqQs0ZlZWVyu93Bbd57772Q56v71FrdNt9lt9tlt9sbeS8AAEBrF9YzRpZlKScnR6tWrdKGDRvUo0ePkPX9+/dXdHS0CgsLg8sOHDigkpISeTweSZLH49Hu3btVXl4e3KagoEAOh0MpKSnNsyMAAKBNCOsZo+zsbC1fvlz//d//rY4dOwavCXI6nWrXrp2cTqdGjx6t3NxcxcfHy+FwaNy4cfJ4PBo4cKAkaejQoUpJSVFmZqbmzp0rr9erqVOnKjs7m7NCAACgXsIaRosWLZIk/fSnPw1ZvmTJEj344IOSpHnz5ikiIkIZGRkKBAJKT0/XwoULg9tGRkZq9erVGjt2rDwejzp06KCsrCzNmjWruXYDAAC0ES3qe4zChe8xAgCg6fE9RgAAAK0IYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAARljDaPPmzbrjjjuUlJQkm82m1157LWS9ZVmaPn26EhMT1a5dO6WlpemTTz4J2ebYsWMaNWqUHA6H4uLiNHr0aJ04caIZ9wIAALQVYQ2jkydPql+/flqwYMEZ18+dO1fPPfecFi9erO3bt6tDhw5KT0/XqVOngtuMGjVKe/fuVUFBgVavXq3NmzdrzJgxzbULAACgDbFZlmWFewhJstlsWrVqlUaMGCHp27NFSUlJ+s1vfqPHH39ckuTz+eRyubR06VLdd9992r9/v1JSUvT+++9rwIABkqR169bptttu02effaakpKQzvlYgEFAgEAje9/v9Sk5Ols/nk8PhaNodBQDgPNV9yppGfb7awNcqnT+yUX9+t9hrjA4dOiSv16u0tLTgMqfTqdTUVBUVFUmSioqKFBcXF4wiSUpLS1NERIS2b99+1ufOz8+X0+kM3pKTk5tuRwAAQKvRYsPI6/VKklwuV8hyl8sVXOf1epWQkBCyPioqSvHx8cFtziQvL08+ny94Ky0tbeTpAQBAaxQV7gHCwW63y263h3sMAADQwrTYM0Zut1uSVFZWFrK8rKwsuM7tdqu8vDxkfXV1tY4dOxbcBgAA4Fy12DDq0aOH3G63CgsLg8v8fr+2b98uj8cjSfJ4PKqoqFBxcXFwmw0bNqi2tlapqanNPjMAAGjdwvpW2okTJ/Tpp58G7x86dEi7du1SfHy8unbtqvHjx+vpp5/WZZddph49emjatGlKSkoKfnKtd+/eGjZsmB555BEtXrxYVVVVysnJ0X333XfWT6QBAACcTVjD6IMPPtDgwYOD93NzcyVJWVlZWrp0qSZNmqSTJ09qzJgxqqio0I033qh169YpNjY2+JiXX35ZOTk5uvnmmxUREaGMjAw999xzzb4vAACg9Wsx32MUTn6/X06nk+8xAgCgCfE9RgAAAK0IYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAABGVLgHAAAArVv3KWvCPUKj4YwRAACAQRgBAAAYhBEAAIBBGAEAABiEEQAAgEEYAQAAGIQRAACAQRgBAAAYhBEAAIBBGAEAABiEEQAAgEEYAQAAGIQRAACAQRgBAAAYhBEAAIBBGAEAABiEEQAAgEEYAQAAGIQRAACAQRgBAAAYhBEAAIBBGAEAABiEEQAAgEEYAQAAGIQRAACAQRgBAAAYhBEAAIBBGAEAABiEEQAAgBEV7gEAAEDL1H3KmnCP0Ow4YwQAAGAQRgAAAEabCaMFCxaoe/fuio2NVWpqqt57771wjwQAAFqZNhFG//Vf/6Xc3FzNmDFDO3bsUL9+/ZSenq7y8vJwjwYAAFqRNhFGzz77rB555BE99NBDSklJ0eLFi9W+fXu99NJL4R4NAAC0Iq3+U2mVlZUqLi5WXl5ecFlERITS0tJUVFR0xscEAgEFAoHgfZ/PJ0ny+/1NOywAAA1w1Yz157TdnifTG/X5WrrawNeSJMuyGu05W30YffHFF6qpqZHL5QpZ7nK59NFHH53xMfn5+XryySdPW56cnNwkMwIA0Byc88M9QXh8+eWXcjqdjfJcrT6MGiIvL0+5ubnB+xUVFerWrZtKSkoa7S8WDeP3+5WcnKzS0lI5HI5wj3Ne41i0HByLloNj0bL4fD517dpV8fHxjfacrT6MunTposjISJWVlYUsLysrk9vtPuNj7Ha77Hb7acudTif/obcQDoeDY9FCcCxaDo5Fy8GxaFkiIhrvkulWf/F1TEyM+vfvr8LCwuCy2tpaFRYWyuPxhHEyAADQ2rT6M0aSlJubq6ysLA0YMEDXX3+95s+fr5MnT+qhhx4K92gAAKAVaRNhdO+99+ro0aOaPn26vF6vrr76aq1bt+60C7LPxm63a8aMGWd8ew3Ni2PRcnAsWg6ORcvBsWhZmuJ42KzG/IwbAABAK9bqrzECAABoLIQRAACAQRgBAAAYhBEAAIBx3oTRggUL1L17d8XGxio1NVXvvffe926/cuVK9erVS7GxserTp4/Wrl3bTJO2ffU5Fi+88IJuuukmderUSZ06dVJaWtoPHjucu/r+u6izYsUK2Ww2jRgxomkHPI/U91hUVFQoOztbiYmJstvtuvzyy/n/VCOp77GYP3++rrjiCrVr107JycmaMGGCTp061UzTtl2bN2/WHXfcoaSkJNlsNr322ms/+JiNGzfq2muvld1u16WXXqqlS5fW/4Wt88CKFSusmJgY66WXXrL27t1rPfLII1ZcXJxVVlZ2xu23bNliRUZGWnPnzrX27dtnTZ061YqOjrZ2797dzJO3PfU9Fj//+c+tBQsWWDt37rT2799vPfjgg5bT6bQ+++yzZp687anvsahz6NAh66KLLrJuuukm66677mqeYdu4+h6LQCBgDRgwwLrtttusv/71r9ahQ4esjRs3Wrt27Wrmydue+h6Ll19+2bLb7dbLL79sHTp0yFq/fr2VmJhoTZgwoZknb3vWrl1rPfHEE9Zf/vIXS5K1atWq793+4MGDVvv27a3c3Fxr37591h//+EcrMjLSWrduXb1e97wIo+uvv97Kzs4O3q+pqbGSkpKs/Pz8M24/cuRIa/jw4SHLUlNTrUcffbRJ5zwf1PdYfFd1dbXVsWNHa9myZU014nmjIceiurrauuGGG6x///d/t7KysgijRlLfY7Fo0SKrZ8+eVmVlZXONeN6o77HIzs62hgwZErIsNzfXGjRoUJPOeb45lzCaNGmSdeWVV4Ysu/fee6309PR6vVabfyutsrJSxcXFSktLCy6LiIhQWlqaioqKzviYoqKikO0lKT09/azb49w05Fh819dff62qqqpG/YWB56OGHotZs2YpISFBo0ePbo4xzwsNORavv/66PB6PsrOz5XK5dNVVV2n27NmqqalprrHbpIYcixtuuEHFxcXBt9sOHjyotWvX6rbbbmuWmfH/Gutnd5v45uvv88UXX6impua0b8F2uVz66KOPzvgYr9d7xu29Xm+TzXk+aMix+K7JkycrKSnptP/4UT8NORZ//etf9eKLL2rXrl3NMOH5oyHH4uDBg9qwYYNGjRqltWvX6tNPP9WvfvUrVVVVacaMGc0xdpvUkGPx85//XF988YVuvPFGWZal6upq/fKXv9Rvf/vb5hgZ/+BsP7v9fr+++eYbtWvX7pyep82fMULbMWfOHK1YsUKrVq1SbGxsuMc5rxw/flyZmZl64YUX1KVLl3CPc96rra1VQkKCnn/+efXv31/33nuvnnjiCS1evDjco513Nm7cqNmzZ2vhwoXasWOH/vKXv2jNmjV66qmnwj0aGqjNnzHq0qWLIiMjVVZWFrK8rKxMbrf7jI9xu9312h7npiHHos4zzzyjOXPm6O2331bfvn2bcszzQn2Pxd/+9jf9/e9/1x133BFcVltbK0mKiorSgQMHdMkllzTt0G1UQ/5dJCYmKjo6WpGRkcFlvXv3ltfrVWVlpWJiYpp05raqIcdi2rRpyszM1MMPPyxJ6tOnj06ePKkxY8boiSeeUEQE5x+ay9l+djscjnM+WySdB2eMYmJi1L9/fxUWFgaX1dbWqrCwUB6P54yP8Xg8IdtLUkFBwVm3x7lpyLGQpLlz5+qpp57SunXrNGDAgOYYtc2r77Ho1auXdu/erV27dgVvd955pwYPHqxdu3YpOTm5OcdvUxry72LQoEH69NNPg3EqSR9//LESExOJoh+hIcfi66+/Pi1+6oLV4leRNqtG+9ldv+vCW6cVK1ZYdrvdWrp0qbVv3z5rzJgxVlxcnOX1ei3LsqzMzExrypQpwe23bNliRUVFWc8884y1f/9+a8aMGXxcv5HU91jMmTPHiomJsV599VXr888/D96OHz8erl1oM+p7LL6LT6U1nvoei5KSEqtjx45WTk6OdeDAAWv16tVWQkKC9fTTT4drF9qM+h6LGTNmWB07drT+/Oc/WwcPHrTeeust65JLLrFGjhwZrl1oM44fP27t3LnT2rlzpyXJevbZZ62dO3da//u//2tZlmVNmTLFyszMDG5f93H9iRMnWvv377cWLFjAx/W/zx//+Eera9euVkxMjHX99ddb27ZtC677yU9+YmVlZYVs/8orr1iXX365FRMTY1155ZXWmjVrmnnitqs+x6Jbt26WpNNuM2bMaP7B26D6/rv4R4RR46rvsdi6dauVmppq2e12q2fPnta//Mu/WNXV1c08ddtUn2NRVVVlzZw507rkkkus2NhYKzk52frVr35lffXVV80/eBvzzjvvnPH//3V//1lZWdZPfvKT0x5z9dVXWzExMVbPnj2tJUuW1Pt1bZbFuT4AAADpPLjGCAAA4FwRRgAAAAZhBAAAYBBGAAAABmEEAABgEEYAAAAGYQQAAGAQRgAAAAZhBACGzWbTa6+9Fu4xAIQRYQSgxXvwwQdls9lks9kUHR2tHj16aNKkSTp16lS4RwPQxkSFewAAOBfDhg3TkiVLVFVVpeLiYmVlZclms+lf//Vfwz0agDaEM0YAWgW73S63263k5GSNGDFCaWlpKigokCR9+eWXuv/++3XRRRepffv26tOnj/785z+HPP6nP/2pHnvsMU2aNEnx8fFyu92aOXPm977mjBkzlJiYqA8//LCpdgtAC0MYAWh19uzZo61btyomJkaSdOrUKfXv319r1qzRnj17NGbMGGVmZuq9994LedyyZcvUoUMHbd++XXPnztWsWbOCcfWPLMvSuHHj9B//8R9699131bdv32bZLwDhZ7Msywr3EADwfR588EH96U9/UmxsrKqrqxUIBBQREaFXXnlFGRkZZ3zM7bffrl69eumZZ56R9O0Zo5qaGr377rvBba6//noNGTJEc+bMkfTtxdcrV67UqlWrtHPnThUUFOiiiy5q+h0E0GJwjRGAVmHw4MFatGiRTp48qXnz5ikqKioYRTU1NZo9e7ZeeeUVHT58WJWVlQoEAmrfvn3Ic3z3zE9iYqLKy8tDlk2YMEF2u13btm1Tly5dmnanALQ4vJUGoFXo0KGDLr30UvXr108vvfSStm/frhdffFGS9Lvf/U5/+MMfNHnyZL3zzjvatWuX0tPTVVlZGfIc0dHRIfdtNptqa2tDlt1yyy06fPiw1q9f37Q7BKBFIowAtDoRERH67W9/q6lTp+qbb77Rli1bdNddd+mBBx5Qv3791LNnT3388ccNeu4777xTy5cv18MPP6wVK1Y08uQAWjrCCECr9LOf/UyRkZFasGCBLrvsMhUUFGjr1q3av3+/Hn30UZWVlTX4uf/5n/9Z//mf/6mHHnpIr776aiNODaCl4xojAK1SVFSUcnJyNHfuXO3cuVMHDx5Uenq62rdvrzFjxmjEiBHy+XwNfv577rlHtbW1yszMVEREhO6+++5GnB5AS8Wn0gAAAAzeSgMAADAIIwAAAIMwAgAAMAgjAAAAgzACAAAwCCMAAACDMAIAADAIIwAAAIMwAgAAMAgjAAAAgzACAAAw/g9wjoUOqDo0pwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "active_ranks_dist = beta(100, 2) #models rank of true positives\n",
    "\n",
    "plt.hist(active_ranks_dist.rvs(1000),bins=np.linspace(0, 1));\n",
    "plt.xlim(0, 1)\n",
    "plt.xlabel('Rank')\n",
    "plt.ylabel('Count')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b49eb3f2",
   "metadata": {},
   "source": [
    "because we're modelling ranks, and not scores, then the background distribution of decoys can be modelled with a uniform distribution:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5998288f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Count')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAGwCAYAAABCV9SaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkN0lEQVR4nO3dfVhUdf7/8dcgMGgBhiY3NSpqea9taka3aiRqd6ZtdudiW9kN2ibXrkVZaLXhWqu2LelVW1q7matdaXemq5R2I1qRrFFqedNqKpj2lVHMAeH8/ujnXE2iwTjDmZnP83Fdc13OmTOH93hQnteZwxmHZVmWAAAAIlyU3QMAAAA0BaIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEaItnuAYKurq9OuXbsUHx8vh8Nh9zgAAKABLMvSgQMHlJaWpqiowByjifjo2bVrl1wul91jAAAAP+zYsUNnnnlmQLYV8dETHx8v6ae/tISEBJunAQAADeF2u+Vyubw/xwMh4qPn6FtaCQkJRA8AAGEmkKemcCIzAAAwAtEDAACMQPQAAAAjED0AAMAIRA8AADAC0QMAAIxA9AAAACMQPQAAwAhEDwAAMALRAwAAjED0AAAAIxA9AADACEQPAAAwAtEDAACMQPQAAAAjRNs9AADYpf0D7zRovW+nXhHkSQA0BY70AAAAIxA9AADACEQPAAAwAtEDAACMQPQAAAAjED0AAMAIRA8AADAC0QMAAIxA9AAAACMQPQAAwAhEDwAAMALRAwAAjED0AAAAIxA9AADACEQPAAAwAtEDAACMQPQAAAAjED0AAMAIRA8AADAC0QMAAIwQbfcAABDq2j/wToPW+3bqFUGeBMDJ4EgPAAAwAtEDAACMQPQAAAAjED0AAMAIRA8AADAC0QMAAIxga/TMmjVLvXr1UkJCghISEpSRkaF3333X+/jhw4eVk5OjVq1a6dRTT9XIkSNVUVFh48QAACBc2Ro9Z555pqZOnaqSkhJ99tlnGjRokK655hp9+eWXkqQJEyborbfe0sKFC7Vq1Srt2rVLI0aMsHNkAAAQphyWZVl2D/FzSUlJevLJJ3Xdddfp9NNP17x583TddddJkjZu3KiuXbuquLhY559/foO253a7lZiYqMrKSiUkJARzdABhpqEXHWwoLk4IBE4wfn6HzDk9tbW1mj9/vqqqqpSRkaGSkhLV1NQoMzPTu06XLl3Utm1bFRcXH3c7Ho9Hbrfb5wYAAGB79HzxxRc69dRT5XQ6ddddd2nRokXq1q2bysvLFRsbq5YtW/qsn5ycrPLy8uNur6CgQImJid6by+UK8isAAADhwPbo6dy5s0pLS7V27Vrdfffdys7O1ldffeX39vLy8lRZWem97dixI4DTAgCAcGX7B47GxsaqU6dOkqQ+ffro008/1dNPP61Ro0apurpa+/fv9znaU1FRoZSUlONuz+l0yul0BntsAAAQZmw/0vNLdXV18ng86tOnj2JiYlRUVOR9bNOmTdq+fbsyMjJsnBAAAIQjW4/05OXlaejQoWrbtq0OHDigefPmaeXKlVq2bJkSExN12223KTc3V0lJSUpISND48eOVkZHR4N/cAgAAOMrW6NmzZ49+97vfaffu3UpMTFSvXr20bNkyXX755ZKkGTNmKCoqSiNHjpTH41FWVpaeffZZO0cGAABhKuSu0xNoXKcHwPFwnR4gdEX0dXoAAACCiegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGiLZ7AAAAgqn9A+80aL1vp14R5ElgN470AAAAIxA9AADACEQPAAAwAtEDAACMQPQAAAAjED0AAMAIRA8AADAC0QMAAIxA9AAAACMQPQAAwAh8DAWAsMHHCQA4GRzpAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEbgiswAGqWhV0WWuDIygqcx34fAURzpAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARbI2egoIC9evXT/Hx8WrTpo2GDx+uTZs2+awzYMAAORwOn9tdd91l08QAACBc2Ro9q1atUk5OjtasWaPly5erpqZGgwcPVlVVlc96d9xxh3bv3u29TZs2zaaJAQBAuLL14oRLly71uT937ly1adNGJSUluuSSS7zLW7RooZSUlKYeDwAARJCQOqensrJSkpSUlOSz/JVXXlHr1q3Vo0cP5eXl6dChQ8fdhsfjkdvt9rkBAACEzMdQ1NXV6b777tOFF16oHj16eJffdNNNateundLS0rR+/Xrdf//92rRpk15//fV6t1NQUKApU6Y01dgIsoZeap6POwAA/JqQiZ6cnByVlZXpo48+8lk+duxY75979uyp1NRUXXbZZdqyZYs6dux4zHby8vKUm5vrve92u+VyuYI3OAAACAshET3jxo3T22+/rQ8++EBnnnnmCdft37+/JGnz5s31Ro/T6ZTT6QzKnAAAIHzZGj2WZWn8+PFatGiRVq5cqfT09F99TmlpqSQpNTU1yNMBAIBIYmv05OTkaN68eXrjjTcUHx+v8vJySVJiYqKaN2+uLVu2aN68eRo2bJhatWql9evXa8KECbrkkkvUq1cvO0cHAABhxtbomTVrlqSfLkD4c3PmzNGYMWMUGxurFStWaObMmaqqqpLL5dLIkSM1adIkG6YFAADhzPa3t07E5XJp1apVTTQNAACIZCF1nR4AAIBgIXoAAIARiB4AAGAEogcAABghJC5OCMB+Df3ID5iBj4AJPeyTk8eRHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAErsgMAEAE4crNx8eRHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiBj6EIM4G+vDiXK286Jv5d2/WaG/p1ATRMpPz/xZEeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGMGYj6Hokb9MUc4WJ1wn1C+fjeOLlEukI7yZ+H0YSa+Zj/lpOg35u6nzHAr41+VIDwAAMALRAwAAjED0AAAAIxA9AADACEQPAAAwAtEDAACMYGv0FBQUqF+/foqPj1ebNm00fPhwbdq0yWedw4cPKycnR61atdKpp56qkSNHqqKiwqaJAQBAuLI1elatWqWcnBytWbNGy5cvV01NjQYPHqyqqirvOhMmTNBbb72lhQsXatWqVdq1a5dGjBhh49QAACAc2XpxwqVLl/rcnzt3rtq0aaOSkhJdcsklqqys1AsvvKB58+Zp0KBBkqQ5c+aoa9euWrNmjc4///xjtunxeOTxeLz33W53cF8EAAAICyF1RebKykpJUlJSkiSppKRENTU1yszM9K7TpUsXtW3bVsXFxfVGT0FBgaZMmdI0A+MYXIE09DR0nwBomEj5NxUpr6MxQuZE5rq6Ot1333268MIL1aNHD0lSeXm5YmNj1bJlS591k5OTVV5eXu928vLyVFlZ6b3t2LEj2KMDAIAwEDJHenJyclRWVqaPPvropLbjdDrldDoDNBUAAIgUIXGkZ9y4cXr77bf1/vvv68wzz/QuT0lJUXV1tfbv3++zfkVFhVJSUpp4SgAAEM5sjR7LsjRu3DgtWrRI7733ntLT030e79Onj2JiYlRUVORdtmnTJm3fvl0ZGRlNPS4AAAhjtr69lZOTo3nz5umNN95QfHy89zydxMRENW/eXImJibrtttuUm5urpKQkJSQkaPz48crIyKj3JGYAAIDjsTV6Zs2aJUkaMGCAz/I5c+ZozJgxkqQZM2YoKipKI0eOlMfjUVZWlp599tkmnhQAAIQ7W6PHsqxfXScuLk6FhYUqLCxsgokAAECkCokTmQEAAIKN6AEAAEYgegAAgBH8OqenQ4cO+vTTT9WqVSuf5fv379e5556rrVu3BmQ4+M/Ey4vbhb9rhIJQ/z7kI2rMEOrfh34d6fn2229VW1t7zHKPx6OdO3ee9FAAAACB1qgjPW+++ab3z8uWLVNiYqL3fm1trYqKitS+ffuADQcAABAojYqe4cOHS5IcDoeys7N9HouJiVH79u3117/+NWDDAQAABEqjoqeurk6SlJ6erk8//VStW7cOylAAAACB5teJzNu2bQv0HAAAAEHl9xWZi4qKVFRUpD179niPAB314osvnvRgAAAAgeRX9EyZMkWPPvqo+vbtq9TUVDkcjkDPBQAAEFB+Rc/s2bM1d+5cjR49OtDzAAAABIVf1+mprq7WBRdcEOhZAAAAgsav6Ln99ts1b968QM8CAAAQNH69vXX48GE999xzWrFihXr16qWYmBifx6dPnx6Q4UJVoC+nHuqX7Q4GE18zAMBefkXP+vXrdc4550iSysrKfB7jpGYAABCK/Iqe999/P9BzAAAABJVf5/QAAACEG7+O9AwcOPCEb2O99957fg8EAAAQDH5Fz9HzeY6qqalRaWmpysrKjvkgUgAAgFDgV/TMmDGj3uWTJ0/WwYMHT2ogAACAYAjoOT233HILn7sFAABCUkCjp7i4WHFxcYHcJAAAQED49fbWiBEjfO5blqXdu3frs88+08MPPxyQwQAAAALJr+hJTEz0uR8VFaXOnTvr0Ucf1eDBgwMyGGA3rhqNYGnM91ZDr+wO4Nf5FT1z5swJ9BwAAABB5Vf0HFVSUqINGzZIkrp3767f/OY3ARkKAAAg0PyKnj179uiGG27QypUr1bJlS0nS/v37NXDgQM2fP1+nn356IGcEAAA4aX799tb48eN14MABffnll/rhhx/0ww8/qKysTG63W/fee2+gZwQAADhpfh3pWbp0qVasWKGuXbt6l3Xr1k2FhYWcyAwAAEKSX0d66urqFBMTc8zymJgY1dXVnfRQAAAAgeZX9AwaNEh/+MMftGvXLu+ynTt3asKECbrssssCNhwAAECg+BU9f//73+V2u9W+fXt17NhRHTt2VHp6utxut5555plAzwgAAHDS/Dqnx+Vy6fPPP9eKFSu0ceNGSVLXrl2VmZkZ0OEAAAACpVFHet577z1169ZNbrdbDodDl19+ucaPH6/x48erX79+6t69uz788MNgzQoAAOC3Rh3pmTlzpu644w4lJCQc81hiYqLuvPNOTZ8+XRdffHHABgRM09CPKODjCcxg2sehmPZ6JTNfs10adaTnv//9r4YMGXLcxwcPHqySkpKTHgoAACDQGhU9FRUV9f6q+lHR0dH6/vvvT3ooAACAQGtU9JxxxhkqKys77uPr169XamrqSQ8FAAAQaI2KnmHDhunhhx/W4cOHj3nsxx9/VH5+vq688sqADQcAABAojTqRedKkSXr99dd19tlna9y4cercubMkaePGjSosLFRtba0eeuihoAwKAABwMhoVPcnJyVq9erXuvvtu5eXlybIsSZLD4VBWVpYKCwuVnJwclEEBAABORqOvyNyuXTstWbJEe/fu1dq1a7VmzRrt3btXS5YsUXp6eqO29cEHH+iqq65SWlqaHA6HFi9e7PP4mDFj5HA4fG4n+u0xAACA4/HrisySdNppp6lfv34n9cWrqqrUu3dv/f73v9eIESPqXWfIkCGaM2eO977T6TyprwkAAMzkd/QEwtChQzV06NATruN0OpWSktJEEwEAgEjl1weONqWVK1eqTZs26ty5s+6++27t27fvhOt7PB653W6fGwAAgK1Hen7NkCFDNGLECKWnp2vLli168MEHNXToUBUXF6tZs2b1PqegoEBTpkzx6+txKXAAACJXSEfPDTfc4P1zz5491atXL3Xs2FErV67UZZddVu9z8vLylJub673vdrvlcrmCPisAAAhtIf/21s916NBBrVu31ubNm4+7jtPpVEJCgs8NAAAgrKLnu+++0759+/ioCwAA0Gi2vr118OBBn6M227ZtU2lpqZKSkpSUlKQpU6Zo5MiRSklJ0ZYtWzRx4kR16tRJWVlZNk4NAADCka3R89lnn2ngwIHe+0fPxcnOztasWbO0fv16vfTSS9q/f7/S0tI0ePBgPfbYY1yrBwAANJqt0TNgwADvR1nUZ9myZU04DQAAiGRhdU4PAACAv4geAABgBKIHAAAYIaQvThjuIukKz5HyWiLldUiR9VoAoClwpAcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABG4GMoANiOj9QA0BQ40gMAAIxA9AAAACMQPQAAwAhEDwAAMALRAwAAjED0AAAAIxA9AADACEQPAAAwAtEDAACMQPQAAAAjED0AAMAIRA8AADAC0QMAAIxA9AAAACMQPQAAwAhEDwAAMALRAwAAjED0AAAAIxA9AADACEQPAAAwAtEDAACMQPQAAAAjED0AAMAIRA8AADAC0QMAAIxA9AAAACMQPQAAwAi2Rs8HH3ygq666SmlpaXI4HFq8eLHP45Zl6ZFHHlFqaqqaN2+uzMxMffPNN/YMCwAAwpqt0VNVVaXevXursLCw3senTZumv/3tb5o9e7bWrl2rU045RVlZWTp8+HATTwoAAMJdtJ1ffOjQoRo6dGi9j1mWpZkzZ2rSpEm65pprJEkvv/yykpOTtXjxYt1www1NOSoAAAhzIXtOz7Zt21ReXq7MzEzvssTERPXv31/FxcXHfZ7H45Hb7fa5AQAAhGz0lJeXS5KSk5N9licnJ3sfq09BQYESExO9N5fLFdQ5AQBAeAjZ6PFXXl6eKisrvbcdO3bYPRIAAAgBIRs9KSkpkqSKigqf5RUVFd7H6uN0OpWQkOBzAwAACNnoSU9PV0pKioqKirzL3G631q5dq4yMDBsnAwAA4cjW3946ePCgNm/e7L2/bds2lZaWKikpSW3bttV9992nxx9/XGeddZbS09P18MMPKy0tTcOHD7dvaAAAEJZsjZ7PPvtMAwcO9N7Pzc2VJGVnZ2vu3LmaOHGiqqqqNHbsWO3fv18XXXSRli5dqri4OLtGBgAAYcphWZZl9xDB5Ha7f/otrvsWKMrZwu5xAABAA9R5DmnHzOtVWVkZsPNzQ/acHgAAgEAiegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABghJCOnsmTJ8vhcPjcunTpYvdYAAAgDEXbPcCv6d69u1asWOG9Hx0d8iMDAIAQFPIFER0drZSUlAav7/F45PF4vPfdbncwxgIAAGEmpN/ekqRvvvlGaWlp6tChg26++WZt3779hOsXFBQoMTHRe3O5XE00KQAACGUOy7Isu4c4nnfffVcHDx5U586dtXv3bk2ZMkU7d+5UWVmZ4uPj631OfUd6XC6XXPctUJSzRVONDgAATkKd55B2zLxelZWVSkhICMg2Q/rtraFDh3r/3KtXL/Xv31/t2rXTggULdNttt9X7HKfTKafT2VQjAgCAMBHyb2/9XMuWLXX22Wdr8+bNdo8CAADCTFhFz8GDB7VlyxalpqbaPQoAAAgzIR09f/zjH7Vq1Sp9++23Wr16ta699lo1a9ZMN954o92jAQCAMBPS5/R89913uvHGG7Vv3z6dfvrpuuiii7RmzRqdfvrpdo8GAADCTEhHz/z58+0eAQAARIiQfnsLAAAgUIgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGAEogcAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIARiB4AAGCEsIiewsJCtW/fXnFxcerfv78++eQTu0cCAABhJuSj59///rdyc3OVn5+vzz//XL1791ZWVpb27Nlj92gAACCMhHz0TJ8+XXfccYduvfVWdevWTbNnz1aLFi304osv2j0aAAAII9F2D3Ai1dXVKikpUV5enndZVFSUMjMzVVxcXO9zPB6PPB6P935lZaUkqc5zKLjDAgCAgDn6c9uyrIBtM6SjZ+/evaqtrVVycrLP8uTkZG3cuLHe5xQUFGjKlCnHLN85a0wwRgQAAEG0b98+JSYmBmRbIR09/sjLy1Nubq73/v79+9WuXTtt3749YH9p8I/b7ZbL5dKOHTuUkJBg9zhGY1+EFvZH6GBfhI7Kykq1bdtWSUlJAdtmSEdP69at1axZM1VUVPgsr6ioUEpKSr3PcTqdcjqdxyxPTEzkGzhEJCQksC9CBPsitLA/Qgf7InRERQXu9OOQPpE5NjZWffr0UVFRkXdZXV2dioqKlJGRYeNkAAAg3IT0kR5Jys3NVXZ2tvr27avzzjtPM2fOVFVVlW699Va7RwMAAGEk5KNn1KhR+v777/XII4+ovLxc55xzjpYuXXrMyc3H43Q6lZ+fX+9bXmha7IvQwb4ILeyP0MG+CB3B2BcOK5C/CwYAABCiQvqcHgAAgEAhegAAgBGIHgAAYASiBwAAGCEioqewsFDt27dXXFyc+vfvr08++eSE6y9cuFBdunRRXFycevbsqSVLljTRpJGvMfvi+eef18UXX6zTTjtNp512mjIzM39136HhGvvv4qj58+fL4XBo+PDhwR3QII3dF/v371dOTo5SU1PldDp19tln8/9UADV2f8ycOVOdO3dW8+bN5XK5NGHCBB0+fLiJpo1MH3zwga666iqlpaXJ4XBo8eLFv/qclStX6txzz5XT6VSnTp00d+7cxn9hK8zNnz/fio2NtV588UXryy+/tO644w6rZcuWVkVFRb3rf/zxx1azZs2sadOmWV999ZU1adIkKyYmxvriiy+aePLI09h9cdNNN1mFhYXWunXrrA0bNlhjxoyxEhMTre+++66JJ488jd0XR23bts0644wzrIsvvti65pprmmbYCNfYfeHxeKy+fftaw4YNsz766CNr27Zt1sqVK63S0tImnjwyNXZ/vPLKK5bT6bReeeUVa9u2bdayZcus1NRUa8KECU08eWRZsmSJ9dBDD1mvv/66JclatGjRCdffunWr1aJFCys3N9f66quvrGeeecZq1qyZtXTp0kZ93bCPnvPOO8/Kycnx3q+trbXS0tKsgoKCete//vrrrSuuuMJnWf/+/a0777wzqHOaoLH74peOHDlixcfHWy+99FKwRjSGP/viyJEj1gUXXGD94x//sLKzs4meAGnsvpg1a5bVoUMHq7q6uqlGNEpj90dOTo41aNAgn2W5ubnWhRdeGNQ5TdKQ6Jk4caLVvXt3n2WjRo2ysrKyGvW1wvrtrerqapWUlCgzM9O7LCoqSpmZmSouLq73OcXFxT7rS1JWVtZx10fD+LMvfunQoUOqqakJ6IfLmcjfffHoo4+qTZs2uu2225piTCP4sy/efPNNZWRkKCcnR8nJyerRo4eeeOIJ1dbWNtXYEcuf/XHBBReopKTE+xbY1q1btWTJEg0bNqxJZsZPAvWzO+SvyHwie/fuVW1t7TFXZ05OTtbGjRvrfU55eXm965eXlwdtThP4sy9+6f7771daWtox39hoHH/2xUcffaQXXnhBpaWlTTChOfzZF1u3btV7772nm2++WUuWLNHmzZt1zz33qKamRvn5+U0xdsTyZ3/cdNNN2rt3ry666CJZlqUjR47orrvu0oMPPtgUI+P/O97PbrfbrR9//FHNmzdv0HbC+kgPIsfUqVM1f/58LVq0SHFxcXaPY5QDBw5o9OjRev7559W6dWu7xzFeXV2d2rRpo+eee059+vTRqFGj9NBDD2n27Nl2j2aklStX6oknntCzzz6rzz//XK+//rreeecdPfbYY3aPBj+E9ZGe1q1bq1mzZqqoqPBZXlFRoZSUlHqfk5KS0qj10TD+7IujnnrqKU2dOlUrVqxQr169gjmmERq7L7Zs2aJvv/1WV111lXdZXV2dJCk6OlqbNm1Sx44dgzt0hPLn30VqaqpiYmLUrFkz77KuXbuqvLxc1dXVio2NDerMkcyf/fHwww9r9OjRuv322yVJPXv2VFVVlcaOHauHHnpIUVEcO2gKx/vZnZCQ0OCjPFKYH+mJjY1Vnz59VFRU5F1WV1enoqIiZWRk1PucjIwMn/Ulafny5cddHw3jz76QpGnTpumxxx7T0qVL1bdv36YYNeI1dl906dJFX3zxhUpLS723q6++WgMHDlRpaalcLldTjh9R/Pl3ceGFF2rz5s3e8JSkr7/+WqmpqQTPSfJnfxw6dOiYsDkapBYfXdlkAvazu3HnWIee+fPnW06n05o7d6711VdfWWPHjrVatmxplZeXW5ZlWaNHj7YeeOAB7/off/yxFR0dbT311FPWhg0brPz8fH5lPUAauy+mTp1qxcbGWq+99pq1e/du7+3AgQN2vYSI0dh98Uv89lbgNHZfbN++3YqPj7fGjRtnbdq0yXr77betNm3aWI8//rhdLyGiNHZ/5OfnW/Hx8darr75qbd261frPf/5jdezY0br++uvtegkR4cCBA9a6deusdevWWZKs6dOnW+vWrbP+97//WZZlWQ888IA1evRo7/pHf2X9T3/6k7VhwwarsLDQzF9ZtyzLeuaZZ6y2bdtasbGx1nnnnWetWbPG+9ill15qZWdn+6y/YMEC6+yzz7ZiY2Ot7t27W++8804TTxy5GrMv2rVrZ0k65pafn9/0g0egxv67+DmiJ7Aauy9Wr15t9e/f33I6nVaHDh2sP//5z9aRI0eaeOrI1Zj9UVNTY02ePNnq2LGjFRcXZ7lcLuuee+6x/u///q/pB48g77//fr3//x/9u8/OzrYuvfTSY55zzjnnWLGxsVaHDh2sOXPmNPrrOiyL43MAACDyhfU5PQAAAA1F9AAAACMQPQAAwAhEDwAAMALRAwAAjED0AAAAIxA9AADACEQPAAAwAtEDwAgOh0OLFy+2ewwANiJ6ANhqzJgxcjgccjgciomJUXp6uiZOnKjDhw/bPRqACBNt9wAAMGTIEM2ZM0c1NTUqKSlRdna2HA6H/vKXv9g9GoAIwpEeALZzOp1KSUmRy+XS8OHDlZmZqeXLl0uS9u3bpxtvvFFnnHGGWrRooZ49e+rVV1/1ef6AAQN07733auLEiUpKSlJKSoomT558wq+Zn5+v1NRUrV+/PlgvC0CIIXoAhJSysjKtXr1asbGxkqTDhw+rT58+euedd1RWVqaxY8dq9OjR+uSTT3ye99JLL+mUU07R2rVrNW3aND366KPecPo5y7I0fvx4vfzyy/rwww/Vq1evJnldAOzHp6wDsNWYMWP0r3/9S3FxcTpy5Ig8Ho+ioqK0YMECjRw5st7nXHnllerSpYueeuopST8d6amtrdWHH37oXee8887ToEGDNHXqVEk/nci8cOFCLVq0SOvWrdPy5ct1xhlnBP8FAggZnNMDwHYDBw7UrFmzVFVVpRkzZig6OtobPLW1tXriiSe0YMEC7dy5U9XV1fJ4PGrRooXPNn55xCY1NVV79uzxWTZhwgQ5nU6tWbNGrVu3Du6LAhByeHsLgO1OOeUUderUSb1799aLL76otWvX6oUXXpAkPfnkk3r66ad1//336/3331dpaamysrJUXV3ts42YmBif+w6HQ3V1dT7LLr/8cu3cuVPLli0L7gsCEJKIHgAhJSoqSg8++KAmTZqkH3/8UR9//LGuueYa3XLLLerdu7c6dOigr7/+2q9tX3311Zo3b55uv/12zZ8/P8CTAwh1RA+AkPPb3/5WzZo1U2Fhoc466ywtX75cq1ev1oYNG3TnnXeqoqLC721fe+21+uc//6lbb71Vr732WgCnBhDqOKcHQMiJjo7WuHHjNG3aNK1bt05bt25VVlaWWrRoobFjx2r48OGqrKz0e/vXXXed6urqNHr0aEVFRWnEiBEBnB5AqOK3twAAgBF4ewsAABiB6AEAAEYgegAAgBGIHgAAYASiBwAAGIHoAQAARiB6AACAEYgeAABgBKIHAAAYgegBAABGIHoAAIAR/h/7g3ZWnI2QrwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "decoy_ranks_dist = uniform(0, 1)\n",
    "plt.hist(decoy_ranks_dist.rvs(1000),bins=np.linspace(0, 1));\n",
    "plt.xlim(0, 1)\n",
    "plt.xlabel('Rank')\n",
    "plt.ylabel('Count')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28906ad9",
   "metadata": {},
   "source": [
    "now let's move to a benchmark that one might perform to evaluate scoring functions. in this case, we'll assess performance using the AUROC of 100 actives among 10,000 decoys:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f79eaaaf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AUROC is: 0.9819749999999999\n"
     ]
    }
   ],
   "source": [
    "n_pos = 100\n",
    "n_neg = 10_000\n",
    "pos_ranks = active_ranks_dist.rvs(n_pos)\n",
    "neg_ranks = decoy_ranks_dist.rvs(n_neg)\n",
    "\n",
    "preds = np.concatenate([pos_ranks, neg_ranks])\n",
    "labels = np.concatenate([np.ones(n_pos), np.zeros(n_neg)])\n",
    "\n",
    "print('AUROC is:', roc_auc_score(labels, preds))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34c05fe2",
   "metadata": {},
   "source": [
    "phenomenal AUROC score! now what happens when we brute-force dock and score a library of 40 billion compounds using this scoring function? \n",
    "\n",
    "for the simulation purposes, we'll assume there's a 'background' hit rate of 0.1%, which is a reasonable hit rate from HTS. in practice, HTS libraries are curated for diversity and bio-like compounds, so an ultra-large library would have a lower hit rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e024bb4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "library_size = 40e9\n",
    "pos_rate = 0.001\n",
    "n_pos = int(library_size*pos_rate) #number of actives in the whole 40 billion cmpd library.\n",
    "\n",
    "pos_ranks = active_ranks_dist.rvs(n_pos) #these are the ranks the actives achieve by our scoring function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a43a18a",
   "metadata": {},
   "source": [
    "now let's ask what the _best_ rank of any of the true positives is, by figuring out where the highest scoring positive sits within the 40 billion compounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "915bda8b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best rank is: 39999758564\n"
     ]
    }
   ],
   "source": [
    "best_rank = pos_ranks.max()\n",
    "print('Best rank is:', int(best_rank * library_size))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "273843b7",
   "metadata": {},
   "source": [
    "seems pretty good - it's very close to 40 billion, so it's right near the top. but consider that we'll only test on the order of 1e3 molecules, selected from the very, very, _very_ top of the list. unfortunately, even though the best-scored active is very close from the top, it's not _that_ close - it still woulnd't be found unless we tested on the order of 100k molecules. this is how many inactives are ranked better than the best active:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "25c46131",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of negatives before first positive is found: 241435\n"
     ]
    }
   ],
   "source": [
    "print('Number of negatives before first positive is found:', int(library_size - (best_rank*library_size)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7843bcae",
   "metadata": {},
   "source": [
    "# oh,\n",
    "So AUROC on the scale of 10s or 100s of thousands of decoys is a relatively weak benchmark. a severe benchmark would use millions or billions of decoys - but of course, this requires random decoy selection, so now you can't trust that every high-scoring decoy is actually a true negative!\n",
    "\n",
    "the problem comes from changing the size of the slice that we can synthesise and test. taking the top-1000 compounds from a library of a million in-stock compounds would likely find a hit. but the same slice out of a billion compounds is, proportionally, 1000X more extreme, so you'd want an even more stringent model to ensure there are some hits in the slice. \n",
    "\n",
    "one way to address this might be to estimate the change in AUROC with library size, showing that performance increases as the number of decoys increases, all the way up to 40 billion to show that some actives are still in the top 100/1000/etc. \n",
    "\n",
    "the gold standard, of course, is actual synthesis and testing of high-scoring compounds - each prospective experiment is often, unfortunately, a rich source of high-scoring decoy data for benchmarking. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d42fff",
   "metadata": {},
   "source": [
    "### appendix: a million-scale library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8cab509e",
   "metadata": {},
   "outputs": [],
   "source": [
    "library_size = 1e6\n",
    "pos_rate = 0.001\n",
    "n_pos = int(library_size*pos_rate) #number of actives in the whole 40 billion cmpd library.\n",
    "\n",
    "pos_ranks = active_ranks_dist.rvs(n_pos) #these are the ranks the actives achieve by our scoring function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e53ed10b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best rank is: 462\n"
     ]
    }
   ],
   "source": [
    "best_rank = pos_ranks.max()\n",
    "print('Best rank is:', int(library_size - (best_rank * library_size)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7098b795",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = (pos_ranks * library_size).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "49989aba",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of hits in top 1000: 6\n"
     ]
    }
   ],
   "source": [
    "n_in_top_1k = (r[(-r).argsort()]>999000).sum()\n",
    "print('Number of hits in top 1000:', n_in_top_1k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26e94e83",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}