bayesian-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian analysis of performance metrics \n",
    "\n",
    "This notebook reevaluates the performance of a binary classifier by using the Bayes rule. As an example, we'll use the performance of two lie detectors that are being proposed to perform massive screening of migrants and/or refugees.\n",
    "\n",
    "Although the scientific consensus is that lie detectors do not exist [0], in this post we will not discuss this issue neither we won't discuss assumptions of the methods, experimental validation, etc. Rather we will try to analyze the actual feasibility of the proposed solutions to detect rare events through massive screening of populations using tools with low specificity. \n",
    "\n",
    "## Summary of a test\n",
    "\n",
    "A test is a method to detect events such as cancer, spam, deception, etc. There can be two outcomes of the test for the binary case: *positive* (the test says there is a disease, deception, etc.) and *negative* (the test says there is no disease, lie, etc.). There are four results of a test: \n",
    "- **True positive**: the test says the person lies, and the person lied. \n",
    "- **False positive**: the test says person lies, and actually the person didn't lie.\n",
    "- **True negative**: the test says the person doesn't lie, and they didn't. \n",
    "- **False negative**: the test says the person doesn't lie, and they lied.\n",
    "\n",
    "These outcomes can be represented in the [confusion matrix](https://en.wikipedia.org/wiki/Confusion_matrix), for instance, the reported performance of iBorderCtrl is (see Table IV in [5]): \n",
    "\n",
    "| Test Result\\ Ground Truth | Lie    | no-lie |\n",
    "| ------------------------- | ------ | ------ |\n",
    "| Test positive             | 73.66% | 24.45% |\n",
    "| Test negative             | 26.34% | 75.55% |\n",
    "\n",
    "How we interpret this performance: \n",
    "\n",
    "- If the person lies, there is a 73.66% of chances of having a positive test (true positive) and 26.34% of having a negative test (false negative). \n",
    "- If the person tells the truth, there is a 24.45% of chance the test will be positive (false positive) and 75.55% the test will be negative (true negative). \n",
    "\n",
    "Some notes on tests and probability (adapted from [1]):\n",
    "- **Tests are not the event**. A test for detecting lying is a separate event from actually lying. \n",
    "- **Tests are flawed**. Tests detect things that don’t exist (false positive) and miss things that do exist (false negative). People often use test results without adjusting for test errors.\n",
    "- **Bayes’ theorem** estimates the real probability of an event from the test results.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes' Theorem\n",
    "\n",
    "\n",
    "![Bayes Theorem summary](pics/Bayes_Theorem.png)\n",
    "\n",
    "Bayes' Theorem summary. Source <a href=\"https://betterexplained.com/articles/an-intuitive-and-short-explanation-of-bayes-theorem/\">Betterexplained</a>\n",
    "\n",
    "The probability of an event is the number of ways it could happen given all possible outcomes. In our case:\n",
    "\n",
    "\\begin{align}\n",
    "\\text{Probability} &= \\frac{\\text{desired event}}{\\text{all posibilities}} \\\\\n",
    "P(\\text{Lie}\\mid\\text{+}) &= \\frac{P(\\text{+}\\mid\\text{Lie}) P(\\text{Lie})}{P(+)} \\\\\n",
    " &= \\frac{P(\\text{+}\\mid\\text{Lie}) P(\\text{Lie})}{P(\\text{+}\\mid\\text{Lie}) P(\\text{Lie}) + P(\\text{+}\\mid\\text{No-lie}) P(\\text{No-lie})} \n",
    "\\end{align}\n",
    "\n",
    "where:\n",
    "- $P(\\text{Lie}\\mid\\text{+})$ is the probability of having a liar ($\\text{Lie}$) given a positive test ($\\text{+}$), formally the *posterior probability*. \n",
    "- $P(\\text{+}\\mid\\text{Lie})$ is the chance of having a positive test ($\\text{+}$) when the person is lying $(\\text{Lie}$). This is the chance of true positive (73.66% in our example).\n",
    "- $P(Lie)$ is the chance of having of persons that would lie in an interview, formally *prior probability*. It can be also referred to as *frequency* or *prevalence*.\n",
    "- $P(No-lie)$ is the chance of having a person that is telling the truth ($P(No-lie)=1-P(Lie)$). \n",
    "- $P(\\text{+}\\mid\\text{No-lie})$ is the chance of having a positive test ($\\text{+}$) when the person is not lying $(\\text{No-lie}$). This is a false positive (24.45% in our example).\n",
    "\n",
    "The posterior probability is also referred to as **positive predictive value** (PPV) and it answers the question: **provided we have a positive test, what is the probability the person is lying?**\n",
    "\n",
    "Analogously, we can use the rule to calculate the **negative predictive value** (NPV): \n",
    "\n",
    "\\begin{align}\n",
    "P(\\text{No-lie}\\mid\\text{-}) &= \\frac{P(\\text{-}\\mid\\text{No-lie}) P(\\text{No-lie})}{P(\\text{-})} \\\\\n",
    " &= \\frac{P(\\text{-}\\mid\\text{No-lie}) P(\\text{No-lie})}{P(\\text{-}\\mid\\text{No-lie}) P(\\text{No-lie}) + P(\\text{-}\\mid\\text{Lie}) P(\\text{Lie})} \n",
    "\\end{align}\n",
    "\n",
    "where $P(\\text{No-lie}) = 1 - P(\\text{Lie})$. It answers the question **provided we have a negative test, what is the probability the person is telling the truth?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formulating hypothesis of scenarios\n",
    "\n",
    "### Correcting the confusion matrix\n",
    "\n",
    "We can use the prevalence or prior probability to update the confusion matrix to estimate the chances of all the four possible results of the test. Assuming $P(Lie)=0.15$:\n",
    "\n",
    "| Test Result\\ Ground Truth | Lie    | No-lie |\n",
    "| ------------------------- | ------ | ------ |\n",
    "| Test positive             | True Pos | False Pos |\n",
    "| Test negative             | False Neg | True Neg |\n",
    "\n",
    "\n",
    "| Test Result\\ Ground Truth | Lie (15%) | No-lie (85%)|\n",
    "| ------------------------- | ------ | ------ |\n",
    "| Test positive             | 0.15x0.7366=<br>0.1105 | 0.85x0.2445=<br>0.2078 |\n",
    "| Test negative             | 0.15x0.2634=<br>0.0395 | 0.85x0.7555=<br>0.6422 |\n",
    "\n",
    "Assuming $P(Lie)=0.01$:\n",
    "\n",
    "| Test Result\\ Ground Truth | Lie (1%) | No-lie (99%)|\n",
    "| ------------------------- | ------ | ------ |\n",
    "| Test positive             | 0.01x0.7366=<br>0.0074 | 0.99x0.2445=<br>0.2421 |\n",
    "| Test negative             | 0.01x0.2634=<br>0.0026 | 0.99x0.7555=<br>0.7479 |\n",
    "\n",
    "In absolute numbers, in a population of 1000 persons with 10 liars (1%), the test will discover 7 liars, but will wrongly label as liars 242 innocent people. In other words, to discover one actual liar the test will wrongly classify 34 people that are actually telling the truth. The test will correctly classify as non-liars 747 people. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Posterior probability\n",
    "$Lie$ is called *hypothesis* whose probability may be affected by data (evidence or test). By varying $P(Lie)$ we can observe the variation of the posterior probability under different hypotheses. For example, we can calculate the *posterior probability* (**positive predictive value**, PPV) with the hypothesis of 15% of liars ($P(Lie)=0.15$) in the population that will be screened:\n",
    "\n",
    "\\begin{align}\n",
    "P(\\text{Liar}\\mid\\text{+}) &= \\frac{0.7366 \\times 0.15}{0.7366 \\times 0.15 + 0.2445 \\times 0.85} \\\\\n",
    "&\\approx 34.71\\%\n",
    "\\end{align}\n",
    "\n",
    "And the expected rate of false positives (**false discovery rate**, FDR): \n",
    "\n",
    "\\begin{align}\n",
    "P(\\text{No-liar}\\mid\\text{+}) &= 1-P(\\text{Liar}\\mid\\text{+})=1-0.3471  \\\\\n",
    "&\\approx 65.29\\%\n",
    "\\end{align}\n",
    "\n",
    "This means that, in a scenario of 15% of liars, 34.71% of positive tests will correspond to actual liars crossing the borders and 65.29% will correspond to false positives. This performance decay that might seem counterintuitive can be explained by the [base rate fallacy](https://en.wikipedia.org/wiki/Base_rate_fallacy).\n",
    "\n",
    "Analogously, we can calculate the **negative predictive value** (NPV) with the Bayes' Theorem:\n",
    "\n",
    "\\begin{align}\n",
    "P(\\text{No-liar}\\mid\\text{-}) &= \\frac{0.7555 \\times 0.85}{0.7555 \\times 0.85 + 0.2445 \\times 0.15} \\\\\n",
    "&\\approx 94.20\\%\\\\\n",
    "\\end{align}\n",
    "\n",
    "And the expected rate of false negatives (**false omission rate**, FOR)): \n",
    "\\begin{align}\n",
    "P(\\text{Liar}\\mid\\text{-}) &= 1-P(\\text{No-liar}\\mid\\text{-})=1-0.9446  \\\\\n",
    "&\\approx 5.79\\%\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python example code\n",
    "\n",
    "We can better observe how the test behaves for different base rates by testing different values of $P(Lie)$.\n",
    "\n",
    "The following code calculates the posterior probability for finding an actual liar given a positive test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rc('text', usetex=True)\n",
    "\n",
    "# E is the test outcome\n",
    "# H is the hypothesis\n",
    "posterior = lambda PrEH, PrH, PrEnH: (PrEH * PrH) / (PrEH * PrH + PrEnH * (1 - PrH))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### iBorderCtlr\n",
    "\n",
    "The [iBorderCtlr](http://iborderctrl.eu/) is an EU funded project to develop a intelligent portable control system for borders. The project features many of the components of smart borders systems (biometrics, document authentication, etc.) but it also performs *risk-scoring* of travelers and *automated deception detection* to evaluate travelers and migrants at the border. \n",
    "\n",
    "We can get $P(\\text{+}\\mid\\text{Lie})$ and $P(\\text{-}\\mid\\text{No-Lie})$ iBorderCtlr last publication [5]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PrHE=0.3471089958060412 - PrnHnE=0.9420406786125556\n"
     ]
    }
   ],
   "source": [
    "# Pr(E|H) Chance of a positive test (E) given there is a liar (H). \n",
    "# This is the chance of a true positive. Also sensitivity\n",
    "PrEH = 0.7366\n",
    "# Chance of a positive test (E) given that the person wasn't a lying. \n",
    "# This is the probability of a false positive, also (1-specificity)\n",
    "PrEnH = 0.2445\n",
    "\n",
    "# Test the example\n",
    "PrH = 0.15\n",
    "PrHE = posterior(PrEH, PrH, PrEnH)\n",
    "PrnHnE = posterior( (1-PrEnH), (1-PrH), (1-PrEH) )\n",
    "\n",
    "print('PrHE={} - PrnHnE={}'.format(PrHE, PrnHnE))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can build the posterior probabilities of true positives and true negatives and annotate a specific hypothesis of 10% of people lying to pass the border:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.050010000000000006"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEFCAYAAAAPCDf9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeVxc1fn48c9hDwQYIGQPSYbsMRoHSOISjWZQ61JbA0Fba2s1jNa2trYNjbZ2dQFtqz9rLUTtt9aNQLXuCxOTuJsAmtUYw2QjCVmAAQJhP78/7sxk2CdhNobzfr3mFebemTvPTOC5Z8495zlCSomiKIoyPAT5OgBFURTFe1TSVxRFGUZU0lcURRlGVNJXFEUZRlTSVxRFGUZU0lcURRlGQnwdwEBGjRolp0yZ4uswFEVRhpSysrLjUsrE7tv9PulPmTKF0tJSX4ehKIoypAgh9vW2XXXvKIqiDCMq6SuKogwjKukriqIMIx5L+kIIQz/7MoQQRiHESk+9vqIoitKTR5K+EMIIrO5jnwFASmkGrP2dHBRFURT38kjStyX0mj52ZwFW288WwOiJGBRFUZSefDFkU0fXE0KCR17l9bsgeiwkLYIJKRAW5ZGXURRFcZf2jk52VjXwxQErn++3ctvFeqaPiXbra/jlOH0hRDaQDZCUlHT6B2iug7J/gezU7geFwLhzIOk8mLRQOxGMHO3GiBVFUU5fTWMrn++vpWyfdttSWcfJtg7H/nOTdAGR9K1AvO1nHVDd/QFSygKgACA1NfX0V3kJCoWMp2H/p7D/E6jaCgfLtNsnf9ceE5+snQSSFmn/JiSDENS99hpH//YI7YcPEzJuHKN//jNir7nmDN+qoiiKprNTsvvYCUeCL99Xi+V4Y4/HTU6I5NxJOuZP0nHhtFFuj8NrSV8IoZNSWoFCINW2WQ+Y3f5iYZEw99vaDaClASo3nToJVJZCTYV2++JZ7TGRo6izTufw6weRre0AtB86xOHf3gugEr+iKKelobmNLw5YKd9npWx/LZ/vr6Whub3LY8JDgjhnko6UyXEYkuIwJOlIGBnu0bg8kvSFEBlAqhAiQ0pZbNu8FkiRUpYLIVJtI3ysUspyT8TQRXg0JF+q3QA62qBqC+z/TDsJ7P8UGo9y9O0gZGvXj0Q2N3P0oVxiv3E5hIR5PFRFUYamqrpmPttTzaa9NZTureWrIw10X412gm4E5yZpST5lchyzx8UQGuzd6VLC39fITU1NlR6vvSMl1Fj48oKr+3oAs79Xp10PmHIhTFkM489VJwFFGaaklOytbmLTnho+21PDxr3VHKg52eUxocGCueNjMSRpCd4wWce42BFei1EIUSalTO2+3S8v5HqdEJCQTMj48bQfOtRjd0h0ELQ1gWWddgMIjVQnAUUZJjo7JV8daWDjnho27q1h454ajjW0dHlMdHgIKVPiWDA1nrQp8cybEEtEaLCPIu6bSvpORv/8Zxz+7b3I5mbHNhERweh7/wiXLIS9H566Hf+q75OAfgmMmw/B6uNVlKGovaOTrQfrtCS/p4ZNe2uo79YfnxAV5kjwC6bGM3tcDMFBwkcRu05173Tj8uidE0d7ngSchcfC1MXaCUC/BBKmad8oFEXxO52dki+r6vmkopqPK6rZuKeGEy1dk/wE3QgWTI13JPrkxCiEH/9N99W9o5K+u9hPAnvehz0boMbSdX/MxFMnAP3Fap6AoviQlJKKYyf4uKKaj3dX8+meaqxNbV0eM3VUFIv0p5L8xLhIH0V7ZlTS97bafVryt6zXbk3dpiOMnqudAJIv0eYJhI/0foyKMowcqGni44rjWqKvqO7RJz9BN4LzkhM4PzmB85ITvHrR1RNU0velzk44uh0q1mkngH0fQ7vTlf6gUJi0AKYthWnpMHae6gpSlEGyNrXy4e7jfLDrOB/uPs5Ba9fRNaNGhnO+U5JPio/06+6a06WSvj9pb4EDG099CzhUfqpkBMDIMTDNqN2SL4ERcb6KVFGGjLaOTj7fb+WDr4/x/tfH2VJp7TJOPnZEKIv08ZyfPIrzkxOYNnpkQCX57lTS92cna7VrAV+XwO610OA0bFQEwcQ07RvAtKXaqKAgtfaNogDsPd7oSPKfVFR3ufgaGixInRzPRTMSWTx91JAZXeMuKukPFVLC0R22E4BZmy3c6XSBKSoRkpfavgVcClGeKVKqKP6oobmNj3ZX88HXx/jg6+Psr2nqsn/a6JEsnj6Ki6YnslAfT2TY8B02rZL+UNXSAJYN2glgtxnqDjjtFNrcgJlXwIxvQOJMdS1ACSjaKJtG1u08yns7j7Jpbw3tnadyli4ylAumjeKi6aNYPD2R8bqhffHVnVTSDwRSwvFdtm8BJdoF4Y7WU/vjpsLMb8CMK2Dy+RAc6rtYFeUMNbd18KmlmnU7j7Luq2NdWvPBQQJDko6Lpidy0YxEzpoQO6y6bE6HSvoekpeXx6ZNm7BarZjNZjIyMgAwmUwYjR5eFKylASreg6/ehq/f6TosNDwWphu1bwDTjV0uBpvNZoqKisjPz3dbKCaTyXE8IQT5+flkZ2cDkJmZyerVq9HpdMTFxbF8+XJqarR1dFavXk1mZia5ubkYDKdWzszLy8NgMKDX6ykuLmblSteXU87JycFisZCWltbr8+z79Xo9ubm5jhgB4uPjHe8jMzMTq9WKwWBwPE7xjEPWk6z76ijrdh7lo93VXWrKx0WGsmTmaC6ZNZqLpycSG6kaM65QtXc8xDmpJCcnU1RU5L0XD4+GOddqt84OrXz0V2/Brrfh2E7Y9l/tJoIpPjqVjO/eDLOuxGg0uvWEVFxcTHp6uuO+wWDokvSd6fV6R1I1m83k5ORgMpkoLCzskvRLSkq6fLbl5eVd9velvFwr2lpUVNQluduZzWbHfpPJhMVioby8nPT0dLKzsykoKMBsNmOxWFi1ahUGg4HMzEyXX19xTWenZMvBOt7dXsV7O4+ys6qhy/6542O4dJaW6M+ZqFOteTdSw0ACRVCwtiBM+h/gjs/gp5/D5Q/A1ItACArXb4F374H/dy48cSGsz4UjO+hR+/UM5OfnO77hgNZaNhqNjgTbl9TUVCwWCxkZGRQXFzu2WywWdDqd435GRobL30rMZrPjBJSWltYjhpKSEsf+lJQUzGYzer0eq1VbtrmiogIAo9HoSPLp6en487fNoaK1vZMNu45xz8tbOe/BtXzr8Y/4x/oKdlY1EBkWzGVzxvDgdfP47O6lvPHTxfzispkYkuJUwnczlfQ9xGKxYDKZyMzMZP369ZhMJkDrtrAnuMzMTNLT07skPOfnp6enYzKZSElJcSSl7s+x3zeZTOTl5WG1WrVtWSYKtobA918jp+5GzAfCMX00DmvnSCw7N5Nzz71kXnoO1txzoOR35N3zU8wlJT1ew2w2k5mZiclkcrSiXWEymQbsEikoKHAkYIPB4Dh+cXGx4/MC7duBq0m3urqa+HhtYTadTudI4nYJCQmOzxK0JG8wGKiuriY5ORnQEr7zt4OysjJSU3t8S1ZcUN/cxqubD/Hj58tJ+VMJ3396I899tp8j9S2Mj43gB+dP4T+3LODze9MpuCmV6xckMSYmwtdhBzTVveNBa9asYc+ePY7+a2d5eXmOLoXMzMwuLWVn+fn55OTkUFpa2qUbIjMzE51O5+i3zszMdLSG7V1MKSkpZGdnk/uXRynfsoP8khJob6FmQyFUPkrW1ErWfPQ12c2PUPJiIysTzOS9kEj6uQvJ/vUDZGZdT3x8vKObozfdu0/s9Ho9Op0Oi8XS4/Emk4mampoufepZWVmOLp7uXTvdn9/XN4jeupO6MxqNjm8mJSUlpKWlOU42FRUVPbqEysvLqampUV07p6GqrpmSL4/w7vYqPrVU09Zx6tvkrLHRXDZ3LJfNGcPc8TEBPTnKX6mk70HLly9Hp9P1mvQrKiqoqamhrKwMq9WK1Wrt0qUBOBKPvXXa/TkhIdp/n9VqdbRuQTuhVFdX90i4AISEw9QLYdp2Mu6/j/SLz8M46Vz0o4ug/iAVm3dT89nHlL3yT6wtOu75w93cd/99WOvqyc/P75HgnVvN3a1atapHa9/5QqmzjIwMcnJyWLVqVa8nEfv70+v1/Sb3hIQEx+dttVpJSOg6j8GevE0mE3q9Hr1eT2FhoeMbh/0bzsqVK7FareTn53v3Os0QdaCmibe2HebNrVV8ceDU70SQgIVT4x2JflL80CpaFohU0veg7kkctGSv1+tJSUkhPj6+zxZ+b7o/p7y8nJKSEjZt2uRIrnl5eej1elauXNlrt1EXwSHoxk2huHoapsdKYCykWH9PfN02MibXAnXweQ5FafGYW+eTn/sbcp/4j3b9wEav1/d6UgMtwZaWlnY5IfXHYDCQk5PjGEnjzP4aA7X0jUYjhYWFGI1GSkpKunQT2TmPzsnNzaWmpoby8nLHc9LS0gBtlI87RzgFmgM1Tby59TBvbj3M5so6x/aI0CAump5I+pwxLJ09hvgotbiQP1FJ3wvsfdL20SKAo4smPz8fnU7nUmuy+3NWrVrluOiZk5NDbm4uRqORFStWUGLrn3eWnp7eI4llZWWRk5Pj6E7JfuR17TXWVaLjBGlxtZRsPQAcINcYAX+drY0WmnsdTFqITqcbsLXfWxLvjf0aSH+JdqCWvsFgoLCwkMzMTPR6fZeLsSUlJY5rHvbPD059runp6ej1ejIyMsjLy8NsNpOSkuJ4H6dzgg5U+6obeXNrFW9uPczWg6cSfWRYMJfOGs1V88axZOZoRoT534pRikaN0x/C7C1VvV7vkbH3gDa658h22P6SNvyzdu+pfTETYO63Sb/PTMmGTz02G9hqtbJixQrVzeIje4838oatRb/9UL1je1RYMEtnj+HKeWO5eIZK9P5GTc4KQOXl5V2SvMlk8uwFRynh0Oe2E8DLUF8JQPGONohKJOO7P4Szl8Po2W592ZycHLKystTFVC86Wt/Ma1sO8+oXB7t03YwMD2Hp7NFcOW8cF89I9Ms1YBWNSvqKe3V2apPBtr8E21/G9MJectMj0EUIGHs2nHM9nJUB0WMG9TJWq5UHHnhAzYj1grqTbbyzrYpXNh/kk4pq7CVuosKCSZ8zhivnjeMileiHDDUj14OEEF36e+1dLs4sFgv5+fl+lbwG1SUUFARJC7Xb5feTv+xj2FqE6c9Pk3/lFqjagjj/x2SkjQfdJIgeR+5Df0Gv13e5tuE8Iqi3Eg01NTU9RuAMxNUyDN33O5eSgN5LMwSa5rYO3tt5lFe+OMi6ncdo7dDWdQgLDmLJzESunT+BpbNHq0QfQFTSdwO9Xj9k+puLi4sdJyi3lWMICoapiykuO0L6XYvg7CjYXIg+rpCiK08AX0Lofvgil8w7drLqTw9jSE3DarWydOlS1q5di06n61Giwbkf311lGHrbHx8fT05ODmvWrHG8vr20hHNpBo/XUvKSjk7JJxXV/O+Lg7yzrYoGWw16IeA8fQLfOnc8V8wdp2rcBCg1I3eYKSws9Nix8/Pzybj+O9ronhueB10SXPmwtghMWyOWDc/D/o8xbPgevPsbdK1VrFq1ioKCgh7HMhqNjlFB7izD0Nt+nU7XYw5Cb6UZhjrLsRM89M5OLsx9jxuf+oziskoaWtqZNyGW31w1m09+vZQXsheRlZakEn4AU0nfAxylENLTeyS03soaeLQcg1MMOTk5mM1mTCYTVqsVi8XiGBdvP659qGL31zijcgwiCBasgFvN8JNyyuOuIS15FDQcho8fg38sxLDrr2x6twia67s81Ww2O7rC3FmGYaD9dr2VZhiK6k628fxn+7nuHx9x6V828Pi6Cg7XNTMpfgQ/XTqdtb+4mNd+ciG3LtYzNlaVPxgOVPeOG1gsFkf/r720QPdSCHZFRUVdyhp4tRxDbm6XET/2vvOsrCzWrFlDdna2owRC97jOpBxD988lLe06qpsmwg+XwebnYet/se7bgnVPMzw8A8tXtWReeQmWI/WsWrWq10Q72DIMruqvNIO/6+iUfLj7OMVllbyzvYrWdq2fPjIsmKvmjSMjZSJpU+IJUoXMhiWV9N2gtz79vkoh5OTkkJOT45ji77NyDE4yMjJIT0/vUmis+2vcc8893HfffY64XSnH0P1zsV/MJinPdgH4AUr/dCfpKR9A+2H0Me0ULSinvHEsmXfeRsZl50PMeMB9ZRgG2m/XV2kGf7bneCNrSg/wUnklR+pbHNvPT05gmWEiV5w1lqhw9Sc/3KnfAA/orxSCPRGazWby8/P9oxwDWleHc3XL3uJyjrv7KKT+yjE4PwZwXBS1NrWS+/x7lJWVQUc1PL8YomMwcBiDroniW6eTce1VYPgeNTXaAjGDLcPgSpkG0NZG6K00g79pae/g7W1VvLjxAJ9YTi2ikxQfSUbKRL597gRV70bpQiV9D+ivFEJeXp5ju321KL8sx9DtNdLS0rrE3d1A5Rjs7InWfoyioiLbNxsdRI+Dn2+E3WtZnfQUS3/zXzLmvKOtCnakBdY9gN5w06DKMPS13z6M1GQykZOT02tpBn+y++gJXty4n/+WV1Lb1AZoNW+uPns8WWmTSJ0cpypYKr1Sk7MCgFfKMbjAnljdpvE4bH4R6wdPsuKZ7RRlRoII1tYBTr0Z9Jdq8wWGiea2Dt7adpgXPjvAxr2nvlXNGRfDDQuTuHb+eGIi1KgbRePVGblCiAzAChiklHn97NdLKXuO13Oikv7AvF6OoQ/2biR3t4pzVq4ka/EMDM0fwpevQac2rpy4KZByM5x7I0SNcutr+pPdR0/w7Kf7ePnzg9Sd1Fr1kWHBXDt/PDcsSGLehFjVqld68FrSF0IY0JJ5sRAiGyiVUpZ324+UslwIYQRqnPd3p5L+0GLvuumtrPSZ6FGGoeEIfP4fKPs/qDugbQsO0+YGpN6iLRkZAAmwo1Oy9ssjPPPJPj7cfdyx/eyJsdywIIlrzhnPSHVRVumHN5N+LlAipTTbknqX1r4t6edKKdOFENmqpa+ckc4O2G2G0qdh1zuA7fd47Nmw8DY4axmEDr1x5zWNrby4aT/Pfbqfg9aTAIwIDeZb507guwuTOGtCrI8jVIYKb9be0QHOwzi6jImztfAtQogKIMcDr68MB0HBMONy7Wbdr7X8y/4NVVvglR9Byb1av3/qLRAzztfRDmhLpZV/f7yP17Yccoyrn5IQyffOm0JGykRiR6i+esU9vP79UAihQ+vPzwdWCyHKpZSWbo/JBrIBkpKSvB2iMtTokmDpvXDRSq3m/2dPQNVWeP8h+PBvMPfbsPB2mJji60i7aG3v5I2th/j3x/scSwwKAZfOGs1N503moumJagKV4nae7t7JQOvfd+7eWQkUSCmtvXX/dKe6d5TTJiXs/wQ+fQJ2vg5SazkzIRUW3a71/wf7ruVsbWrl+Y37+ffHex2TqGIiQshKm8SNiyYzOSHKZ7EpgcOb3TuFgP2F9IDZFoBOStllILftxDA05rYrQ4cQMPl87WbdDxtXQ/m/4WAp/PcWePe3WvJP+QFExHgtrL3HG/nXR3tYU1rJybYOAKaPHskPL5zKt+ZPUCtPKV7hqSGb2YAFpyGZtrNOiu3nlbb98epCruIVrY2wpRA+/Scc/0rbFh6jJf5FtzvKPbiblJLSfbU8+YGFd3ccwf7ntnj6KG5drOei6aPUcEvFI9TKWYoC2opfu0vgo/8H+z7UtgWFass8nv8Tty312N7RyVvbqnjyA4tjucGw4CCunT+eWxZPZdZY733DUIYnlfQVpbuDZVry//LVU/3+0y+D838KUy48o/H+zW0d/Le8kvwNFvbXNAEQFxnKjYsm873zJjM6eugNI1WGJpX0FaUvNRb45HH4/Dlo18bGM94Ai++CmVe5VOrhREs7z3+2jyc/2MPRBu3i7JSESG5drGeZYaLqr1e8TiV9RRlIYzVsWg0bC6DJVrFy9BxY/Att2GdQz8Rd09jK/320h39/ss9RImH2uBh+tCSZK+eNI1gNuVR8RCV9RXFVa5NW6uGjR6H+oLYtPllL/mcvh+BQDtedZPX7e3hh437HSJy0KXH8aMk0lsxMVBdnFZ9TSV9RTld7C2x+AT74K1j3aZuiJ/JGzPXcvfdsGju0Ec9LZibyoyXTWDA1vr+jKYpXqaSvKGeqo53ajc/Tvv4hElv2A1Al4/hw9HeYffVPmTt5rI8DVJSe+kr6w6cYuaKcgcN1J/ntaztZ+HoiC+vu5462n3IwPJmxopaMY48zt+gi+Cxf+1agKEOAaukrSi+q6pr5x/rdvLjxAK0dnQgBV80bx51LpzN99EjY9TasfwAOb9aeEDMRLvqlVtvfhyUeFMVOde8oiguOn2jh7+/t5vnP9ndJ9j9dOp0ZY6K7PlhK2PkGrLsPju7Qtukmw5Jfw7zlEKzq3Su+o5K+ovSjobmN1R/s4ckPLDS1aqNxrjp7HD+9dDozx0b3/+TOTtjxMqx7AKq/1rYlTIMlq2DudcNqSUfFf6ikryi9aG7r4NlP9/H4ut2OBcaXzhrNLy+fyexxp1kqoaMdthbBhgehdq+2bfQcMP5em+mrhnEqXqSSvqI4ae/o5KXygzxi3sWhumZAG2efc8UsUqcMcuhlRxt88TxsyIP6Sm3b5Ash/Q8wscffoKJ4hEr6ioJW9fKd7Ud46J2dVBxrBGDW2Ghyrpjl/klV7S2w6UltMZeTtdq2OdfC0t9BQrL7XkdReqGSvjLsbT5g5c9v7GDTXi0BJ8VH8ovLZnDN2eM9u0LVSSt89Ii2qEt7MwSFQMrNcHEOjEz03Osqw5pK+sqwdch6kofe+YqXP9dKKsRHhXHn0uncsCCJsBAvXmStOwjr79e6fmQnhI3UKnqe/2MIU6tlKe6lkr4y7DS2tJO/oYKCDyw0t3USFhzEzRdO4Y5LphET4cOx9Ed2wNo/aGP9AaLHaRd75y1XI30Ut1FJXxk2OjolxWUHePjdXRyzlTm+6uxx/PqKWUyKj/RxdE72fgTv3gOHPtfuT0iBKx6ESQt8G5cSEFTSV4aFTXtruPeV7Xx5uB6A+ZN0/Pbq2aRM9tNiaJ2dsOVFMP8BTlRp2+Zlai3/2Im+jEwZ4lTSVwLa0fpmHnhrp6PffoJuBDnfmMU1Z48bGmWOW07Ah3+Djx+DjhYIGQEX/kzr8w/zo28nypChkr4SkNo6Ovn3x3t5xPw1J1raCQsJ4vaLk7l9STIRoUNwtarafVByL+z4n3Y/ZgKk/xHOWqYmdymnRSV9JeB8tPs4v3t1O7uPngDAOHsM9149h6SEAGgZ7/0I3v41VG3R7k+9CK78CyTO8G1cypChkr4SMA7XneTPr3/JG1sPA9patL+7Zi6XzBrt48jcrLMDPn8WzL/TJncFhcL5P9GqeaohnsoA3JL0hRAxUsp6t0Y2AJX0FbuOTsm/P97LX979isbWDiJCg/jJpdO5dfFUwkOGYFeOqxqrYe3vofwZ7X7sJG2Uz6yrVJeP0qdBJX0hxHXA9UCslPJyIUShlDLLA3H2oJK+ArDtYB13v7yVLZV1AFw+dwz3XjOXCboRPo7Miw5sgjfuOtXlM/0y+EYexE/1bVyKX+or6bta8NtkS/b/tN2Pc19oitK3ptZ2/layi6c/2ktHp2RcbAR/vPYs0ueM8XVo3jcpDbLXw6an4L0/wdfvgmUDXPwrOP9OCAnzdYTKEODq9L86IcStQJyt1W/1YEyKAsC6nUdJ/+v7rP5gD1JKbr5gCiV3XTw8E75dUDAszIYfl8LZWdrwzvf+DAVLoLLM19EpQ4DLffpCiF8ByUCplPJJj0blRHXvDD/HT7Tw+1e38/oW7ULtnHExPLhsHmdP1Pk4Mj9kWQ+v3Wmr3y9g0e1wyT0QPtLHgSm+5o4+fQABSAAp5UtujbAPKukPH1JK3th6mHtf2U5NYysjQoO5K30GN18whZBgVZOmT61N2sItH/8dZAfEJsHVf4PpRl9HpvjQYPv0k7v9PBXwStJXhodjDS3c+8o23tqmlSK4YFoCD153tn/VyvFXYZHaBK6518GrP9Eu9D63TCvgdsUDEDXK1xEqfuSMxukLIZ6QUt7ugXh6UC39wCal5LUth/ndK9uobWojKiyYu6+azXcWJA2N8gn+pqMdPn0c1t2v1e4fEQ9X/xXmftvXkSleNqiWvm3Ujv3sIAC15psyaMcaWvjt/7bx9natdX/htFE8uGweE+NU6/6MBYfABXfCrKvh9Z/Bnveh6Aew41W48mGISvB1hIqPudqn32UgsJRyzwCPz0Ab4WOQUub1st8A6G3HKu7vWKqlH5je3HqYe17eSm1TGyPDQ7jnqtlcnzZJte7dSUoofQrevRfaGiFqNFzzKMy60teRKV5w2i19IcSDnGrdd9+HlHJVH/sMAFJKsxBCL4QwSCnLuz3MJKU0CSFW9rFfCVANzW384bUdFJdpC4Yvnj6KB5edPbwmWXmLEJB2KyQvhVfugH0fwYs3wDk3aDN6R6jRUMNRf907+Wd4zCygxPazBTACjqRu+xZQAdDbtwAlcJXureHna77gQM1JwkOCuOeq2Xxv0WTVuve0+Knw/dfhs39qK3ZtfkGb1PXNx9QIn2Goz6Q/UBdOP3RAjdP97p2IaeD4RmDso/snG8gGSEpKOsMwFH/R1tHJo+av+cf63XRKOGtCDI9kzWfa6GhfhzZ8BAXBeT+C6enwv9uhcpM2wif1Frjsz6pm/zDi0uBnIcStQohSIUS1EGK3EOLrQb5utb1Lx9by70JKWSClTJVSpiYmJg7ypRRfqjh2gmVPfMzf1+1GArcvSeal2y9QCd9XRk2HH76jrcwVFKr1+RcsgcNbfByY4i2uznjJtF0QWC2lnAas7eexVsC+Np0OqO62vwKt2wfbv2kuxqAMIVJK1mw6wFX/7wO2VNYxQTeCF1csIueKWYSFqIlWPhUUDBf+HFashVEz4PhX8ORSbXJXZ6evo1M8zOXaO7Z/q22zc1P6eWwhtpE5tn/NAEII+1Ujc7f9m1yOVhkSTrS08/PCL1j53y00t3Xy7XMn8NbPFrNQr4YL+pVx50D2Bkj9IXS0aou0P7cMGqp8HZniQa4m/RUAUsqH0A4ViWIAACAASURBVGbkZvf1QKduGyNgdRqZs9a23wJY7d06Aw3ZVIaW7YfquOaxD/nfF4cYERrMw5nn8Les+cREhPo6NKU3YZFayYbrn9cmclW8B/84D3a+6evIFA9xdZz+pVLK97wQTw9qnP7QIKXk2U/38ac3vqS1vZNZY6P5+3cMTButCn8NGQ1V8PJtYFmn3U9bAZffByHhvo1LOSN9jdN3taWfIoR4VwjxhBBivptjU4a4upNt/Oi5cn77ynZa2zv5zsIk/nfHBSrhDzXRY+HGl+Cy+7SLvJtWw1OXQc2ZDuRT/NHpLpc4FcgAsqWU0z0WlRPV0vdv2w7WcftzZRyoOcnI8BAeuG4e15wz3tdhKYN1sFwr32DdB+GxcO3fYc43fR2VchoG1dIXQsQIIVagTdhKA37t5viUIai4rJJlT3zMgZqTnDUhhjd+eqFK+IFiggFM72s1fFrqYM334K0caG/1dWTKILnavfMk2tj6y6SUy6WU//VkUIp/a23v5Df/28ovizbT0t7J9WmTKL7tfCYnRPk6NMWdRugg61mtZENQqDaj9+nLoXafryNTBuGMSit7k+re8S9Vdc3c/lwZn++3EhYcxB+vncv1C9Ss6YBXWaZ199Tth4hY+HYBzLzC11Ep/RjshVxF4ZOKaq5+7AM+329lfGwERbedpxL+cDExBUwbYOaV0FwHL2TB+gfVZK4hSCV9ZUBSSp7+cA83PvUZx0+0csG0BF77yYWcM0lVaRxWIuO18fxLfwcIWP8AvPgd7SSgDBmnlfSFEDGeCkTxT63tnax6aSt/fH0HHZ2S25ck88wPF5IwUo3dHpaEgMV3wXeLIUIHu96C1ZfC0Z2+jkxxkaujd64TQqwBimz3Cz0aleIXqk+0cOOTn/HipgOEhwTx2A3nknPFLIKDVCnkYW+6EbLXw5izoHq3Vrtnx6u+jkpxgastfZOUcjlgn6UR56F4FD+xs6qeax//iI17axgTE07Rbeep4ZhKV/FT4ZZ34axl0HpCG9Zp/gN0dvg6MqUfLhdcE0LcCsTZCq5ZPRiT4mMlO46w7B8fU1l7knMmxvLqjy/k7Imq/17pRVgULHtKm8UrguDDv0LhjdBywteRKX1wKenbWvlxQC0Qb7uvBBgpJf/cUEH2f0ppbO3gm+eMp9B0HmNiInwdmuLPhIDzfwzfe1nr5//qTW08v/WAryNTetHfcokOtoJrD3k6GMV32js6+e0r23lh434AfnX5TH60JFktZai4Tr8Ebl2rDec8sk27wHv98zBJLZnhT1TBNYXGlnZWPFPKCxv3Ex4SxBPfNXDHJdNUwldO36hpcKsZpl4EjUfh/66Crap6uj9xtXvnISnlZUAekO6G5RIVP3G0vpmsgk9Y99Ux4iJDeX7FQr4xb5yvw1KGshFxWrXOlJuhowX+ewusu19N5PITrnbvxABZQCbaRVxVcC0AfH2kgR/8axMHrSeZnBDJ/928gKmjVP0cxQ2CQ7XFWRJnwTurYEOuNrTzW0+o+vw+5lLSB1YDhbbWvhIAPrVUk/1MKfXN7cyfpOOp76eqCVeKewkBi26DhGla3Z5t/4WGI3D9c1oxN8Un+uzeEUI84XR3D7BACPGA/eb50BRPeWvrYW56aiP1ze1cNmcML6xYpBK+4jnTjXDzmzByLOz7EJ6+AuoqfR3VsNVfS9/56ku+pwNRvOOFjfu55+WtdEr4/nmTufeauWqGreJ5487WLvA+lwHHvoQnjVoph7Fn+TqyYafPlr6Ucq3T3Wop5R77DfDvesxKr55YX8Gql7SEf1f6DH7/TZXwFS/STYIfvg2TL4CGw1qL37Le11ENO64O2cztdl+1/IcQKSUPvPkluW/vRAj447Vz+enS6WpIpuJ99pE9c78NrQ3wbAZsVqW8vKnfC7lCiGVoo3aMQgi9fbPHo1Lcpr2jk7tf3sqa0kpCggR/WX4O186f4OuwlOEsNAKWPQ0xE+CTv8PL2XCyVrvoq3hcv0nftizif4UQD0op1TDNIaalvYM7X/iCt7dXEREaxBM3pnDJzNG+DktRICgILr8PosfBu/fA2znQbIWLc7RRP4rHuNq9c67zHVVa2f81t3Vg+k8Zb2+vIiYihGdvWagSvuJ/zv8xfPPvWrG29Q/A26vUJC4PczXpdz/1qtLKfuxkawcrnill/VfHiI8K44XsRaROifd1WIrSO8P3IPPfEBwGnz0Br9wBHe2+jipguZr0LUKIXwoh5gshfgVYPBmUcuYaW9q5+f828sHXxxk1MowXVixi7vhYX4elKP2b8034TiGERsHm56Ho+9DW7OuoApKrtXduA+qA2wCr7b7iZxqa2/j+0xv51FLD6OhwXsw+j5ljo30dlqK4JvlSuOkViIiFna/D88uhtdHXUQUcl9fIlVKullLeJqVc7cmAlDNTd7KNm57eSOm+WsbFRlBoOo9po0f6OixFOT2T0uDmt2DkGNizAZ7LhJYGX0cVUFwqwyCEeNC5BIMqw+Bf6k628b2nPuPz/VYm6EZQmH2eKpymDF1j5mqJP3o87PtIG8vfXO/rqAJGfy397mUYCrrd+iSEyBBCGIUQKwd4XL/7lYGdaGnnB//ayJbKOibFj6DQtIikhEhfh6Uog5OQDDe/ATET4cCn8Ox10Fzn66gCQn9JXwohLhVCXApM7eXWKyGEAUBKaQas9vu9PM4IqCV1BqGptZ0f/muTo4X/wopFTIxTCV8JEPF6LfHHJkHlJnjmW3BSLc89WP0l/WSnWw6QYrtdZrvflyxOLZxuAYyDD1PprrlNG5a5cW8NY2MieH7FQpXwlcATN0VL/LrJcKgcnvkmNNX4Oqohrb+Ca6vtN7QROw/Zbr/mVFLvjQ5w/l9J6P4AIYTB9k1AOQMt7R3c9mwZH+2uZtTIcJ5bsZDJCaoPXwlQuiStNHPcVDi8WSX+QXJ19E6cEOJWIcQUIcQKtNb/YKiZQmeoraOTHz//uWPi1fMrFpKcqEbpKAEudqKW+BOmQdVWeHaZ6uM/Q66O078MbRbu74BY4NJ+Hm7lVFLXAdXOO11p5QshsoUQpUKI0mPHjrkS4rDQ0Sm5a81mSnYcIXZEKM/espAZY9Q4fGWYiBkPN716qqvnOTWO/0y4lPSFENehXXQdL6V8GG35xL4UAvaKnHrAbDuGfX00vW10T4bt5x4XeqWUBVLKVCllamJiootvJbBJKfn9q9t5bfMhRoaH8MwPFzBnfIyvw1IU74qdAN9/TavQeeBTeOF6aDvp66iGFFe7d0xSyuWcKr/QZ+0dKWU5OEbnWO33gbW2/cVSymK0bwNqoUwXPWL+mv98uo+wkCBW35TKOZPUR6cMU3GTtcQ/cgzseR8KvwftLb6OashwNenXCSFuBeJtrf5+x03ZWupmKWWB07aUXh6T7HRSUPrw74/38ujarwkS8NgN53Jeco9r44oyvCQka109kQmwuwSKfwgdbb6OakhwNennoLXua4F4W6tf8YJXvjjI71/bDsAD183j8rljfRyRoviJ0bPge/87Vavn5dtUWWYX9LuIipMHgRVSSjUX2os27DrGL9ZsRkrIuWIWWWlJvg5JUfzLuLPhxpfhmWthWzFExsM38tRCLP1wtaVvBfYIIZ5QtXe8Y/MBK7f9p4z2TsmKxVO57WL9wE9SlOFoYgrc8LxWj39jAbz/sK8j8mun09J/0JOBKKccqGniln9v4mRbB9cZJrDqG7PVIuaK0p+pF8Gyp7Q6/Ov+DFEJkPpDX0fll/qrshkjhPinEOIdYJmUco/zzYsxDit1TW384F8bOX6ilQunjSJ32dkEBamErygDmvNNuOqv2s+v3wXb/+fbePxUf907q4E1wHJAqC4dz2tp7yD7P6VUHGtk5pho/nGjgdBgl5c8UBQl9Wa49DeAhJdWgGWDryPyO/1lFCGlfE9KWSelfIjBl15Q+iGlJKd4C5/t0Va9+tfNacREhPo6LEUZehb/EhbeBh2t8OJ34NAXvo7Ir/SX9KsHuK+40V9LdvG/Lw4RFRbM0z9IY7xuhK9DUpShSQi4/AE4KwNaT2jLLloP+Doqv9HfhdwsIUSq7WeBVjLBaPtZSimnezy6YWJN6QEee283wUGCv3/XwFkT1ELmijIoQUHwrSeg8ag2a/e5TLjlHW1M/zDXZ9KXUqpKmF6waW8N97y8FYA/fHMul8wc7eOIFCVAhITB8v/A05fDsS+1cg3fLda2D2PqKqEPVdY2cdt/ymjrkNx8wRRuXDTZ1yEpSmAZoYPvrIGo0dpC66//HKT0dVQ+pZK+jzS2tLPimTKqG1tZPH0U91w529chKUpgipsM3ymE0Ej44tlhP3lLJX0f6OyU/LJoM18ermfqqCj+foOBEDU0U1E8Z4IBlj0JCG3y1pY1vo7IZ1Sm8YFH137NW9uqiI4IYfVNqcRGqqGZiuJxs66CK2yFBV65A/Z/5tt4fEQlfS97c+vhLmWSp41WSx0qitcsug0WmLQx/IXfHZZDOVXS96JdRxr4xZrNANx95WyWqJE6iuJ9l98P+iXQeEybvNXa5OuIvEolfS9paG7jtv+UcbKtg2+fO4FbLpzq65AUZXgKDoGMf0HcVKjaAq/8aFiN6FFJ3wuk1C7cWo43MmtsNPd/e56qmqkovhQZDze8CGHRsP3lYTWiRyV9L8h/38I7248QHRHCP29MYURYsK9DUhRl9KyuI3q+fM3XEXmFSvoe9vHu4+S9vROAvy6fz5RRUT6OSFEUh5lXgPF32s8vmeDIdt/G4wUq6XvQ4bqT/OSFz+mUcMclyaTPGePrkBRF6e6Cn8G85dDWCIU3QnOdryPyKJX0PaSto5M7niunulFbDOWu9Jm+DklRlN4IAdc8CmPmQY0FXr49oBdYV0nfQx5+5yvK91sZFxvBo9fPJ1itfqUo/issErKegfBY+OoN+OgRX0fkMSrpe8C6r46S/76F4CDBYzecS8LIcF+HpCjKQOL1cF2B9vN7fwLLep+G4ykq6bvZkfpmxwSsu9JnkDpFVahWlCFj5hVw0a9AdkLxD6Gu0tcRuZ1K+m7U0Sn52YtfUGPrx7/9YrXCpKIMOUtWQfKl0FQNa74P7S2+jsitVNJ3o8fX7eYTSzWjRobx16xzCFL9+Ioy9AQFw3VPQuwkOFgK79zj64jcSiV9N9m4p4ZHzLsA+FvWfEZHR/g4IkVRzlhUAiz/NwSFwqbVsONVX0fkNirpu4G1qZU7X9TG49++JJnF0xN9HZKiKIM1IQUu+5P286s/htp9vo3HTVTSHyQpJb/53zYO1zVzbpKOu9Jn+DokRVHcZeFtMPNKbcLWf2+BjjZfRzRoKukP0itfHOL1LYeJDAvmkaz5hKoVsBQlcAgB1z4OMROgchO892dfRzRoKkMNwkHrSX77yjYA7r16DpMTVF0dRQk4kfGw7CkQwdqkrd1mX0c0KB5J+kKIDCGEUQixso/92bZbride3xs6OyW/XLOZhuZ2jLPHkJU2ydchKYriKZPPg0tWaT+/ZIKGKt/GMwhuT/pCCAOAlNIMWO33nfYbAbOUsgDQ2+4POU9/tMcxPPPBZao+vqIEvAvvgqkXQ9NxeNk0ZOvzeKKlnwVYbT9bgO5JXe+0zWK7P6TsrKon7+2vAMhddjajVJkFRQl8QcFw3WqIHKWVaPjsn76O6Ix4IunrgBqn+wnOO6WUBbZWPoABKPVADB7T2t7Jz178gtaOTm5YkMTS2apcsqIMG9Fj4JuPaT+bfw9Hdvg0nDPhswu5tm6fEilleS/7soUQpUKI0mPHjvkgur49vm43O6samJwQyW+umu3rcBRF8bZZV4LhJuhogZdWDLkyDZ5I+lbAXmVMB1T38TijlDKvtx22bwOpUsrUxET/mei041A9j6/bDWjdOlHhIT6OSFEUn7j8AW1h9SPbhtwwTk8k/UJO9dPrATOAEEJnf4AQItue8IfKhdy2jk5+VbyZ9k7JTedNZpE+YeAnKYoSmMJHamWYRRB8/Bjs/dDXEbnM7Unf3l1jS+ZWp+6btU7bc4UQFUKIWne/vqfkb6hg+6F6JsaNIOeKWb4OR1EUX5u0ABb/EpDw8m1DZplFj/RPOF2odd6WYvvXDMR54nU9ZdeRBv7fWtWtoyhKNxev1CZrHSqHN391ahEWP6Zm5A6gvaOTXxVtdozWuWDaKF+HpCiKvwgO1YZxhoyALYWw801fRzQglfQH8NSHe9hcWcf42AjuvlJ16yiK0s2oabD0Xu3n138GTTX9P97HVNLvx/7qJv5aotXIv/+6eURHhPo4IkVR/NJCE0xaCCeOwDt3+zqafqmk3wcpJb95ZRst7Z18a/54lswc3evj8vLyyMzMJD09HSEEmZmZZGZmYjb7riiT2WzGZDK5vP1MOR8rLi4Ok8nkeP9Wq7XP7enp6ZSXd52ekZeXh9lsxmKxkJfX60jePuXk5JCZmdnv87q/b3tM7vw8lGEsKFirxhkSAZtfgF3v+Dqivkkp/fqWkpIivcn66qty1yWXyu2zZst18xfJ277/Z3m0vtml5+r1eg9Hd/qKioo8dlznYxsMBsfPJSUlMjs7u8/tRUVFcuXKlV2OZzQaHT/n5ubKsrIyl+IoKytzHGvlypWyoqKiy/7a2lqZnZ0tdTpdl+fYH5eRkSHLyspkfn6+NBqN0mg0SoPBIHNzc116fUXp4sNHpfxdjJQPz5SyqdanoQClspecKrR9/is1NVWWlnqnUkPda69x+Lf3IpubHds6w8KZeN+fiL3mmgGfn5ycTEVFhSdDPG2ZmZkUFRW5/bjp6emUlJQ47qekpFBWVgaA1WolMzOTkpKSPrc7f1a7d+/miy++YO7cuQC0t7dTV1dHQsLAcyHq6uoIDw8nIiKCpqYmOjo6iI6O7vG4w4cPM27cuF63jxkzhqCgU196q6uriY2NJSQkMEdpRUREMHHiREJDVXel23V2wNOXa7X3z71Ra/37iBCiTEqZ2n17YP5Wn6Gjf3ukS8IHCGpt4ejfHnEp6TuzWCzk5uZSU1PDHXfcwQsvvEB+fj55eXno9XoyMjIc3R0mk4mMjIwez8/MzCQ1NRWz2UxJSQl6vTbnzf48vV5Pfn4+ZrOZ/Px84uPjMZlM6HQ68vPzgVNdOvZY8vPzsVgsrF69Gp1OR15eHgaDAaPR2CUe+zHsxzQYDD3eY18KCgpIT0/vd7vBYKC8vByDwcCePXtYsGABkyZNclQr3bFjB7NnD1zmorKykri4OKKioqivr6e+vp6JEyf2eJyUssvxGhsbsVgszJkzh7Fjxzq219fXEx0d3esxAoGUkurqaiorK5k6daqvwwk89m6efy6Gz5+FOd+G6f41/1T16TtpP3z4tLYPZM2aNaxevZqkpKQe+/Ly8hyt5cLCwl6fb0/qubm5jiSel5dHVlYWJSUlpKenU1BQQFFREatWrSI/P79Lcs7NzSU1NZX8/Hx0OseEaLKyslizZg0AJSUlGI3GHvH0dUzQTkj2E5DzNns/eXV1NStXrux3e1ZWluN9h4SEdEn4zlpaWjh27Fivt8GIiopi9uzZ1NfX09JyqnZKVVVVwCZ8ACEECQkJNHdr3ChulDjzVO39138GLSd8G083qqXvJGTs2F4TfEgv3QKuWL58OTqdjpqankO4KioqqKmpoaysDKvVitVq7ZKYAeLjtRJGGRkZjqS/adMmsrOzAa21nJOTQ25uLjk5OVitVsfj+pORkUF6ejpGo9GRvLvHc88993Dfffc5jumc5O0XabvH2ttr97U9IyODnJwcVq1aRWhoaI+Eb+9aCQ8Pp7/6SyEhIXR0dABat9DpdMmEhIQQFxdHbW0tY8eOpb293eXnDmVq7QcvOO8nsO0lqNoC6x+Ay+/zdUQOqqXvZNuV36U5uGs/p4iIYPTPf3ZGx+uexAFHP3ZKSgpZWVnk5+dTUlLS62PtJwuz2exIumlpaY6RQWazmbS0NPR6PUVFReTk5LiU9O2xFRcXO0avdI9n/vz5fR5Tr9f3eiI7XfaTVlRUz2Um7Ql4oJZ+dHQ09fX1ADQ0NPTan9+dc+xNTU2EhYUB8Nhjj5Genk5xcTEFBQU9RvaYzWasVivl5eUkJydjMpmwWCyO/fYuOVcUFxe79Dg7++sqQ0RwCFzzqFab59Mn4PBmX0fkoJK+zdH6Zn7TMIFH52fQnjgGhCBk/HjG/emPp92f351er6e0tLRLksjOzqawsJD09PQ+E4U9idhb8wArV66ksLCQlJQUSkpKWLlypaNrJjc3l6ysrB7HSU9P75KcAEeCt3fddI+nv2PqdLpeW/uny2QysWbNGiIiIvp8jL2l39sNcJwwKioqCA4OdtzftWuX4xj79u2jpaXF8W9MTAz79u1znIDt36jmzp3LxRdfTEZGBtnZ2eh0OgoKtGn19ver0+kwGAwYDAZMJlOXb0B6vZ7Vq1e79N77O2n29tnqdLoe/4eKn5tggAUmkB3w2p3aRV5/0NuQHn+6eWvI5s9f/FxOznld3vJ/m7zyegOpqKhwDHv0R85DLAdrx44dXe63tbXJ3bt3u+34rsrNzZX5+fmO+/bhpfZ9zuxDPZ1VVFS4PNTU+XW6KyoqkrW1PYf7VVRUuG0IbvfPXPGQ5nop/zJbG8b5yRNefWn6GLKpWvpA2b4aXvr8IGEhQdx79RxfhzMkmEym0+6icFVVVVWXETXesmnTJlJTtRFuVqsVi8XiGFXl6lDcnJwcx8/2CWf2bwuDpdfruwyTVYaA8Gi48iHt5/f+BHWVvo0HdSGXjk7Jva9sB8B0kZ6khEgfR6Sxj9zxVxkZGZhMJoxGY6/XI87UlF+/4bZjOdv74FUDPqa8vByLxUJNTQ1Wq/W0E6xer3d8Fjk5OY6L5c4ngsFyx7UUxctmXQWzroadr8NbOXD9cz4NZ9gn/Rc37Wf7oXrGx0bwoyXTzugYQogu4+xzc3N7HdJoH37pLWazmaKioh4nj762ny77800mk+Pnvj6LuLg4li9f7kha9nkCAAcOHCAyMtLjo2fsrfWQkBAmT54MaOP1q6qqHBeDMzIy2LdvH+3t7ezbt8/xuIEUFxd3ed8Wi6XLRd+BWCwWxwX6srIyLBaL4/Oxj9aCU9cflCHmyofAskFL/Dvf0E4EPjKsk761qZWH3/kKgHuumsOIsOAzOo599Iy/MRqNGI3axBDnpOS8fbCKi4u7TMTq67Nw/uZisVhYunSpY6aucxLecncix44dIzk5maqqKqKjo3sd3dNdY2MjtbW1TJw4kcrKShITEwkPD3fsr6mpISYmhsRE7fj19fWOC7pz5szBbDZjNBppbGxk7NixhIeHU1FRQWNjo0uv3/0iq31UlcFg6NEA6I1er3ck9+LiYrd/g1J8LGY8LP0tvLUS3vo16C+BMN/0KgzrPv2/vLuL2qY2ztMncOU87/che1NfE8AGKz8/v8ds4oHo9XqMRmOv1wRiYmIc4+7j4uI4fvy4S8dsaGggJiYGgMjISEfL3S48PNxxXPtkrJqaGuLj47vMno6KinKcLFpaWrqcOOzKy8spLy+nsLCQ4uJix8Qz+3aLxcLKlSsxm82YzWa8VUZE8XOpt8CYeVC3Hz56xGdhDNuW/ldVDTz32T6CgwR/uHauWyesWK1WVqxY4ag14/z1vHvJBIPBELDlGPpjn/3b/f3W19czYcIEQEvUjY2NLh2vvb2d4GDtm1pISAhNTU1d9kdFRVFbW8vWrVuJi4sjJiaGyspKOjo6kFLy+OOPM2PGDOBUiYbExETHZK+UlBTHsQwGQ5cLu87vwXm7ffaxu1gsll7LWyhDRHAIXPWwVpvnw0fgnBsg3vulMIZtS//+N7+kU8J3FyYxY8zAE3r6Y0/M9jH1Op2OoqIiSkpKevSbdy9vEIjlGJw/C1e0t7dTUVHBjh07aG9v77U7ZbDlGOwnj3nz5jmOBxAcHExycjIxMTGObwe9lWhYvny5W0cr9dc3n5GR0WvXTnl5+Wl/q1L8TNIiOPt66GiBt1f5JIRh2dJ/f9cxNuw6RnR4CHcunT7o4/XWj52Xl0d1dXWPvt6cnJwuJRMCrRyDq9c3ysvLSUtLA7SWeXJysqOF7ZwQ3VWOoba21tH9ExMTQ21tLeHh4Y7um7CwMJqamhyP6V6iQafTOSaluaOv/XSTd28nWGWISv+DdjF311ta3f0Zl3v15Ydd0u/olNz/5pcA3HHpNBJG9uyzHSx7Jc2VK1f2aB3ak6K9KyYlJYX4+Ph+k0B/5RgyMjJ6lGOwH9uVBUJ6K8fQPR7nYzqPPjrTcgxWq5XCwkLHhVy7qKgoIiMjHX3t0LUcQ/d+ervExESio6Mdib2hoYFRo7quZRweHu5I6vX19URGRhIVFcWxY8eIiYmhtbWVsLCwLq/d1NTUpayDuy5+nwmV8ANI9FhY8mt49x5tCOfUiyG071np7jbskn5x2QF2VjUwQTeCH5w/xSOvYTQaWbFiRa/jvPPy8hzbc3NzHX369q6X3lrJ9i4Ti8XC2rVrAa2/ODMzkwceeMCR7Lsfu7v09PQe3wCysrLIyclx9D9nZ2d3iSctLa3PY55OOQZ7tU37SaKvbwOTJ09m165dPbo/Bmrp2/vsKyoqCA8P71KOYcaMGSQmJlJRUcGuXbsIDw/vcvyKigpCQkIcRdfsQzZDQkLUEEnFMxaaoPwZOP4VfPJ3uOiXXnvpYbWISmNLO5c8vJ6jDS08ev18rp0/wS3H9ST7yBJ/najVfTGVM/Hll1/2WjvfnoBdGeeunJ6+PnPFiyzr4ZlrIWQE/KQUYt1b0ruvRVSG1YXcgvctHG1o4ZxJOr55znhfhxMQArEcg6J4hX4JzLkW2k/C2j967WWHTdI/Ut9MwfvaRdXfXDV7yNQUHwrlGEpKStxSddOZvS/flYlRijJkpf8RgsNgSyEcLBv48W4wbJL+30p2cbKtgyvmjiVtiuqnKZkcBwAACG5JREFUdafuQ0HdISQkxOsrWBUUFJCcnOxSPf2UlBTHkFT7fAxXLpybzWbHWHtP1t+3x6Vq8Pu5uCmw6Hbt53fuAS90tw+LpG85doKiskqCgwQrr5jp63AUP5WamorRaHSpnn5ubi5ms9lRIyc3N9elukrO5RXcVX+/L6oG/xCx+BcQmQD7P4EvX/X4yw2LpP+Xkl10dEoyUyaiTxzp63AUP2U2m7vMvLVarY7ROwUFBT2GbNrnRZwpi8Xi8aRsMBg8ds1FcZOIWLjkbu3nknuhvaX/xw9SwA/Z3Hawjje2HCYsJIg7jYOfiKV42O9jPXTcugEfsmnTJlat0mZJulJP32g0kpubS3l5eZfurYKCAlJTUyktLe1SgqM3OTk5jtFP9lIYFotlwOe5Sq/Xk5ubq2by+jvDD+CzAm0I58YCOP8nHnupgG/pP/yuVkXzpkWTGRc7wsfRKP7MXizNXijNlaGo+fn5rFixwnHfPjHPXl2zvwVUutfft9c+cnXBFlepGvxDQHDIqcXTNzwEjdUeeymPtPSFEBmAFTBIKfNOd7+7bNxTw/qvjhEVFsyPLjmzWvmKl7nQIvcEe5/9mVYMtffpb9q0yXEM+8irsrIyampqSEtL67MI2+nW37fr7aRij8lOTTAbIqanQ/JSqFgLGx48teKWm7k96QshDABSSrMQQi+EMEgpy13d7y5SSvLe3gnArYv1xEeFufsllABSWlp6xmUWVq1axdSpU8nNzSUtLc1RJ8disfSb6J2dbv19O3d1Ayl+4rI/wz/XQenTsPA2SHD/xERPdO9kobXiASxA97+kgfa7xfqvjlG6r5a4yFBuXez98qXK0OFcT98VZrOZnJwcx0pX9tE7oJXHKC8vx2w2U15e3iPh22vuO99U/X3FYcwcmP8d6GzX1tT1AE907+gA57+ehNPcP2idnZI824pYd1wyjeiIUHe/hBJAXFlw3HlUj9Fo7FEszrnFbU/0vX1z6F6L35P190HV4B+SltwNW4th+8vaBd0JKQM/5zT45YVcIUS2EKJUCFHqar10Z/tqmjha38zYmAhuXOTaGqeK0h9319M/XWfaL69q8A9BsRO0rp2gEDjo/sl1nmjpWwH7b6gO6H4ZeqD9SCkLgALQCq6dbgBTR0WxYeUl7D3eSEToma17qyjO3F1P/3SdSeJWNfiHsAt/DoabPNKn74mkXwjYK7vpATOAEEInpbT2td/dRoaHcNYED435VoYlX9bTPxMq4Q9hI3TazQPc3r1jH4kjhDACVqeROWsH2K8MY/5e4juQqM96ePPIOH1b90z3bSn97VeGr4iICKqrq0lISBgy1U+HKikl1dXVRER4b6Umxb8EfBkGxf9NnDiRyspKlxc5VwYnIiLC6xVMFf+hkr7ic6GhoUydquZSKIo3+OWQTUVRFMUzVNJXFEUZRlTSVxRFGUaEvw/fEkIcA/ad4dNHAcfdGM5QoN7z8KDec+Ab7PudLKVM7L7R75P+YAghSqWUqQM/MnCo9zw8qPcc+Dz1flX3jqIoyjCikr6iKMowEuhJfzjO/FXveXhQ7znweeT9BnSfvqIoitJVoLf0lQAihMgQQhiFEP2uNjLQfkXxR/alZPvY59LvvisCJukP9KG480PzFy6852zbLdfbsbmb89rKgLWvPxBb9dY0b8bmSS78HxtsjwmYlVJO4285YBYItv3eru5jn0u/+64KiKQ/0Ifi7g/NH7jwno2A2VbRVG+7P5R5ZW1lf+Li761JSlmM9n88HH6vDYDFtt8SCO8ZHO+3r0Wa3fq7HxBJHz9ZjN3LBnpPeqdtFtv9oWzAtZWFEAbbH0+g6Pf/2Na6rwCQUuYFyNoUrvyt2r+56gPkPQ/EreuKB0rS9/li7D7Q73uSUhY4rVtgAEq9FZgPndlCsv5roN/bNCDB1sUTKN2WA/1el6O18Cvou2Ws9CNQkr7SB9vX35IAaBH1u7ZyALbyXVXttBpdwPTr90UIoUP7XcgHVgshhvo3WFcMuK746QiUpD/oxdiHIFffk1FKmeedkDyqkFNdVF3WXrZvc7qgGRD92wz8f1yB1gWC7d9AuIA90HvOBh6w/U5nAgF7onP63e71d/9MBUrSHyghuPVD8xMDvWeEENn2hD/UL+S6sPZyse2CZjxasggEA/0fm7vt3+TV6DxjwN9rO/vFXu+F5jm2xkpqt29rHllXPGAmZ9mGb1nQLu4U2LaV2dfm7W3/UNffe7b9ghSh9XvGA5nDtPtjSHPx97rGtj8QvtG58p5X2vbHB8rfsjcFTNJXFEVRBhYo3TuKoiiKC1TSVxRFGUZU0lcURRlGVNJXFBtbnaIK29DPbCFEfrf9RiGEzjYZqkIIke88TlwIoRf/v727u2kYBqMwfCwxADOEDSLYoGyA1AlgBNiBDVgBugmQDWCGbPBxkc84qlLq0kj58ftctVUT5erEcZyTEN4O7PtyJctIsXCEPpB8qOsr2vmqkDaWesUlg2YWl8w1kl7MLK6Tl3++H9qxmbVafhUGVoDQB5KNpM/e934lwMOxJa8+6v8r2JsSnprFvBH6QHIj7yjykX3lD3xJ0lXmPn5rrEMIj/0KYL8SuB3xeIGTEfpAUivVUG/M7KSA9lBvJcnfYdD41UH/hLG2UjgszMXUBwDMQW/Ofnfsvwe2v9vbtpIUb95+9X6nGRKTIvSBzrXO62Tan8t/V/eyjyaE8D20ATAFQh/F8xuwT8os7/LRey1p69tuld7iVIcQKjN79jn9OJ1D7xFmge4dIIM3lp5V7uUniPq/U0jAGLiRC+R5HWG5JYGPyRH6QAZ/uKod6nXP4aN85vYxOaZ3AKAgjPQBoCCEPgAUhNAHgIIQ+gBQEEIfAApC6ANAQX4Aa65Z/wYd5t8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's test different hypotheses\n",
    "vPrH = np.arange(0.00001, 1.0, 0.0001)\n",
    "vfunc = np.vectorize(posterior)\n",
    "\n",
    "# posterior of true positives / positive predictive value\n",
    "vPrHE = vfunc(PrEH, vPrH, PrEnH)\n",
    "# posterior of false positives\n",
    "vPrnHE = 1 - vPrHE \n",
    "\n",
    "# posterior of true negatives / false predictive value\n",
    "# 1-PrEnH = specificity \n",
    "# 1-PrH = prevalence of negative class\n",
    "# 1-PrEH = 1-sensitivity\n",
    "vPrnHnE = vfunc( (1-PrEnH), (1-vPrH), (1-PrEH) )\n",
    "# posterior of false negatives\n",
    "vPrHnE = 1 - vPrnHnE \n",
    "\n",
    "# Plot the two posteriors\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(vPrH, vPrHE, '-', lw=2, color='tab:blue')\n",
    "ax.plot(vPrH, vPrnHnE, '-', lw=2, color='tab:orange')\n",
    "plt.xlabel(r'$\\displaystyle P(\\textrm{Lie})$')       \n",
    "plt.ylabel(r\"Predictive value\") \n",
    "plt.legend([r\"$\\displaystyle P(\\textrm{Lie}\\mid +)$\", r\"$\\displaystyle P(\\textrm{No-lie}\\mid -)$\"])\n",
    "\n",
    "#plt.title(\"Bayesian evaluation of performance\")\n",
    "\n",
    "i = 500\n",
    "\n",
    "ax.plot(vPrH[i], vPrHE[i], 'o', lw=2, color='tab:red')\n",
    "ax.plot(vPrH[i], vPrnHnE[i], 'o', lw=2, color='tab:red')\n",
    "\n",
    "ax.annotate('True positives (PPV)={:.3f} \\n False positives (FDR)={:.3f}'.format(vPrHE[i],vPrnHE[i]),\n",
    "            xy=(vPrH[i], vPrHE[i]),\n",
    "            xycoords='data',\n",
    "            xytext=(vPrH[i] + 0.07, vPrHE[i] + 0.05),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "ax.annotate('True negatives (NPV)={:.3f} \\n False negatives (FOR)={:.3f}'.format(vPrnHnE[i], vPrHnE[i]),\n",
    "            xy=(vPrH[i], vPrnHnE[i]),\n",
    "            xycoords='data',\n",
    "            xytext=(vPrH[i] + 0.07, vPrnHnE[i] - 0.10),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "vPrH[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate the distribution of each classification outcome:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEFCAYAAAAPCDf9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eVRjWX7n+bla2BchsUWwi4x9FyLTVXZtWZDptj1eIcNdaZ/u6nZCzbE9Hnu6gwpn95x2z0mXwd1lz9hzbEi7vfSk7cig5ozdtquyUFbW4kpXZkjEkrFHIvYINgmB2EG688eThFiDiEAI0P2c806g956efhDwe/f9fvd+v0JKiUKhUCgSA128A1AoFArFzqGSvkKhUCQQKukrFApFAqGSvkKhUCQQKukrFApFAqGSvkKhUCQQhngH8Dhyc3NleXl5vMNQKBSKPYXL5RqTUuat3r/rk355eTlOpzPeYSgUCsWeQgjRu95+Vd5RKBSKBEIlfYVCoUggVNJXKBSKBEIlfYVCoUggVNJXKBSKBEIlfYVCoUgg9m3Svz56nXveeyjpaIVCoVhm18/TfxomFya5P34fgIGpAewFdrKTs+MclUKhUMSffTnSz0rK4kcO/giphlS8c14cvQ5ueW4RlMF4h6ZQKBRxZV8mfYADGQd4uexlKrMrCRLktuc2jl4HnllPvENTKBSKuLFvkz6AUW/EVmDjs8WfJcOYwcTCBO/1v8f10essBZfiHZ5CoVDsOPs66YfJS8ujtqyWIzlHALg/fp+O3g5GZkbiHJlCoVCsZXhyjm/eGmJuMbDt106IpA9g0Bk4nXeaF0tfJDspm6nFKb4z8B2cQ04WA4vxDk+hUCiYWwzwA7eHd++MMDa1wL0h/7Z/xr6cvbMZ5hQzNWU13PXe5Y7nDt2T3TyafkRVQRUHMw7GOzyFQpGg9IxN4+odZ34piF4HJw5mc/xA1rZ/TsKM9KPRCR3HLcepKavBnGJmLjDH9x9+nx88+gF/8d//gvLycnQ6HeXl5bz11lvxDlehUOxjpuaXeO/uCO93eZhfClKQlcy/OHWAk0XZ6HRi2z8v4Ub60WQnZ/NiyYt87PuYm2M3+au/+iva/kMb83PzAPT29tLQ0ADAq6++Gs9QFQrFPiMYlNwb9vPRwARLQUmSQce5UhOVeRkx/Vyx21es2u12uRMmKtOL01RWVDI8OLzmWFlZGT09PTGPQaFQJAbe6QU+7Pbgndb6ieWWNGxlOaQY9dv2GUIIl5TSvnp/Qo/0o0k3pjPycP3ZPH19fTscjUKh2I8sBoLcGJjg/rAfKSE9WU91uZmDptQdiyEmSV8IUQf4AJuUsmXVMRvgAtyhXQ4pZWMs4nhSSktL6e1d6zCWfzAf/4KfzKTMOESlUCj2A4O+WZw9XqbnAwgBRw9kcrooG4N+Z1ur2/5poaSOlNIB+MKvozBLKYWUshKoB5q3O4an5Y033iAtLW3FvuSUZF759Vfo6O3gnveeknJQKBRPxNxigO9/PMZ37o0yPR/AnG7k5ROF2EpzdjzhQ2xm75xHG+WDNpqviT4YuhmEsUop3ewSXn31Vdra2igrK0MIQVlZGa1trfziq79IQAa4MXaDb/V9C9+c7/EXUygUCc/HI1P8j+sP6fXMYNAJzpWaeOl4Ieb0pLjFFIvyjgnwRr22rHeSEKIGcKx3LJ68+uqr687UKckswTXsYnx+nHf73uWo+ShHzUfR67av8aJQKPYHE7OLXOn2MuLXZgIeMKVQXW4mIzn+bdR4ztOvlVKuO2QWQjQIIZxCCOfo6OhOx7UuhemFvFT+Es+ZntME3Ly36ejrYGx2LN6hKRSKXUIwKLk5OMHXP3rEiH+eFKOOH37OwueO5O+KhA+xGen7AHPoaxOwkazl6lp/BCllG9AG2pTNbY3uGTDqjJzLP0dJZgnOISf+BT/v9b/HIdMhTuaexKDbHf+pCoVi5xnxz3Gle5yJWW0apjUvnXOlJpINu6saEIssdQkIzw21EirhCCFM4ZG9EMLKct1/z5GbmkttWS23Pbe5N36PB74HPJx6SFVBFQXpBfEOT6FQ7CALS0GuD/h4MDwFQGaKgecrzBRkpcQ5svXZ9qQvpewUQthDNXuflLIzdOhdoCrqVO/ad+8d9Do9p/JOUZxZjHPYiW/ex3cHv0t5Vjln8s6QpI9fo0ahUOwM/d4ZnL1eZheC6AQcP5jFiYPZ6GMgn7BdqBW520BQBrk/fp/bntsEZIAUfQq2AhtFGUXxDk2hUMSAmYUlnD3jDIzPApCbkcTzFWZMabtnsKdW5MYQndBx1HyUgxkHcQ458cx5eP/h+xRnFHMu/xwpht35mKdQKJ4MKSUPRqa41u9jKSAx6AXnSkw8l5+BELt3dB+NSvrbSFZSFp8r+Rxdvi4+GvuIgakBRmZGOJN3hvLs8niHp1AongHfzAIfdHvxTC0AUGJOpaosh7SkvZVG91a0ewAhBM/lPKeN+oedDM8Mc2X4Cv3+fqoKqkgzpj3+IgqFYtewFAhy8+Ekdx9NEpSQlqSnqiyHEvPe/FtWST9GpBnT+HTxp+mZ6OH66HWGZoZ4p+cdTuWeotJUuWceBRWKRGZoYo4Pe7xMzWme2ocLMjhdbCLJsHetSFTSjzHl2eUUpBdwbeQaA1MDXB29qo36C6vIStp+VxyFQvHszC0GuNrno3tsGgBTmpHqcjN5mclxjuzZUUl/B0g1pPKJg59gwD/A1ZGrjM2N4eh1cNxynMM5h9GJvTtqUCj2GxvZFsbCxSoeqKS/gxRnFpOfls/10ev0TPZozV7/AFUFVeSk5MQ7PIUiofHPLeLsGefRxBwABVnJVFeYyUoxxjmy7UUl/R0mSZ9EdWE1JZkldA53rhBwO2Y+pgTcFIodJhiU3B3yc3Nw2bbQVmrCGmPbwnihkn6cKEwvpLa8lltjt3jge8Ad7x0GpgawF9jJTc2Nd3gKRULgmZrnw24v4zOxsy3cbaikH0eMOiNn889GpBzCAm7PmZ7jZO5JjLr99VipUOwWdoNtYbxQHcRnpKWlhfr6empraxFCUF9fT319PQ7H1q0CclNzqS2t5Zj5GDp0fOz7mG/2fJOh6aGYxOxwOGhsjJ1DZVNTE/X19bS0tDz+5ChWxxT+WcYyVkXiMeib5R8/esS9IT+g2Rb++KkDCZHwQWnvbCuVlZV0dXU90zV8cz6cw07G58cBtk3Arb29nbq6ume6xlbo7Ozk0qVLNDc309TURGNjI1arNXK8ra2Ny5cvA+D1ejl//jwNDQ00NTXx9ttvMz4+HonX6/XS0NBAW1sbVquVmpqadT9TodgKswsBXL3j9HlnADCnG3m+whJXF6tYspH2jhrp7zJMKSZeLH2R07mn0Qs9PZM9vNPzDgP+gWe67qVLl7Ypws1xOBzU1tYCUF1dveaJp6GhgY6ODjo6OrDb7dTV1WEymWhtbV1xc7Barfh8mvr2s95IFYmNlJKPR6b4+xsP6fNqtoW2svjbFsYLlfRjhNvtprGxkfr6er797W9HShQtLS20t7cDRMpC4ddhdEKH0WfkD7/0h/zlb/0lv/ZTv4bjnoP3H77Pz9b97Ir3hK/R2NhIS0sLPp8vsq+trQ3Qyi3hko7P58PtdkdKMOHE2tLSEknQ0XE5HI5IiaWzs5PH4fF4MJs1Dx2TybRhwnY4HJhMphWJPhqbzYbH46GyshJAjfIVT8XE7CLv3hnhw24viwHJAVMKP3b6AEcL98+8+ydFNXJjyNtvv013dzde71rrgJaWFmpra2loaKC+vn7d0otep+fv/p+/40v/65foudmD45KDwrOF/Ndf/a9c+KULmEwmqquruXDhAvX19bS2tgJEyidVVVU0NDTQ3NxMZ2dn5Hg4nvPnz/P2229HRt8XLlxYE5fZbObixYvYbJrRmdvt3rBf0dDQsOWfTXNzMx0dHRseD99gurq6aGpqwu12b3iDUChWEwhK7jya5ObgBEEJKUYdVWU5lFnS4x1a3FFJP4a88sormEymdZN+V1cXXq8Xl8uFz+fD5/NhMplWnGO1WhFCUFlUSVFKEZ2POun6qIuGmw2MD4+ziDbNzOfzRUbXoN1QPB4Pbrd70/jq6uqora2lpqYmklBXx/X666/zxhtv4PP5IiWYzZK7xWKJfL8+nw+LxbLmnPDTxWZcunQpUiYKP3VcuHDhse9TKEb8c3zY7WVyVtPLqcxL5+wutC2MFyrpx5DVSRy0pGq1WqmqqsJsNm+5uZpsSOYnPv0TBFICFP1wEfOBeW7dvsXf/uPfcuXKFZqbmwEt4VutVi5cuLCmbLRRjO3t7ZHy03pxXb58GYfDQWtrK42NjZuO9Gtqarh06RI1NTV0dHSsO/PG4XBEnhw2orKyks7Ozsh1qqurH/u9KBKbhaUg1/p9fDyybFv4QoWZ/F1qWxgvYpL0hRB1aB64Ninlmnl7Qggbmn8uUsrHZ6Y9jtVqxel00tjYGBl9h8snra2tmEymSElmM8Lv8f6lF5kqefnfvEyXu4tRzyi/8e9/g6/+7lepqanhtddeW7d0UltbGynxhDl//jxNTU2RUfTquKqrqyPXam5ufuxI32azcenSJerr67FarZHkXltbG7mO2+2O1OrDhH82jY2NNDU1ReKora3FarXuyMwjxd6lzzODq29v2RbGi22fshlO6FLKdiFEA+CM8skNn9MqpWwUQlwAHKuPR7OXpmzuND/xMz/BT/0vP4WpyMTN929y5707/M2f/40ScFMkDNPzSzh7xxmMsi18ocJCdppa2LiTdonngfAw0w3UAJGkHnoK6AJY7ylAsXX+83/8z/zRH/8RnlkP/kU/NedrcPQ6qC6sVgJuin2NlJL7w1NcH9BsC416wdk9ZlsYL2KR9E1AdOdydSevGiJPBDUq8T89NpuNN9veBGB4ehjXsIuJhQne7XuXIzlHOGY5hkGn2jaK/cV6toX2MjOpSapRuxXiVQfwhEs6oZH/CoQQDUIIpxDCOTo6uvPR7UEK0gt4qfwlDuccBuDu+F0cvQ5GZ9TPT7E/WApojdpv3BzCM7VAWpKeTx3K5VOH8lTCfwJikfR9QHj+oAnwrDrehVb2IfTvmmkZUso2KaVdSmnPy8uLQYj7E4POwJm8M3yu5HNkJWXhX/Tz7YFv0zncyWJgMd7hKRRPzdDEHP94c4jbDzWf2sMFGfzYqQN71qc2nsTi2f8SEG4eWAEHgBDCJKX0hV7XRR2/EoMYEhpLqoWa0hrujt/lrucuXRNdPJp+hC3fxoGMA/EOT6HYMuvZFj5fYSY3Y+/bFsaLbU/6UspOIYRdCFED+KJm5rwLVEkp3UIIX7iskwhTNuOBXqfnhOUExRnFOIeceOe9/NPDf6Iss4wz+WdI1qs/GsXupntsms4o28KTRdkcS2D5hO1CqWwmAEEZ5MH4A255bhGQAZL1yZzLP0dJZkm8Q1Mo1uCfW+RKj5ehiXlg/9oWxpqdnLKp2GXohI4j5iMUZRThGnYxMjvCDx79gP7Jfs4VnCPVkBg64ordTaLZFsYLlfQTiIykDD5d/Gm6J7q5MXaDwelBRnpGOJ13Gmu2EjNTxI81toW5adhK97dtYbxQST/BEEJgNVkpTC+kc6STR9OPcA276J/sp6qgiowkNapS7ByabaGPe0OaXk56sp7nK8wcyFZPn7FCJf0EJc2Yxo8U/Qh9k31cG73GyOwIHb0dnMg9wSHTIbWqURFzBn2zOHu8TM8HEAKOFmZyqigbg17JiMQSlfQTnNKsUgrSC7g2co0+fx/XR6/T7+/HXmAnOzk73uEp9iFrbQuTeKHCTE4CuljFA5X0FSTrk3nhwAuUZpbiGnHhnfPi6HVwzHKMo+ajSsBNsS1IKekaneJqn4/FgMSgE5wuyeZIQaZ6stxBVNJXRDiQcYCXU1/mxtgN3BNubnluMeAfwF5ox5xifvwFFIoNmJhd5MNuL6N+bRrmQVMK1eVm0pNVCtpp1E9csQKj3khVQRUlmSURAbdv9X2LQzmHOGE5oQTcFE9EICi5/XCSWw+XbQvtZWZKLUo+IV6ov2DFuuSn5VNbVsstzy0ejD/g/vh9Hk49pKqgivy0/HiHp9gDKNvC3YlK+ooNCQu4lWSW4BxyMrEwwXcGvoM128rp3NMY9WqFpGItq20Ls1INPF+ubAt3CyrpKx6LOcVMTVkNd713ueO5g3vCHRFwO5hxMN7hKXYRq20LTxzM5vjBLGVbuItQSV+xJXRCx3HLcYoyinAOO/HOefn+w+9TmlnK2fyzSsAtwZmeX+JKj5eHvjkA8jKTeb7crGwLdyEq6SueiOzkbF4seZEHvgfcGrtFn7+P4ZlhzuadpTSrNN7hKXaYiG1hv4+loGZbeK7URGWesi3craikr3hihBAczjnMwfSDdI50MjwzzAdDH9Dn78OWbyPNqGZmJALj05ptoXdasy0sNadRVZajXKx2OSrpK56asICbe8LNjdEbPJp+xDs973Am7wwV2RVqpLdPWQoE+WhwgrtDfqSEtCQ99vIcinPUzX4voJK+4pmxZls5kH6Aq8NXGZwexDXios/fR1VBFZlJmfEOT7GNPJqY5UrPOFNz2jTMI4UZnC42YVR6OXsGlfQV20KqIZVPFn2Sfn8/10auMTo7qgm4WU5wKOeQknLY48wtBujsG6dnTNPLUbaFe5eY/CUKIeqEEDVCiAsbHG8O/dsQi89XxI+SzBJeKn+JsswyAjLAjbEbvNf3HhPzE/EOTfGUuEen+Icbj+gZm0GvgzMl2fzoiUKV8Pco2570hRA2ACmlA/CFX6+iQQjRBbi3+/MV8SdZn8zzB57nU0WfIs2QhndeE3C7NXaLQDAQ7/AUW8Q/t8i37g7zA7eX+aUghdnJ/NipA5w4mK18avcwsSjvnAc6Ql+7gRqgc9U59aGbgmIfU5heyEtlL/HR2Ed0TXRx23ubgakB7AV2LKmWeIen2IBgUHJnaJKbgxMEgpBs0HFO2RbuG2KR9E2AN+r1en/dttDMDpuUsiUGMSh2CUa9EVuBLSLgNrkwybf6v8Uh0yFO5p5UAm67jLGQbaFP2RbuW+LyFxdO9EKIWiFEzepRf6jW3wBQWqoW/OwH8tLyqCmr4Y7nDvfG7/HA9yAi4FaQXhDv8BKe1baFGSkGqstzlG3hPiQWSd8HhMXXTYAn+qAQog5AStkeOrbGkVtK2Qa0Adjtdrn6+OLiIgMDA8zNzW1v5LuAlJQUiouLMRr33/J1g87AqbxTFGcW4xx24pv38d3B71KRVcGZvDNKwC1ODIzP4OwZZ2ZBsy08diCLkwezlG3hPiUWSf8SYA99bQUcAEIIk5TSh1bnDzdwK4HWJ/2AgYEBMjMzKS8v31cLgKSUeDweBgYGqKioiHc4MSMnJYfPl36ee9573PbcpnuyWxNwK7BRlFEU7/ASBmVbmJhse9KXUnYKIexCiBrAJ6UMN3HfBapCxxuEEF6gK+r4lpmbm9t3CR80eQOLxcLo6Gi8Q4k5OqHjmOUYRZlFOIeceOY8vP/wfYozijmXf44Ug5LhjRXr2RaeKTFxuEDp5SQCManph8ozq/dVbXb8Sdmvv5z79fvaiKykLD5X8jk+9n3MzbGbDEwNMDIzwtn8s5RllcU7vH3HxMwiH/Yo28JERv1PPwUtLS1cuXIFn8+Hw+Ggrq4OgMbGRmpqasjJyeGVV17B69UmMb355puYTKYN9yc6QggO5RziYMZBXMMuhmeG+XDoQ/omNSkHJeD27ASCklsPJ7j9cFLZFiY4Qso1fdJdhd1ul06nc8W+O3fucOzYsThFtJLKykq6urpW7KuqqsLlcgHgcDi4fPkyra2ta/a3trZy+fLlNdfcTd9fPOiZ6OHa6DUWg4sYhIHTeaexZlsT7ilouxiZnOPDnmXbwufyMzhTkq1sC/c5QgiXlNK+er9qz8cYu92O27124XFNTQ0+ny8OEe1+yrPLebn8ZYoyiliSS3SOdPLt/m/jX/DHO7Q9xfxSgA/cHhx3RpicXSIr1UDN8XyerzCrhJ/A7Pnyzl990BeT637hhe1ZH9DW1kZtbe2a/Q6Hg+bm5m35jP1IqiGVTx78JAP+Aa6OXGVsboyO3g6OW45zOOewEnB7DL2eaVy948wtKttCxUr2fNLfjbjdbhobG/F6vVit1khyd7vd1NfX43a7uXjxIjU1NXGOdPdTnFlMflo+10ev0zPZw0djHzHg16QcTCmqH7KadW0LK8xkp6o1EAqNPZ/0t2tEvp2YzWZaW9cuP7BarVy+fJnOzk7q6+sjDWDF5iTpk6gurI5IOYzPj+Poc3DUfJRj5mPodapUIaXk3rCfG/0TyrZQsSnqGTkO2Gw2bDYb7e3t8Q5lT1GYXshL5S9xyHQIieSO9w4dfR2MzY7FO7S44p1e4J1bw3T2aj61peY0fuL0QZ7Lz1QJX7GGPT/S36u8+eabfP7zn1ej/SfEqDNyNv9sRMrBv+Dnvf73OGQ6xIncExh1iVPGWG1bmJ6sp6pM2RYqNkcl/Wdk9XTNjfYBkemaACaTacVrxZORm5pLbWktd7x3uOu9ywPfAwanBqkqqKIwvTDe4cWcRxOzfNjtZXpe8ydQtoWKrfJESV8IkSWlnIxVMArFk6DX6TmZe5LiDG3UPz4/zvcGv0d5Vjln8s6QpN9/GjJziwE6e8fp8Wh6OTlpRqqVbaHiCdhS0hdC/Czw80A28LIQ4pKU8nxMI1MotogpxcSLpS9yf/w+tz236ZnsYWh6iHP55yjOLI53eNuGe3SKzj4fC0tB9Do4VWTiaGGmcrFSPBFbHek3SilfFkL8ceh1TqwCUiieBp3QcdR8VJNyGHIxNjfGPz/6Z4r8RZzLP0eqYe/qwvvnFvmw28vwpKaXU5idTHW5mcyUxOlfKLaPrSb9CSHELwE5oVG/Wkqq2JVkJWXx2ZLP4p5wc2P0BoNTg4zOjHIm7wzl2eXxDu+JWM+20FaWQ0VuerxDU+xhtpT0pZSvCCH+PTAOmKWUr8Q2LIXi6RFCUGmq5ED6AVzDLoZmhrgyfIU+vybglm7c/UlztW1hRW4650pNyrZQ8cxsuZErpfxdIUQWqIauYm+QZkzjU8Wfoneyl2sj1xieGeabPd/kZO5JnjM9tyvnsC8sabaF94eXbQufLzdTmK38BRTbw1YbuW+jNXF9gAAkkLCN3MdJKwshaG1tpaGhAYD6+nolrxxHyrLKKEgr4OrIVQamBrg2eo1+fz/2QjtZSVnxDi9Cv3cGV6+yLVTElicZ6b8cy0D2EhcuXIh8XVlZuUYe2WazrUj60Vit1ohEg8Ph4LXXXltXXlmxvaQYUvjEwU8wODVI53AnnjkPHT2agNsR85G4CrjNLgRw9nrp984CYMnQbAtNaftvyqki/mz1N71DCPEVIcTPhrfNThZC1AkhaoQQFx5z3qbH9ypms5mamhocDsem5yl55Z2nKKOIl8tfpiKrgiBBbnpu8m7fu4zPje94LFJKHgz7+fsbD+n3zmLQC6rKcnjpeIFK+IqYseUpm2iG55bHnSiEsAFIKR1CCKsQwraeD27IQ7f6SYJdj8v3YzNKrj9c/0zvb2xsjJR7NkLJK8eHJH0S9kJ7RMDNN+/j3b53OZJzhGOWYxh0sV+oPjGzyAfdHsamFgAoyknFXpajbAsVMWerI32nlPJ3pZRvhrdNzj3P8pRON5CQ+sFWqxWTybTGQCUsr1xVVYXP58Nms8UpQkVBesEKAbe743dx9DpiKuAWCEpuDPj4+s1HjE0tkJqk41OHcvnM4TyV8BU7wlZ/y6xCiCtAxLdQSvk/b3CuCfBGvV7zdBAa/TuEEI1bjnQDnnVEHksuXry4ZiSv5JV3FwadgbP5ZynJLME57GRyYZL3+t+jMruSU3mntlXAbWRyjg+6vfjnlm0Lz5aYSDKoRq1i53iS8s52Yt7soBCiAWgAKC3dfXr5W8Vms+F0OjGb13670fLKKvHHH0uqhZrSGu5673LXe5euiS4eTT/aFgG3+aUA1/p8dI1OA5CVauD5CjP5mWoapmLn2dIQQ0rZjVam+R3g50KvN8LHclI3AZ7og+FR/mM+r01KaZdS2vPy8rYS4q7l4sWLGzZ033zzTb7yla/scESKjdDr9JzIPUFNWQ05yTnMLM3wvcHv8eGjD5kPzD/VNXs90/zDjUd0jU6jE3C6OJt/cfKASviKuCGklI8/SYg/AjqAq2jJv2YjwbVQI9cupWwLzc5xSCk7hRAmKaVPCBE9rL0IvLZeozeM3W6XTqdzxb47d+5w7Nixx8a9V9nv399eICiDPBh/wC3PLQIyQLI+mXP55yjJLNnS+6dCtoWPlG2hIk4IIVxSSvvq/Vuu6Usp/9/Q129uVosPJXh7aHaOLyqhvwtUSSnbQwE1oD0JKBS7Dp3QccR8hKKMIpzDTkZnR/nBox/QP9nPuYKNBdyCQc228KOBaNvCHCrz0nflCmBF4vEkgms/CziAWrRZORsipWxbZ1/VOuesOU+h2E1kJGXwmeLP0D3RzfXR6wxODzLSM8LpvNNYs60rzvVOL/BhtwfvtKaXU2ZJw1aaQ2qS0stR7B6eVHDt54EPleCaIpEQQmA1WSlML6RzpJNH049wDbsY8A9gy7eRok/jxuAE96JsC+3lZopMe1fOWbF/2ar2zjvRMgzKREWRiKQZ0/iRoh+hb7KPa6OagNvlO//AwkwRGaIEnU5wpDCT08XZyrZQsWvZanlndTFSmagoEpbSrFKyk3L52s1/4o6nG/CRnzZE/YnPUG5WfxqK3c1Wk75bCPHv2GJNX6HYz3SNTnG1z0dK4CiHMszIlG5y0iWuse8yI49x1Hw0rgJuCsVmbLWm/yUhxGvAl9AkGb4U27B2P0KIFYuqmpubsVqtG0onK0nlvc/k3CJXomwLD2SnYC8/Q4rxDDfGbuCecHPLc4sB/wD2QjvmlE3XICoUcWGr8/RXmKYIIcqllD2xDCzMbp2nX1lZSVdX15r9VVVVuFwuQBNUu3z5Mq2trWv2t7a2biipvBu+P8UywaDk9qNJbj1cti2sKsuhfJVt4cjMCK5hF1OLUwgEh3IOccJyYkcE3BSK1Ww0T3+rz6CrpexkUP0AACAASURBVCBbnz2k/Y/dbl8juAZKUnkvMeqf5xu3hrgxoCX8itx0fvz0gTUJHyA/LZ/asloO5xwG4P74fTp6OxiZGdnpsBWKDdl0CCKE+Dk01cwaIYSVtQ3d+OP8s9hc1/7FZ75EW1sbtbW1a/YrSeXdz9PaFhp0Bs7knaE4oxjXsIuJhQm+M/AdrNlWTueexqhXK3IV8WXTpC+l/BrwNSHE70gpv7xDMe0JwhLJoClnhpO42+2msbERr9e7Zn99fT1ut5uLFy9uqrOviC/RtoU6AUefwrbQkmqhpkwTcLvjuYN7wq0JuOVXcSDjQAyjVyg2Z6s1/Z9D88WNECXLEFP2Wk3/cbX+sKTyeueE2Q3fXyIys7CEs2ecgfHttS2cmJ/AOezEO6c18UszSzmbf5ZkffIzx6xQbMSz1vStQGVoe4ntl1pOGKIllRW7g2XbwkcMjGu2hfby7bMtzE7O5sWSFzmTdwaDMNDn7+Odnnfom+zbhugViidjq1M2fzf6dUh1U/GUvPnmm3z+859XOvq7AN/MAh92e2NuWyiE4HDOYQ6mH8Q17GJkdoQPhj6gz9+HLd9GmjFtWz9PodiIrcow/DHL5R0BVG1yekKwUXlmo/3h6ZoAJpNpxWvFzhMISm4OTnDn0SRBCalJOuxlZkrMsU2+GUkZfKbkM7gn3NwYvcGj6Ud8s/ebnM49TUV2hVLiVMScrQ5nmtGSvoCIqYpCsScZnpzjwyjbwkMFGZwp3lnbQmu2lcK0Qq6OXOXh9ENcIy76/f1UFVSRkZSxY3EoEo/HTdnMQqvfVwPZaEYqSg5ZsSeZXwpwtc+HO2RbmJ1qpLoiJ24uVmnGNH646Ifpn+zn6uhVRmZH6Ojt4LjlOIdzDqtRvyImbJj0hRAVwGXgt4F2NMMTO+ASQtRJKa/vTIgKxbPTMzaNq3ec+aUgOgEni7I5fiALnS7+ibUkq4T89Hyuj1yn19/LjbEbESmH7OTseIen2GdsNtL/Y+DzUsqJqH1XhRCXgUvAy+u/TaHYPay2LczPTKZ6F9oWJuuTef7A85RkluAaceGd9+LodXDUfJSj5qPodcqIRbE9bJb0xaqED0DI5zb+wyOFYhP2qm3hgYwDvJz6Mh+NfUTXRBe3vbcZmBrAXmDHkmqJd3iKfcBmSd8thDgrpbwWvVMI8SKPkVYOmZ/7AJuUsmWd4+HlqLVSyqYnjFmh2JT1bAurynJIMe6N0bJRb8RWYNNG/cMuJhcm+Vb/tzhkOsTJ3JNKwE3xTGz42xOSU35bCCGBK4AFbaqmjHbRWo0QwhZ6v0MIYRVC2KLM0cPHa6WUTUKIptXH9wobSSsLIWhtbaWhoQGA+vp6Ja+8QywFgvvKtjAvLY+ashpue25zf/w+D3wPNCmHgiry0/LjHZ5ij/I47Z1XQg1dG9ANtG1huuZ5tFk+oD0R1ACRpB5K8OHX1r2Y8EHT21lPGtlms61I+qvf09qqCZQ6HA5ee+21DeWVFU/GQ98sV3q8TM8HEIJ9Y1to0Bk4nXea4sxinEPOiIBbRVYFZ/LOKAE3xRPz2L8IKWW3lPJroW0r8/NNgDfq9bqFSCHEBfahnIPZbKampgaHw7HpeUpeeXuYWwzw/Y/H+Pa9UabnA5jTjbx8opCqspw9n/CjMaeYqSmr4YTlBDp0dE92842ebzA4NRjv0BR7jLgVB6WULUKIy0IIp5RyRfYTQjQADQClpaWbXmf80tsxiS/n/CtP/d7GxkYaGxs3VdJU8srPTti2cGEpiEEnOFWczZGCzF0xDTMW6ISO45bjFGUU4Rp24Znz8P7D9ynJLOFs3llSDPFZb6DYW8RiKOQDwj5xJsATfVAIYQvX/dHKP2vqIFLKNimlXUppz8vLi0GIz05YKrm+vp6mppW9aKvVislkWmOgEn5PVVUVPp8Pm82G4smZnFvk3TvDfOD2srAU5EB2Cj92+gDHdsm8+1iTnZzN50o+x9m8sxiEgX5/P+/0vEPvZG+8Q1PsAWIx0r+EtogLNHVOB4AQwhQa0UfX+E1oTeKn5llG5M/CRjX9MBcvXlwzkg+/JyyvrATXnoyt2hYmAkJodowHMzQBt+GZYT4c+pB+f78ScFNsyraP9MON2dC0TF9Uo/bd0L9tgDU0rRMp5b7UGLbZbDidznXr9kpe+ckZ9c/z9ZvLtoXWvI1tCxOJdGM6ny7+NPYCO0adkUfTj3in5x26fF1sxStDkXjEpKYvpVyjzyOlrAr962NZv2dfZ72LFy9G3LVWo+SVt8bCUpDrAz4eRNkWvlBhpiBL1a+jqciuoDC9kKvDVxmcHqRzpJO+yT7shXYykzLjHZ5iF7El56x4sluds2LJfv/+tkq/dwZnr5fZBU0v59iBLE48oW1hItLv7+fayDXmAnPohZ4TlhMcyjmETqifWyKxkXOWWtqn2HXEyrYwUSjJLKEgrYDro9fpmezhxtgN+v392AvsmFLUYsBERyV9xa5BSsmDkSmu9ftYCkgMesHZEhOH8jN2tV7ObiRJn0R1YXVEymF8fpx3+95VAm4KlfQVuwPfzAIfdHvxhGwLi3NSsZfnkJakfkWfhcL0Ql4qf4mbYzf52Pcxt7236Z/SRv25qbnxDk8RhQwGCUxMEBgbY8njYWnMQ/oLz2PY5mnre/YvSkq5L0d/u73Hst0sBYLceji547aFiYRRZ+Rc/jlKMktwDjvxL/h5r/89DpkOcSL3BEadknKIB8GFBQKh5L40NkrA60UuBVacs+TxqKQPkJKSgsfjwWKx7KvEL6XE4/GQkpIYM1OGJ+f4oNvLVBxtCxOJ3NRcaktrue25zb3xezzwPeDh1ENsBTYK0wvjHd6+JzA5ydKYh4BHG8kHJibXnKPLSMdgycWQa8GQm4sue/tNdPZk0i8uLmZgYIDR0dF4h7LtpKSkUFxcHO8wYsrcomZb2D22bFv4fIWZvMzkOEe2/9Hr9JzKO6UJuA078c37+N7g9yjPKudM3hmS9KpZvh3IpSWWvN7QSF5L8nJ+YcU5Qq9DbzZjsFjQhxK9bgcGfHtyyuaWkBL20VPAfiHatlCvgxMHd49tYaIRlEHuj9/ntuc2ARkgRZ/CufxzFGfu70FHLAhOT0fq8EueMQI+HwRX5lZdagp6izaCN+TmojeZEPrYNdQTb8rmza+BIRnSciE9D9ItkGJSN4I44Z9bxNkzzqMJzbawIEuzLcxKUfXkeKETOo6aj2pSDkMuxubG+OdH/0yRv4hz+edINexNH4JYI4NBAuPjK0o1wZnZlScJ0OeYtAQfSvS69N2xenx/Jv3FOZj3a9v0GIze1fbrjStvAmm5kJyx4q1vvfUWr7/+On19fZSWlvLGG2/w6quvxuGb2B8Eg5K7Q35uDmq2hUkGHedKTVTmZTz+zYodISspi8+WfJYuXxcfjX3E4NQgozOjnMk7Q3l2ebzDizvB+fmoGTVjWsM1EFxxjjAaMeTlLpdqzDkI4+4c0Ozf8k5gUUv4M2MwPQrTHliYWnueMU27AaTn8dY/fI+GX2tiZmYmcjgtLY22trZNE/9GLlrRuN1uWltbd1RO2eFwcPny5Yhxy+P2bxdNTU243W6Onz7HD/30Fxmf0WwLyy1p2Mpy+MUv/HxEZTT882hsbMTr9WI2myNx1dfXrzlPEVtmFmdwDjsZnhkGoCCtAHuBPWEE3KSUBCcmIqWagGeMgH9t3tBnZa4o1egyM3fdpJKNyjv7N+mvx8JM6CYQ3kYhsNxcKf+f/j29Q541bysrLaWnd2PZ2srKSrq6ujb96Hgk/Wja29t3ROens7OTv/rrv+Ff/spFfus//iYv/vSrVFgrqC43c9CUSltbG3a7HZvNRn19PRcvXgTAZDJhtVoj+5xO55rzlBT1ztEz0cP10essBBcwCAOnck9RaarcdYntWZELCyyNj2sj+LExljxe5OLiinOEQa81XHPDI3kLuuTdP+kg8Wr665GUBkmlYIoyZpmbDD0JjNE37F33bX39fXD770JloVxt22P9gUuXLu1I0v/a//g6WZU27g35sR47zdj9K/zKT30yopdTU1MTeQqqra3F6XSusJZ0u90RP4LV56mkv3OUZ5drAm4jVxmYGuDq6FX6/f1UFVaRlZQV7/CemsDU1IrFT4EJH6wa9+rS0jDkasndkJeHPjsbods/04gTK+mvR0qWtlkqKS0tpXedEX1pgQVmPNoWniWqN0KaVhYiuATzU5H+gM/n47XXXsPn81FfX78iqTkcDlpbWzGbzTQ2NkZGsj6fj8bGxjWJOWy8YrfbcTgcdHR0RJJh+H1h793V1zaZTCs8eRsbG2lubsbr9dLa2orb7Y6Ys7e0tGCz2aipqVkRT/ga0fGux+xCgM6+cW58PMAPlZ/FnG7k06fK+eB7314hkBZd9nK5XDQ2ao6ZYY+B8GdGG8ZHn6fYOVIMKXzi4CcY8A9wdeQqY3NjdPR0cNxynCPmI7tewE0GAlrD1eNhaXSMgGeM4Nz8ypN0AoM5Z7lUY7GgS9vfpSyV9KN44403aGhoWFPTf6Plq3DkpeWS0PSY1h/wD4F/CHdPH/U//jkQeqwVZTT/p4tc/pPfh/Rcql745Iqkf/ny5RWlipaWFmpra2loaNjQWCWc1Nvb2yMlopaWFs6fP09dXR3t7e20tbXhcrlWXDvs3NXc3ExnZ2fkBuD1ak8058+f5+2336ahoYGOjg4uXLiwJh6z2bzmmtH+v1JKRv3z9Htn+MxPfQGdgKOFmbx0vJBvDW3cyOrs7MTr9Uaua7PZcLlckZ9B+Oaw+jzFzlOcWUx+Wn5EwO2m5yYDUwPYC+zkpOTEO7wIwdlZbdHT2Jg2N358fM20SZGchCE3T1v8ZLGgz8lBGBIrDSbWd/sYws3aDWfvZEatWlycjdwErCUHuPxffh2WQqOIh1dp+cuv45mYwv3gHri/A2OzsDBD07/732j6zdfx+Xy0trbS1dWF1+vF5XLh8/nw+XwrRrmgma0D1NXVRRL3lStXIjcTm81GU1MTzc3NNDU1Ra79OOrq6qitrV1Rclkdz+uvv84bb7wRuabVao187sTsIle6vWT657ECB0wp2I+Wka1fQKcT+Hw+LBbLms8NX2u185jJZKK+vp729nYuXLiw4XmKnSdawK1zuBPfvI93+97lSM4RjluO77iA23o6NcHp6TXn6bOzV5ZqMtSsMZX0V/Hqq69ubYqmMRVMJdpmTIWzX9D6AzNjtPyX38Na+RwXPnWM9ned4HXDwxHwdmGdfJ/Lv/ULOK720Pp7X6Hq1HHMBcXUbWC2Assjc4fDEUnO1dXVOBwO6urqcDgcVFdXR+wYw2WerZRETCYT7e3tkXOrqqowm80rnjhWX/Ob3+zg0cQcD32zSMCgE5Ra0vjCb/wq2T/6EpcuXaKmpoaOjo51Y2hqalpxU4puMrtcLmpra9c9TxF/CtMLqS2v5dbYLR74HnB3/C6DU4PYC2Mr4LYVnRphMIQSfEjGwGxGJKkVxquJSdIPWSH6AJuUsmWd4+F6R6WUsmn18T1LqD9QU/dFXnvtNTpujUFSOpT9MMxcA8MHtPzFP9LxwS0Amn+1HttRE/Vf/l1af++3MeWYufzf/0RrFCcvux2F6/put5t339VcJy9cuEB9fT1f+cpXIsm+paWFjo4O7drrzBKqra1dk0TPnz9PU1MTFy5cAIiUdVpbWzGZTFRXV6+4ZkbeQUo++ZNkzS5yBM228FypiWSDNtKz2WxcunSJ+vp6rFZrpCxTW1tLR0cHLS0tOBwOqqqqAM1drKamZsWUzbq6unXPUy5juwOjzsjZ/LMRKYewgNtzpuc4mXtyWwTcnlanZr/NLooF2z5lUwhhA6xSyvZQcndG+eSGvXPdUkq3EOIy0CqldGx0vW2dsrkbCCxpDeHp0eXpo/P+tecZUiA9D/fYLM1tl2ht+xMwxk+IbbVtYWaKgeeVbWHCEwgGuOu9y13vXYIESTOkUVVQ9UQCbmt0asY8yIXdoVOzl9nJKZvngY7Q126gBuiMOm4NbW2h4ytXMe139AbILNC2MItzUTeB0A1haQ4m+mF4BCYG4Ppfa6P/9NzlVcVpFu16MWa1beHxg1mcOJiNXunlJDx6nZ4TuScoyijCOexkfH48IuB2Ou80yfq189nX6NSMj6+dNpmaounThBJ8rHVqEolYZAwTED3hfUUnb5Vpug24FIMY9hbGlOX+QJh5v9YkLvDQevSU9nQQlpbwdmvnCAGpOVHSEqH1A9s0p3h6fgln7ziDIdvC3Iwknle2hYp1MKWYeLH0RR6MP+CW5xY9kz0MTQ9xJvc0B5cy9qxOzX4kbo3cUBmoI7r0E3WsAWgAKC0tXX04MUjO1DZz6EEoGIQ530ppidlxmPFq29h97TydIbR+IHf5qSDlyRbTSCm5PzzF9YFl28JzJSaeU7aFik3QCR2H0svJnxTcevA+U0N3uT/xfUaM2ZRmlWDUaYMFkZS0PGVyl+vU7EdikfR9gDn0tQlYq2ugUbNekxciTwNtoNX0tz3CvYhOB2lmbeOwti/cH5iJkpWY98PUsLaFMaSEbgCW5ScC4/oKiqttC0vMqVSVKdtCxVo206k5QSqjATMDcgaPcQ5PygjPWe2UV5zdlTo1iUQs/pIvAeHmgRVwAAghTFJKX+jrhnDCF0LUbNbIVWzCRv2B1fpCS3NaX2BiYPm85MwVN4Gl5BxuDs9yN2RbmJakp6osR9kWKiI8qU5NhcVCUVYqnb6bDM0M0ckwg5PXqEqtIt2oyjfxIiaCa6HyjBttFk9baJ9LSlkVmr1zGa3ubwbqE2r2TjyYn1peSRy+IQSXIocnZhfp9szgJ4O5ZDMHDhRz2FpBUoZl2/oDir3HdurU9E32cW30GvOBeQzCwMnckzxnek6N+GOIUtlULCMlzI4zPznCA3c33tFHJC9OkGYUlOemk5kcegDU6ZefBsL/PmF/QLE32LJOTc7T69TMLc1xbfQa/f5+ACwpFqoKqshO3n4fWIVK+opVdI9N0xltW1iYzvHsJXSzoSmj4f7AaiJuZLmP7Q8odi/x1Kl5OPUQ17CLucAcOnQcsxzjqPnorhdw22soaWUFsAXbwqyo/sDS/HJfIFwWWpyFyUFtC5OUsTxbKLJ+QM3G2C3sNp2agxkHyUvN4/rodbonu7nlucWAfwB7oR1zivnxF1A8E2qknyCsZ1toKzVhfVLbwvmplY3imTHNpWw1qabQDSB0I0jNUf2BHWIv6dSMzIzgGnYxtTiFQHA45zDHLccx6NR49FlR5Z0ExjM1z4fd3jW2hSnGbVjhKGVo/UBUWWh2HORKD1F0ekg1rzKiUbXc7WCv69QsBZe45bnFg/EHSCQZxgzsBXby0vLiHdqeRiX9BGQxEOTGwAT3h/1ICenJ+ohtYUwJBrQFY5Gy0KimQLoafdKySX34qSBJTRHdjP2sU+OZ9eAcdjK5oP2uVGZXcir3FEZVKnwqVNJPMAZ9szh7vEzPBxACjhRmcrooe4WL1Y6ytBBlUh/uD8ysPS8pfa2+kCFxZR9W6NSMjRLwrTNtch/p1ARlkDveO9z1aAJuqYZUqvKrOJBxIN6h7TlU0k8QZhcCuHrH6fNqCdWcbuT5Cgvm9F2YOBemVy4i26g/kJIdeiIIlYZSzfuyPyCDQW3a5ON0akz7X6dmYn4C55AT77wm41WaWcrZ/LPrCrgp1kcl/X2OlJKu0Wmu9o2zGJAYdILTJdkczs9Et1fUMKWEuYmV+kIz3rX9AaGL0hcKPQ2kZO8po3qA4Px81IyaMa3hGlj5vSayTo2Ukvvj97ntuc2SXCJZn8y5vHOUZJU8/s0KlfT3MxOzi3zY7WXUry2mOWBKobrcTEbyPpgBEe4PROsLzU2sPU+fFFUWCt0MdlF/YDOdmmj0WZnLi59yc5VODTC1MIVr2MXI7AgAB9MPci7/HGnG3fP/uxtRSX8fEghKbj+c5NbDCYISUow6qspyKLPsv8f9FUT6A+GnAY9WKlpNUvpKkbm03B3rDzypTo02kregS1bli/WQUtI92c2N0RssBhcx6oyczjuNNTux7DieBLU4a58x4p/jw24vk7Oahk5lXjpno2wL9zWGJMg6qG1hwv2B6DUEC9Pa5utbPi8lO6oslKuplm6DqfeWdWryQgk+N3dDnRrFWoQQWLOtFKYVcnXkKg+ntVW9/ZP9VBVUkZGkDM+3ihrp7zEWloJc6/fx8YiyLdyUcH8gbE05PaZ9vW5/wLzKiGbz/oAMBAh4vStKNdutU6PYnP7Jfq6OXo0IuB23HOdwzuGEL4VFo8o7+4A+zwyuPmVb+NQEgzDrXTlbaNa39jx9UlRZyEJQl8mSfzYuOjWKjZkPzHNt5Bp9fu1Jzpxixl5gVwJuIVTS38OsZ1v4QoWF7LTEmMURU5YWVhjVS/8ogfExAhMzLE1MszQxTXBuQbsRJGdoOkNJGejzizDkF8Zcp0bxeB5NPcI14mJ2aRYdOo5ajnLMfCzhBdxUTX8Pstq20KgXnFW2hdtKMAiBaVga07E0Jgh4Dcj5bFjQQ6hiI4xgyExGn52EIcuIIVuPMIxCyjwwA9OzQN629QcUT8aBjAO8nPoyH419RNdEF7c9txn0D1JVUIUl1fL4CyQYaqS/Sxmf1mwLvdPLtoX2MjOpSSqpPAtPrVMzP7lyIdmsV5tOGo3QacJykYVkFs2oXt2gd4zRmVGcw86IgNuhnEOcsJxISAG3HS3vCCHq0LxybRv54AohbOuZoq8m0ZL+UiDIzYeT3Hk0iVS2hc9ETHVqgkFNWC5aX2jd/oAxau1AaNposioFxZKl4BK3Pbe5P34/IuBWVVBFflp+vEPbUXasvCOEsAFIKR1CCOt6yT1kmdgMVG335+9lHk3McqVnnKk5bRrm4YIMThebSDIkdm1yq+yoTo1OFxKKiyofBBZXzhaaHoOFKfA/0rYwxtSV+kLpuZo5jWJbMOgMnM47TXFmMc4hJxMLE3xn4DtUZFVwJu9Mwgu4xeKZ5zzQEfraDdQAK5J+6IbgjcFn70nmFgN09o3TM6bp5ZjSjDxfYSY3QyWCjdiyTk3ODurU6I2QWahtYRZnl0tC4XUEi7Pg69e2MMmZa/WF9IlXkthOzClmaspquOu9yx3PHbonu3k0/YiqgioOZhx8/AX2KbH4rTKhmZ6HUZ2UTXCPTnG1zxexLTxZlM2xwqy9o5ezQ+xZnRpjKphKtC3M3OTySuLwv/N+bfO6tXMi/YEofaHUHNUfeEJ0Qsdxy3GKMopwDjvxznn5/sPvU5JZwtm8s6QYEm99ixpKxAn/3CJXerwMTWhTRNbYFiYw+16nJiVL2yyV2utgMGREMxqlL+TTbgYzHhi9p52nM6zVF1L9gS2RnZzNiyUv8rHvY26O3aTf38/IzAhn885SmlUa7/B2lFgkfR8QNro0AZ4nvYAQogFoACgt3V//IcGg5M7QJDcHJwgEeXrbwn1EwuvU6MKrgs2Qd0TbF1hcNqKJfhrwD2lbGEPKSjeytFwwJt7odSsIoc3mOZB+gM6RToZnhvlg6AP6/H3Y8m0JI+C27bN3Qo1cu5SyTQhxAXBIKTuFECYppS/qvA4pZe3jrrefZu+MhWwLfWHbwtw0bKXbZFu4h1ihUzM6RmByQunUbIXFuZUm9dNjsDS39rzkzFX6QhbVH1iH7oluro9eZzG4iEEYIgJue+JpcQvs2OydUIK3h2bo+KJm7rxLaLZOaEqnXQhRJ6Vs3+4YdhuabaGPe0NaiSI9Wc/zFWYOZMfYtnAXsGWdGrPSqXksxpT1+wMzY8sexSv6A93aOUJo/YAV+kKmfWlE8yRUZFdQmF7I1eGrDE4P0jnSSb9fE3DLTMqMd3gxQy3OijED4zM4e8aZWdBsC48WZnIqnraFMSY4O6stelI6NfEh0h+I1hca1wTootEZVugLaf2B/ZvoHke/v59rI9eYC8yhF3pOWE5wKOfQnpZyUNo7O8xa28IkXqgwk7MbbQufEhkMEpiYWCEpHJxeq2uvz87W5sQrnZr4EFhaoS/E9Jj2JLAaQ8rK3kB6XkL1B+YD89wYvUHPZA8AOck52AvsmFJM8Q3sKVFJf4fQbAu1aZjRtoVHCvbIzJJNCC4shFa3hhY/eb3IpZVSBMJgCCX4kIyB2YxI2j83un3D4lyUUX3ohrBRfyDyRJCXEP2BoekhXMMuZpZmNAE381GOmo+i32O6Sirpx4iWlhauXLmCz+fD4XDwqZd+gvmlIJ//6Vd56aUaqsvNpO+wbaHD4eDy5cu0trZuaf9GPKlOzX/4gz+g5+FDqquruXDhwppzm5qacLvdK443Njbi9Xoxm81bjksRI+b9yw3i8BNBcGnlOUJo/YDIjKG8fdkfWAwucnPsJh/7PgYgMykTe4Gd3NTcOEe2dVTSjyFh28KaHzrN733te6QYddjLzJRa4t+MbG9vp66u7rHnyaUlljzeSIJ/Up2azs5OLl26RHNzM01NTTQ2NmK1LlvZrXfc5/NhMpmwWq3U19dz8eJFnE4nly9fBsDr9XL+/Pl1byCKHSC6PxC+Ccx6N+kP5C4/FaRkxSfmbWZsdgznkBP/olYOO2Q6xMnck3tCwE1JK8eIkck5PuzZvbaFly5dWjfpP5VOTU7OhtMmHQ4HtbXaDNzq6mocDgcNDQ1bPu52u7Fardhstsj+xsbGLd2wFDEiev0Ah7V9gaWQEU3UQrJ5P0wNa1sYQ/JyOSj8VGDce7PVclNzqSmr4Y7nDvfG7/HA94CHUw+pKqiiIL0g3uE9FSrpPyXzSwGu9fnoGtUal1mpBlKNel6waqoTbreb5uZmvF4vv/zLv8xf//Vf09raSktLC1arlbq6Ourr6/H5fOsmN7fbTX19PXa7HYfDQUdHR2TkHH6f1WqltbUVh8NBa2srZrOZxsZGTCZTpFQSTq6//eUvM+ru5s0/oCMC6AAAEKNJREFU/zPc3d38/r/9JbLT0/m//u5vOV1RwWdPn+aL//cfMrmwQMO/+leYS0p48y/+InJN25HDm/48PB4PZrO2Js9kMnHlypUtHe/s7KS+vj4SdxiHwxF5ClDsIvQGyMjXtjCLc6tsKUP6QhMD2hYmKWO5JBR+KtgD4mcGnYFTeac0AbdhJ755H98d/C7lWeWcyTtDkn5v9axU0n8K+jwzOHu9zC1ublv49ttv093djde7VluupaWF2tpaGhoaqK+vX3dEG07q7e3ttLa20tzcTEtLC+fPn6euro729nba2tpwuVxcvHgRm80GwMd37hDw+/k/vvQlrjgctLz0EvKKk7nhYQKTfn76+Rf4/65c4Ze+8C/5bnc3v/nV3+Or/+1P+bFf+IVIPGazecU13W43Dodj3Z9H9Ij9SbHZbLhcrsjPIJzkm5ub6ejoeMy7FbsCYwpkF2lbmPmptfpCC1PaNt6jnRPpD0QtJEvN2bX9gZyUHD5f+nnuee9x23ObnskehqaHsBXYKMooevwFdgkq6T8B0/NLXOnx8tCnzXLIy0zm+XLzhraFr7zyCiaTad2k39XVhdfrxeVy4fP5IvXtaMIj47q6usjI/cqVK5Eka7PZaGpq4rdff50vX7zIuMfL7zU2EJyeYXFgkLk7dwkuLiIDQfRZmRiTjBgPFPKLv/Wf+NGf+zl+sqiI506dwliQj7unB29nZySe119/nTfeeAOfz0draytWq3XT5G6xWCLfp8/nw2KxbPm4yWSivr6e9vZ2Lly4gM+3ji69Ym+RnKFt5grttZRR6wfC/gPjy9vYA+08nV57Aoj2IEjZPZ63OqHjmOUYRZlFOIeceOY8vP/wfYozijmXf25PCLippL8FIraF/T6Wgppt4blSE5V5m9sWrk7ioCV7q9VKVVUVZrN505p1OEk6HI7ICNh+7hzfaG/npz7xCf7+z/6Mk6mp5N67z5/86y/y7Y9u8N/+9u/44o++jC49jZTjxzCYTGT/zE+jS0oize1G/4MfoM/KwmQy0d7eTmNjI8C68Vy+fDlSOmpsbNx0pF9TU8OlS5eoqamho6Mjct0w6x2PbjK7XK5Izd/hcESeMBT7hPCq4NQcyD2k7QsGVukLjWkrjKdGtC2MIXmtEU1SfCdJZCVl8bmSz0UE3AamBhiZGeFM3hnKs8vjGtvjUEl/FW+99Ravv/46fX19lJaW8pv/+29R+vzLEdvCUnMaVWU5T2RbaLVacTqdNDY24nZr0rnhMkpraysmkykyYyUat9tN3c/8DO6PP+bv/+APmPzGOzQ89xz/+qtf5Xd++7cpy8/nz3/9N/iDr3+db9++hTAa+cobb5BZXExSby+pp04hkpN5+cd/fM10yPPnz9PU1BSZGbM6nurq6kh5pbm5+bEjfZvNxqVLl6ivr480ZAFqa2vp6OhY97jVal0xZTN8A3C73VRWVm7556vYo+j0kJGnbWGW5qNmC4XWECzOwOSgtoVJylheSRy+IexwfyAs4HYw4yCuYRfDM8NcGb4SkXLYrQJuaspmFG+99RYNDQ3MzMxE9iWlpPJLX/4dan/y57CX51CcE7v/yGidmgc3bvDVP/1Tvvpv/u3Kk3QCQ47SqVEkEAvTKxeRzYxpKqSrSTVF6Qvl7Xh/oGeih+uj11kILmAQBk7lnqLSVBm3RZlqnv4WKC8vp7e3d83+wqJi+np7MW6zXs5mOjU9w8P8n3/3t/z+r/yK0qlRKKKREuYmVuoLzXhBrjTVQafXHMiipaeTs2JqRDO7NMvVkasMTmlPJbkpuVQVVpGVtPPrFlTS3wI6nY71fh5CCILB4Drv2DpKp0ahiCHh/kCkLDSm3RhWo09aqy8Ug/7AgH+AqyNXIwJuxy3HOZzz/7d3Pk9tJFcc/7YQMggD4oexjY1/YNc63tiJLZlLKj5sSkrlsgdXmeVArqD8BZrlkhtJ0ClX2NpbLgtspXINOHtJaqvWeFIbh/KmHGSvbXZtI2AMNmuEofcwPaPRaKQZgSSkmfepUoHU86Nb03rT8/r1971XVQE3WpxlQ2rlNbqO9yL9fDmvbD+JXEinhiCqSM78wGX1s3cZQ+4BcSPY2QI2vlNfGoGWfH0h/8F+h6dbT6Mn2IOvV77G443HuJ++j2ebzxA5HkFHU8eBjn1QPG/0N9/u4KtHa3ixsY2PfpfAp3/6GNtvswm2g8EgxsfHbY+To1OTTmN3I1/F0He0JcdV42tvr3sRNoKoWfwBoK1XfWlk3uTrC2XeqC/lSXa7pvZcfaHmDvXGUgKBhgAGTgygr7UP917cw/r2Ou48uYNLHZfwftf7hybg5ln3jjlt4RG/D9fPhPDl3/+WE70zPj6O4eHhnH0PqlNDEESNYJwf0G4CW6v58wNMSFIYo4Wa2h3PD+zs7WAxvYiHiroeoRoCbuTeMVAsbeGXFtuXW6eGIIgagTE16qc5BOCi+tnerrpgLCdR/avsE4KGNj9gXEMQaLE8TaOvEdd6rulSDpuZTXzx9AtcDF3Ele4raPRVL9zUUyN9c9rCo01+DJzr0NMWWoVsNh85gj+PjGDwlzezB2JAQyikh0z6u7vha7G+2ARBuIB3mfxENJn8QAw0BrM3AO2pwDQ/sLu3iwdrD/DN2jfg4Aj6g4gcj+BEy4myVrmq0TsiB64CIMw5T5ZabqRcRt+ctvDyyTZc6W2D792OHlFzKRbD05cv8/Y9fewYHnz+edZV09kB1pi9MzPGclayaouZjKRSKV0/p1YoVV+/VKz0852Ux+PxnDqR5j5Rk2S2DDcBcUPYzeRv19SeGy0U7AR8DVDeKlh4sYD17XUAKLuAW9XcO4yxMABwzucZY/2MsbAhObptebn5IbOLhW/X8HTtB4BzdPNtXG/dQ/C759i6n8bu5mt922crK5bHWE6ncfTmTcsyQF1xa7WithYxSh9Eo1FEo9GKnEeW1Us6MzOjG3ezvr65vLOzE5IkYXp6WjfusixDkiRdc1+WZZJoIGqDQBAInAU6zqrvOQe2N0z6Qmuqa+jtK2B1Sd1OzA+Egt34VcsZ/L+xDf99/VQXcLvWcw19rX0Vq3YlfPpDADR5xBSAKAC5hPKywDnHw+V1LC4+BtbX0bGp4Lw/g+NBHwAG7X7M/A3qhGt3N/p6e/FkuTwhm7VKIX39crNfff3JyUkYn+yMBt584yCImoIxdVTf1A50CRmRvT3D/ICWqD4rPOdbUTMV9PFdPNxew0sfx39ePcVy1yX8/NQv0Owvfw6CSsw0hgAYZSW7Siw/MK833uBfn07j0V+mEfz3V+h9+QjhYAbHgw3wBVsQOHsGwfB1tP46hvZbt9D6wQdovnoVf5iYQNAkaeA0ZNOIoigYHBxELBbD1NRUTtn8/LyuH6+NdrVtZ2dn846VSqUQi8UQj8cRiUR0BUrzPtr7eDyOZDJpWQdJkjA/P69nrUqlUpAkSdfnB1TJZ01YzXgOq3oXw6yfv7S0VFK5EVmWceHCBQwNDVmK2BFEzeLzqRpBPT8Bzt8EfnoLuP5b4L3fAKdvqE8JgaNoZg242nQMP9tl6H3+P/gX/4rFf/we3z/5Z9mrVJPRO4yxUQCjwP5G2YFgE3bebKGx0Ydz/adw4lyvI50aLTTTLmTTjJbwBFBdPRMTE7q7JxKJ5IxwZ2ZmcnTqnejqA8Dk5CQkScLCwgJkWc7ZRxNISyQSumiadi5jHSYmJiDLsl6uqXgODQ1henoao6OjmJubQyKRyKtXtfT1rSikuU8QdUlDI9B2Un1pZLbAttLoeZNG28Yyvv3+HjbfptHQUP4w70oYfQVAp/g/BGC1xHJwzqcATAHqRG6pFQj4G3D1ow9xtKMVTU1HStp3eHjY1sibsfLpJ5NJrK6u6qqaGpIkQZIkXafeia6+ZuS6urqgKErePn6hxaMoij56LlYHM7dv30YsFkM0GtXPZT5HNfX1rTBr7hOEqwgEgcAZIHQGTafCuHT5Q6wqj9DVerrsp6qE0f8MgDZj3A9gHgAYYyHOuVKovNx0nzy8rPVaSsREIpHnstFuEJpOvRNdfTPmfWRZxtzcHO7evatHBxWrgxW1oK9vRSHNfYJwO12h8xU5btmNPudcZozdYIxFASiGyJw7ACJFyl1DNBrFyMiIZbq/ZDKZo1MfDodtdfXNmLXvx8bGkEqlEAqFIEkSJiYmitYhFovVpL4+AD3nQDwehyRJiEajlpr7BEHsD08tznIrg4ODuiGudOw9QRD1AckwuJixsbGcRV+FXCUEQRA00icIgnAhhUb6pAhGEAThIcjoEwRBeAgy+gRBEB6i5n36jLEVAPnZyp3RDSBtu5W7oDZ7A2qz+zloe89yzo+ZP6x5o38QGGMLVhMZboba7A2oze6nUu0l9w5BEISHIKNPEAThIdxu9KfsN3Ed1GZvQG12PxVpr6t9+gRBEEQubh/pEy6CMXabMRZljBXVVrYrJ4haREslW6DMUd93gmuMvt2XUs4vrVZw0OZR8aqdbOz7xJhbGYBS6Aci1FsHqlm3SuLgGofFNq6RHy3ht1zebD2HiOi3nxQoc9T3neIKo2/3pZT7S6sFHLQ5CmBeJKTpF+/rmSGoCXiAbG5lV+Ow38Y557NQr7EX+nUYQEqUp9zQZkBv71qB4rL2fVcYfdh/KW40GHZt6jd8lhLv6xnb3MqMsbD48biFotdYjO6XAIBznnRJbgonv1XtybXfJW22o6x5xd1i9A89GfshULRNnPMpMcoHgDAAL0iVdtpvUlfY9dsBAF3CxeMWt6Vdv5ahjvCXUHhkTBTBLUafKIB4/J1zwYioaG5lF47ynbKqXVs3+fULwRgLQe0LkwA+YYzV+xOsE2zzipeCW4z+gZOx1yFO2xTlnCerU6WK8hmyLqqc3MvaZ4YJTVf4t2F/jZegukAg/rphAtuuzaMA/ij69CAA197oDH3bsu/vF7cYfTuDUNYvrUawazMYY6Oawa/3iVzDaNYq9zI457NiQrMTqrFwA3bXeN5UfreqtasMtv1aQ5vsrV7VKocYrNwwPa1pfbtQ39/fudyyOEuEb6WgTu5Mic/ucc4jhcrrnWJtFh1kBqrfsxPAoEfdH3WNw369Jsrd8ETnpM0JUd7plt9yNXGN0ScIgiDscYt7hyAIgnAAGX2CIAgPQUafIAjCQ5DRJwiB0ClaEqGfo4yxSVN5lDEWEouhlhhjk8Y4ccZYP2NspsCxQy4JIyXqHDL6BJFlAape0ayIClE0US8tZJBzroXMyQAmOedanDzE/yNWB+acK6h/KQzCBZDRJ4gsUQD3DO+NkgCjdiGvYtRfzLDLXlg1S9Q2ZPQJIssAhEaRGNn3iwVfAHDB4TF0GWvGWMIoASyeBGJlrC9BlAwZfYLIEkZWhjrKOS/JQAujrgCAyGEgi6cD4w3DbaJwRJ3hP+wKEEQtYPDZz9ptW2D/26Z9+wFok7dLhs9JGZI4VMjoE4TKDRxMk8nsy78LNdmHzBhLWe1AEIcBGX3C84gJWAkOxbvE6D0MYEjsO4RsFqcwY6yfc54UPn3NnUO6R0RNQNo7BOEAoVh6IHEvcYMI79eFRBDlgCZyCcIZ02UItySDTxw6ZPQJwgFicZVipevuBDHKJ98+ceiQe4cgCMJD0EifIAjCQ5DRJwiC8BBk9AmCIDwEGX2CIAgPQUafIAjCQ5DRJwiC8BA/AiGFvoHLdwTsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vPrH = np.arange(0.00001, 1.0, 0.0001)\n",
    "# TP = Pr(H)*Pr(E|H)\n",
    "vtp = vPrH*PrEH\n",
    "# FP = Pr(H')*Pr(E|H') = (1-Pr(H))*(1-Pr(E'|H'))\n",
    "vfp = (1-vPrH)*(PrEnH)\n",
    "# FN = Pr(H)*Pr(E'|H) = Pr(H)*(1-Pr(E|H))\n",
    "vfn = vPrH*(1-PrEH)\n",
    "# TN = Pr(H')*Pr(E'|H')\n",
    "vtn = (1-vPrH)*(1-PrEnH)\n",
    "\n",
    "# Plot the variations of tp,fp,fn,tn\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(vPrH, vtp, '-', lw=2, color='tab:blue', alpha=0.4)\n",
    "ax.plot(vPrH, vtn, '-', lw=2, color='tab:green', alpha=0.4)\n",
    "ax.plot(vPrH, vfp, '-', lw=2, color='tab:orange', alpha=0.4)\n",
    "ax.plot(vPrH, vfn, '-', lw=2, color='tab:red', alpha=0.4)\n",
    "\n",
    "plt.xlabel(r'$\\displaystyle P(\\textrm{Lie})$')       \n",
    "plt.ylabel(r\"Outcome\") \n",
    "plt.legend([r\"TPR\", r\"TNR\", r\"FPR\", r\"FNR\"])\n",
    "\n",
    "\n",
    "i = 500\n",
    "\n",
    "ax.plot(vPrH[i], vtp[i], 'o', lw=2, color='k')\n",
    "ax.plot(vPrH[i], vtn[i], 'o', lw=2, color='k')\n",
    "ax.plot(vPrH[i], vfp[i], 'o', lw=2, color='k')\n",
    "ax.plot(vPrH[i], vfn[i], 'o', lw=2, color='k')\n",
    "\n",
    "point_x = vPrH[i]\n",
    "point_y = vtp[i]\n",
    "ax.annotate('True positives={:.3f}'.format(point_y),\n",
    "            xy=(point_x, point_y),\n",
    "            xycoords='data',\n",
    "            xytext=(point_x + 0.07, point_y + 0.10),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "point_y = vtn[i]\n",
    "ax.annotate('True negatives={:.3f}'.format(point_y),\n",
    "            xy=(point_x, point_y),\n",
    "            xycoords='data',\n",
    "            xytext=(point_x + 0.07, point_y + 0.10),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "point_y = vfp[i]\n",
    "ax.annotate('False positives={:.3f}'.format(point_y),\n",
    "            xy=(point_x, point_y),\n",
    "            xycoords='data',\n",
    "            xytext=(point_x + 0.07, point_y + 0.10),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "point_y = vfn[i]\n",
    "ax.annotate('False negatives={:.3f}'.format(point_y),\n",
    "            xy=(point_x, point_y),\n",
    "            xycoords='data',\n",
    "            xytext=(point_x + 0.07, point_y + 0.05),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "vPrH[i]\n",
    "vtp[i] + vtn[i] + vfp[i] + vfn[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And also the Sensitivity and Specificity so we can see that they are not good performance metrics to evaluate cases with low prevalence, so it's better to use PPV and NPV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.12001, 0.15550000000000008, 'Specificity (TNR)=0.756')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEFCAYAAADwhtBaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dz4/j6Jkf8O8znnF6B7DNVrkX2OzY3cPCBrBPDqXCHnwIkKHW5+yKU0H6mqGc02CBTMn1D7hWSu5rqnLtADVS9pTLguzjHIKWtIkPNhCsaIzTMQyUq5qeQ68Rw/PkIJJF/X5VJba6q74fQGiRfPnqZVU1H70/+L6iqiAiIlrnrV0XgIiI3gwMGEREZIQBg4iIjDBgEBGREQYMIiIywoBBRERG3t51Acr0zW9+Ux89erTrYhARvVGGw+FvVPXB7P5bHTAePXqEwWCw62IQEb1RROTzRfvZJEVEREYYMIiIyAgDBhERGWHAICIiIwwYW/TkyRM8evQIb731Fh49eoQnT57sukhERFtzq0dJvUpPnjyB7/t4+fIlAODzzz+H7/sAgMePH++yaEREW8GAMeMPX/4BP/3NTzc+75MffZIHi8zLly/xyY8+wXfr391W8YiIjNjfsPGNf/aNrebJgDHjS3yJf0z+cePzfv1/f710/3XyIyK6iT9+948ZMMr2FfkKvvfgexuf9yd/+if41fNfLdx/nfyIiG7iG1/dbrAASgoYItIAkABwVLUzc8wBMAQQp7siVW2KSFtVWyLiq2o3Tbto39K8t+EteQt/dv/PNj6v8zedqT4MAHj33XfR+ZvOtfIjInrdbH2UVBoQoKoRgCTbLqioqqjqPgAPQDvd74vIGFeBZG6fQd478/jxY3S7XTx8+BAigocPH6Lb7bLDm4hujTJqGIcAwvR9DMAFMMoOpjf7jK2q/fS9N3Ns0b6Vee/a48ePGSCI6NYqI2BYAC4L23uLEomIC6AYDBwRAaabmmb3GeVNRETbt8sH9+qqmmQbqtpJaxN7aTBZuG8dEfFFZCAig/Pz83JKTkR0B5URMBIAlfS9BeBiSbq8/0FEGmlnNtL09qJ9JnmraldVa6pae/Bgbjp3IiK6pjICxhkmN3ek/0YAICJWlkBEspt/JsZV89Q+gMGSfQvzJiKi8m29D0NVRyJSS5uQElXNOqWfAqgWkl7OnOOLyCWAcXbOkn2L8iYiopKJqu66DKWp1WrKFfeIiDYjIkNVrc3u52y1RERkhAGDiIiMMGAQEZERBgwiIjLCgEFEREYYMIiIyAgDBhERGWHAICIiIwwYRERkhAGDiIiMMGAQEZERBgwiIjLCgEFEREYYMIiIyAgDBhERGWHAICIiIwwYRERkpJSAISINEXFF5GjBMUdEVETG6StI97fTf/1CWj99tQv75tIREVH5th4wRMQBAFWNACTZdkFFVUVV9wF4ALJg4IvIGECc5uMCiFS1C8BOt+fSERHRq1FGDeMQQJK+jwG4xYNpIMnYqprd+D1V3S8ctwvnxun2onRERPQKvF1CnhaAy8L23qJEWQ2isMsREQBwVLWT1izyYwDOFqXbWqmJiGilXXZ611U1q4kgDRIRgL1C81PWxBWq6mhVukJ6X0QGIjI4Pz9/BZdBRHQ3lBEwEgCV9L0F4GJJurxvI+0kb6SbF7hqfgIAN6tJrEkHAFDVrqrWVLX24MGDG1wGEREVlREwznB1I7eRNjuJiJUlEBEbV/0cwKSPImue2gcwSNP5hWDhLktHRETl23rAyJqO0ht8km0DeDqT9HLmnA/T2sNYVUfp+e106O2LZem2XX4iIlpMVHXXZShNrVbTwYCVECKiTYjIUFVrs/v5pPeOtVotVKtVVKtV9Pv9a+cTRRGazabx/nXH1snO63Q68DwP9XodIgLP8+B5HqJo0nJ4//59NJvNfH+SJEv3x3GMTmezgW+tVgue5210Xva51712ojtLVW/tq1qt6ussDEP1fT/ffvHixVby7fV6pZ7T6/UWprdte26f4zj5++L1zu5vNBqqqtput3U4HBqVYzgc6tHRkaqqHh0d6Xg8njoeBIG6rquu66rjOHneWbpGo6HD4XBhOqK7DMBAF9xTWcPYoSRJYFn5WICp9zdxdna2PtENzgmCAI1GY33CGbVaDXE8/4C+67p5zaPRaCAIAqP8oihCvV4HABwcHOS1mozv+wjDEGEYolarodFowHEc2PZkTEYcx7Bte2E6IprHgLFDjUYDURShWq2i2+1OHcuaefr9PuI4Rr1eR7PZRLVaRZIkiKIob1YZjUaI4xitVgutVitvasqaebJmm+ym3Ol0EEXR0nM8z8tv7J1OB6PRdsYWdLvd/AZfFEUR2u3JDDG2bcO03+ni4gKVymQEt2VZGI/HC9NFUQTLsvJAMRqNsL+/j8PDw6kgPZuOiKYxYOzYcDhEu91Gr9eb6heo1+sIw3Dqm38QBHBdF4PBAL1eD8fHxwiCAI5zNV1Xu91GrVZDEARTN8PDw0N8+umnAIAwDOG67tJzms1m/i0/DMOp/LNv5abiOM77DC4uLnB0dJTv9zwvD4DFzyie2+12F7420W6384AEAI7jYDgcIgzDqRrPbDoimlbG1CC0Idd14bouqtUqAGA8HuPy8hLD4RBJkuCLL77Ib9J7e3tIkiSvGSRJYtSE02g0UK/X4bru2hu+67oIgiCv2RRltRRTlUplYfls20av18NoNILneVPNQFmtIWsuWmZvbw+Xl5d5ufb25mehWVZey7LgeR76/T6Ojo42vi6iu4gBY4dGoxFs24ZlWVPfdKvVKiqVSn4TXdTun91woyhCEARGI34sy0K/3zdKe3BwgFarhdPT07nPzW7S2+A4DhzHQb/fz683yz+O47l+iYzv+3BdF2dnZ3BdF2EYLh0lVqy9FD9nOBzmAXE2HRHNY8DYMc/z8htkr9cDMLkZep6XNxEtaibpdDoIwxAAFh6v1+tz3+wPDw/RarXyZqFl5xQ7gmc74i3L2vq38dPTU3zwwQdznc3rahiO4+Ds7Aye58G27fyGnzXnAZOgs7+/n5/jui6azSYuLy/ngnIxHRHN44N7tFDWkb5oxFDxhrxtSZLgo48+yoMnEb16fHCPjHW7XZycnCwdXtpsNm/0kOEqJycnOD4+LiVvIroZBgya4/v+ym/4jUYDYRhuvWkqy499CUSvJzZJERHRFDZJERHRjTBgEBGREQYMIiIywoBBRERGGDCIiMhIKQFDRBoi4orI3CPFIuKIiKZLr45FJEj3t9N//VX5rMqbiIjKs/WAISIOAKhqBCDJtgsqqiqqug/AA5DNa+GLyBhAvCwfg7yJiKgkZdQwDgFkT3TFANziwfRmn7FVNZtZz1PV/cLxRfmszJuIiMpTxuSDFoDidKbzc04DEBEXQDF4OCICAI6qdpbkY5Q3ERFt3y47veuqms8toaqdtHaxlwaTaxERX0QGIjI4Pz/fSkGJiKicgJEAqKTvLQAXS9Ll/Q9pR3Y2090FAHtJPmvzVtWuqtZUtfbgwYObXAcRERWU0SR1BiCbg8RG2uwkIlZWoxCRLCBk4vQFAPsAAgCDRfks2UdERCXbeg1DVUdA3keRZNsAns4kvZw558O0ljFW1dGifFbkTUREJeNstURENIWz1RIR0Y0wYBARkREGDCIiMsKAQURERhgwiIjICAMGEREZYcAgIiIjDBhERGSEAYOIiIwwYBARkREGDCIiMsKAQURERhgwiIjICAMGEREZYcAgIiIjDBhERGSklICRrtHtisjRgmOOiKiIjNNXMHP8aFU6EWmn//pllJ2IiBbb+preIuIAgKpGImKLiDOzlGpFVaWQNimc6wI4WJPOT5dybW677EREtNxGNQwR+bpBskNc3dxjAG7xoKpGhU1bVeNFmaxI56nq/sxxIiIqmVENQ0T+EsC/BfANAD8QkTNVPVyS3AJwWdjeW5KnCyAqbDtpraS5Kh0AR0QAwFHVjkn5iah8v//97/H8+XP87ne/23VRyNC9e/fw3nvv4Z133jFKb9ok1VTVH4jIT9Lt+9cq3bT6TC2hYpIuCxIiUhcRlzUNotfD8+fP8bWvfQ2PHj1C+qWOXmOqiouLCzx//hzvv/++0TmmTVK/FZF/D+B+WttIVqRNcHXztwBcLEnnZG+y2oVBukbaf4E0X3s2sYj4IjIQkcH5+fmKYhLRNv3ud7/D3t4eg8UbQkSwt7e3UY3QKGCo6oeY1CpeYNIZ/eGK5Ge4upHbSJuTRMQqFNTGdNCxC8HAzjrOF6SLcdU8tQ9gsKCsXVWtqWrtwYMHJpdHRFvCYPFm2fT3Zdzprar/CcARgE9XdX5nI6LSvoekMELq6UzSy8I5fVXtY1IzsVakGwH4MA0s45nRV0R0x7VaLVSrVVSrVfT7/WvnE0URms35gZjL9q87tk52XqfTged5qNfrEBF4ngfP8xBFk+/J9+/fR7PZzPcnSbJ0fxzH6HS23M2rqmtfAD4F8PeY1B4+BXBmct6uX9VqVYno1fjZz362088Pw1B938+3X7x4sZV8e71eqef0er2F6W3bntvnOE7+vni9s/sbjYaqqrbbbR0Ohys/f9HvDcBAF9xTN6lh/EBVD1X1Q10+QoqIaCeSJIFlXTVQFN/fxNnZWannBEGARqOxPuGMWq2GOJ5/KsF13bzm0Wg0EATBXJrrMg0YoYiciMhfZq+tlYCIaAsajQaiKEK1WkW32506ljXz9Pt9xHGMer2OZrOJarWKJEkQRRE8z0Oz2cRoNEIcx2i1Wmi1WnlTU9bM02q1ppqDOp0Ooihaeo7nefmNvdPpYDTaTkt6t9tFvV6f2x9FEdrtNgDAtm0MBnNdvddmPKwWk+aohc9UEBEV/df/8ctS8v13f/7tlceHw2F+wxwOhwiCAJ1OB/V6Hb7vw/M8OM5k4GUQBGi1WhgMBuj1ejg+Ps6PZTf4druN0WiUf0u/vJx0qR4eHuLTTz+F7/sIwxBHR0dLz2k2mwiCAO12O0+bieMYtj032HOpOI7RbDZxeXkJ27bzwBDHcR6Yjo+P4brumpyuxzRgDHTS6U1E9FpzXReu66JarQIAxuMxLi8vMRwOkSQJvvjii/wmvbe3hyRJ8ppBkiRGTTiNRgP1eh2u66694buuiyAI8ppNUVZLMVWpVBaWz7Zt9Ho9jEYjeJ431cRVqSx7xG1zpgHDFpFnKAxjVdX/sLVSENGtsq4mUIbRaATbtmFZ1lTbfrVaRaVSyW+ii9r9sxtuFEUIgsBotJNlWej3+0ZpDw4O0Gq1cHp6Ove5Wa1lGxzHgeM46Pf7+fVuM3/TPowmgA8BdAovIqLXiud5qFar8DwPvV4PAOD7Ps7OzlCv1+F53sLzsmardruNw8P5MT31en0u0BweHiIIgrwZa9U5vu/PdcoDk6CzaS1jndPTU5ycnGw1z4xMRlAZJBT5CJOJBJ+p6n8upTRbVqvVdJsdPkS03M9//nN85zvf2XUxXktZR/qi0VD1eh1hGJbyuUmS4KOPPsqD5yKLfm8iMlTV2mxaoxqGiPwtJlNx/AiTaUI2H2dGRHQHdbtdnJycLB0622w2b/SQ4SonJyc4Pj7eWn6mTVK2qv6dqv5CVU8xmZaDiIjW8H1/5Tf8RqOBMAy33jSV5besyew6TDu9f5s+exEBqGMypxMREW3BNh+uy1iWlQ+73ZZNJh/cB/BfALyvqycfJCKiW8h0AaW/V9UfFLZXLaBERES3kGkfxuwcuNtYQImIiN4gpn0YsYj8R7APg4jozjLtw/ghgN8C+CGAF+k2EdFrpYz1MDqdDprNJqIoWrq+xE3WwsiYrokhIlOTK77SdTEWzXk++wLw9ZntRybn7frF9TCIXp3buh5GtraEqeuun7HJmhjF9S8ajUZ+rddZF6OM9TBmx2ZtfwwYEdENlLUexqbPR1xn/YxN1sSoVCpwXTevcSxTxroYK/swROSvABwCcNP1tY0WgE2XUE0AOKramTnmABjiqh8kUtVm4fhRdo6ItFW1JSK+qnbX5U1Er4fe/17+oNpNeP9i8VxQwOSmeHJygmq1imazCd/3AVxN/V2r1RBFEcIwhG3beZNNs9lEo9FAs9nM534KggBBEGBvbw+DwQDNZhPNZhNnZ2dot9tTacMwRBzH+Q05a55qt9uwLAue56HdbsO2bXQ6Hbiue+OH6bLyrJrGvIx1MVbWMFT1v+nkmYuuTlbc+4vsteycNCBAVSMASbZdUFFVUdV9AB4KtZd0HfCDQlpfRMZIg4tB3kR0hw2HQ7TbbfR6vak+Bdu28zUpimtkhGGIs7MzdLtd7O/vIwzDqXmdjo6OUKvVpr6dL0sLTNbCyNJnNZxsPQxgElxmg8Wma2Jk1zM7K2+WVzYBY5IkW33KGzAfJfVsdpU9Vf27JWkPAWQ/xRiTCQvzJabSm33GVtVVPVPeTPqVeRPR62FVTaBss+thAFdrQmRNM9m04tkaGZ999hk+/vhjo/yHw+FGHdyr1sMANm/yyhwfH889yV32uhjGc0lh8qT3PoC/wGS682UsAMUJ2Beu0pfWJqLCtjMTHADAERFXRLIlqozyJqK7ZzQa5Tff2W/e2ZoQURTBtm1Uq9V8evIwDPH9739/bZ9AplqtGqfNZOthZM1kRdddE8NxHAwGg4UBp7guRmYb62IY1TB0ZrW9dPbam6rPBIi58Ffoy6inAYaIaCnP8/IbY3HCv+ISpk+fPs37FrKmo6wJq16vr52Dyff9qbSLJhas1+t5TSY7JwzDhR3xN1kT4/j4eOkaH6enp/jggw+MO9ONLBo6NfsC8BMAf5u+foLJmhjL0rYBuOn7BoCjJenCwnun8L5XOLeRvj8C4JvknaYbABh8+9vfXjiMjIi2b9fDapcZj8dTw213YTgcrhxu67puqZ//4sWLpcODNxlWa9qH0QagSEdJqeovVqQ9A5AtvGEjbXYSEUtVk/S9jclIp4yd7sveO5j0UWT1yn1MhvIOFuVdpJPRVF1gsoCS4fUREZWi2+0iDMOVU5xna2JstTZQsK11MVb2YYjI10XkE0wCRgDgrzBZSGkpVR2l57oAkmwbwNOZpJeFc/o66fyuYNJPkeXzYTqMdqyqoxV5ExEtlI2Q2pV162EA5a2JAWx3XYylS7SKyPsAegB+DOAfMLmR1zBpHmqo6v+68aeXjEu0Er06XKL1zbTJEq2rmqR+AuADVf1tYd8/iEgPk2anHyw+jYjuKlWFiNHzvfQaWFZhWGZVk5TMBIvsAxIYPvFNRHfHvXv3cHFxsfFNiHZDVXFxcYF79+4Zn7OqhhGLyPdU9X8Wd4rIvwanNyeiGe+99x6eP3+O8/PzXReFDN27dw/vvfeecfqlAUNVfygin4qIAniGyUNy1ckhZXMUEU1555138P777++6GFSilcNqVfXDtPPbAfALTOaUWjWkloiIbqm1z2GkAYJBgojojjOdS4qIiO44BgwiIjLCgEFEREYYMIiIyAgDBhERGWHAICIiIwwYRERkhAGDiIiMMGAQEZERBgwiIjLCgEFEREYYMIiIyEgpAUNEGiLiisjRgmOOiKiIjNNXMHP8qPDeT1/twr52dqyMshMR0WJbDxgi4gCAqkYAkmy7oKKqoqr7ADwAxWDgAjgovI9UtQvATrcBwBeRMbiIExHRK1VGDeMQQJK+jwG4xYNpIMnYqrrsxm8Xzo3TbQDwVHV/Jh8iIipZGQHDAnBZ2N5blCirQRS2nWIQUNVuWrsAJgs4DbL3y5q7iIioPLvs9K6ralLYrixKlDZphao6AgBV7aSBZa/QTFVM74vIQEQGXFuYiGh7yggYCa5u/haAiyXp8r6N2drFDFdVO2m6hog00v0XuGqmyqU1k5qq1h48eHCtCyAionlrl2i9hjMAtfS9jbTZSUSsrEYhIjau+jmASae2XXjvqOpIRPxCsHAx6cvI+jz2AUyNsCIiovJsvYaRNR2lN/gk2wbwdCbpZeGcvqr2MamZWIXz2+nQ2xeFvD9MaxnjQt5ERFQyUdVdl6E0tVpNB4PB+oRERJQTkaGq1mb380lvIiIywoBBRERGGDCIiMgIAwYRERlhwCAiIiMMGEREZIQBg4iIjDBgEBGREQYMIiIywoBBRERGGDCIiMgIAwYRERlhwCAiIiMMGEREZIQBg4iIjDBgEBGREQYMIiIyUkrAEJGGiLgicrTgmCMimi69OhaRYOb40ap8VuVNRETleXvbGYqIAwCqGomILSLOzNrbFVWVQtqkcK4L4GBZPlm6FXkTEVFJyqhhHOIqCMQA3OJBVY0Km7aqxhvkszJvIiIqTxkBwwJwWdjeW5QorU1EhW1nJpgsyscobyIi2r5ddnrXVTUpbFd2VhIiIlqrjICR4OrmbwG4WJIu75NYULtYls/avEXEF5GBiAzOz8+vdwVERDRn653eAM4A1NL3NtJmJxGxshqFiNgodHYDsNN92XtnWT5L9uVUtQugCwC1Wk23cUFERFRCDSMbtZT2USSFUUxPZ5JeFs7pq2ofk9qDtSyfFXkTEVHJRPX2fgmv1Wo6GAx2XQwiojeKiAxVtTa7n096ExGREQYMIiIywoBBRERGGDCIiMgIAwYRERlhwCAiIiMMGEREZIQBY8darRaq1Sqq1Sr6/f6184miCM1mEwDQ6XTQbDYRRRE6nc7a9NdV/DzP81Cv1yEi8DwPnuchiiYP4osIut1ufp7neUiSyYP+9+/fR7PZzM9JkgRxHC8t9zKtVgue52183qzRaDRVFiIqUNVb+6pWq/o6C8NQfd/Pt1+8eLGVfBuNxkbpe73exp/R6/UWnmfb9tw+x3HUcZyp8mXXWtwfhmFe9na7rcPh0Kgsw+FQj46OVFX16OhIx+Px1PEgCNR1XXVdVx3H0Xa7rePxWB3H0UajsbQ8RHcVgIEuuKeyhrFDSZLAsqx8u/j+pvlu4uzsbOPPCIIAjUbDKG2lUoHrunmNYxnXdfOyNxoNBEGwMn0miiLU63UAwMHBwdzn+L6PMAwRhiFqtVpebtd10ev10Ov1YFkW+v0+Dg8PjT6T6C5iwNihRqOBKIpQrVanmmziOEa1WkWz2cT+/j7ieLLGVNbskzVdNZtN1Ot11Ot1xHGMVquFTqeDwWCAZrOJ0WiEVqs1lzb7jFarhVarlTdPZTdrz/Pyz+x0OhiNbj5lV7PZRLvdXpkmiqI8jW3bMJ3W5eLiApXKZBJjy7IwHo+X5m9ZFmx7Ms9lv9+H53l509qzZ88wHo/znzMRTStjtlrawHA4zG+Uw+Ew/1Zt2zaCIEC/30cQBNjb20O9Xofv+/A8D5eXl9jf38/TZzf4o6MjhGGIIAjyG323251KW9RutzEajaaONZtNBEGAdruNMAxxdDS9fHocx/lN15Rt27AsKy9nMa8sQB0fH8N15xdRjON4ae3E933jMmTXA0xqPe12G41GY6q/xbIsBEGATqeDKIoWloformLAeA24rgvXdVGtVvN92TfmrGnGtm1cXl5iOBwiSRJ89tln+Pjjj43yHw6HG3Vwu66LIAgQx/HCb9rX7Qw+Pj6eq2XYto1er5d3NhebubKfgW3bKwPD3t4eLi8v87Lt7c0vxDhbZsuy8s86ODhAHMfY39/PA6Ft2xiNRgwYRAUMGDs0Go2WfvPOboBRFMG2bVSrVVQqlfwm1+12EUURHMeZy3dWtVo1Tps5ODhAq9XC6enp3LEseG3KcRwMBoM8EMwecxwH/X4/v8bsM9bVMFzXxdnZGVzXRRiGC4Pj7PUXaw/Pnj3D8fFxXrtwXfdatSii244BY8ey5iUA6PV6+f5iU83Tp09hWRY8z0MQBLAsC71eL++XsCxrZf+A7/tTaYufk6nX63lNJjsnDMOFHfGWZd2oluF53sJjp6en+OCDD+Y609fVMBzHwdnZGTzPg23beWCo1+t5E1RWgyjmmQ2dzYJVxvM8VCqVuaY4oruO62G8huI4RrvdNh4lVIbRaIQ4jpeOhCrejMuQJAk++uijhcGNiMrF9TDIWLfbxcnJycphs81m80YPGq5zcnKC4+Pj0vInos0xYLyGshFSu+L7/tpv9o1GA2EYlvI0dJbnJn0uRFS+UgKGiDRExBWRuUZgEXFEREVknL6CdL+bvtpr0mXHzcdTUimy/pRtW9cnQ0S7sfVObxFxAEBVIxGxRcRR1eKTXxVVlULaJP23rqotEWml23Pp0vN9EWkAuNlESEREtJEyRkkdAsh6Q2MALoA8YKhqcXykrapZQ/iosG/20eJiOm8mDyIiegXKCBgWgOIg/fmnqDBpggIQzew7wkzNYUE6R0QAwFHVm01NSkRExnbZ6V1X1ake0zQANEXEWpZOVTtpDWMvDSZTRMQXkYGIDM7Pz0srPBHRXVNGwEgAZI/yWgAulqTLh8CkHdzZdgzAX5KukfZfIM137lFcVe2qak1Vaw8ePLjmJRAR0awyAsYZrm7kNtLmpGKtQURsXHViA5N+jmKQiZeki3HVPLUP4M17Ko+I6A219YCRdVinzUVJoQP76UzSYj9HF4Cd1R4KHdxT6dK8PkzTjRd0jhMRUUk4NQgREU3h1CBERHQjDBhb9OTJEzx69AhvvfUWHj16hCdPnuy6SEREW8PpzbfkyZMn8H0fL1++BAB8/vnn+ZTcjx8/3mXRiIi2ggFjxu//8CX++09/tfF5f/1JKw8WmZcvX+KvP2nh3e/+q20Vj4jIyJ+/v4d/bv3RVvNkwFjgn/7flxufc/7rxUHm/Ne/ulZ+REQ38Ycvtz+giQFjxttvCf7Nv/zTjc/71re+hf/zy18u3H+d/IiIbuKrb2+/i5oBY4aI4I+++pWNzzv58Y+n+jAA4N1338XJj398rfyIiF43HCW1JY8fP0a328XDhw8hInj48CG63S47vIno1uCDe0RENIUP7hER0Y0wYBARkREGDCIiMsKAQURERhgwiIjIyK0eJSUi5wA+v+bp3wTwmy0W503Aa74b7to137XrBW5+zQ9VdW7J0lsdMG5CRAaLhpXdZrzmu+GuXfNdu16gvGtmkxQRERlhwCAiIiMMGMt1d12AHeA13w137Zrv2vUCJV0z+zCIiMgIaxh0J4hIQ0RcEVrflOIAAAO1SURBVDlak27lcaLXkYg4K44Z/e2bYMDA+h/oNn/grwOD6/XTV/tVl60M2X8mVY0AJMv+c4mIC+DgVZatTAa/ZydN03jVZSvLBv+X/VddtrKkf7enS44Z/e2buvMBY90PdNs/8F0zuF4XQKSqXQB2uv2mOwSQpO9jALfhmlYy/Lttqmofk9/zG/13DRj/X47T4/FtuGYgv97LJYe3+rd/5wMG1v9Ab9vNZt312IV9cbr9prMw/R9qbzaBiDjpf7zbYuXvOa1VjAFAVTuqOnq1xSuFyf/VrNZs35JrXmft3/4mGDDW/0C3+gN/Day8HlXtprULAHAA3JUFRSq7LsCWrfu7PQCwlzZL3YqmVqz/2x5hUrMYY/k3clqBAYMWSqvr4S35FpbgKiBYAC6KB29h7cLURfb7vU39GMuIiIXJ30IA4FREbkPteZ2Vf/ubYsBY/wPd6g/8NWB6Pa6qdl5NkUp3hqumNRtABOQ3EGDShp91/t6K9nys/z2PMWm2QfrvbejsX3fNPoCT9O/aA3Brg2Thb3vh3/51MWCsv5ls9Qf+Glh3vRARPwsWt6HTu/At2gWQFGpNT9Pj/bTzt4LJjeY2WPd7jmaOP3ulpSvH2r/tTNYx/uqKVp70i05tppaY/W0v+9u/3mfxwb3JDRJpB2/Wfp+uaVtddvxNtup60z+sHiZtvBUA3h1trnnjGf5dX6bHb0Vt0uCaj9Ljldvwf/lVY8AgIiIjbJIiIiIjDBhERGSEAYOIiIwwYBBtQTr31jgdnuuLSDBz3BURK31QbiwiQfE5ABGxRaS3JG/rlgz1pTccAwbRdgwwmYOrn46+SbIJ7rJhnaqaDWscAQhUNXsOAun7jxZlrKoJbscULfSGY8Ag2g4XwLCwXZymwl83NDmtbawKCqO78DQ2vd4YMIi24wDpvFtpjcJOHwYEgH3DPPLp5EXkqDgNd1oDqW+xvEQbY8Ag2g4HV9PBu6q60c09DQgJAKTrkIzSWkkx2Ny2CRLpDfP2rgtA9KYr9FH016Vdcn5j5lwbQNbRPS7s5wyrtFMMGEQ3V8PN5hib7bt4hslCPyMRiRedQLQLDBhEN5B2VrdgOJFdWmtwABym5x7iavU3R0RsVe2kfRhZExTn8qLXAueSIipZOvvvjSa6S4OLc91mL6JtYKc3Ufk+3cKQWAYL2jkGDKKSpQ/eJYvWZTCR1i7Yl0E7xyYpIiIywhoGEREZYcAgIiIjDBhERGSEAYOIiIwwYBARkREGDCIiMvL/ARwz28S/lv58AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Also TPR\n",
    "sensitivity = vtp / (vtp + vfn)\n",
    "# Also TNR\n",
    "specificity = vtn / (vtn + vfp)\n",
    "\n",
    "\n",
    "# Plot the variations of tp,fp,fn,tn\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(vPrH, sensitivity, '-', lw=2, color='tab:blue', alpha=0.4)\n",
    "ax.plot(vPrH, specificity, '-', lw=2, color='tab:green', alpha=0.4)\n",
    "\n",
    "plt.xlabel(r'$\\displaystyle P(\\textrm{Lie})$')       \n",
    "plt.ylabel(r\"Outcome\") \n",
    "plt.legend([r\"Sensitivity (TPR)\", r\"Specificity (TNR)\"])\n",
    "\n",
    "\n",
    "i = 500\n",
    "\n",
    "ax.plot(vPrH[i], sensitivity[i], 'o', lw=2, color='k')\n",
    "ax.plot(vPrH[i], specificity[i], 'o', lw=2, color='k')\n",
    "\n",
    "\n",
    "point_x = vPrH[i]\n",
    "point_y = sensitivity[i]\n",
    "ax.annotate('Sensitivity (TPR)={:.3f}'.format(point_y),\n",
    "            xy=(point_x, point_y),\n",
    "            xycoords='data',\n",
    "            xytext=(point_x + 0.07, point_y + 0.05),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "point_y = specificity[i]\n",
    "ax.annotate('Specificity (TNR)={:.3f}'.format(point_y),\n",
    "            xy=(point_x, point_y),\n",
    "            xycoords='data',\n",
    "            xytext=(point_x + 0.07, point_y - 0.6),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thermal analysis to screen refugees \n",
    "\n",
    "Last December the University of Granada published an interview with one of the authors of a system to spot liars based on *the Pinocchio Effect*. The methods use thermography to detect an increase of temperature in the nose of people who are lying while answering a question. The authors propose to use the tool to screen persons at refugee tamps to know \"*what is the real purpose of people trying to cross borders between countries*\" [7]. \n",
    "\n",
    "We can get the TP and FP from the journal publication [6]:\n",
    "\n",
    "- $P(\\text{+}\\mid\\text{Lie})=0.85$\n",
    "- $P(\\text{+}\\mid\\text{Non-lie})=0.25$\n",
    "\n",
    "From the interview we don't have clear what the authors mean with *the real purpose of people* at refugee camps. Let's assume this purpose can be, for instance, to spot terrorists that would be hypothetically hidden in a refugee camp or are trying to enter Europe. So we want to know whether it would be useful to screen all the refugees using a lie detector to spot terrorists. \n",
    "\n",
    "We can calculate the prior of refugees that might be related to terrorism by using the data provided by the conservative think tank 'The Heritage Foundation' which claims, with no sources, that there where 44 refugees or asylum seekers involved in terrorist plots in Europe from 2014 to 2017. Note this is an unrealistic scenario with more people related to terrorism that actually are. In that period, we had 3,915,445 first time asylum applicants in Europe [9]. Then, we have the prior of refugees that during an interview will answer they are not terrorist, but they are: \n",
    "\n",
    "$P(\\text{Lie})=\\frac{44}{3915445}=0.000011$\n",
    "\n",
    "So we can update the relative confusion matrix:\n",
    "\n",
    "| Test Result\\ Ground Truth | Lie (0.0011%) | No-lie (99.9989%)|\n",
    "| ------------------------ | ------ | ------ |\n",
    "| Test positive            | 0.000011x0.85=<br>0.00000935 | 0.999989x0.25=<br>0.24999725 |\n",
    "| Test negative            | 0.000011x0.15=<br>0.00000165 | 0.999989x0.75=<br>0.74999175 |\n",
    "\n",
    "We can build the absolute confusion matrix for the population of ~4 million asylum seekers: \n",
    "\n",
    "| Test Result\\ Ground Truth | Lie | No-lie |\n",
    "| ------------------------ | ------ | ------ |\n",
    "| Test positive            | 36 | 978850 |\n",
    "| Test negative            | 6 | 2936551  |\n",
    "\n",
    "Therefore, to discover one terrorist, ~2719 persons would need to be investigated after been labelled as liars by the system.\n",
    "\n",
    "We can calculate the posterior probabilities of having an actual liar when the thermal analysis test is positive (PPV) and having an person answering the truth when the test is negative (NPV): \n",
    "\n",
    "\\begin{align}\n",
    "P(\\text{Liar}\\mid\\text{+}) &= \\frac{0.85 \\times 0.000011}{0.85 \\times 0.000011 + 0.25 \\times 0.999988} \\\\\n",
    "&\\approx 0.0037\\%\\\\\n",
    "P(\\text{No-liar}\\mid\\text{-}) &= \\frac{0.75 \\times 0.999988}{0.15 \\times 0.000011 + 0.75 \\times 0.999988}  \\\\\n",
    "&\\approx 99.9998\\%\n",
    "\\end{align}\n",
    "\n",
    "In relative numbers, this means that *only 9 terrorists hidden in a population of 1 million will be detected by the lie detector* whereas 249997 innocent people would need to be unnecessarily investigated. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the previous code to see how the model performs under different hypotheses (we highlight the case of 0.5% of terrorists, which does not correspond with real frequency)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00501"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEFCAYAAAAPCDf9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1zb1f348dfh0tIrAUpr722wttXeDOBtXqYN3t2ckuJ0NzcL7n75bsXO3+4XB+6mu7jQOd28jcKcm5eppF6mdioQq9ZW25LeL7aFhl4pt/P745OEAAFCSfiQ8H4+HnmQfD7JJ+8EeOfkfM55H6W1RgghxPCQYHYAQgghBo8kfSGEGEYk6QshxDAiSV8IIYYRSfpCCDGMSNIXQohhJMnsAPoyYcIEPWvWLLPDEEKImFJbW3tQa53ZdfuQT/qzZs2ipqbG7DCEECKmKKW2h9ou3TtCCDGMSNIXQohhRJK+EEIMI1FL+kopWy/78pVSdqXUimg9vxBCiO6ikvSVUnZgVQ/7bABaaxfg7e3DQQghRGRFJen7EnpDD7sLAK/vugewRyMGgMYnn2TzZUvZOP9MNl+2lMYnn4zWUwkhREwwY8imhc4fCBnReJLGH9/E3or16JY2AFr37GHv974PQOp110XjKYUQYkCON7eyce8R1u9uZP3uRpZfbOWMSeMi+hxDcpy+UqoQKASYMWNG/w/QepL9T9SiWzq/PN3UxP4fF5N6+CFInw1pszt+WmZAckokwhdCiD4dbmphw57DrN/dyHu+n3UHjtIetMTJ2TPS4iLpe4F033ULUN/1DlrrMqAMICcnp/+rvGhN64nQL631SDtsfi7EHgXjp/g+CGZ1/kDIyIKU1H6HIYQQAEeaWnh3dyPv7Grk3d2NvLe7kW31x7vdLylBMfe0cSyYMp4FU1M5PyvyHSGDlvSVUhattRcoB3J8m62AK+JPlpxC0uQptO7Z021X0qRMWFYKh7bCoW3QsNW47t0Jh3cbl+2vdj/mmEzImGN8AGScDhPmGD/TZkHSyIi/BCFEbGpqaWPj3sO8s6uRt3d5eWeX0YLvukjhiMQE5k0ex1lTUlkwdTwLp6ZyxqRxpCQnRjW+qCR9pVQ+kKOUytdaV/o2rwGytdZupVSOb4SPV2vtjkYME7/5DfZ+7/vopqaOuFJSmPjtFXBmiD79thZo3NnxIdAQ9KHQ4IFjB4zLjrVdXmwCWGYaHwAZp8OE0zuuj5sCCTIVQoh41dau2bz/CO/s7Ejw7+87TEtb5wyfnKiYP3k8i6alsmiqhQVTU5kzaSzJiYOfH9RQXyM3JydHn2rtncYnn2T/b35L6969JE2ezMRvfuPUTuK2txvfAOq3dL94d4BuD/245NGQnmV8EGTOg8y5xs/0LEgacUqvSQhhnj3eE7h3HGLdDiPBr9/TyPHmtk73UQpOzxzLomkWFk9PZfE0C/Mmj2NkUnRb8F0ppWq11jndtsdz0h8UrSeNbwP1W6B+s+9nHRzcDMcPhn6MSjS6ifwfAv4PhIw5cjJZiCHiZGsb63cf5q0dh3DvOIR7u5d9h5u63W9a2igW+xL8omlGK37sSPPHyPSU9M2PbAgoLS2luroar9eLy+UiPz8fgKKiIuz2PqYRJI2EifOMS1cnDkG9Bw5+AAf8l/eNbqODm4zLxqC5AyrBOEcQ/K0gcy5MmAsjRnc6tMvloqKiAqfTObAXH6SoqChwPKUUTqeTwsJCABwOB6tWrcJisZCWlsayZctoaDBG3q5atQqHw0FJSQk2W8dcu9LSUmw2G1arlcrKSlasCH8CdnFxMR6Ph9zc3JCP8++3Wq2UlJT0a5sQoexrbPIldyPJr999mOa2zt/ix6ckcfaMNM6eYWHxdAuLpqaSMTa2zulJ0odOSSUrK4uKiorIHHhUGkzLNi7BWk4Y3wT8HwIH3jeuN3g6Lh8803F/lQDpWVRuHUv+x6+BSWdhzz4L+9I/RSZOoLKykry8vMBtm83WKekHs1qtgQ8Hl8tFcXExRUVFlJeXd0r6VVVVnd5bt9vdaX9P3G7jNE9FRUWnpO3ncrkC+4uKivB4PHg8nrC2BR9HDF+tbe28t+cwtdsPBRL9nsburfg5E8eSPTMN24w0bDMtWCeMJSFBmRBx5EjSN0PyKJi8yLgEaz1pdA35PwT8HwgHN0P9ZsqfO07++Lc67j9iHEw6Eyad5bssgIlnQsr4fofkdDqpqqoK3E5PT8dms+FyuXr9tpOTk0NJSQn5+fkUFxcHWtMejweLxRK4X35+PiUlJWF9M3G5XIEPoNzcXFwuV6cPn6qqqsD+7OxsXC4XdXV1YW0L9SEm4l9TSxtv7fBSva2B6m0N1G4/1K0vftzIJJbMsPgSfBpLpltIHZVsUsTRI0m/Fx6Ph5KSEhoaGvjyl7/MY489htPppLS0FKvVSn5+Pg6HA6/XS1FRUaBbKPjxRUVFWK1WampqWLNmDRaLpdtj/LetVitZWVkUFhay/EeP4fV6cTi+QOGdn6X467fj2rmaojetlFw3hYbt7+Gs2ovH+yKrrnsdS4qi9LWT2CYnYrdZcTzWiLd1BEWfcWCZtQjnw4+Tnp5OUVFRWK1tMLp7+uriKisrCyRWm80WaM1XVlZSVFQUuJ//PQhHfX096enGVA6LxUJ1dXWn/RkZGXi93sDturq6sLeJ4eFwUwu12w7x5rYG3tzawDu7vN1G1MyeMIbsmWmBlvzpE8eSGOOt+HBI0u/D6tWr2bp1a6D/OlhpaSl5eXkUFhbicDi6JX0/p9NJcXExNTU1uN3uTo+xWCyBfmuHwxFoCfu7mLKzsyksLKTkDw/g3rQL5zNGa7zB44GTv6EgK4PVOz+gcLGFqvK/s+LiJEqf3kzeREVh9ggcZb8kPUWx8sJMbDnnQMOT8N4OmLzYOH+gVI/dHlarFYvFEugm8fN/mDU0NHTqKy8oKAh08XTt2un6eH8XTVfhtMTtdjtOp5P8/HyqqqrIzc0Ne5uITweOnKTal+Df3NrAxn2HO42LVwrmTx7PubPTyZ2VTu7sNCaOG56DJiTp92HZsmVYLJaQSb+uro6GhgZqa2vxer14vd5OXRpAIJn6W51dH5OUZPwKvF5voHULxgdKfX19t4TbyYjR5H/lh+Tl5WEv+hnWi9vhu7+nbvOnaNi3g9rqZrxqB3d+JImfufbjfepJnNe6sKb5xganpMLkxXiPToSGHUY3UnpWp7kFK1eu7HYCND09PWQ3jb+LZ+XKlSE/RPyvz2q19prcMzIyAu+31+slI6PzrET/NxX/tyir1Rr2NhEfDh1r5nVPPWvr6vmfp54t+4922p+cqFg0zULurHTOmZ1G9sz0uOyqORWS9PvQNYmDkeytVivZ2dmkp6f32MIPpetj3G43VVVVVFdXB5Krv/toxYoVVFZW9na4QIyB7pTEZLIvXNo5Lq2p8O7A9fhfcVZVUZI7CfauMyabbf0v1iZNw4YT8Psc4zzBlCXGuYWNT2I7PZuamppOH0i9sdlsFBcX43A4uu3zJ/K+Wvp2u53y8nLsdjtVVVWduon8/B86/lFD/dkmYs+Rphbe3NpgJPm6+m4t+VHJiWTPTOMcX0t+yXQLo0YM7rj4WCFJP0z+Pmn/KBAg0EXjdDqxWCxhjfrp+piVK1cGTnr6T4Ta7XaWL1/e6cSqX15eXrdWdkFBAcXFxYHulK7PkZubGzhWSck9YLOB1nBkH+xdh2Xv23if/i2MTzcmoW17BRqPQfmnAFh51lgcD7nhf3+EeRfRbT55kKKiok7dVD29l7219G02G+Xl5Tgcjk4t9ry8PKqqqnznOhyB9w8Ie5uIDSea26jZbiT5tXX1rN/dSFtQJbIRSQlkz0jjgqwMzs/KYNE0CyOSZPZ7OGRylsn8LVCr1RqVsffh8idUju6H3bWwq8b4udsNJxs731klQOZ8mGqDaTkwLde43UvJCa/Xy/LlyyM3HFbElbZ2zfrdjfx30wFe2XyQt3Ye6nTiNSlBsXi6JZDkbTPSol6jJtbJjNwhyu12d0ry/RldE0n+bqRuXVXt7cYs4901HR8GH66H9tbO9xuZCtPPgRnnwozzYYqt04Sy4uJiCgoKTHltYmja23iCVzYd5OXNB3hty0G8x1sC+5SCBVNSA0k+d1Y6Y4bALNdYIklf9KmoqIiSkpKQ5zE6aTkB+971fRuogZ3V0Lij830SkmDyEphxHl7LAu76+1pKfvO76AUvhrwTzW28sbWe/246yCubD7C5y8nXaWmjuPiMTC6eM4HzrRNIHS0nXgdCyjD0QSnVqZXr73IJ5vF4cDqdQ+qEYCS7hPzH6FqOIdT7UvTT+wPnNpzOf2HNGAk7XyftvFtYdnYaDQ31wH9ZdV01DSc0GRta4N7XjG8BM8+HWRcaQ0Z7EG4ZhuD9PT0m+PWIwaO1ZvP+o7z4/n5e2XyQN7c10NzaUdZgzIhEzs+awMVnTOCiOZnMyhiNUvE/Tt5skvR9rFZrzPQ3V1ZWBhKx3W7vuz5QP48dXI4h1PvicDhYuXIlNpsNr9fL0qVLjYlnC27EOvfnOF+phabDuFY7WX7/36j4/Gz44EXc723C1lAH6x42DpQ6w0j+sy+CWReBZTrQdxmGUPv9k7CCt6Wnp1NcXMzq1asl6Q+SppY2/uep58X397Nm4352e08E9ikFi6alctGcCVw8J5OzZ6TJyVcTSNKPQeXl5f0aJtofXcsxdOVv3fv75v0jY8rKyjq3yFPGY//Mdyh56Hn47L/Jv+ADSn70XZxXLIXtrxmXxh3w9qPGBYx1CWZfhOtFL3mXGiWwQ5VhCFWmwev1hizd4HQ6w54JLE7N3sYTvPD+fl58fz+vbjlIU0tHaz5jzAg+Oncil8zN5MLTJ5A+RkqKm02Sfg/8o038w/66Jh2n09mprEHUyjF0iaG4uBiXyxXof29oaMDpdOLxeAJVMP3VLe12e6fnsFgs3eLuL7fb3W1mq3+IZVculyvQFWadM5eajdvg/C8Zl/Z244TwtleNIaLbXgPvdnhrO/XVTaSf+BdsL8HinUn19gT4tMMoYEfPZRp6K90gIqetXbNu5yFe8LXm3993pNP+BVPHc9nciVw2fxKLpqbGfIGyeCNJ38fj8QQmFPlLC3QtheBXUVER6N6AQS7HUFLSacSPf8JTQUEBq1evprCwMFACoWtc6enpneIO9R6EOo8R/L7k5uZSX995WWP/bOTg+3s8HlauXBmy66nz5Kwk4FJIuwQO76LwwklQ+wdIajCqjXo2QV0rlFqNEUFZl0Hjbmhr7XZcET1NLW28uvkgz2/Yh2vjfhqONQf2jR6RyIWnT+CyeRO5dN5EJo0fnuUNYoUkfZ9Qfdc9lUIoLi6muLgYr9eL0+k0vxwDxlDLvLw87HZ74Lm6Psedd97Jz372s0DcXRN8cHGynt4X/8nsYDU1NYGuFf/93W53tw/AsMswXHSChsWL4awJeB+4h4ymalAHfMNGa8ioO0nDn56DusvxbhtPxqjxMLr30g2i/xqPt/DCBx/y/Hsf8vKmA52qUs5IH81l8yZy2byJnGtNH/RVocSpk6Tfg95KIfgTm7+bZ0iWYwjxHECnuLuOQrJarSFrDHW9DxAouez1eikpKaG2trbT/Ww2W6Dapv/5+12G4YoSqnaNouiOR+CsM4zzAFvWYD/+FOWv1WGf/hRVT56gKGcEjJtM+Z8s2Gdqqp57lqIvfqnP90t0t6+xiaoN+3juvQ953VNPa9As2IVTU7n8zElcseA05kwcKyNtYpQk/R70VgqhtLQ0qKxBSaBPP7bKMXQfdmqxWEK29rvy18PxH6OioiLk2P5Vq1axdOnSbh+Gp1yG4cu/pKqqCtvVpZR//Ys4Xq3BOu0EttmHoWk/5Xt34Lj+aqwZI7F9cBz0FRT93hU4p1JcXCxF10LYdvAYz6zfy3PvfcjbOzt+/4kJiguyMrj8zEnknXUaUy2jTIxSRIpMzhoihlw5hgiLahmG9jajXMTm52HTs7Dvnc77Jy+BM66EM64wrvdSLmK42FF/nKff3cvT7+5h/e7Dge0pyQlcPCeTK846jcvmTSRNRtvELJmRO8QN+XIMAzSoZRgad8Pm52DTc+B5CVqDlsEbO8n4AJj/MZh9MSQNn6S2s+E4z7y7l6ff3cs7uzrqKY0dmYR9/kSuWjiZi+dkSnXKOCFJX4Qt7HIMYfJ6vdx1113mzGRuPm4MCd30rPEhcHh3x76RqTD3Sph/HWQt7bb4fDzYdeg4/3l3H0+9u7dT182YEYnYz5zENQsnc/EZmVK8LA5J0hdCa6Nm0PtPwcYnYf+Gjn3Jo+F0O5z5cZhz+SmtMzxUNBxr5ul39vDEuj3Ubj8U2D56RCJL5xuJ/qNzJdHHO0n6QnR1cAts/LfxAbDH3bE9cQRYPwoLboR518DIcWZFGLYTzW24Nn7IE2/t5uVNBwKjbkYlJ3LZ/Ilcu3AyH507UbpuhhFJ+kL0xruz4xvA9rWA7/8iKcU4Abwg3/gGkDx0Jh61tWv+V1fPE+t28+z6fRw9aUxYS0xQXHj6BD5x9lTyzpwkJYmHKUn6QoTr6H7jG8C7/4Adazu2jxwP866FhTfC7I9CojnJdMOew/zzrV38a90e9h85Gdi+eFoq1589lWsXTSFz3EhTYhNDhyR9IU5F4y5Y/zisr4S9b3dsH50BZ30CFt9srCAW5YlK3uPN/GvdHlbX7OS9PR1DLGdmjObjS6Zy/ZIpWDPHRjUGEVsk6QsxUAe3GMn/3Uqo39yxfcJcWHIzLCqA8ZMj9nRt7ZpXtxxkdc1Oqt77kOY2o3pl6qhkPrZ4Cp+wTeXs6RaZGStCkqQvRKRobUwAe2c1vFMOxw4Y21WCMfRzyc0w9+pT7v/fXn+MytpdVNbuYm+jMcdAKbhoTiaO7GnknTlJRt6IPknSFyIa2lpgiwvWPQIfPAvtvnVeU1KNk7+2T8OUs/s8zMnWNp5dv49H39jBG1s76h/NSB+NI3saN2ZPY4qUQRD9MKhJXymVD3gBm9a6tJf9Vq11WW/HkqQvYsaxeqP7Z90jnfv/Jy+BnM8bQ0BHdu5333bwGI+9uYOK2l2BcsUpyQlcvXAyy3Kmc86sdKlHL07JoCV9pZQNI5lXKqUKgRqttbvLfrTWbqWUHWgI3t+VJH0Rk/atN5L/ukehyTcTdsQ4WFxAy9mfZU1DJo+8sYNXNh8MPOTMyeO55bwZfGzxFMalyKLgYmAGc2H0AsBfscsD2IGuSb0EyMP4cAhdY1eIWHbaArjyLlj6fdjwL6h5AHa+DtV/Jrn6z2S2z2Fi61LGJ13AFYtncct5M1k8LVVOyoqoi0bStwDBRdk7rWbha+F7lFJ1QHEUnl+IIUMnpfC/sXYeSD6dnc013JSwhhsSXyE7YTPZIzbTPmo1CWmfg9TloCJT60iI3gz67BKllAWjP98JrFJKubXWni73KQQKAWbMmDHYIQoxYE0tbTzx1m4eXLstsIZscuIMas/6LmfaJpB77EVU9f0k7F0Hr/4a1t5r1P0570swrds3ciEiJhp9+iVAldba5Tthaw0+mauUWgGUaa29vj79kCd7/aRPX8SSvY0neOh/23nszR0cOm6M5MkcN5JPnTuTm8+d0XmmrNaw80144z7Y8G/QvuUIp+bAeV80PgQSpW9fnJrB7NMvB/xPZAVcvgAsWutOyzL5PhhkKSMR07TWuHcc4i+vbePZ9fto8xU7WzQtlVs/MotrFk5hRFKIhVuUghnnGpfGXfDmKqh90FgL+B9fgOe/ZyT/nFtjouibiA3RGrJZiHESNzAk0/epk+27vsK3P12GbIpY1dauef69ffzpv55ArfrEBMVVC07j1o/MwjYjrf8nZpuPGRO+Xv8THPzA2JaSCrnL4dzbYWxmhF+FiFcyOUuICGlqaeMf7l2s+q+HbfXHAbCMTubmc2bw6fNnMjk1ApOotIbNVfDab40F4cGo+Hn2p+GCr0LazIE/h4hrkvSFGCDv8WYefn07D67dxsGjxkSq6emjWH6RFUf29OjVqt/xhpH8P3jGuK0SjYleF30LJs6PznOKmCdJX4hTtOvQce5/dSvl1Ts53mycbF0wdTxFF2dx1YLTSEocpIXW92+E1+6BdyugvRVQRqXPS4ph4rzBiUHEDEn6QvTT1oPH+MOLW/jnW7sDJ2cvmjOB2y/J4oKsDPMmUnl3wGv3gvuv0NYMKFhwg5H8M+eaE5MYciTpCxGmLfuP8PsXtvDvt/fQro2Ts9cumkzhxVbOmpJqdngdGncbY/zdf+tI/gvzjeQ/YY7Z0QmTSdIXog8b9x7m9y9s4Zn1e9EakhIUN9qm8aVLs5iZMcbs8HrWuAte8SX/9hajxPOim+DS74JlutnRCZNI0heiB+t3N3Lvms08v+FDAEYkJuDImcYXP5rFtLTRJkfXD96d8Mqv4K2HjeSfOBLOLYQLvwWj082OTgwySfpCdLFx72F+9fwmXBuNZD8yKYFPnjODokuskRl2aZaGrfDCT40yz2CM87/wW3BuESTH8OsS/SJJXwifrQeP8ZuqTTz5zh60hlHJiXzqvBksv9jKxHGnttrVkLTnLaj6AWx92bg9fqrR5bP4k5AgK2/FO0n6Ytjb7T3Bva7NVLp30dauGZGYwM3nzuBLl2bFV7IPpjXUvQCuH8C+d41tkxfDlSUw83xzYxNRNZi1d4QYUvYfaeKPL9bx6Bs7aG5rJzFBUZAzna/Z5zA13pcgVApOXwrWS43x/Wt+ZKzq9cCVxnKOeT+C1GlmRykGkbT0Rdw6erIV58t1/PmVrZxoMSZVfWzxFL5hn4M1c2wfj45TzceMCV6v3QOtTZA8Gi78plHaQfr744p074hho7Wtnb9X7+S3rk2Bcgn2+ZP4v8vPYP7k8SZHN0R4dxhVPDc8YdxOnWGs9DX/WnPjEhEj3Tsi7mmteeH9/dz1n/fZsv8oAGfPsHDn1fPJmSVDFjuxzIBlf4Vtr8J/7oAP34XyW2Du1XBVqYzvj2PS0hdx4d1djfzsmQ287jFW6pyRPpriK+dx9cLTZN3ZvrS3QfX9sObH0HwEksfApSvh3C9CorQLY5V074i4tK+xiZJn3+efb+0GIHVUMl9bOodPnTeDkUkyLLFfDu+FZ+/o6PI5bSFcew9MyzY3LnFKpHtHxJWTrW3c/+pWfv/CFo43tzEiMYHPfWQWX/7o6aSOliUGT8n4yUaXz6bn4en/M4Z4/nmpMalr6fdhxBAuRSHCJklfxJw1Gz/kx09tYLtvAZMrzprE/7vmTKanx1DJhKHsjMth1uvwcgms/T288SfY9Bx8/A8w6yNmRycGSLp3RMzwHDjKT57awIsfHADg9Ilj+cF1Z3LRHFlCMGr2rIMnvgT73wOUtPpjiPTpi5h17GQr976wmb+8upWWNs24kUl8I+8MPnP+TJIHawGT4ay1Gf57t1HMTbdB2my4/o8w8wKzIxO9kD59EZOef28fP/z3e+xpbEIpKMiZzneunMuEsSPNDm34SBoBl90J867paPU/cDVc8BW47PvGfhEzJOmLIWmP9wQ//Pd7gXLHC6em8pPrF7BkusXkyIaxKUug8KWOVv/a34HnZbjxfsg8w+zoRJgk6YshpbWtnQfXbuM3VZs41tzG2JFJfPvyM/j0+bNITJDx9qbzt/rn5ME/boN974DzYrjy55B9q1HrRwxp/Ur6SqnxWuvD0QpGDG9v7/Ty3X++y3t7jD+xqxacxg+uO4vTUuO0AmYsm34O3P4qPPMdeOfv8NQ3Ycsa+NjvZMGWIS6spK+UugG4CUgFrlBKlWutC6IamRg2TjS38cvnP+Avr21Fa5hqGcWPP34WS+dPMjs00ZuU8XCD02j1P/VNeP8p2F0L+Q9I2eYhLNyhD0Va62XAVt/ttCjFI4aZNzz1XHXPf7n/1a0kKEXhxVaqvnWxJPxYsjDfaPVPPw+O7IUHrzHG9w/xkYHDVbjdO41KqduANF+r3xvFmMQwcLy5ldJnP+DBtdsAmDtpHL90LGbhtFRzAxOnJm0mfO4po37P2nvh+Tthx/+MoZ0p8jsdSsJq6fta+WnAISDdd1uIU7K27iBX/Pa/PLh2G0kJiq8tncOTX71QEn6sS0yGy38CBY/AyPFGd0/ZRztW7BJDQliTs3ytewAFaACt9eNRjCtAJmfFj2MnW7nrPxt5+PUdAJw5eTx3OxZx1hRJ9nGnwQOrP2Mk/KQUuO5eWCynAQfTQCdnZXW5PhsYlKQv4oN7xyG+Wb6O7fXHSU5UfPWyOXzxo1kyozZepVvhC1XG6J63HoJ/FsKH68H+Q1mU3WRhJX2t9d3Bt5VS90UnHBFvWtra+d0LW/jDi1toa9fMO20cv71pCfNOkxWs4l7yKPj4741JXf8pNvr6D7wPN/5Z+vlNFO6QzT/h69bB6OLp9pVBiK48B47yzfJ1vL2rEaWg6GIr37r8DKlzP9zk3gYTzoDVn4XNz8OqpfDJv8OE082ObFgKt09/dvBtrfXWnu7ru38+xggfm9a6NMR+G2D1Hauyt2NJn37s0VrzyBs7+NnTGznR0saU1BR+tWwJ52dlmB2aMFPDVvj7zbB/A4xMhWUPQtZlZkcVt/rdp6+U+gUdrfuu+9Bar+xhnw1Aa+1SSlmVUjattbvL3Yq01kVKqRU97Bcx6tCxZr5T+TaujfsBuH7JFH708QWkjpKFTYa99Nnwhefh8SL44Gl4xGGc4D37FrMjG1Z6695xnuIxC4Aq33UPYAcCSd33LaAOINS3ABG7qrc18LXH3mJvYxPjU5L46ScW8rHFU8wOSwwlI8dBwcPg+oHRx/+vL0HjLrhkhdTtGSQ9Jv2+unB6YQEagm53/U6fC4FvBPYeun8KgUKAGTNmnGIYYrC0t2vue7mOX1dtoq1dY5th4d5Pns20NFnJSoSQkGCM57fMMEb3vPRzaNwJ1/7GGOsvoiqs8XJKqduUUjVKqXql1Bal1OYBPm+9v0vH1/LvRGtdprXO0VrnZGbKqkhD2YEjJ/nsA29y93Mf0Nauuf2SLMqLzpeEL/p2znK46RFIGmUM63y0AE4eMTuquBfuIGmH74TAKq316cCaXu7rBfxl9ixAfUwAgGUAACAASURBVJf9dRjdPvh+5oYZgxhi1m45yNX3vsIrmw+SPmYED96ayx1XzZOx9yJ8864xyjeMngB1a+CvH4PjDX0/TpyycP87G30/632zc7N7uW85vpE5vp8uAKWUf/ULV5f91WFHK4aE9nbNH17cwi33v8GBIyc5d3Y6//n6RXx07kSzQxOxaFoO3FYFlpmwx22synV4r9lRxa1wk/5yCEzSysLX3x5KULeNHfAGjcxZ49vvAbz+bp2+hmyKoeVwUwtFD9dy93MfAPC1pXN4dPl5TBovNe/FAKRb4fPPQeY8OLARHrgSDm0zO6q4FO44/cu01i8MQjzdyDj9oWPzh0coeqgWz8FjjE9J4p6bzubSedK6FxF0vAEevgH2vAXjJsOnn4CJ88yOKib1NE4/3JZ+tlLqeaXUfUqpJRGOTcSAp9/Zy8f/8Bqeg8eYd9o4nvzqhZLwReSNTofP/BtmXmjU5n/gKtizzuyo4kq4pZXv1lpfDpQCeREYvSNiRGtbOz9/ZiNfftTN8eY2rl8yhX9+6SPMzBhjdmgiXqWMh09Vwpwr4EQD/O3jsPdts6OKG+EO2RyvlFqOMWErF7gjqlGJIaHxRAu3PlhN2X89JCUofnjdmfymYAmjRkjtHBFlyaOMSVxzr4Ymry/xv2N2VHEh3O6dP2OMrb9ca71Ma/2PaAYlzLft4DE+8cfXeGXzQTLGjODR5efxuY/MRsmsSTFYkkaA469wxlVw4hD87WOyIEsEhL1y1mAtmiLMt3bLQaP//oDRf//Elz/CObPT+36gEJGWNAKW/dXX1XPIGMe/b73ZUcU0mUUjOnnkje185i9v0niiBfv8iVR+8QKmp8vsWmGipJFQ8BDMudzXx/8xOLDJ7KhiliR9ARgnbH/47/e485/raW3XFF1ixfnpHMaODHdxNSGiKGkkLHsIspbC8Xp46BPg3Wl2VDGpX0lfKSXLHcWh482tFD1Uy4Nrt5GcqLg7fxErr5pPYoL034shJDnFOLk7/Tw4vAseuh6OHjA7qpgT7uidG5RSq4EK3+3yqEYlBs3Boyf55Ko3WPP+fiyjk3l0+Xk4cqabHZYQoY0YDTeXw6SFUL/FmMjV1Nj340RAuC39Iq31MsBfbjktSvGIQbT14DFuvG8tb+/0Mi1tFP/44gXkzpITtmKIG2WBTz9ulG7Y9w48ehM0Hzc7qpgRdsE1pdRtQJqv4Jo3ijGJQeDecYgb71vL9vrjLJyayuNfuoCszLFmhyVEeMZONEo0jJsMO9bC48uhvc3sqGJC2EM2MVr3h4B0320Ro6o2fMjNq16n4Vgzl5yRyd8Lz2PiOCmYJmJM2kwj8aekwvtPwfP/z+yIYkJYQzN8BdfujnYwIvpW1+zkjn+8Q7uGgpzp/PQTC6T+vYhdE+cZJ3cfugFe/6NRnvm8282OakiTgmvDyJ9f8bCi0kj4X1s6h1/cuFASvoh9sy+Gj//euP7sHfD+0+bGM8RJwbVhQGvNr6s28dOnNwLw/WvP5Ft5Z0hJBRE/Ft8EH/0uoKHyC7C71uyIhiwpuBbn2ts1P3pyA/eu2UyCgl86FvP5C2ebHZYQkXfJClhyC7SegMdultW3ehDud/tVSMG1mNPa1s63K9/mwbXbGJGYwB9vySY/e5rZYQkRHUrBtb81avEf3Qflt0BLk9lRDTk9Jn2l1H1BN7cC5yil7vJfoh+aGIjm1na+8uhbPO7ezegRifzlc7lcueA0s8MSIrr8BdpSZxhdPE99E8JYHXA46W30TvDatc5oByIix0j4bp7f8CHjUpJ48NZzyJ4p8+nEMDFmAtz0CPzlCnj7UThtIZz/JbOjGjJ6bOlrrdcE3azXWm/1XwD56Byimlvb+dIjRsJPHZXMo7edJwlfDD+TF8H1fzSuP38n1L1objxDSLh9+iVdbkvLfwg62drGFx+uxbXRSPiP3HYuC6elmh2WEOY46xNw0bdBt0PF58C7w+yIhoRek75S6kZfobUCpdRzvsvzgIz1G2KaWtq4/aHaoMJp57JgqiR8McxdeqdRh7/JCxW3Qmuz2RGZrtekr7X+h6/kQpnW+grf5XLfmH0xRJxsbeP2h2t58YMDpI02unTOmiIJXwgSEuATTkidDrtrwPUDsyMyXbjdO2cH35DSykNHS5sxSuelDw6Q7lvL9swpsuyBEAGj08HxICQkG6UaNvzL7IhMFW7S79qdI2cGh4C2ds3/rX6bqg0dffjzJ0vCF6KbaTlw+U+M6//6CtTXmRuPicJN+h6l1LeVUkuUUt8BPNEMSvStvV3z3cff5d9v72HsyCT+9vlzJOEL0Ztzb4f5H4OTh6His9B60uyITBFu7Z3bgUbgdsDruy1MorXmJ09voLxmJynJCdz/2RwWT7eYHZYQQ5tSRmG2tNmw711Y82OzIzJF2CUWtdartNa3a61XRTMg0bdfV23igdeM0grOT+dwrjXD7JCEiA0pqXDjn0Elwv9+PyzH74dVhkEp9YvgEgxShsE8Zf+t43cvbCExQfG7m8/mkjMyzQ5JiNgyLQc+6qsZ+cQX4XiDufEMsqiUYVBK5WMsqWjTWpf2cr8Vve0XnT3u3sXPn3kfgF85FnPFWVJLR4hTcuG3YIsLdr4BT34dlv3N6P4ZBnrr3tFKqcuUUpcBs0NcQlJK2QC01i7A678d4n52jDLNIgwvbzrAisp3AKMe/vVnTzU5IiFiWGIS3FAGI8bBxn/DWw+bHdGg6S3pZwVdioFs3+Vy3+2eFNCxcLoHsA88zOHtnV1evvhwLa3tmqJLrFIPX4hISJsF1/zSuP6fYji0zcxoBk1vBddW+S8YI3bu9l3uoCOph2IBgjvJup1lVErZfN8ERB+2HTzGrQ9Uc7y5jRvOnkrxFfPMDkmI+LGoAM68HlqOwb+/NizKMIc7eidNKXWbUmqWbwWtrAE+b/oAHz8sHDhyks8+8Cb1x5q5+IxMSvIXkZAwPPodhRgUSsHVv4TRGbD1Zah90OyIoi7ccfqXY8zC/QGQClzWy929dCR1C1AfvDOcVr5SqlApVaOUqjlw4EA4IcadE81t3PbXarbXH2fRtFTuu8Umi5gLEQ1jM+Hqu43rz38PvDvNjSfKwl0j9waMk65TtNa/xFg+sSflgNV33Qq4fMfwzx6yKqXyfSN8rKFO9Gqty7TWOVrrnMzM4Tcksb1d863V63h7VyPT00fxl8/lMmZkbwOthBADctYNMO9aaD4CT8Z3N0+4TcciX7VNf/mFHmvvaK3dEBid4/XfBtb49ldqrSsxvg3INNIQ7n7+A/6zfh/jUpJ44HO5TBg70uyQhIhvSsE1v4ZRaVD3QlyP5gk36TcqpW4D0n2t/t5O5Ppb6i6tdVnQtuwQ98kK+lAQwOqandz3Uh2JCYr7bsnm9InjzA5JiOFh3CS4yjdt6Lk74eh+c+OJknCTfjFG6/4QkO5r9YsIe91Tz53/fBeAH3/8LC6cM8HkiIQYZhY64PQ8ONkIz/8/s6OJinCT/i8Ap6/2zp+jGdBwtfXgMW5/uJaWNs0XLpzNLefONDskIYYfpYyTukkp8E45eF4yO6KICzfpe4GtSqn7pPZO5B1uauELf63Ge7wF+/yJfPfq+WaHJMTwlT4bLv6Ocf2pb0FLk7nxRFh/Wvo5QClQ5ruICGhv13yrfB2eA8eYO2kc99x0NokyFl8Ic13wNZgwFxrq4LV7zI4monqrsjleKfUnpdRzwI1a663Bl0GMMa7ds2Yzro37SR2VTNlnsmVophBDQdIIuPbXxvVXfhVXK2311tJfBawGlgFKunQir2rDh9yzZjMJCu795NnMzBhjdkhCCL9ZF8Lim6HtpFGbJ070lvSV1voFrXWj1vpuBl56QQTZsv8o3yxfB8B3rpgndfGFGIryfgwjx8OWKthcZXY0EdFb0q/v47Y4RYebWih8qIajJ1u5ZtFkbr/E2veDhBCDb2wmXLLCuP7sSmhrMTeeCOgt6Rcopap9lxrf7c1KqS1Kqc2DFWC80Vrz7dVv4zlwjHmnjePu/EWoYbJ4gxAx6ZwiSM+C+s3wZuyvFttbaeV0rXWu75Ljuz1Ha3261nrOYAYZT+5/dSvPb/iQcSlJOD+dzegRcuJWiCEtaQRc8TPj+su/gGOx3ekhZRsHkXvHIX7xH2O5w7vzF8uJWyFixRlXgvVSaGqEF39mdjQDIkl/kBw61sxXH32L1nbN5z8ymysXyPq2QsQMpeDKu0AlQu0D8OEGsyM6ZZL0B0F7u+b/Kt5mt/cEi6dbuOMqWf1KiJgzcT7k3Aq6Hdb82OxoTpkk/UFQ9oqHF943JmD94eazGZEkb7sQMemSYkgeA5v+A9v/Z3Y0p0SyT5TVbGvg7uc+AODXyxYzLW20yREJIU7Z2IlwwVeM664fxORiK5L0o+hwUwtf//s62to1hRdbWTp/ktkhCSEG6vyvwOgJsPMN+OA/ZkfTb5L0o+j7T6xnt/cEC6em8u3L55odjhAiElLGd1ThXPMjaG8zN55+kqQfJU+8tZsn1u1hVHIi99y0RPrxhYgnObeCZSYceB/efszsaPpFMlEU7Gw4zveeWA/AD647E2vmWJMjEkJEVNJIuMy3stZLv4DWZnPj6QdJ+hHW2tbON8vXceRkK1ecNYmC3OlmhySEiIYF+ZA5Dxp3wrpHzI4mbJL0I+yPL9VRs/0Qk8aP5Bc3SF0dIeJWQoIxhBOMmvsx0tqXpB9B7+zycs8aoxbdrxxLSBszwuSIhBBRdeb1Mdfal6QfISdb2/i/1W/T1m4sbH7hnAlmhySEiLYYbO1L0o+Q37o2s3n/UawTxvCdK2R4phDDRoy19iXpR8BbOw7hfLmOBAV3OxaTkpxodkhCiMESY619SfoD1NTSxrcr3qZdw/KLrGTPTDM7JCHEYDvzesicb7T2311tdjS9kqQ/QL+u2kTdgWNkZY7hm3lnmB2OEMIMCQlw4TeM66/+FtrbzY2nF5L0B6B2ewOrXvGQoOBXy5ZIt44Qw9mCGyF1urGs4gfPmB1NjyTpn6KTrW2sqHwHraHokiyWTLeYHZIQwkyJyUYxNoDXfjtkK3BK0j9F971UR92BY1gzx/D1pbJksBACsH0aRqXDrmrYvtbsaEKSpH8Ktuw/wh9frAPg559YKN06QgjDiDFwTqFx/dXfmBtLDyTp91N7u+a7j6+nua2dm3Knc541w+yQhBBDyTmFkDQKtlTBvvVmR9NNVJK+UipfKWVXSq3oYX+h71ISjeePpvKanby5rYEJY0ey8qr5ZocjhBhqxmSA7TPG9bW/MzeWECKe9JVSNgCttQvw+m8H7bcDLq11GWD13Y4J+w838fNnNgJGyeTU0ckmRySEGJLO/xKoBFj/DzjyodnRdBKNln4B4PVd9wBdk7o1aJvHdzsm/OipDRxpauXSuZlcu2iy2eEIIYaqtFkw92pob4HaB8yOppNoJH0L0BB0u1Ont9a6zNfKB7ABNVGIIeJe3nSAp9/Zy6jkRH5y/QIpmSyE6N25txs/q++H1pPmxhLEtBO5vm6fKq21O8S+QqVUjVKq5sCBAyZE19nJ1jZ++O/3APiGfQ7T0kabHJEQYsibdSFMWgDH9sN7/zQ7moBoJH0vkO67bgHqe7ifXWtdGmqH79tAjtY6JzMzMwoh9s/9r25l60Gj1MKtH5ltdjhCiFigFJxbZFx//b4hM1krGkm/nI5+eivgAlBKBaasKqUK/Ql/qJ/I3dt4gt+t2QLAjz62QBY4F0KEb6HDmKy1dx3sfNPsaIAoJH1/d40vmXuDum/WBG0vUUrVKaUORfr5I+2nT2/kREsbVy88TRZGEUL0T/IoyLnVuP7GfebG4pMUjYMGnagN3pbt++kCYqL+8GtbDgZO3t55zZlmhyOEiEW5txmVNzc+aQzfHDfJ1HCkr6IHLW3t/MB38vYrl53OVMsokyMSQsSk8VNg7lXQ3gpvPWR2NJL0e/LXtdvYsv8oszJGc9tFcvJWCDEA/i4e919Nr7UvST+EQ8eauXfNZgC+f92ZjEySgmpCiAGwXgaWmeDdAXUvmBqKJH2gtLQUh8NBXl4eSikuvPw66v7+E05v8XDp3ImmxORyuSgqKgp7+6kKPlZaWhpFRUU4HA4cDgder7fH7Xl5ebjdnadYlJaW4nK58Hg8lJaGHI3bo+LiYhwOR4+P62l/JN8LIaImIQGyP2tcN3uGrtZ6SF+ys7P1YJoxc7bOWvm0nnXHU3rDnsZBfe6eVFRURO24wce22WyB61VVVbqwsLDH7RUVFXrFihWdjme32wPXS0pKdG1tbVhx1NbWBo61YsUKXVdX1+f+Q4cO6cLCQm2xWAL3czqd2m63a7vdrm02my4pKQnr+YUYFIf3af2jdK1/mKZ1456oPx1Qo0Pk1KiM3ollh443k96uuSl3OvMnjzc7HADKy8vJz8+P+HGdTidVVVUh9+Xk5FBS0r0Iqn97fn4+xcXFgft4PB4slo7Vw/Lz8ykpKcHpdPYZh8vlIi8vD4Dc3FxcLheFhYV97nc6ndTUdFTxKCwsDDyuqKgoKu+Z6F1LSwu7du2iqanJ7FCGpqufgJbjULcdUrx93z8MKSkpTJs2jeTk8ApAStIP8r+6ek60tDF6RCLfuvwMPB4PJSUlNDQ08OUvf5nHHnsMp9NJaWkpVquV/Pz8QHdHqCTj8XhwOBzk5OTgcrmoqqrCajXmrfkfZ7VacTqduFwunE4n6enpFBUVYbFYAgnT36Xjj8XpdOLxeFi1ahUWi4XS0lJsNht2u71TPP5j+I9ps9m6veaelJWVBRJtT9ttNhtutxubzUZlZWWnrhar1dopIfemvr6e9HRjErfFYqG6urpf+7tyuVxYLJbAey0Gz65duxg3bhyzZs2S+lShNB2GhjpISIZJ84xZuwOgtaa+vp5du3Yxe3Z4A04k6fu0t2t++vQGAL54SRYTx6Vw9ACsXr2arVu30tDQ0O0xpaWl5OXlUVhYiMPhCNmy9Cf1yspKnE4nJSUllJaWUlBQQH5+PpWVlZSVlVFbW8vKlSsDidnj8QBQUlKC2+0OfAD44ygoKGD16tUUFhZSVVXFihUrusWTnp7e6ZjBPB5Pt6To8XgoKiqioaEBq9XaqRUfantBQQHl5eXYbLZADKF4PB5cLlfIfcEt+kgpKSnp8RuMiK6mpiZJ+L0ZOQ4SR0BbM5w8AikD601QSpGRkUF/apRJ0vd5/K3dvLfnMIkJitsu6kiGy5Ytw2KxhEz6dXV1NDQ0UFtbi9frxev1duriAAIt1Pz8/EDirq6uDiQ7m80W6CYpLi7G6/WG1SWSn59PXl4edrs9kLy7xnPnnXfys5/9LHDM4CTvP0nbNdZQz93Tdn8Xz8qVK0O2qv2v3Wq19prcMzIyAu+v1+slIyOjX/uDhXpdYnBJwu+FUjA6A47sheMNA076xiH7937L6B2gqaWNXz3/AQCWUcmMGtExRLNrEgcjuQJkZ2dTUFAQ6BsPdV9/snK5XIHE6O+X9m/Pzc3FarVSUVFBcXFxWEnfH1twt0rXeJYsWdLjMa1Wa8gPsv7yf2g5HI5u+/zH93g8lJWVhbwA2O32QMu8qqoKu71zOaa+9gdzuVz96sYS8amsrIysrKzAN+muo7xcLhderxe3201WVhZFRUWBb9fQ0TUbjsrKyn7F5j2pcL+7EZq8xoStQSYtfeDh17ezt7GJ+ZPHs2lk6LfE30cd/Mfh70ZxOp1YLBYqKiq6Pc7/x+PxeFizZg0AK1aswOFwcNdddwWSfWlpaSCxhTqBmpeX1y1xFxQUUFxcHOhW6RpPbm5uj8e0WCwRaRX7h3L29kHVV0vfZrNRXl6Ow+HAarUGknZeXh5VVVU97vf/LoqKiiguLsZqteLxeMjKyhrw6xKxLScnB7vdHuhyLS4upqysjMLCwsDfvcViwWazYbPZKCoq6vRt1Wq1smrVqrCeq7fGU6hv/5YJk/Dsrse2UMMJL4wZ5JpeoYb0DKVLtIdsHj7RrJf86Dk9s/gpvWbjvogeu66uLjDscSgKHmIZaYcOHdL5+flRO74YmjZs2GB2CFprY8iw0+kM3PYPM/bvC5afn99teHFdXV3YQ46Dn6eriooKfejQoW7b69bX6Apnqdb7PwjrOfoS6n2nhyGbw75758+vbOXQ8RZyZqaZNhHLLEVFRf3+ahquu+66i5UrV0bl2EL0pbq6mpycHMBobXs8nkCr398925fi4uLAdf/EQ3+X5EBZ5y+h6pU3oOUYtAzu8NZh3b1Tf/Qkf37F6KopvmpexE9A+UfuDFX5+fkUFRVht9tDno84Vf6vz9K3PrzNuuPpqBx32y+u6fM+brcbj8dDQ0MDXq+336O5rFZr4H+iuLg4MGgi+INgQBISaTh8wrh+ogGSp0TmuGEY1kn/jy/Vcay5jUvnZnLO7IxOQy5LSkpCDmn0D7scLC6Xi4qKim4fHj1t7y//44uKigLXlVIh34u0tDSWLVsW6MP0zxMAQu7zeDxUVlb2OJQzlOLiYjweD7m5uSEf598fPHTUf8IteJSRf4hpTyOPRPzyNzpOdXJeZWVlp8d6PJ5OJ337EjxEuba2ttPExeBzW+mZvhLLxxtg3OQBj9kP17BN+ru9J3jof9sB+PYVc/mv74TqUGO32wOjVYL/GIO3D1RlZWWniVjWHt6L4G8uHo+HpUuXUltb222fy+Vi+fLlgWP4J3D1xV/Lxz/iqOtcAv8/UkVFReAkrtvtDsxNKCsrw+VykZ6eHjix63A4wn5+EVnhtMijoaamZkD/G8GjeIDA6DqbzRbWhL/ggQuVlZU9f5NOSOoYs9981BjDPwiGbZ/+Pa5NNLe1c93iKZw1JdXscMJSXl4eleM6nc5+t4qsVit2uz3kOQG73d6ptRVuSztUuYVgVVVVgf3Z2dmBYbD+5/L31Qb/c4aahCbiV/As+nC43W7cbjfl5eVUVlbicDior68PbPd4PKxYsQKXy4XL5Qp7lnlYlIJRvvWkTgzeIoLDsqW/7eAx/uHeTWKC4lt5Z4S8j9frZfny5Xi9XhwOR7daMF3LG8RrOYbe+IdUdn29Lpcr0PUSyXIMGRkZnYaZ1tXVUVhYSHl5OVlZWeTn5wdaeG63G4fDEXj9YniwWq199t9nZ2cHrttstk4ndoP/loO396eLMhwej8dowIxKg6MfGkM3U6eBin47fFi29H//4hba2jU3nD2V2RPGAB2J2eFwUFxcHBh3X1VV1a2lWlFRwcqVK3E6ndhstkD5g6qqqh5b4/6kHlyEzF+Owd+CLSsr63Zsv5KSEnJycgJj8P385RigY+JS13h6Oqb/dYc6dxH8XoTD/5js7Gy8Xm+PpR96m6TVl66TtDIyMgJdQv5/UP9Xc5vNRm1tLVVVVd2+rovhbdmyZREdteZvqISSn58fstHhdruND5jkUZCUArrNKMswCIZdS397/TH++ZbRyv/KZacHtofqxy4tLaW+vr5b0iguLu5UMiHeyjH01KffldvtJjc3t9Nj/C3s4BZTpMoxBE/KslqtWK1WysvLA10+eXl5nU4cWywWHA5Hv08mi/hmsVgCkxMj8S2wv12j3Rpao9KMsgwnDkFK9Luah13S/4OvlZ+fPY2ZGWN6vJ+/kuaKFSu6tQr8Cc7fFZOdnU16enqvv/zeyjHk5+d3K8fgP3Y4i4SEKsfQNZ7gYwaPPjrVcgxer5fy8vLAiVw//wzH4JPOweUYeiu8ZrfbKS8vD7ToQ712/wehw+EIdHO53e7AY3Jzczs9d21tbchqoWJ4i9QgiFPR7RyTP+k3NUJ7GyREd6W+YZX0dzYc53FfX/5XLj291/va7XaWL18esn+wa8kEf59+vJdjCK62CfT4bWDVqlUsXbq024fgQMsx+M+vWCyWwMQv/2vNy8sLlLv2n8vwD9mUuvpiSEsaCcmjjTr7Jw93nNyNEmXM1h26cnJydKTOmN/xj3f4e/VObrBN5dfLlkTkmOHwjygYquPF/Uk1GvwnxIficFgReRs3bmT+/PlmhxF7ju6Hw7uN7p30/o82C/W+K6VqtdY5Xe87bE7k7mw4TmXtLhIUfPWyOWaHM6RIOQYhTDbKd26h6XDUK28Om6T/x5e20NquuX5Jx4idwRIL5Rj83SeRJOUYhAhT4ggYMRbQRt9+FA2LpL/HeyLQyg8esSM6dB0KGgkWi2VQS1YI4defevrZ2dmBocn+80bhDKAInkwYkfr7gdZ+96TvjzUShkXSv//VrbS0aa5eOBlr5lizwxFCRFlwPf3CwkIsFktgPkjXevolJSW4XK5AjZySkpKwGivB5RUiUn/fP1yz6bAxiieIxWKJ2HyTuE/63uPNPPbmDgBuv0QW1xBiOHC5XJ1m3nq93sB8kbKysm5DNv3zY06Vx+MZeFJOHAHJYwBtjOLpwj8UeqDieshm45NP4vn53VQcOsCR8RlM/0gxTL3O7LCEGB5+GKWJRj/su8+7uro6MIAgnHr6drudkpIS3G53p27OsrIycnJyqKmp6XW4MRiTNv2j4PwlUTweT5+P62RUqlFj/0Rjt6Gb/sqyAx2CHLct/cYnn2Tv977PqEMHSABSD9ez93vfp/HJJ80OTQgRZf5iaf5CaeEMSXY6nSxfvjxw2z9B01/Ar7dyIV3r7/trYIW7YEtAiu8D52Qj6PZuuyOxrnVUWvpKqXzAC9i01qX93R8J+3/zW3RT5xVpdFMT+3/zW1Kvk9a+EFEXRos8Gk61nr6/cqy/T7+6ujpwDP8IvNraWhoaGnpc7wH6X3/fL/ChcngftDfDmJewzj2rU1dUb3V+whXxpK+UsgForV1KKatSyqa1doe7P1Ja9+7t13YhRHwYSD39lStXMnv2bEpKSsjNzQ3UyeltYZ+u+lt/3y/QDXRkn1GWYXQ6WGae0uvoTTS6dwowWvEAHqDru9/X/ohoycgMuT1p8uRoPJ0QZXqWCQAABeBJREFUYgjobz19l8tFcXFxoCZU8DDjFStW4Ha7cblcuN3ubgnfX3M/+BKR+vv+UTwnGiEKFROi0b1jAYLf8Yx+7h8wrTV/X3wtjpceIqWtJbBdpaQw8ZvfiPTTCSGGiP7W07fb7d2KBgafePUn+lDfHLrW4o9Y/f2kFEgcCW0nO62oFajBP0BD8kSuUqpQKVWjlKo5cOBAvx+/o+E4z0xazIPnf9Jo2StF0pQpTP7Jj6U/X4hhLtL19Purz355pYxRPCho7TgvGajBP0DRaOl7Af+rsgD1/dyP1roMKAOj4Fp/A5iZMYa1d1yG58C5zJkqdV+EEB0iXU+/v8JK3GMmwdjTAmWWI7nsZzSSfjngr+xmBVwASimL1trb0/5IGz0iiQVTY2PtWyHE4DKznn5YEjun5kiu8xzx7h3/SByllB3wBo3MWdPHfiFEHBjq5drjTX/f76iM0/d1z3Tdlt3bfiFE7EtJSaG+vp6MjAyUUmaHE/e01tTX15OSkhL2Y+K6DIMQYnBNmzaNXbt2cSoDMMSpSUlJYdq0aWHfX5K+ECJikpOTmT17ttlhiF4MySGbQgghokOSvhBCDCOS9IUQYhhRQ314lVLqALD9FB8+ATgYwXBigbzm4UFec/wb6OudqbXuVoRsyCf9gVBK1Witc/q+Z/yQ1zw8yGuOf9F6vdK9I4QQw4gkfSGEGEbiPekPx5m/8pqHB3nN8S8qrzeu+/SFEEJ0Fu8tfRFHlFL5Sim7UqrXFSr62i/EUORfSraHfWH97YcjbpJ+X29KJN+0oSKM11zou5QMdmyRFry2MuDt6R/EV701dzBji6Ywfsc2330GvrrGENGP/+XCUPtjke/vdlUP+8L62w9XXCT9vt6USL9pQ0EYr9kOuHwVTa2+27FsUNZWHkrC/Lst0lpXYvyOh8PftQ3w+PZ74uE1Q+D19rSwb0T/9uMi6TNEFmMfZH29JmvQNo/vdizrc21lpZTN988TL3r9Hfta93UAWuvSOFmbIpz/Vf83V2ucvOa+RHRd8XhJ+qYvxm6CXl+T1rosaN0CG1AzWIGZqI/FR2NOX3+3uUCGr4snXrot+/q7dmO08OvouWUsehEvSV/0wPf1tyoOWkS9rq0ch638cNUHrUYXN/36PVFKWTD+FpzAKqVUrH+DDUef64r3R7wk/QEvxh6Dwn1Ndq116eCEFFXldHRRdVp72b8t6IRmXPRv0/fvuA6jCwTfz3g4gd3Xay4E7vL9TTuAuP2gC/rbDvm3f6riJen3lRAi+qYNEX29ZpRShf6EH+sncsNYe7nSd0IzHSNZxIO+fseuLvurBzW66Ojz79rPf7J38EKLHl9jJafLt7WorCseN5OzfMO3PBgnd8p822r9a/OG2h/renvNvj+QCox+z3TAMUy7P2JamH/XDb798fCNLpzXvMK3Pz1e/pcHU9wkfSGEEH2Ll+4dIYQQYZCkL4QQw4gkfSGEGEYk6Qvh46tTVOcb+lmolHJ22W9XSll8k6HqlFLO4HHiSimrUqqih2Nb4mQYqYhxkvSF6FCDUa+o0jcqxOsv6uUfMqi19g+ZcwNOrbV/nDy+68tDHVhr7SX2S2GIOCBJX4gOdqA26HZwSYDCvoa8+lr9vSV293CYNSuGNkn6QnTIxVejyNeyt/omfAFkhXmMQBlrpdSK4BLAvm8CeRGMV4h+k6QvRAcbHWWo7VrrfiVoX1L3AvjWMHD7vh0Ef2DEW1E4EWOSzA5AiKEgqM++sq/79vD4/C6PtQL+k7d1QdulMqQwlSR9IQw5DKwmU9e+/GqMxT7cSilPqAcIYQZJ+mLY852ALSbM4l2+1rsNKPA9toCOVZxsSimr1rrU16fv786RukdiSJDaO0KEwVexdEDFvXwfELZT7UISIhLkRK4Q4VkdgeGWkvCF6STpCxEG3+Qqb6i67uHwtfKlb1+YTrp3hBBiGJGWvhBCDCOS9IUQYhiRpC+EEMOIJH0hhBhGJOkLIcQwIklfCCGGkf8PmnS+2Io+uU4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Pr(E|H) Chance of a positive test (E) given there is a liar (H). This is the chance of a true positive.\n",
    "PrEH = 0.85\n",
    "# Chance of a positive test (E) given that the person wasn't a lying. This is the probability of a  false positive\n",
    "PrEnH = 0.25\n",
    "\n",
    "# posterior of true positives / positive predictive value\n",
    "vPrHE = vfunc(PrEH, vPrH, PrEnH)\n",
    "# posterior of false positives\n",
    "vPrnHE = 1 - vPrHE \n",
    "# posterior of true negatives / false predictive value\n",
    "vPrnHnE = vfunc( (1-PrEnH), (1-vPrH), (1-PrEH) )\n",
    "# posterior of false negatives\n",
    "vPrHnE = 1 - vPrnHnE \n",
    "\n",
    "# Plot the two posteriors\n",
    "plt.clf()\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(vPrH, vPrHE, '-', lw=2, color='tab:blue')\n",
    "ax.plot(vPrH, vPrnHnE, '-', lw=2, color='tab:orange')\n",
    "plt.xlabel(r'$\\displaystyle P(\\textrm{Lie})$')       \n",
    "plt.ylabel(r\"Predictive value\") \n",
    "plt.legend([r\"$\\displaystyle P(\\textrm{Lie}\\mid +)$\", r\"$\\displaystyle P(\\textrm{No-lie}\\mid -)$\"],loc='lower right')\n",
    "\n",
    "i = 50\n",
    "\n",
    "ax.plot(vPrH[i], vPrHE[i], 'o', lw=2, color='tab:red')\n",
    "ax.plot(vPrH[i], vPrnHnE[i], 'o', lw=2, color='tab:red')\n",
    "\n",
    "ax.annotate('True positives (PPV)={:.3f} \\n False positives (FDR)={:.3f}'.format(vPrHE[i],vPrnHE[i]),\n",
    "            xy=(vPrH[i], vPrHE[i]),\n",
    "            xycoords='data',\n",
    "            xytext=(vPrH[i] + 0.07, vPrHE[i] + 0.17),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "ax.annotate('True negatives (NPV)={:.3f} \\n False negatives (FOR)={:.3f}'.format(vPrnHnE[i], vPrHnE[i]),\n",
    "            xy=(vPrH[i], vPrnHnE[i]),\n",
    "            xycoords='data',\n",
    "            xytext=(vPrH[i] + 0.07, vPrnHnE[i] - 0.07),\n",
    "            textcoords='axes fraction',\n",
    "            horizontalalignment='left',\n",
    "            verticalalignment='top')\n",
    "\n",
    "vPrH[i]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decision tree of events\n",
    "\n",
    "To sum up all the concepts we present the decision tree of events of the example of iBorderCtrl to relate the corrected test outcomes and the posterior probabilities.\n",
    "\n",
    "![Decision tree of events](pics/decision-tree.png)\n",
    "\n",
    "\n",
    "## Conclusions\n",
    "\n",
    "- Massive screening of entire populations to look for rare events through a flawed test is problematic \n",
    "- Performance analysis should be corrected with the prior probability of the event of interest\n",
    "- None of the proposed tests will be able to detect illegalised migrants or terrorists with the so-called lie detectors, even if they really existed \n",
    "- If we assume that lie detectors work, terrorists would likely adopt countermeasures, which would decrease the true discovery rate \n",
    "- There is no trade-off between human rights, freedom and privacy and security in the massive screening of a population with flawed tests. No security improvement can be achieved by these means. By forcing migrants and asylum seekers to pass through a test, intended to discover criminals in the public imaginary, these collectives are indeed being criminalised. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[0] National Research Council. 2003. The Polygraph and Lie Detection. Washington, DC: The National Academies Press. https://doi.org/10.17226/10420.\n",
    "\n",
    "[1] BetterExplained, [An Intuitive (and Short) Explanation of Bayes’ Theorem](https://betterexplained.com/articles/an-intuitive-and-short-explanation-of-bayes-theorem/)\n",
    "\n",
    "[2] Countbayesie.com, [Han Solo and Bayesian Priors](https://www.countbayesie.com/blog/2015/2/18/hans-solo-and-bayesian-priors)\n",
    "\n",
    "[3] Wikipedia, [Base rate fallacy](https://en.wikipedia.org/wiki/Bayesian_inference)\n",
    "\n",
    "[4] Wikipedia, [Bayesian inference](https://en.wikipedia.org/wiki/Base_rate_fallacy)\n",
    "\n",
    "[4] Vera Wilde, [Refugee screening: a brief introduction (and a request for equipment)](http://www.scq.ubc.ca/refugee-screening-a-brief-introduction-and-a-request-for-equipment/)\n",
    "\n",
    "[5] OrShea, James, Keeley Crockett, Wasiq Khan, Philippos Kindynis, Athos Antoniades, and Georgios Boultadakis. 2018. ‘Intelligent Deception Detection through Machine Based Interviewing’. In 2018 International Joint Conference on Neural Networks (IJCNN), 1–8. Rio de Janeiro, Brazil: IEEE. https://doi.org/10.1109/IJCNN.2018.8489392.\n",
    "\n",
    "[6] Moliné, A., E. Dominguez, E. Salazar‐López, G. Gálvez‐García, J. Fernández‐Gómez, J. De la Fuente, O. Iborra, F. J. Tornay, and E. Gómez Milán. 2018. ‘The Mental Nose and the Pinocchio Effect: Thermography, Planning, Anxiety, and Lies’. Journal of Investigative Psychology and Offender Profiling 15 (2): 234–48. https://doi.org/10.1002/jip.1505.\n",
    "\n",
    "[7] ‘Científicos de La UGR Investigan El Denominado “Efecto Pinocho” | Canal UGR’. n.d. Accessed 3 June 2019. https://perma.cc/799W-364A.\n",
    "\n",
    "[8] Simcox, Robin. n.d. ‘The Asylum–Terror Nexus: How Europe Should Respond’, no. 3314: 12. https://www.heritage.org/sites/default/files/2018-06/BG3314.pdf\n",
    "\n",
    "[9] https://ec.europa.eu/eurostat/en/web/products-datasets/-/MIGR_ASYAPPCTZA\n",
    "\n",
    "[10] Positive and negative predictive values https://en.wikipedia.org/wiki/Positive_and_negative_predictive_values\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}