part_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "80e6d90a-38d8-4a75-9fb9-4637e47ec536",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Statistics for data science and measurement\n",
    "## Babak Moghadas 1 and Brian Caffo 1, 2\n",
    "### 1 Department of Biostatistics \n",
    "### 2 Department of Biomedical Engineering\n",
    "### Bloomberg School of Public Health\n",
    "### Johns Hopkins University"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae7dbfc9-a6b9-485b-9396-844dd473d820",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## About these slides\n",
    "\n",
    "+ Goal is the least that you need to know\n",
    "+ All slides are created in Jupyter notebooks - free to use and open source\n",
    "+ [Github repo](https://github.com/smart-stats/odsc_2023)\n",
    "+ [pyglide](https://pypi.org/project/pyglide/) provides the interactivity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be399fab-4e3d-46dd-b2ab-5329ccad4fa9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Part 3 validation of ML algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d839961-ae68-4030-a72f-20dda2a10e16",
   "metadata": {},
   "source": [
    "+ In this component, we'll discuss validating the output of ML algorithms.\n",
    "+ Typically, we split the data into:\n",
    "   + **Training data**: data used to train the model.\n",
    "   + **validation data**: data used to choose hyperparameters, such as layers in a neural network.\n",
    "   + **testing data**: data used only for final evaluation of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd3e8e58-b026-4921-b4e0-96c5f96f17ba",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "### Binary outcomes (two-class classification)\n",
    "\n",
    "Consider the following definitions for the results of a diagnostic test $\\hat Y\\in \\{0,1\\}$, where the\n",
    "actual disease state is $Y\\in \\{0,1\\}$. Assume 1 represents having/testing for the disease.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e1aab6d-ec2d-4aab-bd2e-0faec20b06d3",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### Measures conditional on disease status:\n",
    "1. **Sensitivity** $P(\\hat Y =1 ~|~ Y=1)$, probability that the prediction is positive given the disease is present.\n",
    "   + Also called the true positive rate.  \n",
    "   + $P(\\hat Y =0 ~|~ Y=1)$, one minus the sensitivity, is the false positive or miss rate.\n",
    "   + Also called recall and hit rate\n",
    "2. **Specificity** $P(\\hat Y=0 ~|~ Y=0)$, probability that the prediction is negative given the disease is absent.\n",
    "    + Also called the true negative rate. \n",
    "    + $P(\\hat Y=1~|~ Y=0)$, one minus the specificity is the false negative rate or fall out.\n",
    "    + Also called selectivity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e2233e0-be31-4d47-9a53-9ae735f202e0",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### Measures conditional on the data\n",
    "\n",
    "3. **PPV**, positive predictive value $P(Y=1 ~|~ \\hat Y=1)$.\n",
    "    + Also called the precision\n",
    "4. **NPV**, negative predictive value $P(Y=0 ~|~ \\hat Y =0)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47d93968-7ed1-4473-a1d3-a96128b06d2b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### Other measures\n",
    "\n",
    "5. **DLR+**, diagnostic likelihood ratio of a positive prediction $P(\\hat Y=1 ~|~ Y=1) / P(\\hat Y=1 ~|~ Y= 0)$ is also the sensitivity over one minus the specificity, i.e. the true positive rate divided by the false positive rate.\n",
    "6. **DLR-**, diagnostic likelihood ratio of a negative prediction $P(\\hat Y=0 ~|~ Y=1) / P(\\hat Y=0 ~|~ Y=0)$ is one minus the sensitivity divided by the specificity or the false negative rate divided by the true negative rate. \n",
    "7. The **disease prevalence** is $P(Y=1)$ (and, generally less discussed, the prediction disease prevalance $P(\\hat Y =1)$.)\n",
    "8. The accuracy is $P(\\hat Y = y) = P(\\hat Y = 1 ~|~ Y = 1) P(Y = 1) + P(\\hat Y = 0 ~|~ Y = 0) P(Y = 0)$, which is the sensitivity times the prevalence plus the specificity times one minus the prevalance.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ab3e7b6-93cd-4a83-8a66-2559457558c1",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "##### Role of sampling in estimation\n",
    "\n",
    "+ If you have a cross-sectional sample, then all of these quatities are directly estimable. \n",
    "+ If the data were sampled by disease status, or data was augmented to force outcome balance, then $Y$ should be conditioned and the sensitivity, specificity and DLR+/- are directly estimable. \n",
    "    + You can obtain the PPV and NPV using Bayes' rule given a disease prevalance.  \n",
    "+ In a rare setting where the design fixes the prediction result, the NPV and PPV would be directly estimable and one would have to use the prevalance of a positive prediction and Bayes' rule to obtain the sensitivity and specificity. \n",
    "+ Accuracy depends on the prevalence, so imbalance in the outcome categories bounds a reasonable accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdcaadc4-5f0d-4786-9997-17b91fd51654",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Basic example\n",
    "\n",
    "A study comparing the efficacy of HIV tests, reports on an experiment which concluded that HIV antibody tests have a sensitivity of 99.7% and a specificity of 98.5%. \n",
    "\n",
    "Suppose that a subject, from a population with a .1% prevalence of HIV, receives a positive test result. What is the positive predictive value?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dd82c05-4920-419e-a6e3-429dcf6d5154",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "+ Mathematically, we want $P(Y=1 | \\hat Y=1)$ \n",
    "+ We have:\n",
    "  + the sensitivity, \n",
    "  $$P(\\hat Y=1 | Y=1) = .997$$\n",
    "  + the specificity, \n",
    "  $$P(\\hat Y=0 | Y=0) =0.985$$\n",
    "  + the prevalence \n",
    "  $$P(Y=1) =0.001$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f3b6efa-068e-471f-99fd-a2401c2770a7",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### Bayes's rule\n",
    "$$\n",
    "\\begin{align*} \n",
    "\\small P(Y=1 ~|~ \\hat Y =1) & = \\frac{\\small P(\\hat Y =1 ~|~ Y=1)\\small P(Y=1)}{\\small P(\\hat Y=1~|~Y=1)\\small P(Y=1) + \\small P(\\hat Y=1 ~|~ Y=0)\\small P(Y=0)}\\\\\n",
    "               & = \\frac{\\small P(\\hat Y=1|Y=1)\\small P(Y=1)}{\\small P(\\hat Y=1|Y=1)\\small P(Y=1) + (1 - \\small P(\\hat Y=0 ~|~ Y = 0))(1 - \\small P(Y=1))} \\\\\n",
    "               & = \\frac{\\small .997\\times .001}{\\small .997 \\times .001 + .015 \\times .999}\\ = \\small .062 \n",
    "\\end{align*}\n",
    " \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45c483b3-4026-4409-946a-958dffc537b1",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### DLR\n",
    "\n",
    "$$\n",
    "P(Y = 1 ~|~ \\hat Y = 1) = \\frac{P(\\hat Y=1~|~ Y = 1)P(Y=1)}{P(\\hat Y=1)} \n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$ P(Y=0 ~|~ \\hat Y=1) = \\frac{P(\\hat Y=1 ~|~ Y=0)P(Y=0)}{P(\\hat Y=1)}. $$\n",
    "\n",
    "\n",
    "$$ \n",
    "\\frac{P(Y = 1 ~|~ \\hat Y=1)}{P(Y=0 ~|~ \\hat Y=1)} = \\frac{P(\\hat Y = 1 ~|~ Y=1)}{P(\\hat Y = 1 ~|~ Y=0)}\\times \\frac{P(Y=1)}{P(Y=0)} \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "620b9329-fdf8-4f03-bf61-0d470dee6ad8",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### DLR continued\n",
    "\n",
    "+ In other words, the post test odds of disease is the pretest odds of disease times the $DLR_+$. \n",
    "+ Similarly, $DLR_-$ relates the decrease in the odds of the disease after a negative test result to the odds of disease prior to the test. \n",
    "+ DLRs are the factors by which you multiply your pretest odds to get your post test odds. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d25c600-a7b8-41ab-9768-795852e6d74c",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### DLR+ example\n",
    "\n",
    "+ Suppose a subject has a positive HIV test, $DLR_+ = .997 / (1 - .985) = 66$. \n",
    "+ The result of the positive test is that the odds of disease is now 66 times the pretest odds. \n",
    "+ Or, equivalently, the hypothesis of disease is 66 times more supported by the data than the hypothesis of no disease\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bee6b665-f4f7-4d3f-a450-98614b42efe2",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### DRL- example\n",
    "\n",
    "+ Suppose instead that a subject has a negative test result. \n",
    "+ Then $DLR_- = (1 - .997) / .985 =.003$ \n",
    "+ Therefore, the post-test odds of disease is now 0.3% of the pretest odds given the negative test. \n",
    "+ Or, the hypothesis of disease is supported $.003$ times that of the hypothesis of absence of disease given the negative test result\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e9c4b06-7d26-48c1-a47d-9aa60880b2d5",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "#### ROC curves\n",
    "\n",
    "+ The measures covered above require a binary diagnostic value. Typically our ML algorithms give a range.\n",
    "+ Often 50% (estimated more likely than not) is not a reasonable cutoff.\n",
    "+ ROC curves consider the sensitivity and 1-specificity; i.e. the true positive and false positive rate, for all possible cutoffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8c42740b-8d03-4bc4-a652-1a55e248497b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>FLAIR</th>\n",
       "      <th>GOLD_Lesions</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.143692</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.652552</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.036099</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.037692</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.580589</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      FLAIR  GOLD_Lesions\n",
       "0  1.143692             0\n",
       "1  1.652552             0\n",
       "2  1.036099             0\n",
       "3  1.037692             0\n",
       "4  1.580589             0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "dat = pd.read_csv(\"https://raw.githubusercontent.com/bcaffo\"\\\n",
    "    \"/ds4bme_intro/master/data/oasis.csv\")\n",
    "dat[ ['FLAIR', 'GOLD_Lesions'] ].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d2f3e94b-641e-4897-8500-ce454d5e03fb",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4lElEQVR4nO3deXRUdZ4//Pe9t5bUXqlKVVJJgECAGECUBhQVRQEnIAi0G92oz/xsG+ccHe3jPDOtTTtsOvZgz3lsRdTpttW28ac0LiwRERUQUAFZZAskEEIWklSSqlRq3+/zRyQSyVJJqu6t5fM6x9NQuXXv53ZC3nW/K8PzPA9CCCGkB6zYBRBCCEleFBKEEEJ6RSFBCCGkVxQShBBCekUhQQghpFcUEoQQQnpFIUEIIaRXErELiLf2dg+i0fhO/TAa1bDZ3HE9ZzLKhPuke0wfmXCfQtwjyzLIzlb1+vW0C4lolI97SFw6bybIhPuke0wfmXCfYt8jNTcRQgjpFYUEIYSQXlFIEEII6RWFBCGEkF5RSBBCCOkVhQQhhJBeUUgQQgjpVdrNkyDCiEZ57DnWiIOnrWBZBjeMz8MNE/LAMozYpRFC4ohCggxYJBrFax+fRIvDh6lX5YLneXx6oA6HKlvw6KIJkEo4sUskhMQJNTeRAfvwq2q0uwNYfNtojCnUYewwPX45czT8wQj++slp0I64hKQPCgkyIHVWF/Ydb8a8aSPAcT/++HAcizuuH4ELzS7sO9EkYoWEkHiikCAD8sHuakwblwtllvSKr0klLOZePxwbd1XD7QuJUB0hJN4oJEjMGlrcqLW6MLHY2OsxudlKjCnUofybC8IVRghJGAoJErMvDzfgmmIjJFzfPzbTxuVh7/FGOD1BgSojhCQKhQSJSSgcxcEzLbh6VO9PEZdolFKUDNPjy8MNAlRGCEkkCgkSkxPnbcjNVkCjlMV0/OSxZuz+/iLCkWiCKyOEJBKFBInJoTMtGF2gi/l4oy4LBm0Wjp5tS2BVhJBEo5Ag/YpGeZw4bxtQSADAhCID9h5rTFBVhBAhCBISa9aswcyZM1FSUoKqqqoej4lEIli1ahVmz56N22+/HRs3bhSiNBKD801OqBVSaFWxNTVdMmaYDucudqCDOrAJSVmChMSsWbPw7rvvoqCgoNdjtm7dirq6OuzYsQMbNmzA2rVr0dBAHZ/JoOKCHcNzNQN+n0zCoThfh8OVLQmoihAiBEFCYsqUKbBYLH0es23bNtx7771gWRYGgwGzZ8/G9u3bhSiP9OPkeTtG5KoH9d6xw3Q4UGGNc0WEEKEkTZ9EU1MT8vPzu/5usVjQ3NwsYkUE6Bz6Wmd1ocA0uJAYadGirsUNl5eanAhJRWm3CqzROLhfZv0xmQbe3JKKfnqfp2vsMBuUyB3C/Y8dpsd5qwezr+t/joUQMuF7mQn3CGTGfYp9j0kTEhaLBY2NjZg4cSKAK58sYmWzuRGNxncVUpNJg9ZWV1zPmYx6us/vTjUiN1sBh8M76PMOM6nw1ZF6XDMye6glDlkmfC8z4R6BzLhPIe6RZZk+P1wnTXPTnDlzsHHjRkSjUdjtdnzxxRcoKysTu6yMd7ahAxaDckjnGGnR4kxtO02sIyQFCRISzz33HG655RY0NzfjoYcewrx58wAAS5cuxYkTJwAACxcuRGFhIf7pn/4J9913Hx577DEMGzZMiPJIHy40OWExqoZ0DlWWFAZtFs42dMSpKkKIUBg+zXaIoeamwfvpfXZ4gvjd/36Lx++6GswQtyXde7wReo0c9946eqhlDkkmfC8z4R6BzLhPam4iSa22ufMpYqgBAQBFeVqcOG+PQ1WEECFRSJBe1VldMOuz4nKu/BwVWh0+GgpLSIqhkCC9utDsRo5eEZdzcSyDYSY1KusccTkfIUQYFBKkV/UtbpjjFBIAUGhS4dQFanIiJJVQSJAeBUIRONwBGLTxaW4CgOG5GpypbY/b+QghiUchQXrUZPPAoJWDY4feaX2JWa9AhydI25oSkkIoJEiPLrZ6kBPHpwigc6jdMLMaVfWOuJ6XEJI4FBKkRw2tbhh18Q0JACjIUeFMHTU5EZIqKCRIjxpaPTDG+UkCAApMalTSkwQhKYNCgvSo2eZNSEjkZSvQ2u6D1x+O+7kJIfFHIUGuEApH0OEJQKeWx/3cHMfCYlShupHWcSIkFVBIkCtY7T5ka+I7suly+TlKnKUmJ0JSAoUEuUKT3QuDJv5NTZfk56hohBMhKYJCglyh2e6FXiNL2PkLclSotboQidL+EoQkOwoJcoVmmwfZ6sQ9SWTJJNCq5Gho8STsGoSQ+KCQIFdotvsS+iQBdD5NnLtIndeEJDsKCXKFVocvoU8SAGAxKqlfgpAUQCFBuvH6wwiEIlArJAm9Tn4ODYMlJBVQSJBuWh0+GDTyuOxG1xeDRg5fIIwOdyCh1yGEDA2FBOmm1eGDPgGT6H6KYRgU5KhQ3ehM+LUIIYNHIUG6aXH4oFMnttP6kjyDCucaqMmJkGRGIUG6sdq90KmECYl8oxJnLzoEuRYhZHAoJEg3LQI1NwGAxahCvdWNcIQm1RGSrCgkSDetDp9gTxJyGQe9Wo6GVrcg1yOEDByFBOkSiUTR4Q5CK1BIAD8Mhb1IndeEJCsKCdLF1uGHSiGBhBPuxyLPoMTZBodg1yOEDAyFBOlibfdCpxKmP+KS/BwlPUkQksQoJEiXFgFHNl1i1GbB4w/B6QkKel1CSGwoJEiXZpsHGqVU0GtemlRHi/0RkpwoJEiXZpvwzU1A52J/NKmOkOQkWEjU1NRg8eLFKCsrw+LFi3HhwoUrjrHZbHjkkUdw5513Ys6cOVi5ciXC4bBQJWa81nYftCphnyQAIN+ookl1hCQpwUJixYoVWLJkCT777DMsWbIEy5cvv+KY119/HcXFxdi6dSu2bt2KU6dOYceOHUKVmPFaHV5olcL2SQCAJYcm1RGSrAQJCZvNhoqKCsyfPx8AMH/+fFRUVMBut3c7jmEYeDweRKNRBINBhEIh5ObmClFixovyPOxOv6BzJC6RSzlka+SotboEvzYhpG+J3TTgB01NTcjNzQXHcQAAjuNgNpvR1NQEg8HQddyjjz6Kxx9/HNOnT4fP58P999+PyZMnD+haRqM6rrVfYjJpEnLeZGHr8EGZJUVOgv7/68/IAh2aHQFMuybx/z+n+/cSyIx7BDLjPsW+R0FCIlbbt29HSUkJ/va3v8Hj8WDp0qXYvn075syZE/M5bDY3olE+rnWZTBq0tqb3p9xzDR3Qa+RwOLyiXN+oluPomWbcNM6c0OtkwvcyE+4RyIz7FOIeWZbp88O1IM1NFosFVqsVkUgEABCJRNDS0gKLxdLtuPXr12PBggVgWRYajQYzZ87EgQMHhCgx49mcfsEW9utJgUmFsw0d4Pn4BjwhZGgECQmj0YjS0lKUl5cDAMrLy1FaWtqtqQkACgsLsWfPHgBAMBjEt99+izFjxghRYsazu8Tpj7jk0iS+1g6/aDUQQq4k2OimlStXYv369SgrK8P69euxatUqAMDSpUtx4sQJAMCyZctw+PBh3HnnnVi0aBGKiopw3333CVViRmvr8IsyR+IShmFQaFLjbL1DtBoIIVdi+DR7vqc+icH508ZjmFCcg0KjUrQajlS1wh8M41fzxiXsGpnwvcyEewQy4z4zpk+CJD+70w+9RrwnCQAoNKlQWU8zrwlJJhQSBADQ7gpAJ2LHNQDk6BRweYPocAdErYMQ8iMKCYJgKAJ/MAKVQvglOS7HsgyGmdWopH4JQpIGhQRBuysArUoGlmHELgUFOSqcqW0XuwxCyA8oJAjsroAoazb1ZJhZgzN1DrHLIIT8gEKCwO70Qy3wPhK9yc1WwOEOoIM2ISIkKVBIELS7AlArkmOFlq5+iTpqciIkGVBIENg6/NAokqO5CQCGmdQ4VWPv/0BCSMJRSBDYXH7Bty3ty4g8DSou0JMEIcmAQoKg3RWAJkk6rgEgR5eFQCiCFodP7FIIyXgUEgQOVwBqkedIXI5hGBTlaVBxgZqcCBEbhUSGC4V/mEiXlRwd15cMz1XjeLVN7DIIyXgUEhmu3R2ERikFkwQT6S5XlKdFZV07IlHa95oQMVFIZDhHkvVHXKJWSKFTyVF90Sl2KYRkNAqJDGd3+aFJov6IyxXlaXD8PDU5ESImCokM53AFRV/YrzcjLVocO9cmdhmEZDQKiQxnc/qSamTT5QpyVLA7/Wh30dLhhIiFQiLDtTuTa/jr5ViWwSiLFser6WmCELFQSGS4dncgqWZb/9RIixZHqlrFLoOQjEUhkeEc7mDSdlwDwKh8LarqOxAIRsQuhZCMFHNIfPnllwiHw4mshQgsyvNweoJJ29wEAFkyCfJzlDhJC/4RIoqYQ+Kll17C9OnTsXr1ahw7diyRNRGBuDxBZMk4cFxyP1COytfhSFWL2GUQkpFi/u2wZcsWvP3225DL5Xj88cdRVlaGV199FQ0NDYmsjyRQZ39E8k2k+6kxBTocq7YhHKHZ14QIbUAfIa+66io89dRT+Oqrr7BixQps374dt99+O+6//35s2bIFUVpCIaV0rv6avE1Nl2hVMujVclTVO8QuhZCMM+BV3erq6rBlyxZs2bIFDMPgiSeegMViwbvvvosdO3bglVdeSUSdJAEcrgBUWckfEgAwpkCL7860YFyRQexSCMkoMYfEu+++i82bN6O2thZz587FCy+8gGuvvbbr62VlZbjxxhsTUSNJEHsSbVvanzGFemzYdQ4PlpWATbLFCAlJZzH/htizZw8eeughzJo1CzLZle3YCoUCa9eujWtxJLHszgB0quTvkwAAgzYLCrkE5xo6MHaYXuxyCMkYMfdJXHfddZg7d+4VAfHWW291/Xn69Onxq4wknMOdvLOtezKmUIfvzljFLoOQjBJzSKxbt67H11977bW4FUOElSod15eUDNPjcGUrojwvdimEZIx+m5u+/fZbAEAkEsH+/fvBX/YPtKGhASqVKqYL1dTU4Omnn4bD4YBer8eaNWtQVFR0xXHbtm3Da6+9Bp7nwTAM3nrrLeTk5MR4O2QgOjyp9SRh1GZBwrGoaXSiuEAndjmEZIR+Q+L3v/89ACAYDGLZsmVdrzMMA5PJhGeeeSamC61YsQJLlizBwoULsXnzZixfvhzvvPNOt2NOnDiBV155BX/7299gMpngcrl67P8gQxcIRRAKR5El48QuJWYMw2DsMD0OnWmhkCBEIP2GxM6dOwEAv/3tb/HCCy8M6iI2mw0VFRVd/Rfz58/Hs88+C7vdDoPhxyGNb7/9Nn71q1/BZDIBADQazaCuR/p3aUe6ZNu2tD9jC/XY8k0N7ps5OuVqJyQVxdwnMdiAAICmpibk5uaC4zo/tXIcB7PZjKampm7HVVdXo76+Hvfffz9+/vOf49VXX+3WvEXiJ9X6Iy4x6bMAAHVWt8iVEJIZ+nySmDt3Lj799FMAwIwZM3r95LZ79+64FBOJRFBZWYm33noLwWAQv/71r5Gfn49FixbFfA6jUR2XWn7KZEqvp5pT9R0w6BTQ65XdXv/p35PRxNEmnKxzYMrV+YN6f7p9L3uSCfcIZMZ9in2PfYbEs88+2/XnP/7xj4O+iMVigdVqRSQSAcdxiEQiaGlpgcVi6XZcfn4+5syZA5lMBplMhlmzZuH48eMDCgmbzY1oNL5PHyaTBq2trrieU2x1jQ5IWQYOh7frNb1e2e3vyWpYjgrbD9bhjqmFA25ySsfv5U9lwj0CmXGfQtwjyzJ9frjuMySmTJnS9efrrrtu0EUYjUaUlpaivLwcCxcuRHl5OUpLS7v1RwCdfRVfffUVFi5ciHA4jP3796OsrGzQ1yW9szv9UGWJP9ua87VD3nYaEm8bmEgIkSwtgvqRCOqLALbn+vIMCgRDEVxs86DQlJgnR0JIp5j7JN566y2cPn0aAPD999/j1ltvxaxZs3D06NGY3r9y5UqsX78eZWVlWL9+PVatWgUAWLp0KU6cOAEAmDdvHoxGI+644w4sWrQIo0ePxj333DPQeyIxaHeKuwKs1FkP46HXYP7mj1BYT4AJ+TpfdzVBX/Eh8navhLp6B5jwlftbMwyDMYU6HKmkHesISTSGj7FneMaMGSgvL4dGo8GDDz6IWbNmQaVS4R//+Ac2btyY6DpjRs1NsXn2b99h2rg8DDP/+ElckOamaASa6u1Q1X8Lb/5k+HPGAeyVw3A5nw3KpiOQeFrgGP8LBHJKun29rsWFfcebsPrh6wd0+XT8Xv5UJtwjkBn3mfTNTZdzuVzQaDRwu92orKzE22+/DY7jsGbNmrgUSoTlcAcFH93EhAMwfP8m2JAf7ePvBS/tfSJmRGGEa9TtkHbUIfvEu/AMuwmu4tsBpvPhtzBHDbsrgDaHDzl6hVC3QEjGibm5yWKx4MiRI9i2bRumTJkCjuPgdru7hrWS1CHGtqVMOADjodcAhkPH2Hl9BsTlQrrhcJTeDYX1GPQn/i8Q7dzrmmUZjC7Q4ei5tkSWTUjGizkkfvvb3+KJJ57A66+/jkcffRQAsGvXLlx99dUJK44khtsbglzKQSLUtqXRCAzH3gYvVcJVdFvX00DMb5ep4Bh7JyTeVhiOvdMVFMX5WhyupG1NCUmkmJubZsyYgX379nV7bc6cOZgzZ07ciyKJ1e4KQKMS7ilCV7kZTMgL5+g7gMHOkuakcI6+A9pz26A/tQGOCb9EUZ4W2/bXweMPpczmSYSkmgGNgXS5XKipqYHH4+n2+g033BDXokhitQu4RHhW8/fIsh6HY9y9PXZQDwjLwVk8B/qqrdCc3QbX2HkYnqvGiWobpo3Pi0/BhJBuYg6Jjz76CKtXr4ZSqURWVlbX6wzD4Msvv0xIcSQxHK4ANIrED3/lfO3QV2yEc8w88BJ5nE4qRcfoudCf/hBhlRkjLSPx/bk2CglCEiTmkHjxxRfx0ksvYcaMGYmshwig3SXARDqeh/7ke/CZJyKsMsf31FIFnGPmQn9mE8Zd/Wv8+bgbkWgUHCtQHwshGSTmf1WRSIR2nksTNmfim5uUFw+ACzjhs0xKyPkjCiNcI27GsNN/h1kVRfVFZ0KuQ0imizkkli5ditdeew3RaDSR9RABtLsCUCdwjgQb9EBbVQ5X0YwBj2QaiKBhNIK6EbhPvhfHq2koLCGJEHObw9tvv422tja88cYb0Ov13b4Wr1VgiTDaXQFoEvgkoTn7CYLZoxBRmhJ2jUs8hdNgcHwApmo3cOvohF+PkEwTc0gMZRVYklw6PAGoE7Ruk8TVCIX1GNon/DIh578Cy8E7+p8w/eQHcFycBX3BSGGuS0iGiDkkhrIKLEkegWDntqWKBG1bqqvcCl/ez8BLsvo/OE6iymxUySagZNf/gv/lc2A48Ve3JSRdxNxgHAwG8eKLL2LWrFmYPHkyAGDfvn1Yv359wooj8dfuDkCrSsy2pTL7OUjcTfCZJ8T93P3x55TCFWIROLxJ8GsTks5iDonnn38eVVVV+J//+Z+uXzBjxozBe++9l7DiSPy1J2qOBM9DV7UV3oKpQ580NwgWowo7XFchdHo3Iq01gl+fkHQV83P5F198gR07dkCpVIL9YTx6bm4urFZrwooj8dfu8idkZJPcVgk26EHAMCbu546FUiFFVKqEd/hNYHf+L5T3PAuGo6U6CBmqmJ8kpFIpIpFIt9fsdvsVI51Icmt3BaCO90Q6nof23KfwWiYndMhrf3INSlQFc8EodAgc3ixaHYSkk5j/Rc+ZMwdPPfUU6uvrAQAtLS1YvXo15s2bl7DiSPzZnP64T6ST26vABt0IGIrjet6Bys1WoMbqhnTCbIRO70Kk7YKo9RCSDmIOiSeffBKFhYVYsGABnE4nysrKYDKZ8NhjjyWyPhJndmf8h79qqr+AN2+SqE8RAGDOVqDZ7kFYooT0qhnw734DfDQsak2EpLqY2x3q6uowatQo/Mu//AsikQhmz56NkpKS/t9Ikkq7K4BriuP3JCHtqIPE2wLnqNlxO+ega5FwMKizcLHVgxGFExBpqkTw2KeQT7pT7NIISVn9hgTP81i2bBk2bdqEvLw8mM1mWK1WrFu3DgsXLsTzzz+fkOGUJDEcrgA0cXyS0Jz/Aj7zRFFGNPXElK3AhSYniiwaSCfcjsC+dyAZORmcPl/s0ghJSf22D2zYsAEHDx7Ehg0bsGvXLmzYsAG7d+/G+++/j0OHDuH9998Xok4SB+FIFN5AGEp5fDquOW8b5PZz8JlK43K+eMjLVuJCc+fG8axSB+mYmzqbnXhac4yQweg3JDZv3oxnnnkGEydO7Pb6xIkTsWzZMmzeTKNIUkW7KwCNUgqWjc+Tn7p2D/ymUoBL/N4UsTJq5Wh3B+APdo7E44omAeEQgidpzxNCBqPfkKiursbUqVN7/NrUqVNRXV0d96JIYrTHsamJCfmgbPwOPnNy7XHOcixMuizUWd0AOjfFkk4sQ/DIx4i6WkWujpDU029IRCIRqNXqHr+mVqtp6fAUYnf6oYnTRDrlxQMI6oYjKuv5Z0NM5mwFapt/3F+CVRshGTkV/q/eBM/zIlZGSOrpt3E6HA5j//79vf7j+ukEO5K82l1x2myIj0JduweukTOHfq4EMOuVOFTV0u01yaipCHzzLlzHvgQKrhepMkJST78hYTQasWzZsl6/bjAY4loQSZzWDl9c9pHIaj0NnpMhrMqNQ1Xxl62Rw+0LweMLQfXD/TIsB9nEObB/+Xco7hoDVk0/t4TEot+Q2LlzpxB1EAHYnQGMsmiHfB5V7Z7OlV6TdOgzwzIwZytRZ3WjtCi763VWa4Zy7BT49rwJxdz/l4ZuExID2jk+g3T2SQyt45rztELmakDAkNy7wJl0Waj9YSjs5dTjb0bU1YZQ5R4RqiIk9VBIZBC7KwCtamjNTar6b+A3lgBscm/sk5utQF3LlSHBcBLIrpmLwIF/IOqifbEJ6Q+FRIYIhCIIhiJDm0gXCULZeBA+0/j4FZYgerUcvkAEbm/oiq+xWnPnaKddf6ZJdoT0Q7CQqKmpweLFi1FWVobFixfjwoULvR57/vx5XHPNNVizZo1Q5aW9dtfQd6RTWI8hrDQhmqWLY2UJwjAwZytQZ73yaQIAJMVTwQd9CB7fLnBhhKQWwUJixYoVWLJkCT777DMsWbIEy5cv7/G4SCSCFStWYPZs8ReMSye2Dj90qqH1R6jqvobfNC5OFSWeSa/AhV5CgmFYSK+9A8HvP0GkrVbgyghJHYKEhM1mQ0VFBebPnw8AmD9/PioqKmC326849s9//jNuvfVWFBUVCVFaxrANsdNa4m6CxGtDUDcijlUlVm62AvU/zLzuCavUQzpuJnxfvgo+FBCwMkJShyC9j01NTcjNzQXHda4UynEczGYzmpqaus2zOHPmDPbt24d33nkHr7766qCuZTQmZgawyaRJyHmF4gtFkZOthF6v7PO43r4uqz6IaP4EqLV9vz+ZqFQyBL9vBDgOeo286/Vu96ifAoezATj0HkwLHhehysRI9Z/XWGXCfYp9j0kzRCUUCuE///M/8Yc//KErTAbDZnMjGo3v0gsmkwatrT03W6SKuqYO6FRyOBzeXo/R65U9fz0ShOXCAbSPuwdRd2p94jbrFTh5thUTRnV+GOnpHvkxt8Lz9d9xcd82yEpuFqPMuEqHn9dYZMJ9CnGPLMv0+eFakJCwWCywWq2IRCLgOA6RSAQtLS2wWCxdx7S2tqKurg6PPPIIAMDpdILnebjdbjz77LNClJnW2jr8GG4e3CcShfU4QiozovKhT8QTmknfOV/iUkj0hJHIIJu0AIH974PLGQHOOFzACglJboKEhNFoRGlpKcrLy7Fw4UKUl5ejtLS0W1NTfn4+Dhw40PX3tWvXwuv14qmnnhKixLRnc/qhGeQcCVX9N/DnJM+eEQNh1iux90Rjv8exWlNn/8SOl6G6axUYuUqA6ghJfoKNblq5ciXWr1+PsrIyrF+/HqtWrQIALF26FCdOnBCqjIwUjfLocAehHUTHtcTTAomnBUF9UfwLE4BWJUU40nn//ZEUjANrGgXf56+Aj9LClYQAAvZJFBcXY+PGjVe8/pe//KXH4x9/PH06EcXmcAegzJJAwg38M4Gy4VsEjGOTZnvSAbtsvsTVamO/h0uvmoHgoY8Q+Ob/Imv6gwIUSEhyoxnXGaCtww+9Wt7/gT8VDUN58buUmhvRE7NegZrL9pfoC8OykE2aj3D9cQRooh0hFBKZoK3DB+0gJtJltVYgkqVHJEsf/6IEdGm+RKz7DTHSLMim3o3QsW0InT+Y2OIISXIUEhmgzeGHdhA70nV2WF+VgIqEpVZIAQZod/ljfg+r1EE29W74976DcMPJBFZHSHKjkMgALQ7fgGdbs34HZB21CGQXJ6gqATEMcrOVPS4d3hdWa4Zs8kL4vnwN4cbTCSqOkORGIZEBWh0+6NUDCwnVxQOde0Zw8dkTW2xmvQIXBhgSAMAZCiGbtAC+z1+hoCAZiUIiAwy445qPQtlwIGXnRvQkN1uBOuvgZuNzOcN/DIra7+NfHCFJjEIizYUjUTg9wQE1N8ns1eBZCcJKUwIrE5YiS4osGYemNs+g3s/lDId8yl3wf/UGgmdoVzuSOSgk0tylLUs5NvZ9JLo6rNNsD2hzthLnGhyDfj+bnQ/ZtF8gePhj+A9+QBsWkYxAIZHmWh3+AfVHMEEPstpOd06gSzN5BgXO1juGdA5WbYT8xvsRqT8O32cvgQ/2vmAiIemAQiLNtTh8A+qPUDYdRlA/ArwkK4FVicOsV6C+xYVwZGhPAIxcBdn19wGsBJ4PV9CmRSStUUikOavdG/uOdDwPVf23aTE3oidSKQejNgsNrYPrl7gcw3KQTZgNyehp8H7yAgLHPqXmJ5KWKCTSXLPd223Dnb6w9lowYR9CmsIEVyUeS44aNY2xLdERC0nBOMhvegDhs9/Au+UPiDpb4nZuQpIBhUSaa3X4kB1jn4Tk/F4E0rDD+nKFJjVqmuIXEkDnNqiyG34B1jgcno9WInC0HHw0HNdrECIWCok0FuX5mOdIMOEAJPWH0mpuRE9M2Qo4vUF4fKG4npdhWEhHTYH8pgcRrj0K78ZnaPIdSQsUEmnM4QogS8ZBJu1/mW9F8/eI6goRlaX3ZjssyyAvW4mapsRsCcmq9JBNvRuS0dPg3/m/8O54GVFna0KuRYgQKCTSWLPdC4M2tlFKqvqvEbaMT3BFySHXoMT5xo6EnZ9hGHCWEshn/AqMXA3PRyvg//Z98IGhd5gTIjQKiTRmtXuRHUNTk8TVCM7vQNQwUoCqxGcxKFHT7EKiByMxnBTSMTcg65aHEO1ohnvDUwgc/wx8hPorSOqgkEhjTTZvTBPpumZYs5nx46BUSKGQS9BkF+aTPZOlhmxiGeTX3YdwzXfwbHgaoeoD4GPd4IIQEWXGb4UM1WjzIFvTd3MTEwlC2XQk7Tusf8piUKL6YuKanHrCak2QT70b0gmzETi0Cd5NqxGxnhO0BkIGikIijTXbvTBq+25uUjQfRUidh6hcI1BVySHfqMI5gUPiEi5nBOTTHwRXMA6+HS/D98WriLrtotRCSH8oJNJUMBSB0xPsd/irqm5fyu9hPRg5uiw4PUG4vfEdChsrhmEgKZwA+YyHAU4CzwfPIHB0K/iIOPUQ0hsKiTTVbPfCoMkC28fqr1JnPbiAE0HdcAErSw4My8BiVOFcAkc5xVSHRAZpyc2d8yvqjsPzwTMIN54RtSZCLkchkaaabF4YdX33R6jqvu7si2Ay88fAYlQOeVXYeGFVesim/BzSMTfC/+Vr8O1+g4bMkqSQmb8dMsDFVjey+1iziQn5oLAeg8+UWR3Wl7MYVWho9SAYiohdCoAf5lfkjYX8lofAB7zwbPw9wnXHxS6LZDgKiTRV3+JGTh9PEsqLBxHUDQcvTe8Z1n2RSTnk6LIGtfd1IjFSOWRX3w7pxDnw73kTvj1vgg8FxC6LZCgKiTR1sc0Dk07R8xd5Huq6vfCZMmOGdV/yjSpU1TnELqNHXM4IyG/+P+Dd7fB8+J+ItF0QuySSgSgk0lAgFEGHJ9hrc5PcVgmeYRFWWwSuLPkUmlSobnQiEk3OiW2MVA7ZtXdAUnw9vJ/8EYETO2gSHhEUhUQaamzzwKjtfWSTuvYr+M0T0npJ8FgpsqTQqqSosyZXk9NPSQrGQX7j/Qif3gXf52vBB31il0QyBIVEGqqzumDS99zUxHlaIXPUwp+Ge1gPVoFJhcpah9hl9ItVZUN2wxKABzwfLkfEflHskkgGECwkampqsHjxYpSVlWHx4sW4cOHCFcesW7cO8+bNw4IFC3DXXXdh7969QpWXVmqbXTD10mmtrt0Dv6kUYCUCV5W8Ck1qVDU4EE3SJqfLMZwEsqtvh6T4Oni3Po/Q+e/ELomkOcFCYsWKFViyZAk+++wzLFmyBMuXL7/imIkTJ+KDDz7Ali1b8Pzzz+PJJ5+E3+8XqsS0Udfihjn7yicJJuSDsukQfOarRagqeakVMqgUUtRZ3WKXEjNJ4QTIr7sHgW/ehf/gB7S/NkkYQULCZrOhoqIC8+fPBwDMnz8fFRUVsNu7r1dz8803Q6Ho/OVWUlICnufhcDiEKDFtRKM8Lrb2HBKqhm8R1I1AVKYWobLkNsyswuna1Fo/idXlQX7TA4jUH4dvx1rwIfpAReJPkJBoampCbm4uOK5zhzSO42A2m9HU1NTrezZt2oThw4cjLy9PiBLTRqPNA7VChizZT5qTohGoa/fAlztRnMKS3HCTBlX1HUk7yqk3jFwF2XX3AnwU3k3PIeq2iV0SSTNJ2TB98OBBvPTSS3jzzTcH/F6jMTGfkk2m1Fgl9VhNO4bnaaDXK7u9zl3YDyh1yMod1uf71TFsUpTqerpHtVoOg1aO1o4ArhppEKGqoeFv/jk8p7+Bd/OzCNz3O5jyR4tdkiBS5d/lUIh9j4KEhMVigdVqRSQSAcdxiEQiaGlpgcVy5Tj9o0eP4j/+4z/w6quvYtSoUQO+ls3mjnsHpMmkQWtrcg+RvORYpRV6lQwOh/fHF3ke5pPb4CqYipC795m7arUc7j6+ng76useCHBUOnGxCXnZsW74mnfxJ4FgVmt57FrKbH4J05GSxK0qoVPp3OVhC3CPLMn1+uBakucloNKK0tBTl5eUAgPLycpSWlsJg6P6J7fjx43jyySfx8ssvY/x4mg08GOcudiA/p/tSG/K2CoCPIqTNvNVeB2KYWYOaZicCSbKW02BweWNhuO1+BPb9DYFj22jiHRkywUY3rVy5EuvXr0dZWRnWr1+PVatWAQCWLl2KEydOAABWrVoFv9+P5cuXY+HChVi4cCEqKyuFKjHl+YNhWNt9yL2805rnoan+HL68a2nyXD/kMg5mvRJVdeIuHz5UUmM+ZDfej9Dp3Qjsewd8NHVDj4hPsD6J4uJibNy48YrX//KXv3T9+cMPPxSqnLRU0+RCbrYCEu7H7JfZz4ELOBEwZEYb9VCNyNPgRE0bri5OvX6Jy7EKLeQ3LEHwyBb4PvsTFLMfAyNN0WY0IiqacZ1GKuvaUWDq3raorf4MXsvPMnbPiIHKNyrR2u5HhycodilDxkjlkE29C2BYeDf/F6KedrFLIimIfnOkkTN17Sg0/dgfIWs/D87bhoBhjIhVpRaOYzE8V4OT59NjKCnDcpBeXQbWPArej1cjYqsXuySSYigk0kQoHEFtswuFOT88SfA8tFXl8OZPBlhO3OJSzEiLBsfP25Aufb4Mw0A6ehokJdPhLV+DcD1tZERiRyGRJs42dMCkV0Au6wwEuf0sOL8DAWOJyJWlHoNGDinHojbJV4YdKEnBOMgmL4R/118QPPWl2OWQFEEhkSZO1tgxPPeHSTc8D23VVnjzp1BfxGAwDEZZtPi+qk3sSuKOMxRCdsMvETz+Kfxf/51GPpF+0W+QNHG82oaivM6QyLIeAxMJUl/EEIzI1eCC1QmPPyx2KXHHqrIhv/F+RFpq4Pv0/wMf9Pb/JpKxKCTSQFuHDw53APlGFRANQ1e1FZ6CaTQvYgikUg7DTGocP5d+TxMAwEizIJt6NyDLguejVYg6msUuiSQpCok0cOycDcX5WrAsA1XdPkRkGoR0fa/RRPpXXKDD9+fakK6rcDMsC9n42ZCMmATP5ucQrjsmdkkkCVFIpIGDp60YXaADG3BBc/5zeIbdKHZJacGgzUKWjEP1xdSegd0fyYhrOju0v/orAoc3094UpBsKiRTncAdQ3+LGSIsW2rOfIGAYi4gitWcLJ5PiAh2+q2wRu4yE4wyFkN/0AMI1h+Db/ifw/tTZgIkkFoVEitt/qhljC/VQuuqQ1XKqc0QTiZvhJjVsTj9a231il5JwTJYGsmmLAWkWPB8tR8R6TuySSBKgkEhhPM9j34lmjBuuhf7k+/AMuxG8JP33gxASy7EYna/DwTPp/zQBdM7Qlo27DdKrZsC7/U8IHN0KPkrNT5mMQiKFnW9ywh8IY5znAKKSLFrEL0FGF+hwrsEBtzckdimC4fLGQj79QYTPfwdv+X8j6moVuyQiEgqJFPbloQZMHxaBpnYP3CNm0JDXBJHJOIzI0+JQhjxNXMIqtJBNWww2Ox+ej1YiWLGTOrUzEIVEimp3BXCy2orprk/hKbweUXn6b+MoppJhehyrtsEfzKwZygzDQlp8PeTTFiN06kt4t/4BEUej2GURAVFIpKhP99fi/pxTgEwJf06p2OWkPZVCioIcFQ5nwEinnrAaE2Q3LgFnHAHvpufg378BfDD9O/MJhURKancF0Hb6IMaGK+EqupWamQRSWpSNw5WtCGTY08QlDMNCMnIysm75P4jaG+DZ8FRnE1Q0/ZYuIT+ikEhBn+86jF8qvoar+HbwEtptTCgapQz5RhW+O52ZTxOXMFkayK69A7LJixCq3AvPht8hVLWPwiJNUUikmNqGVkxq3Ah//s8Q1uSLXU7GGTcyG4erWuFNw4X/BorVWyC//j5Ix89G8OQX8Lz3Hwgc3w4+4BG7NBJHgu1xTYYuEg6j49O1UGpyELJcI3Y5GUmtkKEoT4O9xxtRdt1wsctJClzOcHA5wxFtb0T4wmG4D2+GpOhnkF51C7i8MWBoufqURiGRInieR/XH65AFP9iSudQPIaLxRQZsO1CLn401waRXiF1O0mCz8yHLzgfvdyN8sQL+r/4KhIOQjJwCSdHPOgODk4pdJhkgCokUwPM8Wr94G4ztAgLjF0DC0bdNTDIZh/FFBnx+qAG/nDWG8vonmCw1pMXXQVp8HaLOVkSsZxH45l1EXW3gzMXg8q/q/F9TERi5qv8TElHRb5skx/NRuPf8He7zx9A8vAwFarXYJRF0zsKuaXbiVI0dE0bRgoq9YbUmsFoTMOZG8EEfovZ6RO0NCNccQrSjGYxcBVZvAZtdAFaXB1ZjAqPJAas2gJHSoIxkQCGRxPhIGN7df4XtfBVO6m/FBEuO2CWRHzAsgyklZuw6ehFFFg3UCmpG6Q8jU4DLGwsubyyAzidk3usA72pD1GNHuOEEeJ8TvNcJ3tcBsBwYpQ6sMhuMUg9GnQ1WZQCj1INVZSMsHw4+KgHDUp9HIlFIJKmo3wXfjrVotnnxreRG3DCaRjIlG4M2C8UFWnzybS3uu200NTsNEMMwYFTZgCob3E++xvM8EPKD97vB+13gA57Ovo6OFiDgAe93oSHgQtTv6QwNjQmMLhec3tL1ZMKojWDomzJkFBJJKGI9B98X61CHfHwTnoRbJhWAYemHPRmNG2HArqMN2H+qGTdMyBO7nLTBMAwgU4CRKQCtqcdj9Hol2m3OH54+HOA97Qhbz4E/fxC8qw18KAA2uwCcaQS4nJFgzSPB6gvoyWOAKCSSCB8JIXB4M0Knd+OUfBKOubJx6zX5kHA//ZxFkgXLMrhxfB4+P9wAU7YCowt0YpeUURhOAkZtANRX9gvxQR+izlbwTitC5w8iemQLeL8LXM5wsHljIckbCy53NHWe94NCIkmEG07Cv+/vgEKHXdJb0erhcOu1FkglFBDJTpElxU0T8rDt21rcPaMYBSb6pZMMGJkCXM5wIOfH+Sx80I+o4yKi7U0IHP4Y0fZGsOqczhFXlqvAWcaCVerFKzoJUUiILNxcheChjxF1tiA44kZ8eJqBWiHFLdeYwXH0WJwqjDoFpo3LxYdfVWPRzaMwPJdGoSUjRpbVOfzWXAwA4KMR8B1WROz1CJ76AtE9b4JRaMHllUCS3xkajDono/s2KCREwIeDCF84jODJz8F7HOBGTcUZ9TTsOdyM8UUGjCnU0WS5FJRnVOGG8XnYtPc8brnGgmtGm+jbmOQYlgOTnQ82u3NgCM9HwTtbEbXXI1S1D4H97wGsBFzuaHCWEnDmYrDG4WAyaK6SYHdaU1ODp59+Gg6HA3q9HmvWrEFRUVG3YyKRCJ577jns3bsXDMPgkUcewb333itUiQnFhwKIXKxAqOYQwrVHwOoskAyfiEbGgl3fNyEcseG2SQXQqWn70VSWa1Bi1s8K8c2pZpxrcOL2qcOgU8vELovEiGFYMLpcsLpcYOSPw3Sj9gZELlYgeOoL8G57Z4e4eRS4nCKwOcPB6vPBSNLz+yxYSKxYsQJLlizBwoULsXnzZixfvhzvvPNOt2O2bt2Kuro67NixAw6HA4sWLcINN9yAwsJCocqMC57nwbvbELHVIdJyHpGmSkRtdWD1+eDMo8BM+39Q1RrB0WNt8PjrML7IgKI8DT09pAmNSobZkwtxps6Btz89g5JhevyspHMJD/oWp5ZLw3RZVTYw7GpI0dkSEO2wgu9o7uwQ/74cvMcORmUAq88HaygAp8sDozWD1ZrBKHUpvX4Vw/M8n+iL2Gw2lJWV4cCBA+A4DpFIBNdffz127NgBg+HHUQmPPPII7rrrLsyZMwcAsHr1auTn5+PXv/51zNdqb/cgGo3vLRmNaths7s4N4cMB8CEfEPSBD3jBB9zgfU5EvR2IeuyA246oqw2MVAZGbUJUmQ2/3Ag7k41WVxiNbW7YnQHk6BUYZlIhN1uZNOGgUsng8QTFLiOhhL7HYCiMC80uNLS6wbEsCk0q5BpUMGjk0ChlkMtYxPu7r9Uq4HSm/4ZAyXSfPB8F73Eg6rWD93SA9zkBvxO8zwU+5AOTpQGj0IFRasFkaTv/LleCkavAyBSds8ulWYBU3vlEwsnBsGzX755EYlkG2dm9D7YQ5EmiqakJubm54H4YyslxHMxmM5qamrqFRFNTE/Lzf5w0ZrFY0NzcPKBr9XWzQ2E0UkckGZzrRbhmpixokQn3KfbvntR9BiKEEJJwgoSExWKB1WpFJNK57WMkEkFLSwssFssVxzU2/rjJelNTE/LyaBYrIYSIRZCQMBqNKC0tRXl5OQCgvLwcpaWl3ZqaAGDOnDnYuHEjotEo7HY7vvjiC5SVlQlRIiGEkB4I0nENANXV1Xj66afhdDqh1WqxZs0ajBo1CkuXLsUTTzyBq6++GpFIBKtXr8bXX38NAFi6dCkWL14sRHmEEEJ6IFhIEEIIST3UcU0IIaRXFBKEEEJ6RSFBCCGkVxQShBBCekUh0YeamhosXrwYZWVlWLx4MS5cuCB2SXG3Zs0azJw5EyUlJaiqqhK7nIRob2/H0qVLUVZWhjvvvBP/+q//CrvdLnZZCfHoo49iwYIFWLRoEZYsWYLTp0+LXVLCvPLKK2n7cztz5kzMmTMHCxcuxMKFC7F3717xiuFJrx588EF+06ZNPM/z/KZNm/gHH3xQ5Iri77vvvuMbGxv52267ja+srBS7nIRob2/n9+/f3/X3//7v/+Z/97vfiVhR4jidzq4/f/755/yiRYtErCZxTp48yT/88MP8rbfempY/t8n075GeJHphs9lQUVGB+fPnAwDmz5+PioqKtPsEOmXKlCtmvqcbvV6P66//cQWla6+9ttvM/nSi0Wi6/ux2u9Nys5xgMIjVq1djxYoVaXl/ySZzds4YoFgXJSSpJRqN4r333sPMmTPFLiVhfv/73+Prr78Gz/N44403xC4n7l566SUsWLAAw4YNE7uUhPr3f/938DyPyZMn49/+7d+g1WpFqYOeJEhGefbZZ6FUKvHAAw+IXUrC/Nd//Rd2796NJ598Ei+88ILY5cTV0aNHceLECSxZskTsUhLq3XffxZYtW/Dhhx+C53msXr1atFooJHoR66KEJHWsWbMGtbW1+NOf/gSWTf8f/UWLFuHAgQNob28Xu5S4+e6773D+/HnMmjULM2fORHNzMx5++GHs27dP7NLi6tLvGZlMhiVLluDIkSOi1ZL+/1IGKdZFCUlqePHFF3Hy5EmsW7cOMll6bjPp8XjQ1NTU9fedO3dCp9NBr9eLV1ScPfLII9i3bx927tyJnTt3Ii8vD3/9618xffp0sUuLG6/XC5fLBaBzl8tt27ahtLRUtHpo7aY+9LYoYTp57rnnsGPHDrS1tSE7Oxt6vR6ffPKJ2GXF1dmzZzF//nwUFRUhK6tzm5rCwkKsW7dO5Mriq62tDY8++ih8Ph9YloVOp8NTTz2F8ePHi11awsycOROvv/46xo4dK3YpcVNfX4/HH38ckUgE0WgUxcXFeOaZZ2A2m0Wph0KCEEJIr6i5iRBCSK8oJAghhPSKQoIQQkivKCQIIYT0ikKCEEJIrygkCCGE9IpCgpAYzZw5ExMnTsSkSZO6/jty5AhKSkoQDod7fV99fT2uuuoqrFy58oqvlZSUoLa2FgCwdu1ajB8/HpMmTcKUKVPwi1/8AkePHk3U7RASEwoJQgbg9ddfx9GjR7v+i2WC0+bNm6HT6bBt2zYEg8E+j507dy6OHj2K/fv34/rrr8dvfvObeJVOyKBQSBCSYJs2bcJvfvMbSCQS7Ny5M6b3SCQS3HnnnbBarWm3PD1JLRQShCTQoUOH0NzcjHnz5mHu3LnYtGlTTO8LBoPYtGkT9Hq9aEtEEwLQfhKEDMhjjz3WtcfIddddh2XLlvV5/Mcff4xbbrkFOp0O8+fPxwMPPACbzQaj0djj8du3b8fu3bvh8Xig1Wrx8ssvQyKhf6ZEPPQkQcgArFu3DocOHcKhQ4fw6quv9nms3+/H9u3bceeddwIAJk2aBIvFgq1bt/b6njlz5uDQoUP4+uuvMWbMGJw6dSqu9RMyUBQShCTI559/DrfbjVWrVuGmm27CTTfdBKvVis2bN/f7XoPBgFWrVmHt2rVoaWkRoFpCekbPsYTEQTAY7NqgCgCkUik2bdqEu+++G08++WTX61arFffccw8qKytRUlLS5zmLi4tx880344033ui3WYuQRKEnCULiYNKkSZg4cWLXf5s3b8a3336Lf/7nf4bJZOr6b8KECbj55ptj7sB++OGH8Y9//AM2my2xN0BIL2g/CUIIIb2iJwlCCCG9opAghBDSKwoJQgghvaKQIIQQ0isKCUIIIb2ikCCEENIrCglCCCG9opAghBDSKwoJQgghvfr/AY4860J9lfY+AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set()\n",
    "x0 = dat['FLAIR'][dat['GOLD_Lesions'] == 0]\n",
    "x1 = dat['FLAIR'][dat['GOLD_Lesions'] == 1]\n",
    "sns.kdeplot(x0, shade = True, label = 'Gold Std = 0')\n",
    "sns.kdeplot(x1, shade = True, label = 'Gold Std = 1')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46ad3543-6d3e-42cc-90c9-500df2f13697",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### The ROC curve\n",
    "+ Consider a given FLAIR threshold, say $c$. \n",
    "+ The true positive rate given that threshold is $T(c) = P(X \\geq c ~|~ Y=1)$ \n",
    "+ The false positive rate is $F(c) = P(X \\geq c~|~ Y=0 )$. \n",
    "+ The function $f\\rightarrow T\\{F^{-1}(f)\\}$ is the ROC curve.\n",
    "+ The ROC curve is typically displayed in the plot: $(f, T\\{F^{-1}(f)\\})$\n",
    "+ An empirical estimate of the ROC curve requires empirical estimates of $T$ and $F$ which can be done non-parametrically or parametrically. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f2a6a67-218d-498f-b985-a5bbd1489864",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### The ROC curve satisfies:\n",
    "\n",
    "1. Always starts at the point (0, 0). \n",
    "2. Always ends at the point (1, 1).  \n",
    "3. Is monotonically increasing (always moves laterally or upward).\n",
    "4. A uniformly better ROC curve lies entirely above a worse ROC curve. T\n",
    "5. Is always worse than the theoretical limt of the ROC curve, the upper part of the box at the points (0, 0), (0, 1), (1, 1).  \n",
    "6. If the prediction value is a random uniform number indpendent of the gold standard, the ROC curve is the identity line from (0, 0) to (1, 1). \n",
    "7. The ROC curve is invariant to strictly increasing monotonic transformations of the predictor. \n",
    "8. The ROC curve is an identity line whenever the prediction is independent of the gold standard provided the prediction is continuous. \n",
    "9. The ROC curve for a prediction as a test for one minus the golad standard flips the ROC curve over the identity line."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c774f8e9-0a7f-49f2-a831-353c7b47101e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### Non-parametric estimation\n",
    "\n",
    "We can estimate these probabilities using the fraction of times a FLAIR values is above the threshold in the two gold standard groups\n",
    "\n",
    "$$\n",
    "\\hat T(c)\n",
    "= \\frac{1}{|\\Gamma_1|} \\sum_{i \\in \\Gamma_1} I(x_i \\geq c)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\hat F(c)\n",
    "= \\frac{1}{|\\Gamma_0|} \\sum_{i \\in \\Gamma_0} I(x_i \\geq c)\n",
    "$$\n",
    "\n",
    "where $\\Gamma_1 = \\{i ~|~ y_i = 1\\}$ and $\\Gamma_0 = \\{i ~|~ y_i = 0\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "159bcfeb-c13a-4245-87f9-a1c5c6b02d65",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD7CAYAAABgzo9kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoDklEQVR4nO3deXRUdZ428OfeWysQCAlJqBAwgBgCArIpIgGEQAIEQitLN8iItDDT+mq3dtuNPSPLwdeZzJn5Q3tk5u05shltFRWQSCNNOxIiEHYChJ1gWIosFQJkqe3e3/tHWkZkSVVSVTdV9XzO8ZwOfSt5flR4cvOrW98rCSEEiIgo4sh6ByAiouBgwRMRRSgWPBFRhGLBExFFKBY8EVGEYsETEUUoFjwRUYQy6B3gh65dq4em+X9Zfnx8BzgcdUFI1HZxzdGBa44OLV2zLEvo3Ln9Pf//NlXwmiZaVPDfPzbacM3RgWuODsFYM7doiIgiFAueiChCseCJiCJUswWfl5eHcePGIS0tDadPn77rMaqqYvny5cjMzMSECROwfv36gAclIiL/NFvw48ePxwcffIBu3brd85jNmzejvLwc27Ztw8cff4w//OEPuHTpUkCDEhGRf5ot+GHDhsFms933mC1btmDmzJmQZRlxcXHIzMzE1q1bAxaSiIj8F5DLJO12O5KTk299bLPZcPXq1UB8aiKiiPDN4csoPl7xv38gBPp5j2GoZx9OD1+Ahx4ZEvCv2aaug4+P79DixyYkxAQwSXjgmqMD1xwZDp6pxqWqOvTs1gkd1WvIaNiOZO8lXDZ0h9qha1DWHJCCt9lsuHLlCgYOHAjgzjN6XzkcdS262D8hIQZVVTf9flw445qjA9ccOTxuFd27tMP/eegKXPs/BxQDzKOfQ1raaCQmdmzRmmVZuu+JcUAuk8zOzsb69euhaRpqamqwfft2ZGVlBeJTExFFhHitGk87P4Gr+GMYUh5G+5lvwdR3DCRJCtrXbPYM/s0338S2bdtQXV2N5557DrGxsfjyyy+xcOFCvPzyyxgwYAByc3Nx5MgRTJw4EQDw4osvonv37kELTUSR6ZvDl3HwTDU8blXvKAEjCy+GevZjhnsfXJIFlswXYOg5PKjF/j2pLd10m1s0vuOao0O0rTnvg4O4VFWHlISWvx7XliSpdjzp+iviRA1OGfqiccDTyBje547jWvo8N7dF06ZeZCUi6tmtE16dOUjvGK0iPC649n0Gz7G/QGrfGZaMVzGsx8CQ52DBExEFkPdyKZyFqyFuVsHYbxzMj86EZLLqkoUFT0RBc8e1380or6xD75ROQUwUPMJVD9eej+E5VQipUxKsU1+HwZamayYWPBEFTfHxCpRX1qFHom976j0SO2DM4JQgpwo8z4UDcBW9D9F4A6ZHpsA0JBeSwaR3LBY8EQVXj8QO+N1c39+lGU4vLGsN1+HalQ/v+X2Q47vDmvUrKAmpese6hQVPROQnIQS8Z3bBuftDwOOCafjTMA2aBEluW5XattIQke783Te/H3+2Z8KFVueAc+daqBdLICc9CMvoBVA6+//O/VBgwRPRbfzdN7+fHokd8Fj/pACk0p8QGjyl/wPX3vWAEDCPnAtjv/GQ5LZ73yQWPBHdwd9980in1V6Fs3AV1KunoXTrD8vo+ZBjEvSO1SwWPBHRPQhNhbtkK9wHNgCKCZYxP4fhoVEhGTMQCCx4oih1r732SNw3bwm1+js4C1dBq/4OhtShMI+aB7ldrN6x/MKCJ4pS99prj6R985YQXjfcB7+A+8gWSJYOsGS+CGOv4XrHahEWPFEU41777dSrZ5rO2mvtMDz0BCwjfgbJEr6/zbDgiSjqCY8Trr2fwnP8r5A6xME66dcwdB+gd6xWY8ETUVTzXjrWNBysrgbG/uNhHv60bsPBAo0FT0RRSTjr4NzzEbyniyB36grLtN/D0PXOWe3hjAVPRFHHc34fXN++D+Gsg+mRHJiGTGsTw8ECjQVPFCFaMpo32i6H1Bpq4fo2H96y/ZDjH4B10q+hdHlA71hBw4InihAtGc0bLZdDCiHgPV0E556PAK8LpkdnwDQwu80NBwu0yF4dUZThZY930m5WNQ0Hu3QMSteHYBn9HORYm96xQoIFT0QRSQgNnuN/hWvvp4AkwfzEPBj7PQlJarvDwQKNBU8UZjhioHlq7RW4dqyGWnEGSvcBsIx6FnJMF71jhRwLnijMcMTAvQnNC/eRP8N9YBNgNMMydiEMfUaGzXCwQGPBE4Uh7rXfSa2+AOeOVdAc5TD0Gg7zyGcgtwvPG3gHCgueiMJa03CwTXAf+TMkSwwsE16CsedQvWO1CSx4ohBpya3wjCYFHrd6259xr/1/ea+ehnPHKojrV2FMy4B5xE8hmdvrHavNYMEThUigboXHvXZAuBubhoOV/hVSTBdYJ78GQ0p/vWO1OSx4ohDyd+88ISEGVVU3g5go/HjLS+DcuQai/hqMD09sGg5mNOsdq01iwRNRWBDOOjh3fwjvmV2QY5Nhzf1HKEkP6h2rTWPBE91HS/bN74V75y0jhIC3bB9c3+ZDOOthGjINpsFTISlGvaO1eSx4ovsI1L45wL3zltAaauEqWgfvhYOQu6TCOvk3UOJ76B0rbPhU8GVlZVi8eDFqa2sRGxuLvLw8pKam3naMw+HA66+/DrvdDo/HgxEjRuCf/umfYDDwZwiFN15zHnpCCHhP7YRzz58A1QvzY7NgHJAFSVb0jhZWfBrKsHTpUsyZMwdfffUV5syZgyVLltxxzH/913+hd+/e2Lx5MzZv3ozjx49j27ZtAQ9MRJHNc+0qGrf8G5yFq6DE90D7GStgGjSZ5d4CzZ5eOxwOlJaWYvXq1QCAnJwcrFixAjU1NYiLi7t1nCRJqK+vh6ZpcLvd8Hg8SErir6MUGN8cvoyDZ6rvuCY82LhvHjpC0+A5/hdc2v85BCSYR/0djOljo2o4WKA1W/B2ux1JSUlQlKafnoqiIDExEXa7/baCf+GFF/DSSy9h1KhRaGxsxNy5czF0qH/vJouPb/k/pISEmBY/NlxF05oPnqlG2eXr6NkttG89753SCWMGp+j6dx0Nz7O76iKq/rwSrsunYe09BAmT/x6GjtE1HCwYz3PANsi3bt2KtLQ0rF27FvX19Vi4cCG2bt2K7Oxsnz+Hw1EHTRN+f+1ovFY42tbscavo2a0TXp05SJevr9ffdaQ/z0L1wn3kS7gPboZktMDy5CJ0fXwiqqvrgAhe94+19HmWZem+J8bN/u5js9lQUVEBVW361VhVVVRWVsJmu31gfn5+PqZNmwZZlhETE4Nx48ahuLjY78BEFB3UqjI0bFgO9/4NMPQcinaz3oIxiic/BkOzZ/Dx8fFIT09HQUEBcnNzUVBQgPT09Nu2ZwAgJSUFhYWFGDhwINxuN3bv3o0JEyYELThFpvvNOu+dEt2TASOF8Lrh2r8BnqNbIVk7wTrxlzCkDtY7VkTy6dWLZcuWIT8/H1lZWcjPz8fy5csBAAsXLsTRo0cBAL///e9x4MABTJ06FdOnT0dqaipmzZoVvOQUkb6/7vzHeiR2wJjBKTokokDyXjmJ+k/fgKfkzzCmjUb7mf+X5R5EkhDC/03vIOEevO8idc15HxwEgLtedx6pa76fSFmzcDfCVfwJPCf+B1JMAiyjn4OhW7+7Hhspa/ZHsPbg+S4kIgoqb/lhOHeug2i4BuOALJiHPwXJwOFgocCCp6Dyd5YLrzuPHJrzJly7PoT37G7InbvBOuFFKIm99Y4VVVjwFFT+znLhvJbwJ4SA91wxXLs+gHA3wDQk92/DwVg3oca/cQo6znKJHlr9tabhYN8dgpzQE9YxC6DEddc7VtRiwRNRqwkh4Dm5A649HwOaCvOIn8L48ERIMscM6IkFT0Stot2ohLNwNdQrJ6Akp8My+jnIHRP1jkVgwRNRCwlNg+fYNrj2fQ7ICswZ82HsO4bvRG1DWPBE5De15hKcO1ZBqzoPpccjsGQ8C7l9Z71j0Y+w4InIZ0L1wn24AO5DmyGZ2sEy7h9g6P0Yz9rbKBY8EflErTzfdNZ+7RIMD46AeeRcyJbIH2UczljwRHRfwuuCa9/n8BzbBqldLKxZv4LhgUf0jkU+YMET0T15r5yAc8cqiJtVMKY/CfNjsyCZrHrHIh+x4InoDsLdANeej+E5uQNSxyRYcxbDkNxX71jkJxY8Ed3G+90hOHeuhWi8DuPASTAPm87hYGGKBU9EAACt8QZcuz6A91wx5LgUWLN+CSWhp96xqBVY8ERRrmk42B64vv0AwtMI07CfwDRoCoeDRQA+g0RRTKtzwFm0Dmr5EciJvWAd/XMocd30jkUBwoKngLjfvVQ5373tEUKD58Q3cBV/AggN5sd/BmP/CRwOFmFY8BQQ95r7zvnubY92/WrTcDD7KSjd+sGSMZ/DwSIUC54ChnPf2zahqfAc/Qqu/RsAxQDL6AUwpGVwzEAEY8ETRQHVcRHOwlXQqspgeGAwzKP+jsPBogALnvzCvfbwIlQP3Ic2w33oS0iW9rBkvgBDz+E8a48SLHjyC/faw4dacbZpOFjtFRj6jITl8TmQLPwhHE1Y8OQ37rW3bcLjgmvfp/Ac2w6pfWdYs1+FocdAvWORDljwRBHEe+k4nDtXQ9yshrHfeJgfncHhYFGMBR/F7rWffj/ca2+bhKserj0fwXNqJ6ROSbBOfR0GW5resUhnLPgodq/99PvhXnvb47lwAK6i9yEab8D0yBSYhuRCMpj0jkVtAAs+ynE/PXxpDdfh2pUP7/l9kOO7w5r1KygJqXrHojaEBU8UZoQQ8J7ZBefuDwGPC6bhT8M0aBIkmf+c6Xb8jiAKI1qdA86da6BePAo56UFYRi+A0jlZ71jURvlU8GVlZVi8eDFqa2sRGxuLvLw8pKam3nHcli1b8J//+Z8QQkCSJKxevRpdunQJdGaiqCOEBk/p13Dt/RQQAuaRc2HsPx6SxOFgdG8+FfzSpUsxZ84c5ObmYtOmTViyZAnWrVt32zFHjx7Ff/zHf2Dt2rVISEjAzZs3YTLxhR6i1tJq7U3Dwa6ehtKtPyyj50OOSdA7FoWBZn/8OxwOlJaWIicnBwCQk5OD0tJS1NTU3HbcmjVrsGDBAiQkNH3jxcTEwGzmbb6IWkpoKmp3fY76z96AWnMJljE/h3Xyb1ju5LNmz+DtdjuSkpKgKAoAQFEUJCYmwm63Iy4u7tZx586dQ0pKCubOnYuGhgZMmDABv/jFLzjzgqgF1Orv4Cxchbrq72BIHQrzqHmQ28XqHYvCTMBeZFVVFadOncLq1avhdrvx/PPPIzk5GdOnT/f5c8THt/wNNAkJMS1+bLhq7ZqNJiUgnyeUwilrS2heN2p3rsfN3RuhtOuIxKd/gw59H9c7VshF+vN8N8FYc7MFb7PZUFFRAVVVoSgKVFVFZWUlbDbbbcclJycjOzsbJpMJJpMJ48ePR0lJiV8F73DUQdOE34tISIhBVdVNvx8XzgKxZo9bBYCw+buL9OfZe/UMXDveg3b9KgwPjYJlxE/Robstotd8N5H+PN9NS9csy9J9T4yb3YOPj49Heno6CgoKAAAFBQVIT0+/bXsGaNqbLyoqghACHo8He/bsQd++ff0OTBRthMcJ57f5aPziLQjVA+vk38A69nlOfqRW82mLZtmyZVi8eDFWrlyJjh07Ii8vDwCwcOFCvPzyyxgwYACmTJmCY8eOYfLkyZBlGaNGjcKMGTOCGp58wxnubZf30jE4C1dD1NXA2P9vw8GMFr1jUYSQhBD+74kECbdofOfPmvM+OHjPMn+sfxLGPtIt0PGCIpKeZ+Gsg3PPR/CeLoLcqSvMY34OQ9c+dxwXSWv2Fdfsu+a2aPhO1ijBmTNth+f8Pri+fR/CWQfTIzkwDZnG4WAUFCx4ohDRGmrhKnof3gsHIMc/AOukX0Pp8oDesSiCseCJgkwIAe/pIjh3/wlQ3TA9OhOmgVkcDkZBx+8woiDSblbBWbgG6uXjULo+BMvo5yDH2pp/IFEAsOCJgkAIDZ7jf20aDiZJMD8xD8Z+T3I4GIUUC54owNRrV+AsXAWt4iyU7gNgGfUs5BhOVaXQY8FHEF7vri+heeE+8me4D2wCjGZYxi6Eoc9IzmMi3bDgI8i97rHK+6gGn1p9Ac4d70FzXISh13CYRz4DuV0nvWNRlGPBRxhe7x5awuuG+8BGuEu2QrLEwDLhJRh7DtU7FhEAFjxRi3ntp5rGDFy/CmPaaJhHzIZkbq93LKJbWPBh6JvDl3HwTPWtaZDf4157aAh3I1x718NT+jWkmC6wTn4NhpT+esciugMLPgwVH6/Apao6pCRwrz3UvOUlcO5cA1F/DcaHJ8I8/GlIRt65jNomFnyY6tmtE16dOUjvGFFDOOvg3P0hvGd2QY5NhjX3H6EkPah3LKL7YsET3YcQAt7vh4O5GmAaMg2mwVMhKUa9oxE1iwXfht3vuvbeKbwEL9i0+mtwffs+vBcOQu6SCuuU30KJ7653LCKfseDbsPtd1z5mcIpOqSKfEAKeU4Vw7fkIUL0wPzYLxgFZkGRF72hEfmHBt3H3uq49Gm+KEArajUo4d66BerkUii2taThYp656xyJqERY8EQChafAc/wtc+z4DJBnmUX8HY/pYDgejsMaCD4J77Z37i9e1h4Z67TKcO1ZBqzwHpftAWDKehdwhXu9YRK3Ggg+Ce+2d+4vXtQeXUL1wH/kS7oNfQDJaYXlyEQwPPs7hYBQxWPBBwpkwbZtaeb5ppG/NJRh6PwbzyLmQrR31jkUUUCx4iirC64Jr/0Z4jm6FZO0E68RfwpA6WO9YREHBgqeo4b1ysmk42I0KGPuOaRoOZmqndyyioGHBU8QT7ka4ij+B58T/QIpJgHXKb2Ho1k/vWERBx4KniOYtPwznznUQDddgHJAF8/CnIBk4HIyiAwv+B3h5Y+TQGm/AtftDeM/ugdy5G6wTXoSS2FvvWEQhxYL/AV7eGP6EEPCeK4Zr1wcQ7gaYhk6H6ZEcSAq/1Sn68Lv+R3h5Y/jS6q/BuXMt1PLDkBN6wTpmAZQ4zuyh6MWCp7AnhIDn5A649nwMaCrMI34K48MTIckcM0DRjQVPYU27UQln4WqoV05ASU5vGg7WMVHvWERtAguewpLQNHiOfQXXvg2ArMCcMR/GvmM4ZoDoB3z6HbasrAyzZ89GVlYWZs+ejQsXLtzz2PPnz2PQoEHIy8sLVEai26g1l9CwaQVcez6G0q0f2s96C6b0sSx3oh/xqeCXLl2KOXPm4KuvvsKcOXOwZMmSux6nqiqWLl2KzMzMgIYkAgCheuDavwENny+FuFkNy7h/gDXrl5Dbd9Y7GlGb1GzBOxwOlJaWIicnBwCQk5OD0tJS1NTU3HHsH//4R4wdOxapqakBD0rRTa08j0vvvQb3wU0w9HoU7Wa9BeODI3jWTnQfze7B2+12JCUlQVGablemKAoSExNht9sRFxd367iTJ0+iqKgI69atw8qVK1sUJj6+5defJyTEtPix3zOalIB9rlAIl5ytoXlcuLbjT7i590soHTqj66zfo12foXrHCqloeJ5/jGsOjIC8yOrxePDGG2/gn//5n2/9IGgJh6MOmib8flygbl/ncasAEBa3wouGW/Z5L5c2DQe7WQVj+pPoNmUBHDdU1Ef4un8oGp7nH+OafSfL0n1PjJsteJvNhoqKCqiqCkVRoKoqKisrYbPZbh1TVVWF8vJyLFq0CABw48YNCCFQV1eHFStW+B2aoptw1cNV/DE8JwshdUyCNWcxDMl9IZvbAYiuf/hErdFswcfHxyM9PR0FBQXIzc1FQUEB0tPTb9ueSU5ORnFx8a2P//CHP6ChoQG/+93vgpO6Fe43b4YzZPTnvXAIzqK1EI3XYRw4CeZh0zkcjKiFfLqKZtmyZcjPz0dWVhby8/OxfPlyAMDChQtx9OjRoAYMtO/nzdwNZ8joR2u8gcbtK9G47W1Ilg5oN30JLCNms9yJWkESQvi/6R0kodiDz/vgIACE/byZSNmnFELAe3Y3XLs+hPA4YRoyDaZBk+86HCxS1uwPrjk66LYHTxQsWp0DzqJ1UMuPQE7s3TQcrHM3vWMRRYyILfh77bVzn11/QmjwnPgGruJPAKHB/PgcGPtncjgYUYBFbMHfa7Y799n1pV2/2jQczH4KSrd+sGTM53AwoiCJ2IIHONu9LRGaCnfJV3Af2AAoBlhGL4AhLYPvRCUKoogueGobVEc5nDtWQau+AEPqEJifmMf5MUQhwIKnoBGqB+6DX8B9eAskS3tYMl+AoedwnrUThQgLnoJCrTjbdNZeewWGPk/A8vjPIFn44jZRKLHgKaCExwnXvs/gObYdUvvOsGa/CkOPgXrHIopKLHgKGO+lY3DuXANxsxrGfuNhfnQGJJNV71hEUYsFT60mXPVw7fkInlM7IXVKgnXq6zDY0vSORRT1WPDUKp6yA3AVrYNw3oTpkSkwDcmFZDDpHYuIwIKnFtIarsO1Kx/e8/sgx/eAddIrULqk6h2LiH6ABU9+EULAe2YXnLs/BDwumIY/DdOgSZBkfisRtTVh/6/ym8OXcfBM9a27MX2PM2cCT7tZDWfRWqgXj0JOehCWMQugxCbrHYuI7iHsC774eAUuVdUhJYEzZ4JFCA2e0q/h2vspIATMI+fC2H88JInDwYjasrAveADo2a0TXp05SO8YEUmrtTcNB7t6Gkq3/rCMng85JkHvWETkg4goeAo8oXnhLtkK94GNgMEMy9jnYejzBMcMEIURFjzdQa3+rmnMgOM7GHoOg/mJZyC3i9U7FhH5iQVPtwivu2k42JEtkCwdYMl8EcZew/WORUQtxIInAID36hm4drwH7fpVGB4aBcuIn3I4GFGYY8FHOeFuhGvfp/Ac/xpShzhYJ/8GhpSH9Y5FRAHAgo9i3otHm4aD1dXA2P9vw8GMFr1jEVGAsOCjkHDWwbnnT/Ce/hZyrA2Wab+HoWsfvWMRUYCx4KOM5/w+uL59H8JZB9PgqTANnsrhYEQRigUfJbSGWriK3of3wgHI8Q/AOunXULo8oHcsIgoiFnyEE0LAe7oIzt1/AlQ3TI/OhGlgNiRZ0TsaEQUZCz6CaTer4CxcA/XycShdH4Jl9HOQY216xyKiEGHBRyChafCU/rVpOJgkwfzEPBj7PcnhYERRhgUfYdRrV+AsXAWt4iyU7gNgyZgPuUO83rGISAcs+AghNC/ch7fAffALwGiGZexCGPqM5HAwoijmU8GXlZVh8eLFqK2tRWxsLPLy8pCamnrbMe+++y62bNkCRVFgMBjwyiuvICMjIxiZ6UfU6gtw7ngPmuMiDL0ehXnkXMjtOukdi4h05lPBL126FHPmzEFubi42bdqEJUuWYN26dbcdM3DgQCxYsABWqxUnT57EM888g6KiIlgsfGdksAivG+4DG+Eu2QrJEgPLxJdgTB2qdywiaiOafdXN4XCgtLQUOTk5AICcnByUlpaipqbmtuMyMjJgtVoBAGlpaRBCoLa2NvCJCQDQWH4c9Z+9AfeRLTA+NArtZ73Fciei2zR7Bm+325GUlARFabpuWlEUJCYmwm63Iy4u7q6P2bhxI3r06IGuXbsGNi01DQfbux43S7+GFNMF1smvwZDSX+9YRNQGBfxF1r179+Ltt9/GqlWr/H5sfLz/42mNpqYfPAkJMX4/Ntw0nD2Iqj//P6g3HOj4aA7ixvwMsim6tsCi4Xn+Ma45OgRjzc0WvM1mQ0VFBVRVhaIoUFUVlZWVsNnufMPMoUOH8Nprr2HlypXo1auX32EcjjpomvDrMR63CqNJQVXVTb+/XrgQzjo4d38I75ldkDsno13uP6LLw4P/tmaP3vFCJiEhJqKf57vhmqNDS9csy9J9T4ybLfj4+Hikp6ejoKAAubm5KCgoQHp6+h3bMyUlJXjllVfwzjvvoH9/bhkEghAC3u+Hg7kaYBoyrWk4mGLUOxoRhQGftmiWLVuGxYsXY+XKlejYsSPy8vIAAAsXLsTLL7+MAQMGYPny5XA6nViyZMmtx/3rv/4r0tLSgpM8wmn11+AqWgfvd4cgd0mFdcpvocR31zsWEYURnwq+d+/eWL9+/R1//t///d+3/vdnn30WuFRRTAgBz6lCuPZ8BKhemB+bBeOALA4HIyK/8Z2sbYh2oxLOwtVQr5yAYktrGg7WiVciEVHLsODbAKFp8Bz/C1z7PgMkGeZRz8KYPobDwYioVVjwOlNrLsNZ+B60yvNQegyCZdSzkDvc/f0FRET+YMHrRKheuA9/CfehLyAZrbCM+3sYeo/gcDAiChgWvA7UyvNNI31rLsHQ+7Gm4WDWjnrHIqIIw4IPIeF1wbV/AzxHv4Jk7QTrxF/CkDpY71hEFKFY8CHivXICzsI1EDcqYOw7BuYRsyGZ2ukdi4giGAs+yIS7Aa7iT+A58Q2kmARYp/wWhm799I5FRFGABR9E3vLDcO5cC9FQC+PAbJiH/QSSwax3LCKKEiz4INAab8C1+0N4z+6B3DkF1gkvQUn0f/gaEVFrsOADSAgB77liuHZ9AOFugGnodJgeyYGk8K+ZiEKPzRMgWl0NnEVroZYfgZzQC9YxC6DEpegdi4iiGAu+lYTQ4DlZCNeejwFNhXnET2F8eCIkmWMGiEhfLPhW0K5XNA0Hs5+EkpzeNBysY6LesYiIALDgW0RoGjzHvoJr3wZAVmAe/RyMaaM5ZoCI2hQWvJ/Umotw7lgFraoMSo9HYMl4FnL7znrHIiK6AwveR0L1wH2oAO5DBZDM7WAZ/wsYej3Ks3YiarNY8D5QK881nbVfuwzDg4/DPHIOZEv03fWdiMILC/4+hMcF1/7P4Tm6DVL7zrBm/wqGHo/oHYuIyCcs+HvwXi6Fs3A1xM0qGNOfhPmxWZBMVr1jERH5jAX/I8JVD1fxx/CcLITUMQnWnMUwJPfVOxYRkd9Y8D/gvXAIzqK1EI3XYRo0Gaah0yEZTHrHIiJqERY8/jYc7Nt8eM/vhRyXAmvWL6Ek9NQ7FhFRq0R1wQsh4D27G85dHwAeF0zDnoJp0GQOByOiiBC1TabVOeDcuRbqxRLIib1hGbMASuduesciIgqYqCt4ITR4TnwDV/EngNBgfnwOjP0zORyMiCJOVBW8dv3q34aDnYLSrT8sGfMhd0zQOxYRUVBERcELTYW75Cu4D2wAFAMsoxfAkJbBMQNEFNEivuBVRzmcO96DVv0dDKlDYH5iHoeDEVFUiNiCF6oH7oNfwH14CyRLe1gyX4Ch53CetRNR1IjIglcrzjYNB6u9AkOfJ2B5/GeQLB30jkVEFFI+FXxZWRkWL16M2tpaxMbGIi8vD6mpqbcdo6oq3nzzTezcuROSJGHRokWYOXNmMDLfk/A44dr3GTzHtkPqEAfrpFdh6D4wpBmIiNoKnwp+6dKlmDNnDnJzc7Fp0yYsWbIE69atu+2YzZs3o7y8HNu2bUNtbS2mT5+Oxx9/HCkpobnxtPfSMTh3roG4WQ1jv/EwPzqDw8GIKKo1e/G3w+FAaWkpcnJyAAA5OTkoLS1FTU3Nbcdt2bIFM2fOhCzLiIuLQ2ZmJrZu3Rqc1D9gEB6Mrt+Gxi3/BsgGWKe+DsuoeSx3Iop6zZ7B2+12JCUlQVEUAICiKEhMTITdbkdcXNxtxyUnJ9/62Gaz4erVq36FiY/3f598+kMuJBw5gdiRP0FsxizIUTQcLCEh+m46wjVHB645MNrUi6wORx00Tfj1mF7DR6FL5jg4bnjhuOYC4ApOuDYmISEGVVU39Y4RUlxzdOCafSfL0n1PjJvdorHZbKioqICqqgCaXkytrKyEzWa747grV67c+thut6Nr165+B/aXJMmQzdyOISL6sWYLPj4+Hunp6SgoKAAAFBQUID09/bbtGQDIzs7G+vXroWkaampqsH37dmRlZQUnNRERNcunCVvLli1Dfn4+srKykJ+fj+XLlwMAFi5ciKNHjwIAcnNzkZKSgokTJ2LWrFl48cUX0b179+AlJyKi+5KEEP5tegdRS/bgAe7ZRQuuOTpwzb5r9R48ERGFJxY8EVGEYsETEUWoNnUdvCy3fNJjax4brrjm6MA1R4eWrLm5x7SpF1mJiChwuEVDRBShWPBERBGKBU9EFKFY8EREEYoFT0QUoVjwREQRigVPRBShWPBERBGKBU9EFKHCpuDLysowe/ZsZGVlYfbs2bhw4cIdx6iqiuXLlyMzMxMTJkzA+vXrQx80gHxZ87vvvospU6Zg2rRpeOqpp7Bz587QBw0gX9b8vfPnz2PQoEHIy8sLXcAg8HXNW7ZswdSpU5GTk4OpU6eiuro6tEEDyJc1OxwOLFq0CFOnTkV2djaWLVsGr9cb+rABkJeXh3HjxiEtLQ2nT5++6zFB6S8RJubNmyc2btwohBBi48aNYt68eXccs2HDBrFgwQKhqqpwOBwiIyNDXLx4MdRRA8aXNRcWFoqGhgYhhBAnTpwQQ4cOFY2NjSHNGUi+rFkIIbxer3jmmWfEq6++Kv7lX/4llBEDzpc1l5SUiEmTJonKykohhBA3btwQTqczpDkDyZc1v/nmm7eeW7fbLWbMmCG+/PLLkOYMlH379okrV66IJ598Upw6dequxwSjv8LiDN7hcKC0tBQ5OTkAgJycHJSWlqKmpua247Zs2YKZM2dClmXExcUhMzMTW7du1SNyq/m65oyMDFitTfekTUtLgxACtbW1oY4bEL6uGQD++Mc/YuzYsUhNTQ1xysDydc1r1qzBggULkJCQAACIiYmB2WwOed5A8HXNkiShvr4emqbB7XbD4/EgKSlJj8itNmzYsDvuY/1jweivsCh4u92OpKQkKIoCAFAUBYmJibDb7Xccl5ycfOtjm82Gq1evhjRroPi65h/auHEjevToEZKbnQeDr2s+efIkioqKMH/+fB1SBpavaz537hwuXryIuXPn4ic/+QlWrlwJEaZzAn1d8wsvvICysjKMGjXq1n9Dhw7VI3JIBKO/wqLgqXl79+7F22+/jX//93/XO0pQeTwevPHGG1i+fPmtgogGqqri1KlTWL16Nd5//30UFhZi06ZNescKqq1btyItLQ1FRUUoLCzE/v37w/Y3cr2ERcHbbDZUVFRAVVUATd/slZWVd/zKY7PZcOXKlVsf2+32sD2b9XXNAHDo0CG89tprePfdd9GrV69QRw0YX9ZcVVWF8vJyLFq0COPGjcPatWvxySef4I033tArdqv4+jwnJycjOzsbJpMJHTp0wPjx41FSUqJH5Fbzdc35+fmYNm0aZFlGTEwMxo0bh+LiYj0ih0Qw+issCj4+Ph7p6ekoKCgAABQUFCA9PR1xcXG3HZednY3169dD0zTU1NRg+/btyMrK0iNyq/m65pKSErzyyit455130L9/fz2iBowva05OTkZxcTG+/vprfP3113j22Wcxa9YsrFixQq/YreLr85yTk4OioiIIIeDxeLBnzx707dtXj8it5uuaU1JSUFhYCABwu93YvXs3+vTpE/K8oRKU/mrVS7QhdPbsWTFjxgwxceJEMWPGDHHu3DkhhBDPP/+8KCkpEUI0XVmxZMkSMX78eDF+/Hjx0Ucf6Rm51XxZ81NPPSUee+wxMW3atFv/nTx5Us/YreLLmn/onXfeCfuraHxZs6qq4q233hLZ2dli8uTJ4q233hKqquoZu1V8WfN3330n5s+fL3JycsSkSZPEsmXLhMfj0TN2i61YsUJkZGSI9PR0MXLkSDF58mQhRPD7i3d0IiKKUGGxRUNERP5jwRMRRSgWPBFRhGLBExFFKBY8EVGEYsETEUUoFjwRUYRiwRMRRaj/Dz1xCICtlAIlAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "x = dat['FLAIR']; y = dat['GOLD_Lesions']\n",
    "c = np.concatenate( [[ np.min(x) - 1], np.sort(np.unique(x)) \n",
    "                     , [np.max(x) + 1]])\n",
    "tpr = [np.mean( (x1 >= citer) ) for citer in c]\n",
    "fpr = [np.mean( (x0 >= citer) ) for citer in c]\n",
    "plt.plot(fpr, tpr);\n",
    "plt.plot([0,1], [0,1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74e25a3d-8df7-467e-a124-1e023405f07e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### Binormal estimation\n",
    "\n",
    "+ We could also assume distributional forms for $T$ and $F$ \n",
    "+ For example, suppose $X ~|~ Y=y \\sim N(\\mu_y, \\sigma_y^2)$. \n",
    "+ Then, note if $\\Phi$ is the standard normal distribution function then \n",
    "$$T(c) = 1 - \\Phi\\{ (c - \\mu_1) / \\sigma_1 \\}$$ \n",
    "$$F(c) = 1 - \\Phi\\{ (c - \\mu_0) / \\sigma_0 \\}$$ \n",
    "$$F^{-1}(f) = \\mu_0 + \\sigma_0 \\Phi^{-1}(1-f)$$ \n",
    "+ Thus, the ROC curve is \n",
    "$$\n",
    "T\\{F^{-1}(f)\\} = 1 - \\Phi\\left\\{  \\frac{\\mu_0 -\\mu_1}{\\sigma_1} + \\frac{\\sigma_0}{{\\sigma_1}} \\Phi^{-1}(1-f) \\right\\}\n",
    "$$\n",
    "where $\\mu_y$ and $\\sigma_y$ can be estimated from the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0aa86fc7-bd7b-45fc-9819-d344bb47bc32",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD7CAYAAABgzo9kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA17ElEQVR4nO3deXjV5Z3//+fZz4FsJGRfCGGJYV9EUEGUHQzGuraoVZnitParHZ1xhk6noJe9nGGu3/zRdurM2KsiiFalVmmRIqVWAgoIshN2shOy78nZPp/798fBVGRL4Jyc7f24Lq6roZ/kvG9PeOXO/bk/79uglFIIIYSIOMZgFyCEECIwJOCFECJCScALIUSEkoAXQogIJQEvhBARSgJeCCEilAS8EEJEKHOwC/i65uZOdL3v2/KTkmJobOwIQEWhS8YcHWTM0eF6x2w0Ghg0aOAV//+QCnhdV9cV8F99brSRMUcHGXN0CMSYZYlGCCEilAS8EEJEKAl4IYSIUNcM+FWrVjFr1izy8/M5efLkZa/RNI2XXnqJOXPmMHfuXNavX+/3QoUQQvTNNQN+9uzZvPXWW2RmZl7xmj/+8Y9UVFSwZcsW3n33XX75y19SVVXl10KFEEL0zTUD/uabbyY9Pf2q12zatIkHH3wQo9FIYmIic+bMYfPmzX4rUgghRN/5ZZtkTU0NGRkZPR+np6dz/vx5f3xpIYSICJ8eqGb30dq//YVSjPIeIZe9mMd/l4KJU/z+miG1Dz4pKea6Pzc5OdaPlYQHGXN0kDFHhn2nGqiq72BoZjxxWjMzurYSq6r599xk5pq6uCMAY/ZLwKenp3Pu3DnGjRsHXDqj763Gxo7r2uyfnBxLfX17nz8vnMmYo4OMOXJ43BrZgwfw/0aew7X392Ayc27SPegNO5lQMOS6xmw0Gq46MfbLNskFCxawfv16dF2nqamJrVu3Mn/+fH98aSGEiAhJegP3O9/DtftdzFljGPjgK5wc6MCAgWGDcgLymtecwf/sZz9jy5YtNDQ08OSTT5KQkMBHH33EsmXLePbZZxk7dixFRUUcPHiQefPmAfDDH/6Q7OzsgBQshIhcnx6oZt+pBjxuLdil+I1ReZns2csD7j24DHbsc57GPHQKjc5mPq36jHGDR5HgiKe+w/+/tRhC6dBtWaLpPRlzdIi2Ma96ax9V9R1kJV///bhQkqrVcJfrLySqJk6Yb6J77P3MmDKCFlcrv9z/a9rc7fzLlGcpyMkNyBJNSN1kFUKIoZnxPP/g+GCXcUOUx4Vrz/t4jvwZw8BB2Gc8z805vnuUZ1rKWH30bbq8Xfxg3JMMdiQFrA4JeCGE8CNvdQnO4tWo9noso2Zhu+VBDFYHXZ5u/lS2lb9W7iDJPojnJv2A7NgrP0DqDxLwQoiAuWTv9zVU1HUwLCs+gBUFjnJ14tr1Lp4TxRjiU3Es/jHm9HzcmpvPKnewuewvdHq6uC3jFu4bfjd2sz3gNUnACyECZvfRWirqOshJ6d2aek5KDDMnZgW4Kv/zlH2Ja8ebqO42rBPuxjqpiE7dw/bSrXxa9Rkdnk5GJORx/4jFAZ+1f50EvBAioHJSYviXRyb1+vpwurGsd7Xi+nwd3rN7MCZlY5/3I0otOjtO/I4DdYfxKo0xSTcxd8hdDE8Y2u/1ScALIUQfKaXwnvoc5863weOiZdJCDg6K58sz66nrbsBhtnN75jSmZ0wlIyYtaHVKwAshLtLXdfOr6cvyTLjQOxpxbl9DXc1RjqZncDAhjuq2LzG0GRiRkMf83FlMShmH1WQNdqkS8EKIi/V13fxqclJimDo61Q9VBZ+mezlzcAMHy4opcZipy00CXAyx2Lk/azGTUsaRYAutG8QS8EKIS/R13TxSdXq6ONF8mqM1+zlSX0KHUWGMszE8Poc7UicwdvAoBjsSg13mFUnACyHEBV7dS2lrOceaTnG86RQV7VUoFHZNJ9+pMS5rKmNHL2agdWCwS+0VCXghotSV1tojcd38SnSlU9V+jtMtZznRfJqTLWdxa26MBiO5jlTmdJsY3lhHbto4Bsz+LsYBCcEuuU8k4IWIUldaa4+kdfNv0nSNivZqTrec5VTLWc60lOHUnACkOAYzLW0yN8XnMaT8GKZDWzDYY7BN/3ssef4/jKM/SMALEcUifa3do3kob6/yBXrzWc62lePW3ACkDkhhcup4RiTkMWJQHgm2eLTzp3AWv47eUoN55O3Yp30Hgz18f5uRgBdCRIwWVytnW8spvfCnor0aTflaD2cMTGNa2s2MGJTH8IShxFn/doKS8jhxfrYOz9G/YIhJxLHwHzFnjw3WMPxGAl4IEZY0XaO6o4azreWcbS2jtK2CJmczABajmZzYLGZlz2BofA7D4ocSc4Ubo96qI77mYB1NWEbPxjblfgxWR38OJWAk4IUQYaHV1UZZWyVlbRWUtpZT3laJW/cAkGCLJy9+CHdlTycvfghZMRmYjVePN+XswLnrHbwnd2CMT8N+z79iThvRH0PpNxLwQoiQ0+11UtFWRXlbJeXtlZS1VdLiagXAaDCSHZPJ7RlTGRqfQ158LoPsCX36+p6ze3B99ibK2YF1QiHWSfdgMAf/yVN/k4AXIkJcT2veUNgO6dG9nOuoobzNF+RVXdWca6tF4TvdLdmRxPCEoeTG5TAkLpusmAysJst1vZbe1YLrs3V4S/diTBqCY+E/Yho8xJ/DCSkS8EJEiOtpzdvf2yF1pVPf1UDZ12bm1e3n8F64ERpriWFk8lAmJo1jSFw2OXFZxFhu/KEipRTekztw7noHvC6stzyAddwCDNdYxgl3kT06IaJMKG171JVOQ3cjFe3VVLRXUdFWRWX7uZ595zaTlZzYLO7Mns6QuGxy47IZZEsgJSXOr+2C9fZ6nNvXoFUdwZQ2EvsdT2JMSPfb1w9lEvBCiBumlKK+u5HK9irK26uobKumsqOabq8vzM1GM5kD07k5bQK5sdkMicsmbWAKRoMxgDXpeI7+BdcXvwODAdvtj2EZdReGAL5mqJGAFyLMBLvFwNfD3Dc7r6ayvepvYW4wkRmTweTUCeTEZpITm03GwFRMRlPAa/uK1nIO17bVaLWnMGWPxT79cYyxg/vt9UOFBLwQYaY/WwwopWjobvItsVwI9Mr2arq93YAvzDNi0pmcMp6cuCxyYrNIH5h6zS2KgaJ0L+6Df8L95Qaw2LDfuQzziNswGAxBqSfYJOCFCEOBWGtXStHobKK8rYrKr9bNvxbmJoOJzJg0JqeMIyc2i+y4TDIGpgUtzL9JayjDue119MYKzHlTsN32KMYBodWfvb+FxjsjhOhXvjBv7rn5WdHuC/Wub4T5pJRxvmWWuKyQCvOvU1437n0bcB/8EwZ7LPa5z2AZOjnYZYWE0Hu3hIhQ13MUnsVqwuPWLvq7vq61K6Vocjb7bn62V1/YzVJNp7cL8IV5xsBUJqaMJTs2iyGxWaTHpGEJwTD/Ju/5kzi3vY5qPY8lfwa2ad/GYAuPXu39IfTfQSEihL+OwrvaWrsvzFv+tpvlwlJLp8cX5kaDkcyBaYxPHnNhzTyTjJj0sAjzr1Publxf/A5PyV8wxA7GsegFzFmjg11WyAmvd1WIMNfXtfPk5Ngr7glXStHsavHtZPnaMkuHpxPwhXn6wFTGDx7tm5lfWGaxXOdToKHCW3EI5/Y3UJ3NWMbM8zUHs9iCXVZIkoAXIgwopWhxtfbc+Pwq0L8Z5mMHj+pZM88cmB72Yf51ytmBc+fbeE99jjEhA0fRTzClDg92WSFNAl6Iq7iedfMr6e3yjFKKVncb5W1VNJ6v51jtGSrbqmn3dAC+ME8bkMKYwQXkxPq2JmbGpF93f5ZQp5TCW7oH12frUM5OrJPuwTpxMYYIHa8/ScALcRX+WjeHK6+dt7haL8zIq3vWztvdvjA3GAykD0hlVFJ+zz7zrJh0rKbI63x4OXpXC64da/GW7cM4OBfHon/ClJQT7LLCRq8CvrS0lOXLl9PS0kJCQgKrVq0iNzf3omsaGxv58Y9/TE1NDR6Ph2nTpvFv//ZvmM3yM0SEN3/uOW91tXG4oaQn0Cvaq2hz+9bYDRhIG5jCqMR838w8LpPxuSNpb3b75bXDiVIK74ntOHf9FjQvtqkPYRk7H0M/Pg0bCXqVvitXrmTJkiUUFRWxYcMGVqxYwdq1ay+65n//938ZNmwYr732Gh6PhyVLlrBlyxYWLVoUkMKFCHWdni5f18S2igv7zatpdbcBvjBPHZjCTYkjepZZsmIzsH1jZm4322gnugLe03ye7k2volUfxZSe72sOFp8W7LLC0jUDvrGxkZKSElavXg1AYWEhL7/8Mk1NTSQmJvZcZzAY6OzsRNd13G43Ho+H1NTIPJld9L9PD1Sz71TDJXvCA623yzMezUNVx7meE4fK2yqp724ELoT5gGRGDhpOTlzmhWWWDOxm2fnxdUrX8Rz9M1V7f4/CgG36d7EU3BlVzcH87ZoBX1NTQ2pqKiaT71cjk8lESkoKNTU1FwX8008/zTPPPMP06dPp7u7mkUceYfLkvj1NlpR0/eucycmx174owkTTmPedaqC0upWhmf376PmwrHhmTsy66L+1rnTOtddyurGM001lnG4so7ylCk3pACQ6EhiemMvcpBkMT8wlLzGHAZbrP+MzGt5nd30l9X96FVf1SRzDJpG86O8xx0VXc7BAvM9+WyDfvHkz+fn5rFmzhs7OTpYtW8bmzZtZsGBBr79GY2MHuq76/NpX2yscqaJtzB63xtDMeJ5/cHy/v3arq42tJbt6ZublbVU9Pc3tJhs5cdnMzpnZ09M8wXbxD6HOFi+dXN97Fenvs9K8uA9+hHvfHzFY7Njveoq0W+fR0NABETzub7re99loNFx1YnzNgE9PT6e2thZN0zCZTGiaRl1dHenpFzfMX7duHa+88gpGo5HY2FhmzZrF7t27+xTwQgSb0+uior2q5/i4sraKi84CzYxJZ0raxJ4wTx2QHNCe5pFMqy/1NQdrqsQ8bCq22x7B6IiL2s6PgXDNgE9KSqKgoICNGzdSVFTExo0bKSgouGh5BiArK4vi4mLGjRuH2+1m586dzJ07N2CFi8h0tV7nw7L8uzzzVcOts61llLZWUNpaRlVHTc9ZoIMvnAXqC/OcGzoLVPyN8rpx7f0Az+HNGBzxOOb9CHPuxGCXFZEMSqlrromcOXOG5cuX09bWRlxcHKtWrSIvL49ly5bx7LPPMnbsWCoqKli5ciUNDQ1omsbUqVP5yU9+0qdtkrJE03uROuZVb+274o3NOVOHMHl40nV/bY/mobKjmrOt5ZS2lnO2tbxni6LNZCU3Loe8+CHkxuWQG5dDjDX4Tasi7X32njuOs3g1qq0Wy00zsU196JLmYJE25t4I1BJNrwK+v0jA916kjnnVW/sALrvvvK9jbnW1cba1vGeGXtle1XO482BHEkPjhpAX7/uTEZMWkkstkfI+K3c3rt3v4Tn2VwyxydjveBJz5qjLXhspY+6LoK3BCxEOdKVT3VHDmdYySi/M0BudzYDvPNAhFw53zosfwtD4IcRZI39nSqjwVhzAuX0tqqsZy9j52Kbch0G2iPYLCXgRUH3t5dLbfede3UtFexWnm0s51XqWsy3lPTtb4q1x5CXkcmfW7QyNzyU7NiMkD6qIdLqzHdfnb+M9vRPjoEwcc3+IKWVYsMuKKvJdLwKqr71crtSvxaW5OXT+GF+WH+V0SyllbRV4dC8AaQNSuDl1PMMT8siLzyXRniA7MYJIKYX3zG5cn7+FcndhnVR0oTmYxE1/k//iIuCup5dLp6eLs61lnGo5y+mWUirbq9GVjgED2bEZTM+cxvCEPIbF5xJrvfFGYMI/9M5mX3Ow8v0Yk4fimLkUU2J2sMuKWhLwIiR0ebo51XKGE81nONV8hnOd5wEwG0wMictmbs6dTB4yikRScJjtQa5WfJNSCs/xbbh2vQu6hm3at7GMmYfBGHo3rqOJBLwICpfm5kxLKSebz3Ci+TSV7dUoFBajhWHxuUxKGc/whKHkxmX3HFoRjbsrwoHeVoezeDXauWOYMgp8zcHiUoJdlkACXvQTj+6lrLWCk82nOdF8mrK2SjSlYTKYyI3LYWHubEYOGk5ufE7YnQ8arZSu4zmyBdee34PRhG3GE1humin3P0KI/EsSAaGUorqjhraBx3DazvNC8ft4dM+FNfRMZmXPIH/QcPISci9pkStCn9ZU5WszUH8WU84E7DMexzhwULDLEt8gAS/8psPdyfGmk5Q0neRY00nfU6JxYPHEc3vGLYwcNJwRCUMZYBkQ7FLFdVKaF/eBjbj3/xGDdQD2Wd/HPGyqzNpDlAS8uG6arlHaVsGxppOUNJ7oWUcfaB7ATYkjKEgcySfbXJj1ATw43z8nIong0erO+mbtzVWYh0/zNQezywNjoUwCXvRJs7OFo43HOdZ0kuNNp3FqTgwYGBqfw91D51KQNJKc2Kyex/6L9X1BrljcKOV14drzezxHtmAYkIBj/j9gHjIh2GWJXpCAF1elK53K9moONxzjSEMJlR3nABhkS2By6jgKEvPJHzT8hg60EKHLe+4Yzm2vo9rrsRTc5WsOZpX3OlxIwItLuDU3J5pPc7ihhCMNx2h1t2PAQF78EO4dtogxgwtIG5Ai664RTLm7cO16F8/xbRjiUnEULseccVOwyxJ9JAEvAGh1tXO44SiHG45xovkUHt2L3WSjIHEkYwePYnTSTSHRPlcEnrd8P87ta1DdrVjGLcR2873SHCxMScBHsWZnCwfqj7C/7jBnW8tQKJLsg7gtYypjBxcwIiFPmnRFEb27Ddfnb+E9sxtjYhaO+T/ClDw02GWJGyD/eqNMQ3cTB+oPs7/uMGVtFQBkDExj4dA5TEgeQ8bANFl6iTK+5mC7cH32FsrTjfXmb2Edf7c0B4sA8g5GgYbuRr6sPcj++sNUtlcDkB2byeK8BUxMHkPqQHmsPFrpHY04d6xFqziIMSUPxx1/hykxM9hlCT+RgI9Q7e4Ovqw7yN7zByhtKwcgNy6He4ctYmLKWAY7rv/ou8u52lmqvW0VLPqPUjqeY5/i2v0eKB3brd/BMnquNAeLMBLwEcTpdXGo4Sh7avdzvOkUutLJGJhGUd5CJqdOIMkRuEfJr9T3/Ur93UXw6K3nfc3Bak5gyhyFfcYT0hwsQknAhzld6RxvOsXu819yqP4obt3DIFsCs7PvYEraRDJj0vutluvp+y76j9I1PIc/xrX3AzCZsd+xFHP+DLnnEsEk4MNUXWcjH539lF01e2l2tTDA7OCWtElMSZtEXvyQkDxAWgSP1liJs/h19PpSzEMmYpv+XWkOFgUk4MOIR/NwsOEoO8/t4UTzaQBuShzBt4bfzbjk0f3SZlfW2sOL0jy49/8R9/6PMNgHYp/zNOahU2TWHiUk4MNAbWcdxdU7+eL8Prq83STaB/HA6EWMjRsX0HX1y5G19vCh1Z72NQdrOYd5xG3Yb12CwS4/hKOJBHyI0pXO0cbjfFr5GcebT2E2mBifPIZbM6aQP2g4qSnxQTvdSNbaQ5vyuHDt+R2eI1sxDByEY8HzmHPGBbssEQQS8CGm09PFzpo9FFftpNHZRIItnsV587k9Y6ocLi2uyVt1FOf21aj2BiyjZmO75QFpDhbFJOBDRH1XI1srt7G75ks8uofhCUO5d/gixg8ejcloCshrXmk9/WpkrT00KVcnrl3v4DmxHUN8Ko7FP8acnh/sskSQScAHWUVbFX+u+JT9dYcxGYzckjaJO7On98v2xiutp1+NrLWHHk/Zl7h2vInqbsM64W6sk4owmOUYRCEBHxRKKY43n+LP5Z9yovk0dpOdOTkzuSt7OvG2uH6tRdbTw5fe1Yrr83V4z+7BmJSNY/4/YErODXZZIoRIwPejr4L9o7NbKG2rIN4ay73DFjE9cxoOsz3Y5YkwoZTCe+pznDvfBo8L65T7sY5fiEE6f4pvkO+IfnKy+Qwbz27hTGspg2wJfCf/Pqam39wve9dF5NA7GnFufwOt8jDG1OHY71iKaVBGsMsSIapX6VJaWsry5ctpaWkhISGBVatWkZube8l1mzZt4n/+539QSmEwGFi9ejWDBw/2d81hpbS1nD+c/ZiTzaeJt8bx8Mh7uTXjFgl20SdK6XhKPsH1xe9AKWy3PYJl9GwM8sSyuIpepczKlStZsmQJRUVFbNiwgRUrVrB27dqLrjl8+DD//d//zZo1a0hOTqa9vR2rNXpv9DR2N/HhmU3sqztErCWG+0csZnrGNKwmS7BLE2FGb6nxNQc7fxJT5mjsdzyBMTY52GWJMHDNgG9sbKSkpITVq1cDUFhYyMsvv0xTUxOJiYk9173xxhssXbqU5GTfN15sbGyASg5t3V4nW8r/yieV2zFgYFHuHGbnzMQuR56JPlK6Rsvnv6ez+F0wWbHP/DvMI6dLmwHRa9cM+JqaGlJTUzGZfHuxTSYTKSkp1NTUXBTwZ86cISsri0ceeYSuri7mzp3LD37wg6j5ZtSVzq6avfzhzGbaPR1MTZvM4rz5DLInBLs0EYa0hnKcxa/T0VCOOXcytumPYRyQEOyyRJjx20KwpmmcOHGC1atX43a7+d73vkdGRgb33ntvr79GUtL1P0CTnBy83xgqW8/x2t63OdFwhpsGD+Pxif+PYYlDAv66Nzpmi9Xkl6/Tn8Kp1uuhe920bF9P+84PMQ2II+X+fyLmpluDXVa/i/T3+XICMeZrBnx6ejq1tbVomobJZELTNOrq6khPv/hBnIyMDBYsWIDVasVqtTJ79mwOHTrUp4BvbOxA11WfB5GcHBuUvixuzcOfyraytWIbDrOdRwseYlraZAyaIeD1+GPMHrcGELSeNn0VrPe5v3jPn8K17Tforecxj5yOfdq3iclOj+gxX06kv8+Xc71jNhoNV50YXzPgk5KSKCgoYOPGjRQVFbFx40YKCgouWp4B39r8tm3bKCoqwuv1smvXLubPn9/ngsNFaWs5a0vepa67gWlpN/Ot4XcTYx0Y7LJEGFIeJ64vfofn6F8wxCTiWPRPmLPGBLssEQF6tUTz4osvsnz5cl599VXi4uJYtWoVAMuWLePZZ59l7Nix3H333Rw5coRFixZhNBqZPn06DzzwQECLDwav7uVPpVv5uPyvJNjieXbCU+QnDg92WVclPdxDl7fqCM7i1aiOJiyjLzQHs8hDb8I/DEqpvq+JBEioL9Gc76zjjaNvU9lxjmlpN/PAyHuC9gRqX8a86q19VwzzqaNTuXNCpr/LC4hI+tVdOTtw7noH78kdGOPTsM38O8xpIy65LpLG3Fsy5t674SUa4bP3/H7eOvE+VqOFp8Z+l/HJ4fUrtPScCR2es3twffYmytmBdUIh1kn3SHMwERAS8Nfg0b28f+qPbK/eSV58Ln835hESbPHBLkuEIb2rBdeON/GWfYkxaQiOhf+IaXDgd1uJ6CUBfxWtrnb+7/AblLdVMjvnDoryFgasN7uIXEopvCd34Nz5W9DcWG95EOu4+dIcTAScfIddQXVHDf9zcDWdnk6WjXmMCSljg12SCEN6ez3O4jfQqo9iShuJ/Y4nMSYEvte/ECABf1lHG4/zmyPrcJgdPD/5abJjw+MmpAgdSul4jv7F1xzMYMB2+2NYRt0lzcFEv5KA/4a9tQdYU/IOmQPT+P74J2W9XfSZ1nwOZ/Hr6LWnMWWPxT79cYyx0d1VVQSHBPzXfHZuN789/nuGJeTy/XFPht0hHLLfPbiU7sV98E+4v9wAFhv2O5dhHnFb1PRjEqFHAv6CbVWf897JDxmVlM+yMY9hNYXftrUrnbEq56gGntZQhnPbb9AbKzHnTcF226MYB8hvfyK4JOCB3TVf8t7JDxk7eBTfG/Mo5jDe3SD73fuX8rpxf/kh7kObMdhjsc99BsvQycEuSwhAAp4D9Ud489h75A8azt+NfiSsw130L2/NCV+bgdbzWPLvwDbtYQw26UckQkdUp1lpawWrj75Nblw2T419HEuYnLb06YFq9p1q6OkG+RVZa+8fyt2N64v1eEo+wRA7GMeiFzBnjQ52WUJcImoDvtnZwmuH1xBvjeP7454MqxOXdh+tpaq+g6xkWWvvb96KQzi3v4HqbMYyZh62KfdjsITP946ILlEZ8G7Nzf8dXoNbc/PMhGVh2eZ3aGY8zz84PthlRA3l7MC58228pz7HmJCBo+gnmFJDu4uoEFEZ8L879Qeq2s/x9+MeJyMmLdjliBCmlML7VXMwVxfWSfdgnbgYQ5gs54noFnUBv6/uEJ+d+4J5Q+5i7OBRwS7nqq62r31YlmzBCzS9sxnXZ2/iLduHcXAujrv/GVNSdrDLEqLXoirgm5zNvH38fYbEZVM4dF6wy7mmq+1rnzkxK0hVRT6lFJ4Txbh2vQOaF9vUh7CMnY9BGs2JMBM1Aa+U4p0TH6ApjSdHLQmbrpBX2tcejYci9Ae9rQ7n9jfQqkswpef7moPFyzKeCE9RE/D76w9ztPE4949YTPKApGCXI0KM0nU8R/+Ma8/7YDBim/5dLAV3SnMwEdaiIuC7vd2sP7mB7NhMZmbeFvDXu9LaeV/Jvvb+oTVX49z2OnrdGUzZ47DPeBxjjEwCRPiLioD/c/k22tztfH/cE/2yNHOltfO+kn3tgaU0L+6DH+He9wcMFgf2u57CPPxWaQ4mIkbEB3yLq5VPKrdzc+oEhsT13w4I6QkT2rS6s76Wvk1VmIdNxXbbIxgdccEuSwi/iviA31S6FV3pLM6bH+xSRAhQXheuvR/iObwZgyMex7wfYc6dGOyyhAiIiA74Flcru2r2clvGLQx2yJpqtPOeO+5rDtZWi+Wmmb7mYNYBwS5LiICJ6ID/tPIzdKUzJ+eOYJcigki5u3Htfg/Psb9iiE3Gcfc/Y84M7YfchPCHiA34bq+T7dW7mJgyVmbvUcxbcQDn9rWormYsY+djm3IfhjBqLCfEjYjYgN9zfj9OzcnsPszeZXtj5NC723DtfBvv6V0YB2XimPtDTCnDgl2WEP0qYgN+Z80eMmPSGRLb+50zsr0x/Cml8J7Zjevzt1DuLqyT78U6oRCDKWK/1YW4ooj8rq/uqKGivYoHRtzT5z3Nsr0xfOmdzTi3r0GrOIAxOQ/HzKWYEqVnj4heERnwe87vx2gwMiVVtr9FA6UUnuPbcO16F3QN27RvYxkzD4NR2gyI6BZxAa+U4kD9YfIHDQ/LgzxE3+htdTiLV6OdO4Ypo8DXHCwuJdhlCRESIi7gazprqe9u7NPNVRF+lK7jOfIxrj0fgNGEbcYTWG6aKW0GhPiaXv0OW1paysMPP8z8+fN5+OGHKSsru+K1Z8+eZfz48axatcpfNfbJwfqjGDAwbrAcghyptKYquja8jGvXu5gyRzHwoVewFtwp4S7EN/Qq4FeuXMmSJUv4+OOPWbJkCStWrLjsdZqmsXLlSubMmePXIvviePNJsmMziLdJX5FIozQPrr0f0PX7laj2Buyzvo9j/o8wDhwU7NKECEnXDPjGxkZKSkooLCwEoLCwkJKSEpqami659rXXXuPOO+8kNzfX74X2hlvzUNZawchBchhypNHqzlL1mxdw79uAOe8WBjz0Cpbh02TWLsRVXHMNvqamhtTUVEwmX5tdk8lESkoKNTU1JCYm9lx3/PhxduzYwdq1a3n11Vevq5ikpOvff56cHMvh2uN4lcaU3DEkJ8f2+WtYrKaerxUOwqXOG6F7XDRv+y3tX3yEKWYQaQ/9KwNGTA52Wf0qGt7nb5Ix+4dfbrJ6PB5++tOf8u///u89PwiuR2NjB7qu+vx5Xx1ft7fsCEaDkcGkXtdxdh63BhAWR+FFw5F93uoSX3Ow9nosBXeRefdSGts0OiN83F8XDe/zN8mYe89oNFx1YnzNgE9PT6e2thZN0zCZTGiaRl1dHenp6T3X1NfXU1FRwVNPPQVAW1sbSik6Ojp4+eWX+1z09SpvqyJ9YCp2s73fXlP4n3J14tr9Lp7jxRjiUnEULseccRNG2wAguv7hC3EjrhnwSUlJFBQUsHHjRoqKiti4cSMFBQUXLc9kZGSwe/funo9/+ctf0tXVxb/8y78EpurLUEpR0V7FmMEFV73uav1mpIdM8HnL9uPcsQbV3Ypl3EJsN98rzcGEuE692kXz4osvsm7dOubPn8+6det46aWXAFi2bBmHDx8OaIG91eJqpcPTSU7s1R9N/6rfzOVID5ng0bvb6N76Kt1bfo7BHsOAe1dgn/awhLsQN6BXa/DDhg1j/fr1l/z9r3/968te/8wzz9xYVdehsr0agOzYzGteK/1mQodSCu/pnbg+fxvlcWK9+T6s4xdJczAh/CBi/hXVdtUDkD5QHlMPF3pHI84da9EqDmJMGeZrDjbo2j+ghRC9EzEBX9dVT6w1BofZAVx5rV3W2YNPKR3PsU9x7X4PlI7t1iVYRs+R5mBC+FnEBHxtVz0pjuSej6/U213W2YNLbz3vaw5WcwJT5ijsM56Q5mBCBEjEBHxdVwNjv7GDRtbaQ4fSNdyHPsb95QdgMmO/Yynm/BnyJKoQARQRAe/yumn3dMjZqyFKa6zAue119IYyzLmTsN3+mPSPEaIfRETAN3W3AJBgiw9uIeIiSvPg3vcH3Ac2YbAPxD7nacxDp8isXYh+ElEBLx0kQ4dWe9o3a285h3nE7dhv/Q4Gu9zcFqI/RUbAd7UAMoMPBcrjxLXnfTxHtmIYOAjHgucx54wLdllCRKXICHiZwYcEb9URnNvfQLU3YBk1G9stD2CwOoJdlhBRKyICvtnZitVkxSFNxoJCuTpx7XoHz4ntGOJTcSz+Meb0/GCXJUTUi4iA73B3MtA8INhlRCVP6Ze4dqxFOduxTrgb66QiDGZrsMsSQhAhAd/p7mKARZYC+pPe1Yrr83V4z+7BmJSDY+FzmAbnBrssIcTXRE7AmyXg+4NSCu+pz3HufBs8LqxT7sc6fiEGY0R8KwkRUcL+X+WnB6o5U9eA2RvLqrf29fy99JzxP729AeeONWiVhzGmDsc+cymmhIxglyWEuIKwD/jdR2vxDHZi1RMv+nvpOeM/Sul4Sj7B9cXvQClstz2CZfRsDAZpDiZEKAv7gAcwmL1Mzcvi/hHSd8bf9JYaX3Ow8ycxZY7GfscTGGOTr/2JQoigC/uAVyiUQcNmkpN//EnpXtyHNuP+8kMw27Df+T3MI26XNgNChJEICHgNAKvJEuRKIofWUO5rM9BYjnnozdhufxTjgIRglyWE6KPwD3iDL+AtRgn4G6W8bl9zsIObMNhjsM/5IZa8KcEuSwhxnSIg4L2AzOBvlPf8KVzbfoPeeh7zyOnYp31bmoMJEeYiIOBlBn8jlLsb157f4Tn6CYaYRByL/glz1phglyWE8IOICXirBHyfeSsP+5qDdTRhGX2hOZhF+vkIESkiJuAtJul/0lvK2YFz12/xnvwMY0I69nv+FXPaiGCXJYTws8gJeHlUvlc8Z/fg+uxNlLMD68TFWCculuZgQkSosE9FhQLAZDAFuZLQpne14NrxJt6yLzEmDcGx8B8xDR4S7LKEEAEU9gGPQQfAKI/NX5ZSCu/JHTh3/hY0N9ZbHsQ6bgEGo/xAFCLShX3AfzWDN8oTlpfQ2+txFr+BVn0UU9pI7Hc8iTEhPdhlCSH6SdgHPD0BLzPSryhdx1PyF19zMIMB2+2PYRl1lzQHEyLKhH/AG2QG/3Va8zmcxa+j157GlD0W+4wnMMYkBbssIUQQhH3A/+0ma3TPTpXuxX1gE+59fwCLDfudyzCPuE2agwkRxXoV8KWlpSxfvpyWlhYSEhJYtWoVubm5F13zq1/9ik2bNmEymTCbzTz33HPMmDEjEDV/g9xk1RrKcG77DXpjJea8W7Dd9gjGAfHBLksIEWS9CviVK1eyZMkSioqK2LBhAytWrGDt2rUXXTNu3DiWLl2Kw+Hg+PHjPProo+zYsQO7PbBPRqqeJZroC3jldeP+8kPchzZjsMdin/cMltzJwS5LCBEirpmKjY2NlJSUUFhYCEBhYSElJSU0NTVddN2MGTNwOHznoubn56OUoqWlxf8VXyI6A7674iid7/8U98FNWEZOZ+BDr0i4CyEucs0ZfE1NDampqZhMvl0qJpOJlJQUampqSExMvOznfPjhh+Tk5JCWlubfai8rupZolLsb1xfraS/5BEPsYByLXsCcNTrYZQkhQpDfb7J+8cUX/PznP+f111/v8+cmJfW9Pa3R7LuJmJwUS4Ijts+fH066Tu+j/k//h9bWSNwthSTO/A5Ga3Q1B0tOjuz3+HJkzNEhEGO+ZsCnp6dTW1uLpmmYTCY0TaOuro709EsfmNm/fz8vvPACr776Knl5eX0uprGxA11XffocTfPN4JuauvBYI3MWr5wdOHe+jffU5xgHZTCg6CcMHjOR+vp2wBPs8vpNcnLshTFHDxlzdLjeMRuNhqtOjK8Z8ElJSRQUFLBx40aKiorYuHEjBQUFlyzPHDp0iOeee45f/OIXjB4tSwb+oJTC+1VzMFcX1kn3+JqDyeEmQohe6NUSzYsvvsjy5ct59dVXiYuLY9WqVQAsW7aMZ599lrFjx/LSSy/hdDpZsWJFz+f953/+J/n5+YGpPMLpnc24dqzFW74f4+BcHHf/M6ak7GCXJYQII70K+GHDhrF+/fpL/v7Xv/51z/9+//33/VdVn/RtSSfUKaXwnCjGtesd0LzYpj6EZex8aQ4mhOizsH+S9SsGwv+JTb2tDmfxarRzxzCl5/uag8X3x04kIUQkipiAD2dK1/Ec/TOuPe+DwYht+uNYCmZKczAhxA2JnIAP0wm81lSNs/g36HVnMeWMxz79cYwxl3++QAgh+iJyAj7MKM2L+8BHuPf/AYPFgX3W32MeNk2agwkh/CbsA/6rbpLhtAav1Z31tfRtqsI8bKqvOZgjLthlCSEiTNgH/FfCId6V14Vr7wd4Dn+MwRGPY96PMOdODHZZQogIFTEBH+q8547hLH4D1VaL5aaZ2KY9jME6INhlCSEiWAQFfGjO4ZW7C9fu9/Ac+xRDbDKOu/8Zc+aoYJclhIgC4R/woZnrAHgrDuDcvgbV1YJl3AJsN38Lg9kW7LKEEFEi/AP+glDafKJ3t+Ha+Tbe07swDsrCMfcZTCl9b74mhBA3IgICPnRaFSil8J7Zjevzt1DuLqyT78U6oRCDKQL+Mwshwk4EJU9wp/B6RxPOHWvQKg5iTM7DMXMppsSsoNYkhIhuERTwwaGUjud4Ma5d74KuYZv2bSxj5mEwSpsBIURwRUzAB2P+rrfW+pqD1RzHlFHgaw4WlxKESoQQ4lIRE/D9Sek6niMf49rzARhN2O54Ekv+HdJmQAgRUiIo4PsnXLWmSpzbXkevL8WUMwH7jMcxDhzUL68thBB9EQEB3z+7aJTmwb1/I+79GzHYBmCf/QPMebfIrF0IEbIiIOB9Ahm0Wt0Z36y9uRrz8Fux3bYEoz36Tn0XQoSXsA94m3swA13ZWIz+H4ryuHDt/T2ew1swDByEY8E/YM6Z4PfXEUKIQAj/gPckE9ORhtHPpx95q0twFq9GtddjKbgL29SHMFgdfn0NIYQIpLAPeH9Trk5cu9/Fc7wYQ1wqjsLlmDNuCnZZQgjRZxLwX+Mt249zxxpUdyvW8YuwTr4Xg9ka7LKEEOK6SMBzoTnYZ+vwnv0CY2IWjvk/wpQ8NNhlCSHEDYnqgFdK4T29E+fnb4HHhfXm+7COXyTNwYQQESFqk0zvaMS5fQ1a5SGMKcOwz1yKaVBmsMsSQgi/ibqAV0rHc+xTXLvfA6Vju3UJltFzpDmYECLiRFXA663nLzQHO4EpczT2GU9gjEsOdllCCBEQURHwStdwH/oY95cfgMmM/Y6lmPNnSJsBIUREi/iA1xorcG77DXpDOebcSdhuf0yagwkhokLEBrzSPLj3/QH3gU0Y7AOxz3ka89ApMmsXQkSNiAx4rfa0rzlYyznMI27Hfut3MNhjgl2WEEL0q14FfGlpKcuXL6elpYWEhARWrVpFbm7uRddomsbPfvYztm/fjsFg4KmnnuLBBx8MRM1XpDxOXHvex3NkK4aYRBwLn8ecPa5faxBCiFDRq4BfuXIlS5YsoaioiA0bNrBixQrWrl170TV//OMfqaioYMuWLbS0tHDvvfdy6623kpXVPwdPe6uO4Nz+Bqq9Acuo2dhueUCagwkhoto1N383NjZSUlJCYWEhAIWFhZSUlNDU1HTRdZs2beLBBx/EaDSSmJjInDlz2Lx5c2Cq/hqz8nBH5xa6N/1/YDTjWPxj7NMfk3AXQkS9a87ga2pqSE1NxWQyAWAymUhJSaGmpobExMSLrsvIyOj5OD09nfPnz/epmKSkvq+T3zvSRfLBYyTc9i0SZjyEMYqagyUnR9+hIzLm6CBj9o+Qusna2NiBrvftCL68KdMZPGcWjW1eGptdgCswxYWY5ORY6uvbg11Gv5IxRwcZc+8ZjYarToyvuUSTnp5ObW0tmqYBvpupdXV1pKenX3LduXPnej6uqakhLS2tzwX3lcFgxGiT5RghhPimawZ8UlISBQUFbNy4EYCNGzdSUFBw0fIMwIIFC1i/fj26rtPU1MTWrVuZP39+YKoWQghxTb3qsPXiiy+ybt065s+fz7p163jppZcAWLZsGYcPHwagqKiIrKws5s2bx0MPPcQPf/hDsrOzA1e5EEKIqzIopfq26B1A17MGD7JmFy1kzNFBxtx7N7wGL4QQIjxJwAshRISSgBdCiAgVUvvgjcbr7/R4I58brmTM0UHGHB2uZ8zX+pyQuskqhBDCf2SJRgghIpQEvBBCRCgJeCGEiFAS8EIIEaEk4IUQIkJJwAshRISSgBdCiAglAS+EEBFKAl4IISJU2AR8aWkpDz/8MPPnz+fhhx+mrKzskms0TeOll15izpw5zJ07l/Xr1/d/oX7UmzH/6le/4u677+aee+7hvvvuY/v27f1fqB/1ZsxfOXv2LOPHj2fVqlX9V2AA9HbMmzZtYvHixRQWFrJ48WIaGhr6t1A/6s2YGxsbeeqpp1i8eDELFizgxRdfxOv19n+xfrBq1SpmzZpFfn4+J0+evOw1AckvFSYee+wx9eGHHyqllPrwww/VY489dsk1H3zwgVq6dKnSNE01NjaqGTNmqMrKyv4u1W96M+bi4mLV1dWllFLq2LFjavLkyaq7u7tf6/Sn3oxZKaW8Xq969NFH1fPPP6/+4z/+oz9L9LvejPnQoUNq4cKFqq6uTimlVFtbm3I6nf1apz/1Zsw/+9nPet5bt9utHnjgAfXRRx/1a53+smfPHnXu3Dl11113qRMnTlz2mkDkV1jM4BsbGykpKaGwsBCAwsJCSkpKaGpquui6TZs28eCDD2I0GklMTGTOnDls3rw5GCXfsN6OecaMGTgcvjNp8/PzUUrR0tLS3+X6RW/HDPDaa69x5513kpub289V+ldvx/zGG2+wdOlSkpOTAYiNjcVms/V7vf7Q2zEbDAY6OzvRdR23243H4yE1NTUYJd+wm2+++ZJzrL8pEPkVFgFfU1NDamoqJpMJAJPJREpKCjU1NZdcl5GR0fNxeno658+f79da/aW3Y/66Dz/8kJycnH457DwQejvm48ePs2PHDp544okgVOlfvR3zmTNnqKys5JFHHuFb3/oWr776KipM+wT2dsxPP/00paWlTJ8+vefP5MmTg1FyvwhEfoVFwItr++KLL/j5z3/Of/3XfwW7lIDyeDz89Kc/5aWXXuoJiGigaRonTpxg9erVvPnmmxQXF7Nhw4ZglxVQmzdvJj8/nx07dlBcXMzevXvD9jfyYAmLgE9PT6e2thZN0wDfN3tdXd0lv/Kkp6dz7ty5no9ramrCdjbb2zED7N+/nxdeeIFf/epX5OXl9XepftObMdfX11NRUcFTTz3FrFmzWLNmDe+99x4//elPg1X2Dent+5yRkcGCBQuwWq3ExMQwe/ZsDh06FIySb1hvx7xu3TruuecejEYjsbGxzJo1i927dwej5H4RiPwKi4BPSkqioKCAjRs3ArBx40YKCgpITEy86LoFCxawfv16dF2nqamJrVu3Mn/+/GCUfMN6O+ZDhw7x3HPP8Ytf/ILRo0cHo1S/6c2YMzIy2L17N5988gmffPIJjz/+OA899BAvv/xysMq+Ib19nwsLC9mxYwdKKTweD7t27eKmm24KRsk3rLdjzsrKori4GAC3283OnTsZMWJEv9fbXwKSXzd0i7YfnT59Wj3wwANq3rx56oEHHlBnzpxRSin1ve99Tx06dEgp5dtZsWLFCjV79mw1e/Zs9c477wSz5BvWmzHfd999aurUqeqee+7p+XP8+PFgln1DejPmr/vFL34R9rtoejNmTdPUK6+8ohYsWKAWLVqkXnnlFaVpWjDLviG9GXN5ebl64oknVGFhoVq4cKF68cUXlcfjCWbZ1+3ll19WM2bMUAUFBeq2225TixYtUkoFPr/kRCchhIhQYbFEI4QQou8k4IUQIkJJwAshRISSgBdCiAglAS+EEBFKAl4IISKUBLwQQkQoCXghhIhQ/z/OFBW4YBaHgQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "mu0, mu1 = np.mean(x0), np.mean(x1)\n",
    "s0, s1 = np.std(x0), np.std(x1)\n",
    "c_seq = np.linspace(0, 3, 1000)\n",
    "\n",
    "fpr_binorm = 1-norm.cdf(c_seq, mu0, s0)\n",
    "tpr_binorm = 1-norm.cdf(c_seq, mu1, s1)\n",
    "\n",
    "plt.plot(fpr, tpr)\n",
    "plt.plot([0,1], [0,1]);\n",
    "plt.plot(fpr_binorm, tpr_binorm);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29abb06b-aa8a-45de-9e36-380f7dfe802c",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "#### AUC\n",
    "+ The area under the ROC curve is given by\n",
    "$$\n",
    "AUC = \\int_0^1 T\\{F^{-1}(f)\\}df.\n",
    "$$\n",
    "+ The ideal test has AUC = 1\n",
    "+ A completely uniformative test has AUC = 0.5\n",
    "+ An informatively bad test has ROC < 0.5.  \n",
    "+ The AUC has a nice interpretation. Let $X_1 \\sim 1-T$ and $X_0 \\sim 1-F$ independent random values. Then \n",
    "$$\n",
    "AUC = P(X_1 > X_0)\n",
    "$$\n",
    "+ It can be shown that the Wilcoxon Rank Sum Test is a test of $AUC=0.5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "bb3d0049-519f-471d-a9f5-da83dac4b3da",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7588\n",
      "0.7587999999999999\n"
     ]
    }
   ],
   "source": [
    "## Calculating AUC directly as the faction of instances where X1 > X0\n",
    "auc = np.sum([j >= i for i in x0 for j in x1]) / (len(x0) * len(x1))\n",
    "print(auc)\n",
    "\n",
    "## using sKlearn\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "fpr, tpr, thresholds = roc_curve(y, x)\n",
    "print(auc(fpr, tpr))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8092871f-4811-4cca-a198-bab9a094fc30",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Binormal AUC\n",
    "\n",
    "$$X_1 - X_0 \\sim N(\\mu_0 - \\mu_1, \\sqrt{\\sigma_0^2 + \\sigma_1^2})$$\n",
    "\n",
    "We can just directly calculate \n",
    "\n",
    "$$P(X_1 > X_0) = P(X_0 - X_1 < 0)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "cc6badf1-2b60-4f99-bf78-936dafd09bc9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.766401836491996"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "stats.norm.cdf(0, loc = mu0 - mu1, scale = np.sqrt(s0**2 + s1**2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:.conda-ds4bio]",
   "language": "python",
   "name": "conda-env-.conda-ds4bio-py"
  },
  "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}