Chapter 16 - Inference about Independence.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 16. Inference about Independence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter addresses two questions:\n",
    "\n",
    "1. How do we test if two random variables are independent?\n",
    "2. How do we estimate the strength of dependence between two random variables?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall we write $Y \\text{ ⫫ } Z$ to mean that $Y$ and $Z$ are independent.\n",
    "\n",
    "When $Y$ and $Z$ are not independent, we say they are **dependent** or **associated** or **related**.\n",
    "\n",
    "Note that dependence does not mean causation:\n",
    "- Smoking is related to heart disease, and quitting smoking will reduce the chance of heart disease.\n",
    "- Owning a TV is related to lower starvation, but giving a starving person a TV does not make them not hungry."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16.1 Two Binary Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that both $Y$ and $Z$ are binary.  Consider a data set $(Y_1, Z_1), \\dots, (Y_n, Z_n)$.  Represent the data as a two-by-two table:\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cc|c} \n",
    "      & Y = 0  & Y = 1 & \\\\\n",
    "\\hline\n",
    "Z = 0 & X_{00} & X_{01} & X_{0\\text{·}}\\\\\n",
    "Z = 1 & X_{10} & X_{11} & X_{1\\text{·}}\\\\\n",
    " \\hline\n",
    "      & X_{\\text{·}0} & X_{\\text{·}1} & n = X_{\\text{··}}\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "where $X_{ij}$ represents the number of observations where $(Z_k, Y_k) = (i, j)$.  The dotted subscripts denote sums, e.g. $X_{i\\text{·}} = \\sum_j X_{ij}$.  Denote the corresponding probabilities by:\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cc|c} \n",
    "      & Y = 0  & Y = 1 & \\\\\n",
    "\\hline\n",
    "Z = 0 & p_{00} & p_{01} & p_{0\\text{·}}\\\\\n",
    "Z = 1 & p_{10} & p_{11} & p_{1\\text{·}}\\\\\n",
    " \\hline\n",
    "      & p_{\\text{·}0} & p_{\\text{·}1} & 1\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "where $p_{ij} = \\mathbb{P}(Z = i, Y = j)$.  Let $X = (X_{00}, X_{01}, X_{10}, X_{11})$ denote the vector of counts.  Then $X \\sim \\text{Multinomial}(n, p)$ where $p = (p_{00}, p_{01}, p_{10}, p_{11})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **odds ratio** is defined to be\n",
    "\n",
    "$$ \\psi = \\frac{p_{00} p_{11}}{p_{01} p_{10}}$$\n",
    "\n",
    "The **log odds ratio** is defined to be\n",
    "\n",
    "$$ \\gamma = \\log \\psi$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.2**.  The following statements are equivalent:\n",
    "\n",
    "1. $Y \\text{ ⫫ } Z$\n",
    "2. $\\psi = 1$\n",
    "3. $\\gamma = 0$\n",
    "4. For $i, j \\in \\{ 0, 1 \\}$, $p_{ij} = p_{i\\text{·}} p_{\\text{·}j}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider testing\n",
    "\n",
    "$$\n",
    "H_0: Y \\text{ ⫫ } Z\n",
    "\\quad \\text{versus} \\quad\n",
    "H_1: \\text{not} (Y \\text{ ⫫ } Z)\n",
    "$$\n",
    "\n",
    "First consider the likelihood ratio test.  Under $H_1$, $X \\sim \\text{Multinomial}(n, p)$ and the MLE is $\\hat{p} = X / n$.  Under $H_0$, again $X \\sim \\text{Multinomial}(n, p)$ but $p$ is subjected to the constraint $p_{ij} = p_{i.} p_{.j}$.  This leads to the following test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.3 (Likelihood Ratio Test for Independence in a 2-by-2 table)**. \n",
    "Let\n",
    "\n",
    "$$ T = 2 \\sum_{i=0}^1 \\sum_{j=0}^1 X_{ij} \\log \\left( \\frac{X_{ij} X_{\\text{··}}}{X_{i\\text{·}} X_{\\text{·}j}} \\right)$$\n",
    "\n",
    "Under $H_0$, $T \\leadsto \\chi_1^2$.  Thus, an approximate level $\\alpha$ test is obtained by rejecting $H_0$ when $T > \\chi_{1, \\alpha}^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.4 (Pearson's $\\chi^2$ test for Independence in a 2-by-2 table)**. Let\n",
    "\n",
    "$$ U = \\sum_{i=0}^1 \\sum_{j=0}^1 \\frac{(X_{ij} - E_{ij})^2}{E_{ij}} $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ E_{ij} = \\frac{X_{i\\text{·}} X_{\\text{·}j}}{n}$$\n",
    "\n",
    "Under $H_0$, $U \\leadsto \\chi_1^2$.  Thus, an approximate level $\\alpha$ test is obtained by rejecting $H_0$ when $U > \\chi_{1, \\alpha}^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the intuition for the Pearson test: Under $H_0$, $p_{ij} = p_{i\\text{·}} p_{\\text{·}j}$, so the MLE of $p_{ij}$ is $\\hat{p}_{ij} = \\hat{p}_{i\\text{·}} \\hat{p}_{\\text{·}j} = \\frac{X_{i\\text{·}}}{n} \\frac{X_{\\text{·}j}}{n}$.  Thus, the expected number of observations in the $(i, j)$ cell is $E_{ij} = n \\hat{p}_{ij} = \\frac{X_{i\\text{·}} X_{\\text{·}j}}{n}$.  The statistic $U$ compares the observed and expected counts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.6**. The MLE's of $\\psi$ and $\\gamma$ are\n",
    "\n",
    "$$\n",
    "\\hat{\\psi} = \\frac{X_{00} X_{11}}{X_{01} X_{10}}\n",
    ", \\quad\n",
    "\\hat{\\gamma} = \\log \\hat{\\psi}\n",
    "$$\n",
    "\n",
    "The asymptotic standard errors (computed from the delta method) are\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{\\text{se}}(\\hat{\\psi}) &= \\sqrt{\\frac{1}{X_{00}} + \\frac{1}{X_{01}} + \\frac{1}{X_{10}} + \\frac{1}{X_{11}}}\\\\\n",
    "\\hat{\\text{se}}(\\hat{\\gamma}) &= \\hat{\\psi} \\hat{\\text{se}}(\\hat{\\gamma})\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yet another test of independence is the Wald test for $\\gamma = 0$ given by $W = (\\hat{\\gamma} - 0) / \\hat{\\text{se}}(\\hat{\\gamma})$. \n",
    "\n",
    "A $1 - \\alpha$ confidence interval for $\\gamma$ is $\\hat{\\gamma} \\pm z_{\\alpha/2} \\hat{\\text{se}}(\\hat{\\gamma})$.\n",
    "\n",
    "A $1 - \\alpha$ confidence interval for $\\psi$ can be obtained in two ways.  First, we could use $\\hat{\\psi} \\pm z_{\\alpha/2} \\hat{\\text{se}}(\\hat{\\psi})$.  Second, since $\\psi = e^{\\gamma}$ we could use  $\\exp \\{\\hat{\\gamma} \\pm z_{\\alpha/2} \\hat{\\text{se}}(\\hat{\\gamma})\\}$.  This second method is usually more accurate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16.2 Interpreting the Odds Ratios"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose event $A$ has probability $\\mathbb{P}(A)$.  The odds of $A$ are defined as\n",
    "\n",
    "$$\\text{odds}(A) = \\frac{\\mathbb{P}(A)}{1 - \\mathbb{P}(A)}$$\n",
    "\n",
    "It follows that\n",
    "\n",
    "$$\\mathbb{P}(A) = \\frac{\\text{odds}(A)}{1 + \\text{odds}(A)}$$\n",
    "\n",
    "Let $E$ be the event that someone is exposed to something (smoking, radiation, etc) and let $D$ be the event that they get a disease.  The odds of getting the disease given exposure are:\n",
    "\n",
    "$$\\text{odds}(D | E) = \\frac{\\mathbb{P}(D | E)}{1 - \\mathbb{P}(D | E)}$$\n",
    "\n",
    "and the odds of getting the disease given non-exposure are:\n",
    "\n",
    "$$\\text{odds}(D | E^c) = \\frac{\\mathbb{P}(D | E^c)}{1 - \\mathbb{P}(D | E^c)}$$\n",
    "\n",
    "The **odds ratio** is defined to be\n",
    "\n",
    "$$\\psi = \\frac{\\text{odds}(D | E)}{\\text{odds}(D | E^c)}$$\n",
    "\n",
    "If $\\psi = 1$ then the disease probability is the same for exposed and unexposed; this implies these events are independent.  Recall that the log-odds ratio is defined as $\\gamma = \\log \\psi$.  Independence corresponds to $\\gamma = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider this table of probabilities:\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cc|c} \n",
    "      & D^c    & D      & \\\\\n",
    "\\hline\n",
    "E^c   & p_{00} & p_{01} & p_{0\\text{·}}\\\\\n",
    "E     & p_{10} & p_{11} & p_{1\\text{·}}\\\\\n",
    " \\hline\n",
    "      & p_{\\text{·}0} & p_{\\text{·}1} & 1\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Denote the data by\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cc|c} \n",
    "      & D^c    & D      & \\\\\n",
    "\\hline\n",
    "E^c   & X_{00} & X_{01} & X_{0\\text{·}}\\\\\n",
    "E     & X_{10} & X_{11} & X_{1\\text{·}}\\\\\n",
    " \\hline\n",
    "      & X_{\\text{·}0} & X_{\\text{·}1} & X_{\\text{··}}\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Now\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(D | E) = \\frac{p_{11}}{p_{10} + p_{11}}\n",
    "\\quad \\text{and} \\quad\n",
    "\\mathbb{P}(D | E^c) = \\frac{p_{01}}{p_{00} + p_{01}}\n",
    "$$\n",
    "\n",
    "and so\n",
    "\n",
    "$$\n",
    "\\text{odds}(D | E) = \\frac{p_{11}}{p_{10}}\n",
    "\\quad \\text{and} \\quad\n",
    "\\text{odds}(D | E^c) = \\frac{p_{01}}{p_{00}}\n",
    "$$\n",
    "\n",
    "and therefore\n",
    "\n",
    "$$ \\psi = \\frac{p_{11}p_{00}}{p_{01}p_{10}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate the parameters, we have to consider how the data were collected.  There are three methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Multinomial Sampling**.  We draw a sample from the population and, for each sample, record their exposure and disease status.  In this case, $X = (X_{00}, X_{01}, X_{10}, X_{11}) \\sim \\text{Multinomial}(n, p)$.  We then estimates the probabilities in the table by $\\hat{p}_{ij$ = X_{ij} / n$ and\n",
    "\n",
    "$$ \\hat{\\psi} = \\frac{\\hat{p}_{11} \\hat{p}_{00}}{\\hat{p}_{01} \\hat{p}_{10}} = \\frac{X_{11} X_{00}}{X_{01} X_{10}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Prospective Sampling (Cohort Sampling)**.  We get some exposed and unexposed people and count the number with disease within each group.  Thus,\n",
    "\n",
    "$$\n",
    "X_{01} \\sim \\text{Binomial}(X_{0\\text{·}}, \\mathbb{P}(D | E^c))\n",
    "\\quad \\text{and} \\quad\n",
    "X_{11} \\sim \\text{Binomial}(X_{1\\text{·}}, \\mathbb{P}(D | E))\n",
    "$$\n",
    "\n",
    "In this case we should write small letters $x_{0\\text{·}},  x_{1\\text{·}}$ instead of capital letters $ X_{0\\text{·}},  X_{1\\text{·}}$ since they are fixed and not random, but we'll keep using capital letters for notational simplicity.\n",
    "\n",
    "We can estimate $\\mathbb{P}(D | E))$ and $\\mathbb{P}(D | E^c)$ but we cannot estimate all probabilities in the table.  Still, we can estimate $\\psi$.  Now:\n",
    "\n",
    "$$\\hat{\\mathbb{P}}(D | E) = \\frac{X_{11}}{X_{1\\text{·}}}\n",
    "\\quad \\text{and} \\quad\n",
    "\\hat{\\mathbb{P}}(D | E^c) = \\frac{X_{01}}{X_{0\\text{·}}}\n",
    "$$\n",
    "\n",
    "Thus,\n",
    "\n",
    "$$ \\hat{\\psi} = \\frac{X_{11} X_{00}}{X_{01} X_{10}}$$\n",
    "\n",
    "as before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Case-Control (Retrospective Sampling)**.  Here we get some diseased and non-diseased people and we observe how many are exposed.  This is much more efficient if the disease is rare.  Hence,\n",
    "\n",
    "$$\n",
    "X_{10} \\sim \\text{Binomial}(X_{\\text{·}0}, \\mathbb{P}(E | D^c))\n",
    "\\quad \\text{and} \\quad\n",
    "X_{11} \\sim \\text{Binomial}(X_{\\text{·}1}, \\mathbb{P}(E | D))\n",
    "$$\n",
    "\n",
    "From this data we can estimate $\\mathbb{P}(E | D)$ and $\\mathbb{P}(E | D^c)$.  Surprisingly, we can still estimate $\\psi$.  To understand why, note that\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(E | D) = \\frac{p_{11}}{p_{01} + p_{11}},\n",
    "\\quad 1 - \\mathbb{P}(E | D) = \\frac{p_{01}}{p_{01} + p_{11}},\n",
    "\\quad \\text{odds}(E | D) = \\frac{p_{11}}{p_{01}}\n",
    "$$\n",
    "\n",
    "By a similar argument,\n",
    "\n",
    "$$\\text{odds}(E | D^c) = \\frac{p_{10}}{p_{00}}$$\n",
    "\n",
    "Hence,\n",
    "\n",
    "$$\\frac{\\text{odds}(E | D)}{\\text{odds}(E | D^c)} = \\frac{p_{11} p_{00}}{p_{01} p_{10}} = \\psi$$\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$\\hat{\\psi} = \\frac{X_{11} X_{00}}{X_{01} X_{10}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In all three methods, the estimate of $\\psi$ turns out to be the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is tempting to try to estimate $\\mathbb{P}(D | E) - \\mathbb{P}(D | E^c)$.  In a case-control design, this quantity is not estimable.  To see this, we apply Bayes' theorem to get\n",
    "\n",
    "$$\\mathbb{P}(D | E) - \\mathbb{P}(D | E^c) = \\frac{\\mathbb{P}(E | D) \\mathbb{P}(D))}{\\mathbb{P}(E)} - \\frac{\\mathbb{P}(E^c | D) \\mathbb{P}(D)}{\\mathbb{P}(E^c)}$$\n",
    "\n",
    "Because of the way we obtained the data, $\\mathbb{P}(D)$ is not estimable from the data.\n",
    "\n",
    "However, we can estimate $\\xi = \\mathbb{P}(D | E) / \\mathbb{P}(D | E^c)$, which is called the **relative risk**, under the **rare disease assumption**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.9**.  Let $\\xi = \\mathbb{P}(D | E) / \\mathbb{P}(D | E^c)$.  Then\n",
    "\n",
    "$$ \\frac{\\psi}{\\xi} \\rightarrow 1$$\n",
    "\n",
    "as $\\mathbb{P}(D) \\rightarrow 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, under the rare disease assumption, the relative risk is approximately the same as the odds ratio, which we can estimate.\n",
    "\n",
    "In a randomized experiment, we can interpret a strong association, that is $\\psi \\neq 1$, as a causal relationship.  In an observational (non-randomized) study, the association can be due to other unobserved **confounding** variables.  We'll discuss causation in more detail later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16.3 Two Discrete Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose that $Y \\in \\{ 1, \\dots, I \\}$ and $Z \\in \\{ 1, \\dots, J \\}$ are two discrete variables.  The data can be represented by an $I \\times J$ table of contents:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{array}{c|cccccc|c} \n",
    "       & Y = 1  & Y = 2  & \\cdots & Y = j & \\cdots & Y = J   & \\\\\n",
    "\\hline\n",
    "Z = 1 & X_{11}  & X_{12} & \\cdots & X_{1j} & \\cdots & X_{1J} & X_{1\\text{·}}\\\\\n",
    "Z = 2 & X_{21}  & X_{22} & \\cdots & X_{2j} & \\cdots & X_{2J} & X_{2\\text{·}}\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n",
    "Z = i & X_{i1}  & X_{i2} & \\cdots & X_{ij} & \\cdots & X_{iJ} & X_{i\\text{·}}\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n",
    "Z = I & X_{I1}  & X_{I2} & \\cdots & X_{Ij} & \\cdots & X_{IJ} & X_{I\\text{·}}\\\\\n",
    " \\hline\n",
    "      & X_{\\text{·}1} & X_{\\text{·}1} & \\cdots & X_{\\text{·}j} & \\cdots & X_{\\text{·}J} & n\n",
    "\\end{array}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider testing\n",
    "\n",
    "$$\n",
    "H_0: Y \\text{ ⫫ } Z\n",
    "\\quad \\text{versus} \\quad\n",
    "H_1: \\text{not } H_0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.10**.  Let\n",
    "\n",
    "$$ T = 2 \\sum_{i=1}^I \\sum_{j=1}^J X_{ij} \\log \\left( \\frac{X_{ij} X_{\\text{··}}}{X_{i\\text{·}} X_{\\text{·}j}} \\right) $$\n",
    "\n",
    "The limiting distribution of $T$ under the null hypothesis of independence is $\\chi^2_\\nu$ where $\\nu = (I - 1)(J - 1)$.\n",
    "\n",
    "Pearson's $\\chi^2$ test statistic is\n",
    "\n",
    "$$ U = \\sum_{i=1}^I \\sum_{j=1}^J \\frac{(n_{ij} - E_{ij})^2}{E_{ij}}$$\n",
    "\n",
    "Asymptotically, under $H_0$, $U$ has a $\\chi^2_\\nu$ distribution where $\\nu = (I - 1)(J - 1)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a variety of ways to quantify the strength of dependence between two discrete variables $Y$ and $Z$.  Most of them are not very intuitive.  The one we shall use is not standard but is more interpretable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define\n",
    "\n",
    "$$\\delta(Y, Z) = \\max_{A, B} \\Big|\\; \\mathbb{P}_{Y, Z}(Y \\in A, Z \\in B) - \\mathbb{P}_Y(Y \\in A)\\mathbb{P}_Z(Z \\in B) \\;\\Big|$$\n",
    "\n",
    "where the maximum is over all pairs of events $A$ and $B$.\n",
    "\n",
    "*Note: the original version on the book contains a typo; it is originally stated as*\n",
    "\n",
    "$$\\delta(Y, Z) = \\max_{A, B} \\Big|\\; \\mathbb{P}_{Y, Z}(Y \\in A, Z \\in B) - \\mathbb{P}_Y(Y \\in A) - \\mathbb{P}_Z(Z \\in B) \\;\\Big|$$\n",
    "\n",
    "*For more context, this is a specific case of the \"total variation distance\" $\\delta(P_1, P_2)$ between the probability metrics $P_1(Y, Z) = \\mathbb{P}_{Y, Z}(Y \\in A, Z \\in B)$ and $P_2(Y, Z) = \\mathbb{P}_Y(Y \\in A)\\mathbb{P}_Z(Z \\in B)$.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.12**.  Properties of $\\delta(Y, Z)$:\n",
    "\n",
    "1. $0 \\leq \\delta(Y, Z) \\leq 1$\n",
    "2. $\\delta(Y, Z) = 0$ if and only if $Y \\text{ ⫫ } Z$\n",
    "3. The following identity holds:\n",
    "\n",
    "$$ \\delta(X, Y) = \\frac{1}{2} \\sum_{i=1}^I \\sum_{j=1}^J \\Big|\\; p_{ij} - p_{i\\text{·}} p_{\\text{·}j} \\;\\Big|$$\n",
    "\n",
    "4.  The MLE of $\\delta$ is\n",
    "\n",
    "$$ \\hat{\\delta}(X, Y) = \\frac{1}{2} \\sum_{i=1}^I \\sum_{j=1}^J \\Big|\\; \\hat{p}_{ij} - \\hat{p}_{i\\text{·}} \\hat{p}_{\\text{·}j} \\;\\Big|$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\hat{p}_{ij} = \\frac{X_{ij}}{n},\n",
    "\\quad \\hat{p}_{i\\text{·}} = \\frac{X_{i\\text{·}}}{n},\n",
    "\\quad \\hat{p}_{\\text{·}j} = \\frac{X_{\\text{·}j}}{n}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpretation of $\\delta$ is this: if one person makes probability statements assuming independence and another person makes probability statements without assuming independence, their probability statements may differ by as much as $\\delta$.  Here's a suggested scale for interpreting $\\delta$:\n",
    "\n",
    "| range        | interpretation             |\n",
    "|--------------|----------------------------|\n",
    "| 0 to 0.01    | negligible association     |\n",
    "| 0.01 to 0.05 | non-negligible association |\n",
    "| 0.05 to 0.1  | substantial association    |\n",
    "| over 0.1     | very strong association    |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A confidence interval for $\\delta$ can be obtained by bootstrapping.  The steps are:\n",
    "\n",
    "1. Draw $X^* \\sim \\text{Multinomial}(n, \\hat{p})$;\n",
    "2. Compute $\\hat{p}_ij, \\hat{p}_{i\\text{·}}, \\hat{p}_{\\text{·}j}$;\n",
    "3. Compute $\\delta^*$;\n",
    "4. Repeat.\n",
    "\n",
    "Now we use any of the methods we learned earlier for constructing bootstrap confidence intervals.  However, we should not use a Wald interval in this case.  The reason is that if $Y$ and $Z$ are independent then $\\delta = 0$ and we are on the boundary of the parameter space.  In this case, the Wald method is not valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16.4 Two Continuous Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose that $Y$ and $Z$ are both continuous.  If we assume that the joint distribution of $Y$ and $Z$ is bivariate Normal, then we measure the dependence between $Y$ and $Z$ by means of the correlation coefficient $\\rho$.  Tests, estimates, and confidence intervals for $\\rho$ in the Normal case are given in the previous chapter.  If we do not assume normality, then we need a nonparametric method for asserting dependence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the correlation is\n",
    "\n",
    "$$\\rho = \\frac{\\mathbb{E}((X_1 - \\mu_1)(X_2 - \\mu_2))}{\\sigma_1 \\sigma_2} $$\n",
    "\n",
    "A nonparametric estimator of $\\rho$ is the plug-in estimator which is:\n",
    "\n",
    "$$\\hat{\\rho} = \\frac{\\sum_{i=1}^n (X_{1i} - \\overline{X}_1)(X_{2i} - \\overline{X}_2)}{\\sqrt{\\sum_{i=1}^n (X_{1i} - \\overline{X}_1)^2 \\sum_{i=1}^n (X_{2i} - \\overline{X}_2)^2}} $$\n",
    "\n",
    "which is just the sample correlation.  A confidence interval can be constructed using the bootstrap.  A test for $\\rho = 0$ can be based on the Wald test using the bootstrap to estimate the standard error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plug-in approach is useful for large samples.  For small samples, we measure the correlation using the **Spearman rank correlation coefficient** $\\hat{\\rho}_S$.  We simply replace the data by their ranks, ranking each variable separately, then compute the correlation coefficient of the ranks.\n",
    "\n",
    "To test the null hypothesis that $\\rho_S = 0$, we need the distribution of $\\hat{\\rho}_S$ under the null hypothesis.  This can be obtained by simulation.  We fix the ranks of the first variable as $1, 2, \\dots, n$.  The ranks of the second variable are chosen at random from the set of $n!$ possible orderings, then we compute the correlation.  This is repeated many times, and the resulting distribution $\\mathbb{P}_0$ is the null distribution of $\\hat{\\rho}_S$. The p-value for the test is $\\mathbb{P}_0(|R| > |\\hat{\\rho}_S|)$ where $R$ is drawn from the null distribution $\\mathbb{P}_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16.5 One Continuous Variable and One Discrete"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $Y \\in \\{ 1, \\dots, I \\}$ is discrete and $Z$ is continuous.  Let $F_i(z) = \\mathbb{P}(Z \\leq z | Y = i)$ denote the CDF of $Z$ conditional on $Y = i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.15**.  When $Y \\in \\{ 1, 2, \\dots, I \\}$ is discrete and $Z$ is continuous, then $Y \\text{ ⫫ } Z$ if and only if $F_1 = F_2 = \\dots = F_I$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It follows that to test for independence, we need to test\n",
    "\n",
    "$$ \n",
    "H_0: F_1 = \\dots = F_I\n",
    "\\quad \\text{versus} \\quad\n",
    "H_1: \\text{not } H_0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity, we consider the case where $I = 2$.  To test the null hypothesis that $F_1 = F_2$ we will use the **two sample Kolmogorov-Smirnov test**.  \n",
    "\n",
    "Let $n_k$ denote the number of observations for which $Y_i = k$.  Let\n",
    "\n",
    "$$\\hat{F}_k(z) = \\frac{1}{n_k} \\sum_{i=1}^n I(Z_i \\leq z) I(Y_i = k)$$\n",
    "\n",
    "denote the empirical distribution function of $Z$ given $Y = k$.  Define the test statistic\n",
    "\n",
    "$$ D = \\sup_x | \\hat{F}_1(x) - \\hat{F}_2(x) |$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 16.16**.  Let\n",
    "\n",
    "$$ H(t) = 1 - 2 \\sum_{j=1}^\\infty (-1)^{j-1} e^{-2j^2t^2} $$\n",
    "\n",
    "Under the null hypothesis that $F_1 = F_2$,\n",
    "\n",
    "$$ \\lim_{n \\rightarrow \\infty} \\mathbb{P} \\left( \\sqrt{\\frac{n_1 n_2}{n_1 + n_2}} D \\leq t \\right) = H(t) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It follows from the theorem than an approximate level $\\alpha$ test is obtained by rejecting $H_0$ when\n",
    "\n",
    "$$ \\sqrt{\\frac{n_1 n_2}{n_1 + n_2}} D > H^{-1}(1 - \\alpha) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16.7 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.1**.  Prove Theorem 16.2.\n",
    "\n",
    "The following statements are equivalent:\n",
    "\n",
    "1. $Y \\text{ ⫫ } Z$\n",
    "2. $\\psi = 1$\n",
    "3. $\\gamma = 0$\n",
    "4. For $i, j \\in \\{ 0, 1 \\}$, $p_{ij} = p_{i\\text{·}} p_{\\text{·}j}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "$Y$ and $Z$ are independent if and only if $\\mathbb{P}(Y = j | Z = i) = \\mathbb{P}(Y = j) \\mathbb{P}(Z = i)$ for all $i$ and $j$, so (1) and (4) are equivalent.\n",
    "\n",
    "$\\gamma = \\log \\psi$, so (2) and (3) are equivalent.\n",
    "\n",
    "Finally,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&\\psi = 1\\\\\n",
    "&\\Leftrightarrow \\frac{\\text{odds}(Y | Z = 1)}{\\text{odds}(Y | Z = 0)} = 1 \\\\\n",
    "&\\Leftrightarrow \\text{odds}(Y | Z = 1) = \\text{odds}(Y | Z = 0) \\\\\n",
    "&\\Leftrightarrow \\mathbb{P}(Y | Z = 1) = \\mathbb{P}(Y | Z = 0) \\\\\n",
    "&\\Leftrightarrow Y \\text{ ⫫ } Z\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "so (1) and (2) are equivalent.\n",
    "\n",
    "Therefore, all statements are equivalent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.2**.  Prove Theorem 16.3.\n",
    "\n",
    "Let\n",
    "\n",
    "$$ T = 2 \\sum_{i=0}^1 \\sum_{j=0}^1 X_{ij} \\log \\left( \\frac{X_{ij} X_{\\text{··}}}{X_{i\\text{·}} X_{\\text{·}j}} \\right)$$\n",
    "\n",
    "Under $H_0$, $T \\leadsto \\chi_1^2$.  Thus, an approximate level $\\alpha$ test is obtained by rejecting $H_0$ when $T > \\chi_{1, \\alpha}^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**. Let $E_{ij} = \\frac{X_{i\\text{·}} X_{\\text{·}j}}{X_\\text{··}}$.  We can rewrite the test statistic as\n",
    "\n",
    "$$ T = 2 \\sum_{i, j} X_{ij} \\log \\frac{X_{ij}}{E_{ij}} $$\n",
    "\n",
    "The interpretation of $E_{ij}$ is that it is the expected count in cell $(i, j)$ under $H_0$:\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{E}_{H_0}(Y = i, Z = j) &= n \\mathbb{P}_{H_0}(Y = i, Z = j) \\\\\n",
    "&= n \\mathbb{P}(Y = i) \\mathbb{P}(Z = j) \\\\\n",
    "&= n \\hat{p}_{i\\text{·}} \\hat{p}_{\\text{·}j} \\\\\n",
    "&= X_{\\text{··}} \\frac{X_{i\\text{·}}}{X_{\\text{··}}} \\frac{X_{\\text{·}j}}{X_{\\text{··}}} \\\\\n",
    "&= E_{ij}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Now, the result will follow by applying the log-likelihood test between $H_0$ and $H_1$.  The test statistic is:\n",
    "\n",
    "$$ 2 \\log \\frac{\\mathcal{L}(\\tilde{p} | X)}{\\mathcal{L}(\\hat{p} | X)} \n",
    "= 2 \\log \\frac{\\prod_{i, j} \\tilde{p}_{ij}^{X_{ij}}}{\\prod_{i, j} \\hat{p}_{ij}^{X_{ij}}}$$\n",
    "\n",
    "where $\\tilde{p}$ is the MLE under $H_1$ and $\\hat{p}$ is the MLE under $H_0$.  But we have:\n",
    "\n",
    "$$\n",
    "\\tilde{p}_{ij} = \\frac{X_{ij}}{n}\n",
    "\\quad \\text{and} \\quad\n",
    "\\hat{p}_{ij} = \\frac{E_{ij}}{n}\n",
    "$$\n",
    "\n",
    "Substituting the MLEs on the log-likelihood ratio, we get\n",
    "\n",
    "$$ 2 \\log \\frac{\\mathcal{L}(\\tilde{p} | X)}{\\mathcal{L}(\\hat{p} | X)} \n",
    "= 2 \\log \\prod_{i, j} \\left( \\frac{X_{ij}}{E_{ij}} \\right)^{X_{ij}}\n",
    "= 2 \\sum_{i, j} X_{ij} \\log \\frac{X_{ij}}{E_{ij}}\n",
    "$$\n",
    "\n",
    "which is the desired result.\n",
    "\n",
    "*Reference: \"G-test.\" Wikipedia: The Free Encyclopedia. Wikimedia Foundation, Inc. 22 July 2004. Web. 17 Mar. 2020, en.wikipedia.org/w/index.php?title=G-test&oldid=914538756*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.3**.  Prove Theorem 16.9.\n",
    "\n",
    "Let $\\xi = \\mathbb{P}(D | E) / \\mathbb{P}(D | E^c)$.  Then\n",
    "\n",
    "$$ \\frac{\\psi}{\\xi} \\rightarrow 1$$\n",
    "\n",
    "as $\\mathbb{P}(D) \\rightarrow 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  We have:\n",
    "\n",
    "$$ \\mathbb{P}(D | E) = \\frac{\\mathbb{P}(D \\land E)}{\\mathbb{P}(E)} = \\frac{p_{11}}{p_{01} + p_{11}}\n",
    "\\quad\n",
    "\\mathbb{P}(D | E^c) = \\frac{\\mathbb{P}(D \\land E^c)}{\\mathbb{P}(E^c)} = \\frac{p_{10}}{p_{00} + p_{10}}\n",
    "$$\n",
    "\n",
    "and so\n",
    "\n",
    "$$\n",
    "\\xi = \\frac{\\mathbb{P}(D | E)}{\\mathbb{P}(D | E^c)}\n",
    "= \\frac{p_{11} (p_{00} + p_{10})}{p_{10} (p_{01} + p_{11})}\n",
    "$$\n",
    "\n",
    "Since\n",
    "\n",
    "$$ \\psi = \\frac{p_{11}p_{00}}{p_{01}p_{10}}$$\n",
    "\n",
    "we have:\n",
    "\n",
    "$$ \\frac{\\psi}{\\xi} = \\frac{p_{11}p_{00}}{p_{01}p_{10}} \\frac{p_{10} (p_{01} + p_{11})}{p_{11} (p_{00} + p_{10})} \n",
    "= \\frac{p_{00}}{p_{01}} \\frac{p_{\\text{·}1}}{p_{\\text{·}0}}\n",
    "$$ \n",
    "\n",
    "As $\\mathbb{P}(D) = p_{10} + p_{11} \\rightarrow 0$ and the probabilites are non-negative, $p_{10} \\rightarrow 0$, $p_{11} \\rightarrow 0$, so\n",
    "\n",
    "$$ \\frac{\\psi}{\\xi} \\rightarrow  \\frac{p_{00}}{p_{01}} \\frac{p_{01}}{p_{00}} = 1 $$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.4**.  Prove equation (16.14).\n",
    "\n",
    "$$ \\delta(X, Y) = \\max_{A, B} \\Big|\\; \\mathbb{P}_{X, Y}(X \\in A, Y \\in B) - \\mathbb{P}_X(X \\in A) \\mathbb{P}_Y(Y \\in B) \\;\\Big| $$\n",
    "\n",
    "is equivalent to\n",
    "\n",
    "$$ \\delta(X, Y) = \\frac{1}{2} \\sum_{i=1}^I \\sum_{j=1}^J \\Big|\\; p_{ij} - p_{i\\text{·}} p_{\\text{·}j} \\;\\Big|$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "As noted when the definition was introduced, this definition of $\\delta(X, Y)$ is a particular case of the total variation distance between metrics $P_1(X, Y) = \\mathbb{P}_{X, Y}(X \\in A, Y \\in B)$ and $P_2(X, Y) = \\mathbb{P}_X(X \\in A) \\mathbb{P}_Y(Y \\in B)$, where\n",
    "\n",
    "$$ \\delta(P, Q) = \\sup_{A \\in \\mathcal{F}} |P(A) - Q(A)| $$\n",
    "\n",
    "and $P, Q$ are two probability measures on a sigma-algebra $\\mathcal{F}$ of subsets from the sample space $\\Omega$.\n",
    "\n",
    "(Note that $P_1(X = i, Y = j) = p_{ij}$ and $P_2(X = i, Y = j) =  p_{i\\text{·}} p_{\\text{·}j}$).\n",
    "\n",
    "It is a known property of the total variation distance in countable spaces that\n",
    "\n",
    "$$ \\delta(P, Q) = \\max_{A \\in \\mathcal{F}} |P(A) - Q(A)| = \\frac{1}{2} \\Vert P - Q \\Vert_1 = \\frac{1}{2} \\sum_{\\omega \\in \\Omega} \\vert P(\\omega) - Q(\\omega) \\vert$$\n",
    "\n",
    "Let's demonstrate it, following the steps in the reference below.\n",
    "\n",
    "Let $\\mathcal{F}^+ \\subset \\mathcal{F}$ be the set of states $\\omega$ such that $P(\\omega) \\geq Q(\\omega)$, and let $\\mathcal{F}^- = \\mathcal{F} - S^+$ be the set of states $\\omega$ such that $P(\\omega) < Q(\\omega)$.  We have:\n",
    "\n",
    "$$ \\max_{A \\in \\mathcal{F}} P(A) - Q(A) = P(\\mathcal{F}^+) - Q(\\mathcal{F}^+)$$\n",
    "\n",
    "and\n",
    "\n",
    "$$ \\max_{A \\in \\mathcal{F}} Q(A) - P(A) = Q(\\mathcal{F}^-) - P(\\mathcal{F}^-)$$\n",
    "\n",
    "But since $P(\\Omega) = Q(\\Omega) = 1$, we have\n",
    "\n",
    "$$ P(\\mathcal{F}^+) + P(\\mathcal{F}^-) = Q(\\mathcal{F}^+) + Q(\\mathcal{F}^-) = 1 $$\n",
    "\n",
    "which implies that\n",
    "\n",
    "$$ P(\\mathcal{F}^+) - Q(\\mathcal{F}^+) = Q(\\mathcal{F}^-) - P(\\mathcal{F}^-)$$\n",
    "\n",
    "hence\n",
    "\n",
    "$$ \\max_{A \\in \\mathcal{F}} P(A) - Q(A) = \\vert P(\\mathcal{F}^+) - Q(\\mathcal{F}^+) \\vert = \\vert Q(\\mathcal{F}^-) - P(\\mathcal{F}^-) \\vert $$\n",
    "\n",
    "Finally, since\n",
    "\n",
    "$$ \n",
    "\\vert P(\\mathcal{F}^+) - Q(\\mathcal{F}^+) \\vert + \\vert Q(\\mathcal{F}^-) - P(\\mathcal{F}^-) \\vert = \\sum_{\\omega \\in \\Omega} | P(\\omega) - Q(\\omega) | = \\Vert P - Q \\Vert_1\n",
    "$$\n",
    "\n",
    "we have that\n",
    "\n",
    "$$ \\max_{A \\in \\mathcal{F}} |P(A) - Q(A)| = \\frac{1}{2} \\Vert P - Q \\Vert_1 $$\n",
    "\n",
    "which completes the proof.\n",
    "\n",
    "*Reference: Upfal, E., and M. Mitzenmacher. \"Probability and computing.\" (2005).  Lemma 11.1, pg. 272-273.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.5**.  The New York Times (January 8, 2003, page A12) reported the following data on death sentencing and rance, from a study in Maryland:\n",
    "\n",
    "|  | Death Sentence | No Death Sentence |\n",
    "|--|----------------|-------------------|\n",
    "| Black Victim | 14 | 641 |\n",
    "| White Victim | 62 | 594 |\n",
    "\n",
    "Analyse the data using the tools from this Chapter.  Interpret the results.  Explain why, based only on this information, you can't make causal conclusions.  (The authors of the study did use much more information in their full report)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "The statistic for the log-likelihood test for independence is is:\n",
    "\n",
    "$$ T = 2 \\sum_{i, j} X_{ij} \\log \\frac{X_{ij} X_{\\text{··}}}{X_{i\\text{·}} X_{\\text{·}j}} $$\n",
    "\n",
    "The statistic for Pearson's $\\chi^2$ test is:\n",
    "\n",
    "$$ U = \\sum_{i, j} \\frac{(X_{ij} - E_{ij})^2}{E_{ij}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "X_ij = np.array([[14, 641], [62, 594]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_idot = X_ij.sum(axis = 1).reshape(2, 1)\n",
    "X_dotj = X_ij.sum(axis = 0).reshape(1, 2)\n",
    "n = X_ij.sum()\n",
    "E_ij = X_idot @ X_dotj / n\n",
    "\n",
    "# Log-likelihood test statistic\n",
    "T = 2 * (X_ij * np.log(X_ij / E_ij)).sum()\n",
    "\n",
    "# Pearson's \\chi^2 test statistic\n",
    "U = ((X_ij - E_ij)**2 / E_ij).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T test statistic: \t\t34.534\n",
      "Pearson chi^2 statistic: \t32.104\n",
      "p-value log likelihood:\t\t 0.000\n",
      "p-value Pearson: \t\t 0.000\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import chi2\n",
    "\n",
    "p_value_log_likelihood = 1 - chi2.cdf(T, 1)\n",
    "p_value_pearson = 1 - chi2.cdf(U, 1)\n",
    "\n",
    "print('T test statistic: \\t\\t%.3f'% T)\n",
    "print('Pearson chi^2 statistic: \\t%.3f'% U)\n",
    "\n",
    "print('p-value log likelihood:\\t\\t %.3f' % p_value_log_likelihood)\n",
    "print('p-value Pearson: \\t\\t %.3f' % p_value_pearson)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both tests indicate that these variables are correlated.  Note that this is not sufficient information, in itself, to imply causation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta(X, Y): 0.037\n"
     ]
    }
   ],
   "source": [
    "from itertools import product\n",
    "\n",
    "p = X_ij / n\n",
    "p_idot = X_idot / n\n",
    "p_dotj = X_dotj / n\n",
    "\n",
    "delta = sum([abs(p[i, j] - p_idot[i, :] * p_dotj[:, j]) for (i, j) in product(range(2), range(2))]) / 2\n",
    "\n",
    "print('delta(X, Y): %.3f' % delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\delta(X, Y)$ is between 0.01 and 0.05, suggesting a non-negligible association."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.6**.  Analyse the data on the variables Age and Financial Status from:\n",
    "\n",
    "http://lib.stat.cmu.edu/DASL/Datafiles/montanadat.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv('data/montana.csv')\n",
    "\n",
    "# Select wanted columns and remove missing data\n",
    "data = data[['AGE', 'FIN']].replace('*', np.nan).dropna().astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[21, 16, 34],\n",
       "       [17, 23, 26],\n",
       "       [22, 37, 11]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Count occurrences in each cell\n",
    "X_ij = np.zeros((3, 3)).astype(int)\n",
    "\n",
    "for index, row in data.iterrows():\n",
    "     X_ij[row['AGE'] - 1, row['FIN'] - 1] += 1\n",
    "           \n",
    "X_ij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_idot = X_ij.sum(axis = 1).reshape(3, 1)\n",
    "X_dotj = X_ij.sum(axis = 0).reshape(1, 3)\n",
    "n = X_ij.sum()\n",
    "E_ij = X_idot @ X_dotj / n\n",
    "\n",
    "# Log-likelihood test statistic\n",
    "T = 2 * (X_ij * np.log(X_ij / E_ij)).sum()\n",
    "\n",
    "# Pearson's \\chi^2 test statistic\n",
    "U = ((X_ij - E_ij)**2 / E_ij).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T test statistic: \t\t22.064\n",
      "Pearson chi^2 statistic: \t20.679\n",
      "p-value log likelihood:\t\t 0.000\n",
      "p-value Pearson: \t\t 0.000\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import chi2\n",
    "\n",
    "p_value_log_likelihood = 1 - chi2.cdf(T, 4)\n",
    "p_value_pearson = 1 - chi2.cdf(U, 4)\n",
    "\n",
    "print('T test statistic: \\t\\t%.3f'% T)\n",
    "print('Pearson chi^2 statistic: \\t%.3f'% U)\n",
    "\n",
    "print('p-value log likelihood:\\t\\t %.3f' % p_value_log_likelihood)\n",
    "print('p-value Pearson: \\t\\t %.3f' % p_value_pearson)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both tests indicate correlation between the two sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta(X, Y): 0.128\n"
     ]
    }
   ],
   "source": [
    "from itertools import product\n",
    "\n",
    "p = X_ij / n\n",
    "p_idot = X_idot / n\n",
    "p_dotj = X_dotj / n\n",
    "\n",
    "delta = sum([abs(p[i, j] - p_idot[i, :] * p_dotj[:, j]) for (i, j) in product(range(3), range(3))]) / 2\n",
    "\n",
    "print('delta(X, Y): %.3f' % delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\delta(X, Y) > 0.1$ indicates a strong correlation between these variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.7**.  Estimate the correlation between temperature and latitude using the data from\n",
    "\n",
    "http://lib.stat.cmu.edu/DASL/Datafiles/USTemperatures.html\n",
    "\n",
    "Use the correlation coefficient and Spearman rank correlation.  Provide estimates, tests, and confidence intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv('data/USTemperatures.txt', sep='\\t')\n",
    "filtered_data = data[['lat', 'JanTF']].dropna()\n",
    "X, Y = filtered_data['lat'], filtered_data['JanTF']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import rankdata\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "def get_pearson_correlation(xx, yy):\n",
    "    mu_1 = xx.mean()\n",
    "    mu_2 = yy.mean()\n",
    "    sigma2_1 = ((xx - mu_1)**2).mean()\n",
    "    sigma2_2 = ((yy - mu_2)**2).mean()\n",
    "    \n",
    "    return ((xx - mu_1) * (yy - mu_2)).mean() / math.sqrt(sigma2_1 * sigma2_2)\n",
    "\n",
    "def get_spearman_correlation(xx, yy):\n",
    "    return get_pearson_correlation(rankdata(xx), rankdata(yy))\n",
    "\n",
    "def bootstrap_correlation(xx, yy, corr_fun, B=10000, alpha=0.05):\n",
    "    n = len(xx)\n",
    "    assert len(yy) == n, 'Sequences must have same length'\n",
    "    \n",
    "    t_boot = np.empty(B)\n",
    "    for i in tqdm_notebook(range(B)):\n",
    "        indexes = np.random.randint(0, n, size=n)\n",
    "        xx_selected, yy_selected = xx.iloc[indexes], yy.iloc[indexes]\n",
    "        t_boot[i] = corr_fun(xx_selected, yy_selected)\n",
    "        \n",
    "    confidence_interval = (np.quantile(t_boot, alpha / 2), np.quantile(t_boot, 1 - alpha / 2))\n",
    "    return confidence_interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4786b969ae4d45e79f4a89ddb4d447d4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=10000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "096eeca2baa84523a41e3885079212ae",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=10000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "pearson_rho = get_pearson_correlation(X, Y)\n",
    "pearson_rho_confidence = bootstrap_correlation(X, Y, corr_fun = get_pearson_correlation, B = 10000, alpha=0.05)\n",
    "\n",
    "spearman_rho = get_spearman_correlation(X, Y)\n",
    "spearman_rho_confidence = bootstrap_correlation(X, Y, corr_fun = get_spearman_correlation, B = 10000, alpha=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pearson confidence: \t -0.857\n",
      "95% confidence interval: -0.964, -0.684\n",
      "\n",
      "Spearman confidence: \t -0.834\n",
      "95% confidence interval: -0.953, -0.649\n"
     ]
    }
   ],
   "source": [
    "print('Pearson confidence: \\t %.3f' % pearson_rho)\n",
    "print('95%% confidence interval: %.3f, %.3f' % pearson_rho_confidence)\n",
    "print()\n",
    "print('Spearman confidence: \\t %.3f' % spearman_rho)\n",
    "print('95%% confidence interval: %.3f, %.3f' % spearman_rho_confidence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 16.7.8**.  Test whether calcium intake and drop in blood pressure are associated.  Use the data in\n",
    "\n",
    "http://lib.stat.cmu.edu/DASL/Datafiles/Calcium.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "We intend to use Theorem 16.16 here:\n",
    "\n",
    "Let\n",
    "\n",
    "$$ H(t) = 1 - 2 \\sum_{j=1}^\\infty (-1)^{j-1} e^{-2j^2t^2} $$\n",
    "\n",
    "Under the null hypothesis that $F_1 = F_2$,\n",
    "\n",
    "$$ \\lim_{n \\rightarrow \\infty} \\mathbb{P} \\left( \\sqrt{\\frac{n_1 n_2}{n_1 + n_2}} D \\leq t \\right) = H(t) $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ D = \\sup_x | \\hat{F}_1(x) - \\hat{F}_2(x) |$$\n",
    "\n",
    "It follows from the theorem than an approximate level $\\alpha$ test is obtained by rejecting $H_0$ when\n",
    "\n",
    "$$ \\sqrt{\\frac{n_1 n_2}{n_1 + n_2}} D > H^{-1}(1 - \\alpha) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = pd.read_csv('data/calcium.txt', sep='\\t')\n",
    "X, Y = data['Treatment'], data['Decrease']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAHSCAYAAADmLK3fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAbVklEQVR4nO3df5Bd5X3f8c+X1WKJ2EwCEsZBJtJQeTBYIPCi4qn5YVNjgRtk1fZUjsfBrjOUjmndzDQpcSY2k0ymKXGmqWsSDaSM4xla7AZjNDGEmgSPGgfHWnWAgiiuQgwsxFgSHTsEy0bw9A8tyrKs2Cvprnfl5/Wa2dE95z577nP3cHXfHJ29p1prAQCA3hw13xMAAID5IIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOjSovl64KVLl7YVK1bM18MDANCJbdu27WqtLZu+ft5CeMWKFRkfH5+vhwcAoBNV9ehM650aAQBAl4QwAABdEsIAAHRp3s4RBgBgds8991wmJiayZ8+e+Z7Kgrd48eIsX748o6OjA40XwgAAC9jExERe85rXZMWKFamq+Z7OgtVay+7duzMxMZGVK1cO9D1OjQAAWMD27NmT448/XgTPoqpy/PHHH9SRcyEMALDAieDBHOzPSQgDAPCKRkZGsmbNmrzpTW/K+973vjz77LNJkle/+tVDfZxrrrkmn/rUp4a6zVcihAEAjiRVw/0awJIlS3LvvffmgQceyNFHH51NmzbN8ZP80RDCAAAM7LzzzsuOHTtesu6ZZ57JRRddlLPPPjurV6/Obbfdtv++z33ucznjjDNy5pln5oMf/GCSZOfOnXnPe96Tc845J+ecc06+9rWv7R9/33335e1vf3tWrVqVG264Icm+X4T7pV/6pbzpTW/K6tWr8/nPf34oz2WgT42oqnVJ/lOSkSR/0Fr7rRnGXJjkd5OMJtnVWrtgKDMEAGBB2Lt3b+64446sW7fuJesXL16cW2+9Nccee2x27dqVc889N5dddlm2b9+e3/zN38zXvva1LF26NE8//XSS5GMf+1h+8Rd/MW9961vz2GOP5Z3vfGceeuihJMn999+fr3/96/m7v/u7nHXWWXnXu96Ve+65J/fee2/uu+++7Nq1K+ecc07OP//8vO51rzus5zNrCFfVSJLrkrwjyUSSrVW1ubW2fcqYn0zye0nWtdYeq6oTDmtWAAAsGN///vezZs2aJPuOCH/kIx95yf2ttXz84x/Pli1bctRRR+WJJ57IU089lT/7sz/Le9/73ixdujRJctxxxyVJ7rrrrmzfvj8l873vfS9/+7d/myRZv359lixZkiVLluRtb3tbvvGNb+TP//zP8/73vz8jIyN57WtfmwsuuCBbt27NZZdddljPa5AjwmuT7GitPZIkVXVzkvVJtk8Z83NJvthae2zyh/Gdw5oVAAALxovnCB/ITTfdlJ07d2bbtm0ZHR3NihUrsmfPnrTWZvwkhxdeeCH33HNPlixZ8rL7po+vqrTWDv9JzGCQc4RPSvL4lOWJyXVTvSHJT1XVV6tqW1X9/LAmCADAwvbd7343J5xwQkZHR3P33Xfn0UcfTZJcdNFF+cIXvpDdu3cnyf5TIy6++OJ85jOf2f/9UyP7tttuy549e7J79+589atf3X8axOc///k8//zz2blzZ7Zs2ZK1a9ce9rwHOSI8068TTs/yRUnenOSiJEuS3FNVX2+tffMlG6q6IskVSXLyyScf/GwBAFhwPvCBD+Rnf/ZnMzY2ljVr1uTUU09Nkpx++un51V/91VxwwQUZGRnJWWedlc9+9rP59Kc/nY9+9KM544wzsnfv3px//vn7P4li7dq1ede73pXHHnssv/Zrv5af/umfzoYNG3LPPffkzDPPTFXl2muvzYknnnjY867ZDjVX1VuSXNNae+fk8q8kSWvt308Zc3WSxa21ayaX/0uSP2mt/fcDbXdsbKyNj48f9hMAAPhx9tBDD+WNb3zjfE/jiDHTz6uqtrXWxqaPHeTUiK1JVlXVyqo6OsnGJJunjbktyXlVtaiqjknyD5M8dEizBwCAH4FZT41ore2tqquS3Jl9H592Y2vtwaq6cvL+Ta21h6rqT5Lcn+SF7PuItQfmcuIAAHA4Bvoc4dba7Ulun7Zu07Tl307y28ObGgAAzJ2BQhjo0ICX3cwcfaQNAMw1l1gGAKBLQhgAgC4JYQAAXtG3v/3tbNy4MaecckpOO+20XHrppfnmN795wPGvfvWrX3F7n/jEJ3LXXXcNe5oHzTnCAABHkEF/hWNQs/2qR2stGzZsyOWXX56bb745yb4rwT311FN5wxvecEiP+eu//uuH9H3D5ogwAAAHdPfdd2d0dDRXXnnl/nVr1qzJWWedlYsuuihnn312Vq9endtuu23G77/22muzevXqnHnmmbn66quTJB/60IfyR3/0R0mSFStWZNeuXUmS8fHxXHjhhUmSa665JpdffnkuvvjirFixIl/84hfzy7/8y1m9enXWrVuX55577rCfmyPCAAAc0AMPPJA3v/nNL1u/ePHi3HrrrTn22GOza9eunHvuubnssstSUw5Z33HHHfnSl76Uv/zLv8wxxxyTp59++qAe+6/+6q9y9913Z/v27XnLW96SW265Jddee202bNiQL3/5y3n3u999WM9NCAMAcNBaa/n4xz+eLVu25KijjsoTTzyRp556KieeeOL+MXfddVc+/OEP55hjjkmSHHfccQf1GJdccklGR0ezevXqPP/881m3bl2SZPXq1fnWt7512M9BCAMAcECnn376/tMYprrpppuyc+fObNu2LaOjo1mxYkX27NnzkjGttZccIZ7JokWL8sILLyTJy77/Va96VZLkqKOOyujo6P5tHXXUUdm7d+8hP6cXOUcYAIADevvb354f/OAHueGGG/av27p1ax599NGccMIJGR0dzd13351HH330Zd978cUX58Ybb8yzzz6bJDOeGrFixYps27YtSXLLLbfM0bOYmRAGAOCAqiq33nprvvKVr+SUU07J6aefnmuuuSaXXnppxsfHMzY2lptuuimnnnrqy7533bp1ueyyyzI2NpY1a9bkU5/61MvGfPKTn8zHPvaxnHfeeRkZGflRPKX9qs3T5VHHxsba+Pj4vDw2MACXWAZYEB566KG88Y1vnO9pHDFm+nlV1bbW2tj0sY4IAwDQJSEMAECXhDAAAF0SwgAAC9x8/U7XkeZgf05CGABgAVu8eHF2794thmfRWsvu3buzePHigb/HBTUAABaw5cuXZ2JiIjt37pzvqSx4ixcvzvLlywceL4QBABaw0dHRrFy5cr6n8WPJqREAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0KWBQriq1lXVw1W1o6qunuH+C6vqu1V17+TXJ4Y/VQAAGJ5Fsw2oqpEk1yV5R5KJJFuranNrbfu0of+ztfZP5mCOAAAwdIMcEV6bZEdr7ZHW2g+T3Jxk/dxOCwAA5tYgIXxSksenLE9MrpvuLVV1X1XdUVWnD2V2AAAwR2Y9NSJJzbCuTVv+X0l+prX2TFVdmuRLSVa9bENVVyS5IklOPvnkg5wqAAAMzyBHhCeSvH7K8vIkT04d0Fr7XmvtmcnbtycZraql0zfUWru+tTbWWhtbtmzZYUwbAAAOzyAhvDXJqqpaWVVHJ9mYZPPUAVV1YlXV5O21k9vdPezJAgDAsMx6akRrbW9VXZXkziQjSW5srT1YVVdO3r8pyXuT/Muq2pvk+0k2ttamnz4BAAALRs1Xr46NjbXx8fF5eWxgADXTrwfMwP/zArDAVdW21trY9PWuLAcAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJcGCuGqWldVD1fVjqq6+hXGnVNVz1fVe4c3RQAAGL5ZQ7iqRpJcl+SSJKcleX9VnXaAcf8hyZ3DniQAAAzbIEeE1ybZ0Vp7pLX2wyQ3J1k/w7h/leSWJN8Z4vwAAGBODBLCJyV5fMryxOS6/arqpCQbkmx6pQ1V1RVVNV5V4zt37jzYuQIAwNAMEsI1w7o2bfl3k/y71trzr7Sh1tr1rbWx1trYsmXLBp0jAAAM3aIBxkwkef2U5eVJnpw2ZizJzVWVJEuTXFpVe1trXxrKLAEAYMgGCeGtSVZV1cokTyTZmOTnpg5ora188XZVfTbJH4tgAAAWsllDuLW2t6quyr5PgxhJcmNr7cGqunLy/lc8LxiYRc109tE0bfrZSADA4RrkiHBaa7cnuX3auhkDuLX2ocOfFgAAzC1XlgMAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOiSEAYAoEsDhXBVrauqh6tqR1VdPcP966vq/qq6t6rGq+qtw58qAAAMz6LZBlTVSJLrkrwjyUSSrVW1ubW2fcqwP02yubXWquqMJF9IcupcTBgAAIZhkCPCa5PsaK090lr7YZKbk6yfOqC19kxrrU0u/kSSFgAAWMAGCeGTkjw+ZXlict1LVNWGqvo/Sb6c5J/PtKGqumLy1InxnTt3Hsp8AQBgKAYJ4Zph3cuO+LbWbm2tnZrk3Ul+Y6YNtdaub62NtdbGli1bdnAzBQCAIRokhCeSvH7K8vIkTx5ocGttS5JTqmrpYc4NAADmzCAhvDXJqqpaWVVHJ9mYZPPUAVX1D6qqJm+fneToJLuHPVkAABiWWT81orW2t6quSnJnkpEkN7bWHqyqKyfv35TkPUl+vqqeS/L9JP9syi/PAQDAglPz1atjY2NtfHx8Xh4bFpSa6TT8aebjdTrIvJL5mRsAHISq2tZaG5u+3pXlAADokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAurRovicAPypVs49pbXjbGnR7lYEGDXVug27rSDfM/cSRo5fXwEL9O40jw7D3+ZH6unNEGACALglhAAC6JIQBAOiSEAYAoEtCGACALglhAAC6JIQBAOjSQCFcVeuq6uGq2lFVV89w/weq6v7Jr7+oqjOHP1UAABieWUO4qkaSXJfkkiSnJXl/VZ02bdhfJ7mgtXZGkt9Icv2wJwoAAMM0yBHhtUl2tNYeaa39MMnNSdZPHdBa+4vW2v+bXPx6kuXDnSYAAAzXICF8UpLHpyxPTK47kI8kueNwJgUAAHNt0QBjZrp69IxXi66qt2VfCL/1APdfkeSKJDn55JMHnCIAAAzfIEeEJ5K8fsry8iRPTh9UVWck+YMk61tru2faUGvt+tbaWGttbNmyZYcyXwAAGIpBjghvTbKqqlYmeSLJxiQ/N3VAVZ2c5ItJPtha++bQZzlsNdNB7mnajAe9OVJV5QD/kAGDGeTvjWSwvzuGua1hW8hz48jgPZYjyKwh3FrbW1VXJbkzyUiSG1trD1bVlZP3b0ryiSTHJ/m92vcC2NtaG5u7aQMAwOEZ5IhwWmu3J7l92rpNU27/QpJfGO7UAABg7riyHAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRLCAAB0SQgDANAlIQwAQJeEMAAAXRoohKtqXVU9XFU7qurqGe4/taruqaofVNW/Hf40AQBguBbNNqCqRpJcl+QdSSaSbK2qza217VOGPZ3kXyd595zMEgAAhmyQI8Jrk+xorT3SWvthkpuTrJ86oLX2ndba1iTPzcEcAQBg6AYJ4ZOSPD5leWJyHQAAHLEGCeGaYV07lAerqiuqaryqxnfu3HkomwAAgKEYJIQnkrx+yvLyJE8eyoO11q5vrY211saWLVt2KJsAAIChGCSEtyZZVVUrq+roJBuTbJ7baQEAwNya9VMjWmt7q+qqJHcmGUlyY2vtwaq6cvL+TVV1YpLxJMcmeaGq/k2S01pr35vDuQMAwCGbNYSTpLV2e5Lbp63bNOX2t7PvlAkAADgiuLIcAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdEsIAAHRJCAMA0CUhDABAl4QwAABdGiiEq2pdVT1cVTuq6uoZ7q+q+vTk/fdX1dnDnyoAAAzPrCFcVSNJrktySZLTkry/qk6bNuySJKsmv65I8vtDnicAAAzVIEeE1ybZ0Vp7pLX2wyQ3J1k/bcz6JJ9r+3w9yU9W1euGPFcAABiaQUL4pCSPT1memFx3sGMAAGDBWDTAmJphXTuEMamqK7Lv1IkkeaaqHh7g8edHzfSUXtHSJLvmYCYMTc26jw5+t8/yiEPc3kLd1pB/aPPyOlqwP49h/wd5+P5+/yy8uR2UI3z6r+Qlr6EF+3faj/EOGMAR1Qs/Rn+l/cxMKwcJ4Ykkr5+yvDzJk4cwJq2165NcP8BjHnGqary1Njbf8+DA7KOFzz5a2Oyfhc8+Wvjso4VlkFMjtiZZVVUrq+roJBuTbJ42ZnOSn5/89Ihzk3y3tfY3Q54rAAAMzaxHhFtre6vqqiR3JhlJcmNr7cGqunLy/k1Jbk9yaZIdSZ5N8uG5mzIAABy+QU6NSGvt9uyL3anrNk253ZJ8dLhTO+L8WJ7y8WPGPlr47KOFzf5Z+Oyjhc8+WkBqX8MCAEBfXGIZAIAuCeHDUFXvq6oHq+qFqhqbdt+vTF5y+uGqeud8zZG/V1XXVNUTVXXv5Nel8z0n9pntMu7Mv6r6VlX978nXzvh8z4ekqm6squ9U1QNT1h1XVV+pqv87+edPzecce3eAfeS9aAERwofngST/NMmWqSsnL0G9McnpSdYl+b3JS1Uz//5ja23N5Nftsw9nrg14GXcWhrdNvnZ89NPC8Nnse4+Z6uokf9paW5XkTyeXmT+fzcv3UeK9aMEQwoehtfZQa22mi4KsT3Jza+0HrbW/zr5P01j7o50dHDEGuYw7ME1rbUuSp6etXp/kDydv/2GSd/9IJ8VLHGAfsYAI4bnhktML11VVdf/kP1f5J8OFwevlyNCS/I+q2jZ5lVAWpte++Dn+k3+eMM/zYWbeixYIITyLqrqrqh6Y4euVjlgNdMlphm+W/fX7SU5JsibJ3yT5nXmdLC/yejky/KPW2tnZdwrLR6vq/PmeEByhvBctIAN9jnDPWmv/+BC+baBLTjN8g+6vqrohyR/P8XQYjNfLEaC19uTkn9+pqluz75SWLa/8XcyDp6rqda21v6mq1yX5znxPiJdqrT314m3vRfPPEeG5sTnJxqp6VVWtTLIqyTfmeU7dm3xTeNGG7PtlR+bfIJdxZx5V1U9U1WtevJ3k4nj9LFSbk1w+efvyJLfN41yYgfeihcUR4cNQVRuS/Ocky5J8uaruba29c/IS1F9Isj3J3iQfba09P59zJUlybVWtyb5/dv9Wkn8xv9MhOfBl3Od5WrzUa5PcWlXJvveN/9pa+5P5nRJV9d+SXJhkaVVNJPlkkt9K8oWq+kiSx5K8b/5myAH20YXeixYOV5YDAKBLTo0AAKBLQhgAgC4JYQAAuiSEAQDokhAGAKBLQhgAgC4JYQAAuiSEAQDo0v8Hndbs0FpM3cQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 8))\n",
    "plt.hist(data[data['Treatment'] == 'Placebo']['Decrease'], \n",
    "         density=True, bins=50, color='red', label='Placebo')\n",
    "plt.hist(data[data['Treatment'] == 'Calcium']['Decrease'], \n",
    "         density=True, bins=50, color='blue', label='Calcium')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's calculate D -- creating the empirical CDFs $\\hat{F}_1$ and $\\hat{F}_2$, and then testing the value of $D$ for each possible $x$ in the empirical distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHSCAYAAAAezFYoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3xV1Z3//9dKIHIVEqCtcjHcvBMQUx3Rqki9VsU7tINCb/zq1PKt1vHSdvxNdfqbjrYd9WvVoVZorZpYvBQVL23VMm21Am0GAUkFQYlYJQmgYjBA1u+PA5kYAiQ5J9knyev5eOQRztlr7/0JofadlbU/K8QYkSRJkvS/cpIuQJIkSco2hmRJkiSpEUOyJEmS1IghWZIkSWrEkCxJkiQ1YkiWJEmSGumWdAFNGThwYCwsLEy6DEmSJHViS5YsqYwxDmrqWFaG5MLCQhYvXpx0GZIkSerEQghv7OmYyy0kSZKkRgzJkiRJUiOGZEmSJKmRtNYkhxDuBc4G3o0xHtnE8QDcBpwFfAjMiDH+pTX32rZtGxUVFWzdujWdkqVm6dGjB0OGDKF79+5JlyJJkhKQ7oN7c4E7gF/s4fiZwOidH8cCd+383GIVFRX07duXwsJCUtlbahsxRqqqqqioqGD48OFJlyNJkhKQ1nKLGONCoHovQyYDv4gpLwH9QwgHtOZeW7duZcCAAQZktbkQAgMGDPC3FpIkdWFt3QJuMLCuweuKne+93ZqLGZDVXvy3JkmZ89prcPHF4NxD5zFhy2+4/t1vksuOjFxv3fmzOKnknzJyrUxp65DcVNKITQ4MYSYwE2DYsGFtWVOr5ebmMmbMmPrXjz32GH379uWiiy5i0aJFzJgxgzvuuKNda/rKV77CVVddxeGHH96iY3vzwgsv8MMf/pAnnnhir+NOPvlkfvjDH1JcXMxZZ53FAw88QP/+/Zsce+uttzJz5kx69eq1z6+jT58+fPDBB82ut6ysjPXr13PWWWcBMH/+fFasWMF1113X7GtIktrOT34Cr74K55+fdCXKlCsW3soneYdln/xsRq6XN/STGblOJrV1SK4AhjZ4PQRY39TAGONsYDZAcXFxk0E6aT179qSsrOxj723ZsoWbbrqJZcuWsWzZsnav6Z577mny/R07duzxWFtYsGDBXo/feuutTJs2rcmQnG6tZWVlLF68uD4kn3vuuZx77rmtvp4kKXNqa+H++2HyZCgpSboaZUR1NXzyWbjySibcfHPS1bSZtm4BNx+4LKT8A7A5xtiqpRbZqnfv3pxwwgn06NFjr+OWLFnCSSedxNFHH83pp5/O22+n/hpOPvlkrrzySk488UQOO+wwFi1axAUXXMDo0aP57ne/C8DatWs59NBDmT59OkVFRVx00UV8+OGH9efv2p2wT58+3HDDDRx77LG8+OKLHzv29NNPM378eMaOHcukSZMAePnll5kwYQJHHXUUEyZMoLy8fK9fQ01NDVOnTqWoqIgpU6ZQU1NTf6ywsJDKykq2bNnC5z73OcaOHcuRRx5JaWkpt99+O+vXr2fixIlMnDhxn7UCfOtb32L8+PFMmjSJDRs27Pa1VlZWUlhYSG1tLTfccAOlpaWMGzeO0tJS5s6dyxVXXAHAG2+8waRJkygqKmLSpEm8+eabAMyYMYNZs2YxYcIERowYwbx58/b6tUuSWmfBAqishBkzkq5EGfPII7B9O0ydmnQlbSrdFnAPAicDA0MIFcD/C3QHiDHeDSwg1f5tFakWcF9M5367fPOb0GhCN23jxsGtt+59TE1NDePGjQNg+PDhPProo8269rZt2/jGN77Br3/9awYNGkRpaSnf+c53uPfeewHIy8tj4cKF3HbbbUyePJklS5ZQUFDAyJEjufLKKwEoLy/nZz/7Gccffzxf+tKXuPPOO7n66qs/dp8tW7Zw5JFHcuONN37s/Q0bNvDVr36VhQsXMnz4cKqrU89aHnrooSxcuJBu3brx29/+lm9/+9s8/PDDe/w67rrrLnr16sXSpUtZunQp48eP323M008/zYEHHsiTTz4JwObNm+nXrx8//vGPef755xk4cOBea911bPz48fzoRz/ixhtv5Hvf+94el7Hk5eVx4403snjx4voxc+fOrT9+xRVXcNlllzF9+nTuvfdeZs2axWOPPQbA22+/zR/+8AdWrlzJueeey0UXXbTHr12S1Dpz5sABB8BppyVdiTKmpARGjYKjjkq6kjaVVkiOMX5+H8cj8PV07pFNmlpu0Rzl5eUsW7aMU089FUgtLzjggP9t8rFracCYMWM44ogj6o+NGDGCdevW0b9/f4YOHcrxxx8PwLRp07j99tt3C8m5ublceOGFu93/pZde4sQTT6xvZ1ZQUACkAuz06dN57bXXCCGwbdu2vX4dCxcuZNasWQAUFRVRVFS025gxY8Zw9dVXc+2113L22Wfzmc98pslr7alWgJycHKZMmVL/tV5wwQV7rWtvXnzxRR555BEALr30Uq655pr6Y+eddx45OTkcfvjhvPPOO62+hySpae+8A08+Cd/6FnRr6wWeah/vvAPPPw/f/jZ08ofcO+Q/2X3N+GabGCNHHHEEL774YpPH99tvPyAVDnf9edfr7du3A7t3W2iq+0KPHj3Izc1t8v5Njf+Xf/kXJk6cyKOPPsratWs5+eST9/m17Kvrw8EHH8ySJUtYsGAB119/Paeddho33HBDs2vd2z27detGXV0dQKvbszWsv+HfdernOUlSJt1/P+zYAdOnJ12JMmbePKir6/RLLcBtqdvFIYccwoYNG+pD8rZt21i+fHmLrvHmm2/Wn//ggw9ywgknNPvc4447jt///vesWbMGoH65xebNmxk8eDDw8SUKe3LiiSdy//33A7Bs2TKWLl2625j169fTq1cvpk2bxtVXX81f/pLaYLFv3768//77zaq3rq6ufo3wAw88UP+1FhYWsmTJEoCPrSHe27UnTJhAyc4nRe6///4W/b1JklovxtRSi2OOgRY2WVI2Ky2FI45IfXRyhuQMKCws5KqrrmLu3LkMGTKEFStWfOx4Xl4e8+bN49prr2Xs2LGMGzeOP/3pTy26x2GHHcbPf/5zioqKqK6u5vLLL2/2uYMGDWL27NlccMEFjB07tn4pwzXXXMP111/P8ccfz44d++5zePnll/PBBx9QVFTEzTffzDHHHLPbmFdeeYVjjjmGcePG8f3vf7/+4cOZM2dy5pln1j+4tze9e/dm+fLlHH300Tz33HP1M9FXX301d911FxMmTKCysrJ+/MSJE1mxYkX9g3sN3X777cyZM4eioiLuu+8+brvttn3eX5KUvr/+FZYtgy9m5GkkZYWKCvjv/+4Ss8gAIRt/zVxcXBwbdjoAePXVVznssMMSqihZa9eu5eyzz06kxVxX1pX/zUlSur7xDfjpT+HttyE/P+lqlBH/+Z9w1VXwt7/B6NFJV5MRIYQlMcbipo45kyxJkjLqo4/ggQfgvPMMyJ1KSQmMH99pAvK+GJI7gMLCQmeRJUkdxhNPpPabcKlFJ/L66/Dyy11mqQUYkiVJUobNmQMHHgifzcyOxcoGDz2U+nzJJcnW0Y4MyZIkKWPefhuefhouuwya2elTHUFJCRx3HBx0UNKVtBtDsiRJyphf/jLVG9ltqDuRlSvhf/6nSy21AEOyJEnKkBhh7tzUhOMhhyRdjTKmtDS1u95FFyVdSbvqkDvuJSU3N5cxY8bUv37sscd47bXXuO6666itrSUvL49bbrmFU045pd1q+spXvsJVV13F4U10at/bsb154YUX+OEPf8gTTzyx13Enn3wyP/zhDykuLuass87igQceoH///k2OvfXWW5k5cya9evXa59fRp08fPvjgg2bXW1ZWxvr16znrrLMAmD9/PitWrOC6665r9jUkpa+8PNXRIAs7i6qdbN4MK1bAf/1X0pUoY2JMLbU46aTUQvMuxJDcAj179qSsrOxj723cuJHHH3+cAw88kGXLlnH66afz1ltvtVtN99xzT5Pv79ixY4/H2sKCBQv2evzWW29l2rRpTYbkdGstKytj8eLF9SH53HPP5dxzz2319SS1zlVXwYIFqQkndV2DB8POPavUGbzySmq5xTe/mXQl7c7lFmk66qijOHDnT1ZHHHEEW7du5aOPPtpt3JIlSzjppJM4+uijOf3003n77beB1GzslVdeyYknnshhhx3GokWLuOCCCxg9enT9bnVr167l0EMPZfr06RQVFXHRRRfx4Ycf1p+/a+OVPn36cMMNN3Dsscfy4osvfuzY008/zfjx4xk7diyTJk0C4OWXX2bChAkcddRRTJgwgfLy8r1+rTU1NUydOpWioiKmTJlCTU1N/bHCwkIqKyvZsmULn/vc5xg7dixHHnkkpaWl3H777axfv56JEyfW77i3t1oBvvWtbzF+/HgmTZrEhg0bdvtaKysrKSwspLa2lhtuuIHS0tL6Hffmzp3LFVdcAcAbb7zBpEmTKCoqYtKkSbz55psAzJgxg1mzZjFhwgRGjBjxsW2uJbXc+vWph7W+/W2oq/OjK39UVEC/fkn/i9QeTZkChYXN/zjllNQTmBdemFDByemYM8nf/CY0mtFN27hxcOutex1SU1PDuHHjABg+fDiPPvrox44//PDDHHXUUey3334fe3/btm184xvf4Ne//jWDBg2itLSU73znO9x7771AatvqhQsXcttttzF58mSWLFlCQUEBI0eO5MorrwSgvLycn/3sZxx//PF86Utf4s477+Tqq6/+2H22bNnCkUceyY033vix9zds2MBXv/pVFi5cyPDhw6murgbg0EMPZeHChXTr1o3f/va3fPvb3+bhhx/e49d/11130atXL5YuXcrSpUsZP378bmOefvppDjzwQJ588kkANm/eTL9+/fjxj3/M888/z8CBA/da665j48eP50c/+hE33ngj3/ve97jjjjuarCkvL48bb7yRxYsX14+ZO3du/fErrriCyy67jOnTp3Pvvfcya9YsHnvsMQDefvtt/vCHP7By5UrOPfdcLupia62kTPrlL1MByYe1pCy2YkWqldvJJ7esS8WnPw07//+7K+mYITkhTS232GX58uVce+21PPvss7sdKy8vZ9myZZx66qlAannBAQccUH9819KAMWPGcMQRR9QfGzFiBOvWraN///4MHTqU448/HoBp06Zx++237xaSc3NzubCJn/ReeuklTjzxRIYPHw5AQUEBkAqw06dP57XXXiOEwLZt2/b69S9cuJBZs2YBUFRURFFR0W5jxowZw9VXX821117L2WefzWc+85kmr7WnWgFycnKYsvN3ddOmTeOCCy7Ya1178+KLL/LII48AcOmll3LNNdfUHzvvvPPIycnh8MMP55133mn1PaSubtfDWscf32U24pI6ptJSyMmBBx+ET30q6WqyXscMyfuY8W1vFRUVnH/++fziF79g5MiRux2PMXLEEUfw4osvNnn+rpnnnJycj81C5+TksH37dgBCo0V+jV8D9OjRg9wmmlLGGJsc/y//8i9MnDiRRx99lLVr13LyySfv+Yvcy30bOvjgg1myZAkLFizg+uuv57TTTuOGG25odq17u2e3bt2oq6sDYOvWrc06d0/XAj72dx190khqtZdfhldfhZ/+NOlKJO3RrgfwTj7ZgNxMrklO06ZNm/jc5z7Hv//7v9fP9DZ2yCGHsGHDhvqQvG3bNpYvX96i+7z55pv15z/44IOccMIJzT73uOOO4/e//z1r1qwBqF9usXnzZgYPHgx8fInCnpx44oncf//9ACxbtoylS5fuNmb9+vX06tWLadOmcfXVV/OXv/wFgL59+/L+++83q966urr6NcIPPPBA/ddaWFjIkiVLAD62hnhv154wYQIlJSUA3H///S36e5PUPHPnQs+eXWojLqnjKSuDv/3NpypbwJCcpjvuuINVq1Zx0003MW7cOMaNG8e77777sTF5eXnMmzePa6+9lrFjxzJu3Dj+9Kc/teg+hx12GD//+c8pKiqiurqayy+/vNnnDho0iNmzZ3PBBRcwduzY+qUM11xzDddffz3HH388O3bs2Od1Lr/8cj744AOKioq4+eabOeaYY3Yb88orr3DMMccwbtw4vv/979c/fDhz5kzOPPPM+gf39qZ3794sX76co48+mueee65+Jvrqq6/mrrvuYsKECVRWVtaPnzhxIitWrKh/cK+h22+/nTlz5lBUVMR9993Hbbfdts/7S2q+mprUb24vvBD23z/paiTtUWkpdOsGaSxh7GpCNv6aubi4ODbsdADw6quvcthhhyVUUbLWrl3L2WefzbJly5IupUvpyv/mpOYqKYHPfx5++1vY2ThHUraJEUaMgEMPhaeeSrqarBJCWBJjLG7qmDPJkqRWmzsXhg2DZvySSFJSXn4Z1q7tcttKp8uQ3AEUFhY6iywp61RUwLPPwvTpqQfmJWWpkhLIy4Pzzku6kg7F/6xJklrlvvtSv8WdPj3pSiTtUV1dqjfymWe6y0sLdaiQnI3rp9U5+W9N2rtdvZFPPBGa6HwpKVv84Q+pLTFdatFiHSYk9+jRg6qqKsOL2lyMkaqqKnr06JF0KVLWevHFVDcpd9iTslxJSapH49lnJ11Jh9NhNhMZMmQIFRUVbNiwIelS1AX06NGDIUOGJF2GlLXmzoVevcDd3KUstn07zJsH55wDffokXU2H02FCcvfu3eu3VZYkJefDD1OTUxdfDH37Jl2NpD16/nnYsMGlFq3UYUKyJHVFtbWpZQ3Z5De/gfffd6mFlPVKSlI/yZ55ZtKVdEhpheQQwhnAbUAucE+M8QeNjvcDfgkM23mvH8YY56RzT0nqSr7+dbjnnqSr2N3w4amH9iRlqdpaeOSRVNs3n7FplVaH5BBCLvAT4FSgAlgUQpgfY1zRYNjXgRUxxnNCCIOA8hDC/THG2rSqlqQu4IMPUls+n3MOXHZZ0tV83Nix9kaWstqzz8KmTTBlStKVdFjpzCQfA6yKMb4OEEIoASYDDUNyBPqGEALQB6gGtqdxT0nqMubNgy1b4Npr4fjjk65GUodSUgL5+XDqqUlX0mGlE5IHA+savK4Ajm005g5gPrAe6AtMiTHWpXFPSeoy5syB0aNhwoSkK1GXde218POfJ12FWqOyEr74xdROe2qVdEJyaOK9xk2MTwfKgFOAkcBvQgj/HWN8b7eLhTATmAkwbNiwNMqSpI5v9WpYuBC+/30ITf3XVmprW7bAHXfA4YfD0UcnXY1aKjcXrrwy6So6tHRCcgUwtMHrIaRmjBv6IvCDmNoBZFUIYQ1wKPBy44vFGGcDswGKi4vdMURSl/aLX6TC8aWXJl2JuqzHH0/1+/vRj3xKU11SOo9dLAJGhxCGhxDygKmkllY09CYwCSCE8EngEOD1NO4pSZ1eXV3qN9ynngpDh+57vNQmSkvhwAPhhBOSrkRKRKtDcoxxO3AF8AzwKvBQjHF5COFrIYSv7Rx2EzAhhPAK8Dvg2hhjZbpFS1Jn9sIL8MYb9iFWgjZvhgUL4JJLbGOiLiutPskxxgXAgkbv3d3gz+uB09K5hyR1NXPnQr9+qfamUiIeeyzVZ9f2YerC/PFQkrLIe++lWr9NnQo9eyZdjbqs0lI46CA4tnHTKqnrMCRLUhb51a+gpsalFkpQVVVq7/EpU2ytoi7NkCxJWWTuXDj0UCfwlKBHHoHt21O/zpC6MEOyJGWJ116DP/whNYvsBJ4SU1KS2sVm3LikK5ESZUiWpCzx85+nGgnYG1mJ+fvfU+1Vpk71JzV1eYZkScoCO3akQvLpp6da00qJmDcv1ajbpRaSIVmSssFzz0FFhQ/sKWGlpXDkkamtqKUuLq0+yZKk3dXWwjnnwLp1zT9nwwbIz4dzz227uqS9WrcutSj+3/4t6UqkrGBIlqQMe/JJePbZ1NKJ/fdv/nnnnAM9erRdXdJe/epXqc9uICIBhmRJyrg5c+CAA+CJJ6Cb/5VVR1FSAkcfDaNGJV2JlBVckyxJGfTOO7BgQapDhQFZHcbq1bBokbPIUgOGZEnKoPvvT3Wq8AE8dSgPPZT6fMklydYhZRHnOSQpQ2JMLbU49lg47LCkq1GX9dZbcMIJsGlT88/54AM47jg46KC2q0vqYAzJkpQhf/kLLFsGd92VdCXq0h54ANauhcsvh+7dm3/eF77QZiVJHZEhWZIyZO5c2G8/92FQwkpL4dOfhjvvTLoSqUNzTbIkZcBHH6Um8M4/H/r3T7oadVmvvQZLlviTmpQBhmRJyoDHH4fqavjiF5OuRF1aaWnq88UXJ1uH1AkYkiUpA+bMgcGDYdKkpCtRl1Zamnpob+jQpCuROjxDsiSl6e234emn4bLLIDc36WrUZS1blvpwqYWUEYZkSUrTL38JdXX2RlbCSkshJwcuuijpSqROwZAsSWnY1Rt5wgQ4+OCkq1GXFWMqJE+cCJ/8ZNLVSJ2CIVmS0rBoEbz6qg/sKWF//Wuqs4XbSksZY0iWpDTMnQs9e9pMQAkrLYVu3eCCC5KuROo03ExEUiKqq+Huu6G2NulK0vPgg6lc0q9f0pWoy4oRSkrgtNNgwICkq5E6DUOypEQ88AB85ztJV5G+vDz4p39Kugp1aS+9BG++CTfdlHQlUqdiSJaUiDVrUssUtmyBEJKuRurASktT+6FPnpx0JVKn4ppkSYlYuxYKCw3IUlp27ICHHoKzznLNj5RhaYXkEMIZIYTyEMKqEMJ1exhzcgihLISwPITw+3TuJ6nz2BWSJaXhv/87tZuNXS2kjGv1cosQQi7wE+BUoAJYFEKYH2Nc0WBMf+BO4IwY45shhE+kW7CkzmHtWjj22KSrkDq40lLo1QvOPjvpSqROJ52Z5GOAVTHG12OMtUAJ0HhB1BeAR2KMbwLEGN9N436SOon33kt1t3AmWUrD9u0wbx6ccw707p10NVKnk86De4OBdQ1eVwCN54UOBrqHEF4A+gK3xRh/kcY9JXUCb7yR+mxIVpv79rfhP/8z6SraRozw0UcwdWrSlUidUjohuanHbWIT1z8amAT0BF4MIbwUY/zbbhcLYSYwE2DYsGFplCUp261dm/psSFabqq2Fu+6CMWNS2zV3Rv37u9RCaiPphOQKYGiD10OA9U2MqYwxbgG2hBAWAmOB3UJyjHE2MBuguLi4cdiW1ImsWZP6bEhWm3r2Wdi0Cf71X1PdHySpBdJZk7wIGB1CGB5CyAOmAvMbjfk18JkQQrcQQi9SyzFeTeOekjqBtWtTzxoNGpR0JerUSkogPx8++9mkK5HUAbV6JjnGuD2EcAXwDJAL3BtjXB5C+NrO43fHGF8NITwNLAXqgHtijMsyUbikjsseyWpzNTXw61+n1uvm5SVdjaQOKK0d92KMC4AFjd67u9HrW4Bb0rmPpM7FHslqcwsWwAcf2D9YUqu5456kdmdIVpsrLYVPfAJOPjnpSiR1UIZkSe1q82bYuNGQrDb0wQfwxBNw0UXQLa1fmErqwgzJktqVPZLV5h5/PLUm2f7BktJgSJbUruyRrDZXUgKDB8PxxyddiaQOzJAsqV0ZktWmNm2Cp5+GSy6BHP8vTlLr+V8QSe1qzZpUj+SBA5OuRJ3SY4+ldtpzqYWkNBmSJbUreySrTZWUwPDh8OlPJ12JpA7OkCypXa1dm8owUsZVVsJvf5vqjexPYZLSZEiW1K7skaw28/DDsGOHG4hIyghDsqR2s2lT6sOQrDZRWgqHHAJjxyZdiaROwJAsqd3YI1lt5u234YUXXGohKWMMyZLaje3f1GbmzYMYXWohKWPcr1NSu8n6kHzvvTBzJtTVJV2JWipGGDMGDj886UokdRKGZEntZs0a6N0bBgxIupI9uPvuVIL/wheSrkStcfbZSVcgqRMxJEtqN1ndI3n1ali0CG6+Gf75n5OuRpKUMNckS2o3Wd0j+aGHUp8vuSTZOiRJWcGQLKndZHWP5JISOO44OOigpCuRJGUBQ7KkdrFpE2zenKUh+dVXYelSmDo16UokSVnCkCypXWR1Z4vS0tRC6YsvTroSSVKWMCRLahdZG5JjTC21OOkkOOCApKuRJGUJQ7KkdpG1Ifl//gfKy11qIUn6GEOypHaxdi306QMFBUlX0khpKeTmwoUXJl2JJCmLGJIltYs1a7KwR/KupRaf/SwMHJh0NZKkLGJIltQusrL926JFqcJcaiFJasSQLKnNxZilG4mUlEBeHpx3XtKVSJKyjCFZUpvbtAneey/LZpLr6lK77J1xBvTvn3Q1kqQsY0iW1OaysrPFH/8Ib73lUgtJUpPSCskhhDNCCOUhhFUhhOv2Mu7TIYQdIYSL0rmfpI4pK0NySQn07AnnnJN0JZKkLNTqkBxCyAV+ApwJHA58PoRw+B7G/QfwTGvvJaljy7qQvH07zJsHZ5+d6ksnSVIj3dI49xhgVYzxdYAQQgkwGVjRaNw3gIeBT6dxL0ntra4uFSYzYN1qKOgD+b2B2oxcMj3PPw/vvutSC0nSHqUTkgcD6xq8rgCObTgghDAYOB84BUOy1HHU1sJhh8Hrr2fkcj/e+UGPjFwuM/r0gTPPTLoKSVKWSickN7UlQGz0+lbg2hjjjrCPHQRCCDOBmQDDhg1LoyxJafvtb1MB+f/5fyAD/3u8/XbIz4dLL81AbZlSXJxakyxJUhPSCckVwNAGr4cA6xuNKQZKdgbkgcBZIYTtMcbHGl8sxjgbmA1QXFzcOGxLak+lpam2aLfdBvvtl9alYoTv/gC+OAUu/XaG6pMkqY2lE5IXAaNDCMOBt4CpwBcaDogx1m8dEEKYCzzRVECWlEW2boVHH4WLL047IANs3Ajvv59FD+1JktQMrQ7JMcbtIYQrSHWtyAXujTEuDyF8befxuzNUo6T29NRTqVQ7ZUpGLpd1nS0kSWqGdGaSiTEuABY0eq/JcBxjnJHOvSS1k9JSGDgQTjklI5czJEuSOqK0QrKk7PLuu/Czn7W+c1v32i1c9cjjLB13GU/9e2b+8/DnP6c+G5IlSR2JIVnqRH76U/jud1t//iU8QR4f8q1FU1m4KHN1HXFE6jlASZI6CkOy1In89a8wciSUl7fu/JwLSxZ7O28AACAASURBVIiLDuS5NSeknjTIkJwc2EcXSEmSsoohWepEysrgqKMgtzUBd/NmeGoBXH45uXkZTMiSJHVAOUkXICkz3nsPVq+GceNaeYFf/zq1055bNUuSZEiWOov/+Z/U51aH5JISOOggOPbYfY+VJKmTMyRLnURZWerzUUe14uSqKvjNb1K9kV08LEmSIVnqLMrKYNAgOOCAVpz8yCOpvnEutZAkCTAkS51GWVlqqUWrJoJLSmD06DTWakiS1LkYkqVOYNs2WLaslRn373+HF15IzSK71EKSJMCQLHUKr76aakzRqpA8bx7U1aXWI0uSJMA+yVKn8LGH9rZsgR07mn/ygw/CkUemtsWTJEmAIVnqFMrKoGdPOGTx/XD4tJZf4KabMl+UJEkdmCFZ6gTKymDMGMiZ8zMoLIRvfKP5J3fvDtOnt1ltkiR1RIZkqYOLMRWSv/y5t+H+F+CGG+Cqq5IuS5KkDs0H96QObt062LgRJn/0q1Ri9gE8SZLSZkiWOri//jX1eezKUigqgsMOS7YgSZI6AUOy1MGVlcEw3qTvK39yxzxJkjLEkCx1cGVl8PVBD6VeuNRCkqSMMCRLHVxZGVy4vQQ+/WkYMSLpciRJ6hQMyVIHtmkT5K5dxciNS5xFliQpgwzJUgdWVgZTKE29uOSSZIuRJKkTMSRLHVhZGUylhNpjT4ChQ5MuR5KkTsOQLHVg7z6/nDEsI2+aSy0kScokQ7LUgR3051J2kAMXXZR0KZIkdSqGZKmDqv0ocso7Jaw56GT41KeSLkeSpE7FkCx1UK8/UsZoXmPjaW4gIklSpqUVkkMIZ4QQykMIq0II1zVx/B9DCEt3fvwphDA2nftJ+l/b7ithG93I//IFSZciSVKn0+qQHELIBX4CnAkcDnw+hHB4o2FrgJNijEXATcDs1t5PUgMxcuAfSnku91SGFw9IuhpJkjqdbmmcewywKsb4OkAIoQSYDKzYNSDG+KcG418ChqRxP6nzqq2FzZubP76sjAHvv8GfR93I6bltV5YkSV1VOiF5MLCuwesK4Ni9jP8y8FQa95M6pxihuBheeaVFp21lPzadOLmNipIkqWtLJySHJt6LTQ4MYSKpkHzCHi8WwkxgJsCwYcPSKEvqYJYsSQXkr34VxjZv2X5lJVz8r4fz+WP7tXFxkiR1TemE5Aqg4RZfQ4D1jQeFEIqAe4AzY4xVe7pYjHE2O9csFxcXNxm2pU6ptBS6d4cf/AAKCpp1ysJH4AXg5qPatDJJkrqsdELyImB0CGE48BYwFfhCwwEhhGHAI8ClMca/pXEvKeusXAmXXQYffdT6a4RYx1MrS/lbj9OZNbF5ARlgwwbIyYEjj2z9vSVJ0p61OiTHGLeHEK4AngFygXtjjMtDCF/befxu4AZgAHBnCAFge4yxOP2ypeT93/8LS5fCmWe2/hqHVL3EAdvWUVr0/zFi6L7H7zJiBBx9NPTs2fp7S5KkPQsxZt/KhuLi4rh48eKky5D2aOtWOOAAOOssuP/+NC40axb89Kfwzjuw//4Zq0+SJO1bCGHJniZw3XFPaoX582HTJpgxI42L7NgBv/pVKmkbkCVJyiqGZKkV5syBoUPhlFPSuMjChfD3v8NUt5WWJCnbGJKlFnrrLXj22dRDe7npbORRUgK9e8PnPpex2iRJUmYYkqUW+uUvoa4uzaUW27bBww/DuedCr16ZKk2SJGWIIVlqgRhTSy1OOAFGjUrjQr/7HVRVudRCkqQsZUiWWuDPf4by8jRnkSG1gUi/fnD66ZkoS5IkZZghWWqBuXNTqyMuuSSNi3z0ETzyCJx/Puy3X6ZKkyRJGWRIlpqppib1rN2FF0Lfvmlc6Omn4b33YMqUjNUmSZIyy5AsNdNjj8HmzfDFL6Z5odJSGDAAJk3KSF2SJCnzDMlSM82dCwcdBCedlMZFPvwwtRPJhRdC9+6ZKk2SJGWYIVlqhnXr4De/genTISed/9U8+SRs2WJXC0mSsly3pAuQOoL77ku1f5s+vdGBd99NPYjXkgt96lNw4okZrU+SJGWWIVnahxhTSy1OOglGjGhw4Kmn4KyzWn7BWbPS3KpPkiS1NUNyO6mrg//6L6iuTroStVRVFbz2GnznO40OzJ0LAwfCD37Q/Ivl5qZ22ZMkSVnNkNxOnn0W/umfkq5CrXXggaln7ept2QJPPJFaf/HlLydWlyRJahuG5HYyZ06q69cbb0BeXtLVqKVycxs9sPf446lOFfY6liSpUzIkt4ONG1M9dr/2NejdO+lqlBGlpanp5RNOSLoSSZLUBmwB1w5KSqC2FmbMSLoSZcTmzbBgQWpvah/AkySpUzIkt4M5c2DsWDjqqKQrUUY89ljqpx6XWkiS1GkZktvY8uWwaJGzyJ1KaWlq671jj026EkmS1EYMyW1s7lzo1g3+8R+TrkQZUVWV2npvyhQIIelqJElSGzEkt6Ht21MbrJ19NgwalHQ1yohHHkl9Y91WWpKkTs2Q3IaeeQbeecelFp1KSQkcfDCMG5d0JZIkqQ0ZktvQnDmpGeTW7FysLPT3v8MLL7jUQpKkLsCQ3EaqqmD+fJg2Dbp3T7oaZcS8ean9xV1qIUlSp2dIbiMPPgjbtrnUolMpKYEjj4TDD0+6EkmS1MYMyW1kzpxUX+SioqQrUUasWwd//KOzyJIkdRFpheQQwhkhhPIQwqoQwnVNHA8hhNt3Hl8aQhifzv06iqVL4S9/gS9+MelKlDEPPZT67AYikiR1Ca0OySGEXOAnwJnA4cDnQwiNfw99JjB658dM4K7W3q8jmTs3tQ75859PuhJlTEkJHH00jBqVdCWSJKkddEvj3GOAVTHG1wFCCCXAZGBFgzGTgV/EGCPwUgihfwjhgBjj22nct028vfgtajfXpH2dujr448/hqxNh4CZgU/q1KWFvvw2LF8MttyRdiSRJaifphOTBwLoGryuAxvv0NjVmMJB9Ifm0yxi/8bmMXOvPAM+Smj9X5xACXHJJ0lVIkqR2kk5IbqpRbGzFmNTAEGaSWpLBsGHD0iirla67nj+uycwi4h49YPx4W+l2KkOHQhL/LiVJUiLSCckVwNAGr4cA61sxBoAY42xgNkBxcXGTQbotjb/ms+19S0mSJGWpdLpbLAJGhxCGhxDygKnA/EZj5gOX7exy8Q/A5mxcjyxJkiQ11OqZ5Bjj9hDCFcAzQC5wb4xxeQjhazuP3w0sAM4CVgEfAjZFkyRJUtZLZ7kFMcYFpIJww/fubvDnCHw9nXtIkiRJ7c0d9yRJkqRGDMmSJElSI4ZkSZIkqRFDsiRJktRISD1bl11CCBuAN5KuI8sMBCqTLkJ+H7KI34vs4PchO/h9yA5+H7JDS74PB8UYBzV1ICtDsnYXQlgcYyxOuo6uzu9D9vB7kR38PmQHvw/Zwe9DdsjU98HlFpIkSVIjhmRJkiSpEUNyxzE76QIE+H3IJn4vsoPfh+zg9yE7+H3IDhn5PrgmWZIkSWrEmWRJkiSpEUOyJEmS1IghWZIkSWrEkCxJkiQ1YkiWJEmSGjEkS5IkSY0YkiVJkqRGDMmSJElSI4ZkSZIkqRFDsiRJktSIIVmSJElqxJAsSZIkNWJIliRJkhoxJEuSJEmNGJIlSZKkRgzJkiRJUiOGZEmSJKkRQ7IkSZLUiCFZkiRJasSQLEmSJDViSJYkSZIaMSRLkiRJjRiSJUmSpEYMyZIkSVIjhmRJkiSpEUOyJEmS1IghWZIkSWrEkCxJkiQ1YkiWJEmSGjEkS5IkSY0YkiVJkqRGDMmSJElSI4ZkSZIkqZFuSRfQlIEDB8bCwsKky5AkSVIntmTJksoY46CmjmVlSC4sLGTx4sVJlyFJkqROLITwxp6OudxCkiRJasSQLEmSJDViSJYkSZIayco1yU3Ztm0bFRUVbN26NelSsl6PHj0YMmQI3bt3T7oUSZKkDqnDhOSKigr69u1LYWEhIYSky8laMUaqqqqoqKhg+PDhSZcjSZLUIXWY5RZbt25lwIABBuR9CCEwYMAAZ9wlSZLS0GFmkgEDcjP59yQpCW+9/BYv3/EyxGTuf+QXjmT0maOTubmkTqdZITmEcAZwG5AL3BNj/MEexn0aeAmYEmOc15JzJUkd26I7F7GsZBn9hvZr93u///b7bF632ZAsKWP2GZJDCLnAT4BTgQpgUQhhfoxxRRPj/gN4pqXndhQnn3wyc+fOpbCwkBkzZvD73/+efv1S/2fwpS99iVmzZvGd73yHX/ziF2zcuJEPPvigxfdo6rpf+cpXuPjii1m9ejW5ubmcc845/OAHqZ81XnjhBebOncvcuXMz9nVKUmtUrqzkoM8cxGW/u6zd7/3Y9MdY89yadr+vpM6rOWuSjwFWxRhfjzHWAiXA5CbGfQN4GHi3Fed2SLfccgtlZWWUlZUxa9YsAM455xxefvnljF/36quvZuXKlfz1r3/lj3/8I0899VTa9UtSpsQYqSqvYsChAxK5f/6ofN6reI9tNdsSub+kzqc5IXkwsK7B64qd79ULIQwGzgfubum5nc0//MM/cMABB2T0mr169WLixIkA5OXlMX78eCoqKjJ6D0lKx5Z3t7B101YGHjIwkfsXjCoAYOPrGxO5v6TOpzlrkpt6CqzxYxm3AtfGGHc0emisOeemBoYwE5gJMGzYsL0W9PQ3n+bvZX/f65iW+tS4T3HGrWe06Jx//ud/5t/+7d8AuO+++xgzZkxGatnbdTdt2sTjjz/O//k//ycj95KkTKgqrwJg4KHJhuTqVdV84ohPJFKDpM6lOSG5Ahja4PUQYH2jMcVAyc6APBA4K4SwvZnnAhBjnA3MBiguLk7o2eiWueWWW7jooova7brbt2/n85//PLNmzWLEiBEZv68ktVblykoABhySzHKL+pnk1c4kS8qM5oTkRcDoEMJw4C1gKvCFhgNijPW7VoQQ5gJPxBgfCyF029e5rdHSGd9sc/rpp/POO+9QXFzMPffc0+zzZs6cyejRo/nmN7/ZhtVJUstVllfSrWe3RDpbAPTM70mP/B5Ur6pO5P6SOp99huQY4/YQwhWkulbkAvfGGJeHEL6283jjdcj7PDczpXdczzzzzL4HNfLd736XzZs3tyhUS1J7qSqvYsDBAwg5yfVpLxhVYEiWlDHN2nEvxrggxnhwjHFkjPH7O9+7u6mAHGOcsatH8p7O7cyuueYahgwZwocffsiQIUP413/917SvWVFRwfe//31WrFjB+PHjGTdunGFZUlapXFmZ2EN7uxiSJWVSh9pxL5vsqS/xzTffzM0335zR6w4ZMoQYO8QybUld0PaPtrNpzSbG/GNmHl5urYJRBSwvXc6O2h3k5uUmWoukjq9ZM8mSJO1J9apqYl3MipnkWBfZtHZTonVI6hwMyS0wY8YM+vfvn3QZH1NYWMh5552XdBmSurBd7d+S6myxS30buNUuuZCUPpdbtMCMGTOSLmE3hYWFFBYWJl2GpC5sV/u3pGeS80fmA7guWVJGdKiZZNflNo9/T5LaU1V5FX0H9yWvT16idfT+RG/y+uQZkiVlRIcJyT169KCqqsoAuA8xRqqqqujRo0fSpUjqIirLKxPbaa+hEAIFowrYuMoNRSSlr8MstxgyZAgVFRVs2LAh6VKyXo8ePRgyZEjSZUjqAmKMVK6sTLyzxS4Fowp4Z+k7SZchqRPoMCG5e/fuDB8+fN8DJUntZsu7W/ho80dZMZMMkD8qn5W/XkndjjpycjvML0slZSH/CyJJarVseWhvl4JRBdRtq+O9de8lXYqkDs6QLElqtWxp/7ZLwcidbeB8eE9SmgzJkqRWqyyvpFvPbvQb2i/pUoAGvZINyZLSZEiWJLVa1coqBhw8gJATki4FgL4H9qVbj26GZElpMyRLklotW9q/7RJyAvkj8w3JktLWrJAcQjgjhFAeQlgVQriuieOTQwhLQwhlIYTFIYQTGhxbG0J4ZdexTBYvSUrO9o+2s2nNpqxZj7xLwagCQ7KktO2zBVwIIRf4CXAqUAEsCiHMjzGuaDDsd8D8GGMMIRQBDwGHNjg+McZYmcG6JUkJq15VTayLWTWTDKmQvPrZ1cS6mDXLQCR1PM2ZST4GWBVjfD3GWAuUAJMbDogxfhD/dyu83oDb4klSJ5dt7d92yR+Zz/aa7bz/9vtJlyKpA2tOSB4MrGvwumLnex8TQjg/hLASeBL4UoNDEXg2hLAkhDAznWIlSdmjvv3bwdm33ALscCEpPc0JyU39rmq3meIY46MxxkOB84CbGhw6PsY4HjgT+HoI4cQmbxLCzJ3rmRe79bQkZb+q8ir2H7I/eX3yki7lYwzJkjKhOSG5Ahja4PUQYP2eBscYFwIjQwgDd75ev/Pzu8CjpJZvNHXe7BhjcYyxeNCgQc0sX5KUlMqVlVn30B5Av6H9yOmeY0iWlJbmhORFwOgQwvAQQh4wFZjfcEAIYVQIIez883ggD6gKIfQOIfTd+X5v4DRgWSa/AElS+4sxZl37t11yuuWQPzyfjas2Jl2KpA5sn90tYozbQwhXAM8AucC9McblIYSv7Tx+N3AhcFkIYRtQA0zZ2enik8CjO/NzN+CBGOPTbfS1SJLayZZ3tvDR5o+yciYZdraBW+1MsqTW22dIBogxLgAWNHrv7gZ//g/gP5o473VgbJo1SpKyTGV5dna22CV/ZD5v/PcbxBjZOVEjSS3ijnuSpBbb1dkiG5dbQGomufb9Wj7c8GHSpUjqoAzJkqQWq1xZSbee3dh/yP5Jl9IkO1xISpchWZLUYlXlVQw8ZGDW7mhnSJaULkOyJKnFsrX92y79C/sTcoIhWVKrGZIlSS2yfet2Nq3dlLXrkQFy83Lpd1A/Nq62DZyk1jEkS5JapHpVNbEuZvVMMkDByAJnkiW1miFZktQi2d7+bZf8UfmGZEmtZkiWJLXIrvZvAw7O8pnkUQXUVNdQU12TdCmSOqBmbSYiSUrG6mdX8/L/fZkYY9Kl1Ht32bvsP2R/8vrkJV3KXtV3uFhdzeCCwQlXI6mjMSRLUhZbfPdi1jy3hoGHZc/Shl4De3Ho+YcmXcY+NWwDN/jThmRJLWNIlqQsVrmykhGnjmDqY1OTLqXDyR+RD2CHC0mt4ppkScpSddvrqF5VnfVdJLJV957d6Tu4rw/vSWoVZ5IlKUttWruJum11Wd2PONsVjCrgb0/8jfvPuj/pUtpEz/yenPPTc+jeq3vSpUidTrNCcgjhDOA2IBe4J8b4g0bHJwM3AXXAduCbMcY/NOdcSVLTKld2jFZr2WzcjHEsunMRH1Z+mHQpGVf7QS2rnlrF+JnjKTypMOlypE5nnyE5hJAL/AQ4FagAFoUQ5scYVzQY9jtgfowxhhCKgIeAQ5t5riSpCbv6EbvcovXGzRjHuBnjki6jTWxcs5HbR9xO9apqQ7LUBpqzJvkYYFWM8fUYYy1QAkxuOCDG+EH83/5EvYHY3HMlSU2rXFlJr4G96DWgV9KlKAv1G9qPnO45PpgotZHmhOTBwLoGryt2vvcxIYTzQwgrgSeBL7XkXEnS7qrKq5xF1h7ldMshf7i7CkptpTkhOTTx3m5d7WOMj8YYDwXOI7U+udnnAoQQZoYQFocQFm/YsKEZZUlS51ZVXuVDe9qr/JGGZKmtNCckVwBDG7weAqzf0+AY40JgZAhhYEvOjTHOjjEWxxiLBw0a1IyyJKnzqtlYw5Z3tziTrL0qGFVA9arqrNqRUeosmhOSFwGjQwjDQwh5wFRgfsMBIYRRIYSw88/jgTygqjnnSpJ2V1VeBeBMsvaqYFQBte/X8uGGzte9Q0raPrtbxBi3hxCuAJ4h1cbt3hjj8hDC13Yevxu4ELgshLANqAGm7HyQr8lz2+hrkaROw/Zvao6GW2/3/kTvhKuROpdm9UmOMS4AFjR67+4Gf/4P4D+ae64kae8qyyvJ6Z5D/+H9ky5FWaw+JK+uZuiEofsYLakl3JZakrJQ1coqCkYWkNs9N+lSlMX6F/Yn5AQf3pPagCFZkrJQZXmlD+1pn3Lzcuk3rB8bV9krWco0Q7IkZZm67XVUr6r2oT01y64OF5Iyy5AsSVlm45qN1G2rcyZZzZI/yl7JUlswJEtSlrH9m1qiYFQBNdU11GysSboUqVMxJEtSlrH9m1piV4eLjatdlyxlkiFZkrJMZXklvQb2omdBz6RLUQfQsFeypMwxJEtSlqkqr3KphZotf0Q+YEiWMs2QLElZpnKl7d/UfN17dqfv4L6GZCnDDMmSlEVqqmv4cMOHziSrRWwDJ2WeIVmSskhleeqhPWeS1RKGZCnzDMmSlEVs/6bWKBhVwJZ3tlD7QW3SpUidhiFZkrJI5cpKcrrnkD88P+lS1IHUd7hY7WyylCnNCskhhDNCCOUhhFUhhOuaOP6PIYSlOz/+FEIY2+DY2hDCKyGEshDC4kwWL0mdTVV5FQUjC8jp5hyGmi9/pB0upEzrtq8BIYRc4CfAqUAFsCiEMD/GuKLBsDXASTHGjSGEM4HZwLENjk+MMVZmsG5J6pQqyytdaqEWKxhpr2Qp05ozVXEMsCrG+HqMsRYoASY3HBBj/FOMcddWPy8BQzJbpiR1fnXb66heVe1De2qx/fbfj96f6G1IljKoOSF5MLCuweuKne/tyZeBpxq8jsCzIYQlIYSZLS9RkrqGjWs2UretzplktUrBqAI2rnJrailT9rncAghNvBebHBjCRFIh+YQGbx8fY1wfQvgE8JsQwsoY48Imzp0JzAQYNmxYM8qSpM6lcqXt39R6BaMKWPP8mqTLkDqN5swkVwBDG7weAqxvPCiEUATcA0yOMVbtej/GuH7n53eBR0kt39hNjHF2jLE4xlg8aNCg5n8FktRJ1Ld/O8SZZLVc/qh83lv3HttqtiVditQpNCckLwJGhxCGhxDygKnA/IYDQgjDgEeAS2OMf2vwfu8QQt9dfwZOA5ZlqnhJ6kwqyyvpNagXPQt6Jl2KOqBdD+9tWrMp4UqkzmGfyy1ijNtDCFcAzwC5wL0xxuUhhK/tPH43cAMwALgzhACwPcZYDHwSeHTne92AB2KMT7fJVyJJHVzVyipnkdVq9b2SV1Uz6HB/IyulqzlrkokxLgAWNHrv7gZ//grwlSbOex0Y2/h9SdLuKssrOWTyIUmXoQ6qYUiWlL5mhWRJyrQNKzbw1DeeYkftjkTuX3hKIRO/NzGRezelprqGDzd86EyyWq1nQU965PcwJEsZ4pZOkhLx+u9eZ81zawg5gdy83Hb92PzmZv5825+JsclGPYmoLE91trD9m9JRMKqAjattAydlgjPJkhJRU10DwGW/u6zdt2B+6daXeObKZ6ipqqHXwF7teu89sf2bMqFgVAFv/fmtpMuQOgVnkiUloqa6hv3236/dAzI0WLu5Ont+LV1VXkVO9xzyh+cnXYo6sPyR+WxauymxZUxSZ+JMsqREbN24NbFWZ/kjU0G0elU1Q44dkvHrx7rIo5c+yqa1zW/FVfVaFQWjChL5oUGdR8GoAmJdZNMbmxgw2t9KSOkwJEtKRE11TXIheXg+hLbrArBxzUZeeeAVPjHmE/T5ZJ9mnfOpsZ/i8IsPb5N61HU07HBhSJbSY0iWlIgkQ3K3Ht3oN7QfG1e1zQNOu3bOO/vusxk6Yeg+RkuZYxs4KXP8vZ6kRCQZkiEVJtoqSPgQnpLS+xO9yeuTZ4cLKQMMyZISUVNdQ4+CHondP39UftuF5PJKeg7oSa8B2dE5Q11HCKFNfwCUuhJDsqR2F2PMipnkDys/ZOvmrRm/dlV5lf2OlZj8kW33A6DUlRiSJbW72vdriTsiPfMTDMkjU2s32+LX0pUrK11qocQUjCpg4+sbqdtRl3QpUodmSJbU7nZtJJL0TDJk/gGnrZu2suWdLc4kKzEFowqo21bHe+veS7oUqUOzu4WkdpcNIblhr+RMqt9e+hBDspKx6wfAB85+gLw+ec0+77irjuOIS45oq7KkDqdZM8khhDNCCOUhhFUhhOuaOP6PIYSlOz/+FEIY29xzJXU92RCS83rn0eeAPhkPybvav7ncQkkZfMxgjpx6JPsP2Z8e/Xs066OqvIql9y1NunQpq+xzJjmEkAv8BDgVqAAWhRDmxxhXNBi2BjgpxrgxhHAmMBs4tpnnSupisiEkQ9u0gassrySnWw75I9xeWsno3qs7Fz54YYvOeejCh9iwYkMbVSR1TM2ZST4GWBVjfD3GWAuUAJMbDogx/inGuOvpl5eAIc09V1LXU7Mxe0Jyph/cq1pZRf7IfHK752b0ulJbyh+V78N+UiPNWZM8GFjX4HUFcOxexn8ZeKqV50rqAnbNJPfIT65PMqTWJb+//n1qt9SS17v5azf3prK80of21OEUjCpgR+0O3qt4j/4H9c/49Tev28yj0x5lW822jF9bTcvplsPpPz6dIf8wZN+D1aTmhOTQxHuxyYEhTCQVkk9oxbkzgZkAw4YNa0ZZkjqqmuoauvXsRvee3ROtY9cDThtf38gnx3wy7evV7aij+rVqRn9udNrXktpTw24vbRGSX//t67yx8A2GnzKc3P38LUt7WP3salY+ttKQnIbmhOQKYGiD10OA9Y0HhRCKgHuAM2OMVS05FyDGOJvUWmaKi4ubDNKSOoekNxLZpWEwyERI3rR2EztqdziTrA6n4f8WRkwakfHrV5VXkdM9h2nPTCOnm91n28Mdh97hpjJpas6/1EXA6BDC8BBCHjAVmN9wQAhhGPAIcGmM8W8tOVdS17O1emt2hOSRme2VXLnS9m/qmPYfvD+5++W2yeY6kArJA0YPMCC3o4KRmX/moqvZ50xyjHF7COEK4BkgF7g3xrg8hPC1ncfvBm4ABgB3hhAAtscYi/d0bht9LZI6iGyZSe7Rvwe9Bv7/7d17kFzleefx79M995Fmpo+u9swISdMSSFBIpmQhDIshhItssGKyqeCqmuvEOAAAGWtJREFUlEkcIlOB8rp2kypi11aSTbK1rrWdTVwEgh1ss+WYpMqLGWdFsJcQy0mBkbgLa0YaCQmNJEvq7hlde+7v/tE98ngY0Onp031Od/8+VSr17T16VGda/dPb73nelsBCstq/SaWymJFYXbrtrFN9KRav038eyymRTHB452Gcc+SzmRTI12YizrkdwI5Zjz064/Z9wH1+x4pIbctmsnhrvLDLAILtcJHqT9GyuIWWRS2BHE+knErREhFgcnySzIEMV3zyisCPLe/NS3qMnRvj/MnzLFi2IOxyKpK+9xCRsovKTDLkOlwENpPcl9YsslSs6ZDsXLCXBQ2/PczU+JTeG2U2c525zI9CsoiUXTaTDb392zQv6XH6ndNMjE4UfSy1f5NK5iU9JrITnDt+LtDjaqv2cCgkF08hWUTKajw7zsTIRGRmkr2kBy4321WMkeERzp84r9kyqVilClXTF7TqvVFeHZd1YHFTSC6CQrKIlFVUtqSeFlQw0GyZVLqL74UDwYaqdH+a1qWtNCei8Z6vFfGGOO0r2tXhoggKySJSViNDI0AVhuTp9m9abiEVqn1FO7G6WOAzj+n+tN4XISnVxZi1QiFZRMoqajPJzYuaaWxvLHr2LN2fJlYXo2NV8LuViZRDrC5Gx8oOhgaCnXlM9aW01CIkCsnFUUgWkbKKWkg2s1zT/SKDQbo/jZf0iNdry12pXEGHqgvpC1xIXVBIDomX9BgZGrn4764URiFZRMoqaiEZggkGmi2TapBIJgJtAze9wY6WW4RDHS6Ko5AsImUVxZCcSCYYPjTM5PjkvMZPTUyRGcgoCEjF85Ieo2dGuZC6EMjxdEFruBI9CUAheb4UkkWkrLKZLLG6GA0LGsIu5SIv6TE1McXpd07Pa/zwoWEmxyY1kywVb3rmMbBdKPtSxBvidKzUWv0wJFbnQ3LAHUtqhUKyiJTV9G57ZhZ2KRcV+5WkZsukWng9wX49P71WP1anuBGG+uZ62rraAr8Ys1bop1ZEyipKW1JPK3b2bHrdpWaSpdJ1rOoACy4kp/q0C2XY1OFi/hSSRaSsohiSFyxfQH1L/fxnkvtStCxuoWVRS8CViZRXXWMd7SvaAwlVk+OTDB0Y0n8eQzZ9MaYUzldINrM7zKzfzAbM7KE5nr/CzF4ws1Ez+4NZzx0yszfN7DUz2x1U4SJSmaIYks2MRM/8P0i0WYJUk6BmHoffHmZqYkohOWRej8f5k+cZPTMadikV55Ih2cziwMPAVmA98CkzWz/rZRngc8CX3+MwNzvnNjrnNhVTrIhUvmwmS1OiKewy3qWYYKD2b1JNggrJ2oUyGkq13Xgt8DOTvBkYcM4ddM6NAU8C22a+wDl30jm3CxgvQY0iUkVGhkYiN5MMuQ+SoQNDTE1OFTQuO5Tl/MnzCslSNbykRzadZWR4pKjj6ILWaAi6Y0kt8ROSO4EjM+4P5h/zywE/NLOXzWx7IcWJSHWZHJ9k9MxoZEPy5NgkZ4+eLWicNkuQanOxt26RM4+pvhSty1pp6ojeN0e1RL2S589PSJ6rT1MhW/Fc75y7htxyjQfM7MY5/xCz7Wa228x2nzp1qoDDi0ilmJ6ZimpIhsKDgWbLpNoEtUtbuj+t90UENC5spHVZq0LyPPgJyYNA94z7XcAxv3+Ac+5Y/veTwFPklm/M9brHnHObnHOblixZ4vfwIlJBorjb3rT5zrak+9PE6mO51lkiVeDiBhRBbNV+hZYhRYHawM2Pn5C8C1hjZqvMrAG4B+j1c3AzazWzhdO3gduAPfMtVkQqW5RDcltXG/GGeMEfJKm+FF6PR7w+XqLKRMqrobWBhR9cWNQGFBdSF8ims5pJjgivRyF5Puou9QLn3ISZPQg8C8SBx51zb5nZ/fnnHzWz5cBuoA2YMrPPk+uEsRh4Kr+zVh3w9865fy7NX0VEoi7KITkWj5FYnSg4GKj9m1SjYmcep5ch6YLWaEgkE5x94izj2XHqm+vDLqdiXDIkAzjndgA7Zj326IzbPye3DGO2M8CGYgoUkeoR5ZAMhQeDqYkp0vvTrL1rbQmrEim/RDLBwDMD8x6vC1qj5WKHi4NDLL1yacjVVA5fIVlEJAhRD8mJZIL9z+zn4fUP+3r91MQUU+PaLEGqj9fjce74OcbOj9HQ2lDw+FR/inhDnI6VWqsfBTMvxlRI9k8hWUTKJpvJgkFje2PYpcxp470bOf/z8wX1Su6+rps1H1tTwqpEym9mb91lVy8reHy6L423xiMW97Wxr5RYUB1Lao1CsoiUTTaTpamjKbIfnMs3LufXv/vrYZchErqZoWo+ITnVn9KMZYQ0J5ppSjQpJBcomp9UIlKVRjLR3G1PRH5ZMRtQTI5PMnRgSMuQIsZLekV1LKlFCskiUjbZoaxCskgFaGpvomVJy7xC8tDBIaYmpnTRXsSoV3LhFJJFpGyymSzNCYVkkUrgJT2GDhQ+85jqU/u3KPKSHqffOc3k2GTYpVQMhWQRKZtsRjPJIpVivhtQXGz/po1EIsVLergpx/Ch4bBLqRgKySJSNtlMliavKewyRMSHRDLB6SOnmRiZKGhcqi9F67JWmjr0Xo+SYtaZ1yqFZBEpCzflGBnShXsilcJLeuBg6G3tQlkN1AaucArJIlIWo2dGcVNOIVmkQsw3VKX6U1qPHEGtS1tpWNCgkFwAhWQRKYuo77YnIr9sPiH5QuoC2XRW65EjyMzU4aJACskiUhYKySKVpdlrpqmjqaAOF6n+XGcLLbeIpvl2LKlV2nFPRMpCIVmkskzPPL72rdc48OwBX2PGzo0Bav8WVYlkgr6n+5iamCJWp3nSS/EVks3sDuCvgDjwDefc/5j1/BXAN4FrgC86577sd6yI1AaFZJHKc8Mf3cDe7+0taExbdxuJ1YkSVSTF8Ho8psanOH3kNIlVOkeXcsmQbGZx4GHgVmAQ2GVmvc65n814WQb4HPBr8xgrIjVAIVmk8qy7ex3r7l4XdhkSkJnrzBWSL83PXPtmYMA5d9A5NwY8CWyb+QLn3Enn3C5gvNCxIlIbskP5kKwd90REQqE2cIXxs9yiEzgy4/4gcK3P4xczVkSqSDaTpWFBA/GGeNiliIjUpIUfXEhdUx3/8oV/4YUvv+B73GU3Xsa2b9beHKefkGxzPOZ8Ht/3WDPbDmwHWLFihc/Di0ilGMloIxERkTBZzLj1f97K0Z8e9T3mxJsneOM7b3DXN+4iFq+ti/38hORBoHvG/S7gmM/j+x7rnHsMeAxg06ZNfkO4iFSIbCZLU0Lb1IqIhGnzg5vhQf+vf+Ubr/CD3/sBZ46coWNlR+kKiyA//yXYBawxs1Vm1gDcA/T6PH4xY0WkimQzWc0ki4hUmFpex3zJmWTn3ISZPQg8S66N2+POubfM7P7884+a2XJgN9AGTJnZ54H1zrkzc40t1V8mykbPjPLELU9wIX0h7FJkHhYsW8CtX76VFddrKdB8ZTNZFq/TBgMiIpVkZkhe/aurQ66mvHz1SXbO7QB2zHrs0Rm3f05uKYWvsbXo2MvHOLb7GMk7krQsaQm7HCnQ4Z2H+eZ/+CYffuDD3PLfb6FxYWPYJVUczSSLiFSe6Yv9NJMsJZPuTwNw59/eSfuK9pCrkUKNnRvjuS8+x0tfe4l9vfu487E7Sd6eDLusiuGcU0gWEalAFjMSPYma3M66ti5TDFGqL0V9Sz1tXW1hlyLz0LCgga1/tZXP/NtnqG+p5zt3fIfv3/v9ixtkyPsbvzDO5NikQrKISAXykp5mkqV00v1pFq1dhMXm6oonlaL7I9189tXPsvPPd/LvX/p39j61t6DgF4vH2Pq1raz52JoSVhk92m1PRKRyJXoSHPjhAdyUq6kco5BcJqn+FF3XzrlsWypMXVMdv/Lnv8L6/7ieXY/sYnJ00vfYvu/3sefJPQrJIiJSMbykx0R2grPHz9LWWTvfiCskl8F4dpzhQ8Ns+PSGsEuRAC3fuJy7/vaugsb8w93/wOCLgyWqKLoUkkVEKtfMDhe1FJK1JrkMMgMZcLDo8kVhlyIh69rSRWZ/hgup2moFODI0Aigki4hUolrtlayQXAapvhQAi69Qj9ha17Ult+Rm8Ke1NZusmWQRkcrV3t1OrD5Wcx0uFJLLYLr926K1mkmudR/c9EEsbjW35EIhWUSkcsXqYiRWJTSTLMFL96dp626jobUh7FIkZPUt9SzfsJyjLx4Nu5SyymayxBvi1DXrMggRkUqU6FFIlhJI9aVYfLmWWkhO55ZOBn86yNTkVNillM30RiJmtdM6SESkmkz3SnbOhV1K2Whap8Scc6T6U2y4V50tJKf7um52/81uUntTLL1qaaDHPvjcQXp/t5epCf8BPF4fZ9u3trHyoysDrWUm7bYnIlLZvKTH2NkxLpy6QOvS1rDLKQuF5BI7d/wcY2fHNJMsF128eO/FwcBD8iuPvcLomVHW3b3O95j9O/bzoz/4Efe9dF/JZnoVkkVEKtvMDhcKyRKIVH+us4Xav8m0RE+C5kXNHHnhCNfcd01gx50YnWD/M/u58jev5BNf/4Tvca8+/iq9v9vLvh/s4/JPXB5YPTNlM1k6LusoybFFRKT0LobkAxm6P9IdcjXl4WtNspndYWb9ZjZgZg/N8byZ2V/nn3/DzK6Z8dwhM3vTzF4zs91BFl8J1P5NZjMzurZ0BX7x3uEfH2bs7BhXbLuioHEbPr0BL+nx/H99HjdVmrVmmkkWEalsHSs7sJjV1MV7lwzJZhYHHga2AuuBT5nZ+lkv2wqsyf/aDjwy6/mbnXMbnXObii+5sqT709S31NfUDjVyaV3XdXHqZ6cYGR4J7Jh9T/dR31LPqltWFTQuVhfjo3/yUU68cYKffe9ngdUzUzaTpclrKsmxRUSk9OINcdpXtDM0UDu9kv3MJG8GBpxzB51zY8CTwLZZr9kGPOFyXgQ6zOwDAddakdL9aRZdvgiL6ap++YXpdclHXwpmNtk5x77effTc3kN9c33B46+65yqWrF/Cv/7xvwbedWNybJLx8+OaSRYRqXDTHS5qhZ+Q3AkcmXF/MP+Y39c44Idm9rKZbZ9voZVK7d9kLp0f7gQjsE1Fjr9ynDODZ+a9pjgWj3HTn95Eam+KPd/dE0hN07JD2khERKQaJJK11SvZT0ieawp09sLF93vN9c65a8gtyXjAzG6c8w8x225mu81s96lTp3yUFX3j2XGGDw/roj15l8a2RpZetTSwkNz/dD8WM9beuXbex1h39zqWbVjGj//0x0yOTwZSF2i3PRGRauElPbKZ7MXJj2rnJyQPAjMvY+wCjvl9jXNu+veTwFPklm+8i3PuMefcJufcpiVLlvirPuIy+zPgdNGezK1rSxeDLw4GcrFc/9P9dF/fTcvilnkfw2LGzX92M5mBDK8/8XrRNU1TSBYRqQ7THS6GDtTGumQ/IXkXsMbMVplZA3AP0DvrNb3Ap/NdLrYAp51zx82s1cwWAphZK3AbEOx3uRGm9m/yfrq2dDEyNEJ6f7qo4wy9PcSJN05w+bbi27etvXMtnZs72fnfdjIxOlH08UAhWUSkWszslVwLLhmSnXMTwIPAs8Be4B+dc2+Z2f1mdn/+ZTuAg8AA8HXg9/OPLwP+zcxeB14C/q9z7p8D/jtEVro/F34WrVVIlnfruu4Xm4oUo7+3HyCQHsdmudnk0++c5tW/e7Xo44FCsohItUisTgC1E5J9bSbinNtBLgjPfOzRGbcd8MAc4w4CNbsfc6ovRVt3Gw2tDWGXIhG0+PLFNLY3MvjCIBvv3Tjv4+zr3cfidYtZtCaY/4ytvnU1K25YwU/+4ids/J2N8+qWMdPFkJxQSBYRqWT1zfUs7FyokCzFS/entR5Z3pPFjK5ru4qaSc4OZTn040N85A8/Elxd+dnkb9/8bb7a+VXi9fGijjd2bgyLGY1tjQFVKCIiYamlNnAKySXinCPVl2LDb9fsRLr40HVdFzv/bCdj58ZoWFD4Nw77d+zHTbqCd9m7lJU3reS2r9xW9HrpaUuvWqpe4SIiVcBLeuz7p31hl1EWCsklcu74OcbOjalHsryvri1duCnH0V1HWXVzYTvlQa6rxYLlC+jcPLt1efGu+8/XBX5MERGpbF7S4/yJ8/Oe3KkkfrpbyDyk+nKdLbTcQt5P57W5cDufJRcToxMMPDPA2rvWapZWRETK4mKHiwPVv+RCIblE1P5N/GhONLP4isUcfbHw7akPPX+IsXNjgXS1EBER8SPRUzsdLhSSSyTdn6a+tZ62zrawS5GIu7ipiCtsU5G+p/uob6ln1S2FL9MQERGZD6+ndnolKySXSKovxaK1i/Q1uFxS55ZOzp88z/Dbw77HOOfY17uPntt7im7RJiIi4ldjWyOtS1sVkmX+1P5N/Oq+LrejeyHrko+/fJyzx84GssueiIhIIbykx9BA9W9Nre4WJTCeHWf48LDav4kvS65cQsOCBnrv6+WZzz3ja8zEyAQWM9Z+fG2JqxMREfllXtLj7effDruMklNILoHM/gw41P5NfInFY3z8kY8X3OFi2YZltCxuKVFVIiIic0skE7z+xOuMZ8eresmfQnIJTHe20HIL8evq37qaq3/r6rDLEBERuaTpi/eG3x5myfolIVdTOlqTXALTPZK9NV7IlYiIiIgE62Kv5Cq/eE8huQTS/WnaV7TT0FrdO9GIiIhI7VFInsHM7jCzfjMbMLOH5njezOyv88+/YWbX+B1bjVJ9KW0iIiIiIlWp2WumKdGkkGxmceBhYCuwHviUma2f9bKtwJr8r+3AIwWMrSrOOdL9aYVkERERqVpe0mPoQHW3gfMzk7wZGHDOHXTOjQFPAttmvWYb8ITLeRHoMLMP+BxbVc4dP8fYuTFdtCciIiJVy0t6VT+T7Ke7RSdwZMb9QeBaH6/p9Dk2Ep7c9iSHdx4u+jhTE1OA2r+JiIhI9fKSHnu+u4cvJb4UyPFu+MINXP+H1wdyrKD4Cclz7avsfL7Gz9jcAcy2k1uqwYoVK3yUFazVt66mfWV7IMdqbGvkshsvC+RYIiIiIlHzoc98iLHzYxcnB4u17OplgRwnSH5C8iDQPeN+F3DM52safIwFwDn3GPAYwKZNm+YM0qW0+cHN5f4jRURERCpSx8oObv/K7WGXUVJ+1iTvAtaY2SozawDuAXpnvaYX+HS+y8UW4LRz7rjPsSIiIiIikXLJmWTn3ISZPQg8C8SBx51zb5nZ/fnnHwV2AB8DBoALwO+839iS/E1ERERERAJizpV9ZcMlbdq0ye3evTvsMkRERESkipnZy865TXM9px33RERERERmUUgWEREREZlFIVlEREREZBaFZBERERGRWSJ54Z6ZnQKK3/6uuiwGUmEXIToPEaJzEQ06D9Gg8xANOg/RUMh5uMw5t2SuJyIZkuXdzGz3e119KeWj8xAdOhfRoPMQDToP0aDzEA1BnQcttxARERERmUUhWURERERkFoXkyvFY2AUIoPMQJToX0aDzEA06D9Gg8xANgZwHrUkWEREREZlFM8kiIiIiIrMoJEecmf2Gmb1lZlNmtmnWc39kZgNm1m9mt4dVY60xsz8xs6Nm9lr+18fCrqmWmNkd+Z/5ATN7KOx6apWZHTKzN/Pvgd1h11NLzOxxMztpZntmPOaZ2Y/MbH/+90SYNdaC9zgP+nwoMzPrNrPnzWxvPi/9p/zjRb8nFJKjbw9wN7Bz5oNmth64B7gSuAP4GzOLl7+8mvWXzrmN+V87wi6mVuR/xh8GtgLrgU/l3wsSjpvz7wG1vCqvb5H7d3+mh4DnnHNrgOfy96W0vsW7zwPo86HcJoD/4pxbB2wBHsh/LhT9nlBIjjjn3F7nXP8cT20DnnTOjTrn3gYGgM3lrU6k7DYDA865g865MeBJcu8FkZrhnNsJZGY9vA34dv72t4FfK2tRNeg9zoOUmXPuuHPulfzts8BeoJMA3hMKyZWrEzgy4/5g/jEpjwfN7I381236WrN89HMfHQ74oZm9bGbbwy5GWOacOw650AAsDbmeWqbPh5CY2UrgQ8BPCeA9oZAcAWb2/8xszxy/3m+GzOZ4TK1KAnKJc/II0ANsBI4DXwm12Nqin/vouN45dw25pS8PmNmNYRckEgH6fAiJmS0Avgd83jl3Johj1gVxECmOc+5X5zFsEOiecb8LOBZMReL3nJjZ14F/KnE58gv6uY8I59yx/O8nzewpckthdr7/KCmhE2b2AefccTP7AHAy7IJqkXPuxPRtfT6Uj5nVkwvI33HO/Z/8w0W/JzSTXLl6gXvMrNHMVgFrgJdCrqkm5N9s0z5J7uJKKY9dwBozW2VmDeQuXu0NuaaaY2atZrZw+jZwG3ofhK0XuDd/+17g6RBrqVn6fCg/MzPg74C9zrmvzniq6PeENhOJODP7JPA1YAkwDLzmnLs9/9wXgc+Qu7Lz8865Z0IrtIaY2f8m91WaAw4Bn51e9ySll2+p9L+AOPC4c+4vQi6p5pjZauCp/N064O91HsrHzL4L3AQsBk4Afwx8H/hHYAXwDvAbzjldVFZC73EebkKfD2VlZjcAPwHeBKbyD3+B3Lrkot4TCskiIiIiIrNouYWIiIiIyCwKySIiIiIisygki4iIiIjMopAsIiIiIjKLQrKIiIiIyCwKySIiIiIisygki4iIiIjMopAsIiIiIjLL/wfUzNk38POtZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f1 = data[data['Treatment'] == 'Placebo']['Decrease'].to_numpy()\n",
    "f2 = data[data['Treatment'] == 'Calcium']['Decrease'].to_numpy()\n",
    "\n",
    "def f1_hat(x):\n",
    "    return sum(f1 <= x) / len(f1)\n",
    "\n",
    "def f2_hat(x):\n",
    "    return sum(f2 <= x) / len(f2)\n",
    "\n",
    "xx = np.linspace(min(data['Decrease']) - 1, max(data['Decrease']) + 1, 100)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(12, 8))\n",
    "ax1.plot(xx, [f1_hat(x) for x in xx], color='blue', label='F1 empirical distribution')\n",
    "ax1.plot(xx, [f2_hat(x) for x in xx], color='red', label='F2 empirical distribution')\n",
    "ax1.legend(loc='upper left')\n",
    "\n",
    "ax2.plot(xx, [abs(f1_hat(x) - f2_hat(x)) for x in xx], color='purple', label='|F1 - F2|')\n",
    "ax2.legend(loc='upper left')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D: 0.409\n",
      "Test statistic: 0.936\n"
     ]
    }
   ],
   "source": [
    "D = max([abs(f1_hat(x) - f2_hat(x)) for x in data['Decrease']])\n",
    "\n",
    "n1 = len(f1)\n",
    "n2 = len(f2)\n",
    "test_statistic = np.sqrt(n1 * n2 / (n1 + n2)) * D\n",
    "\n",
    "print('D: %.3f' % D)\n",
    "print('Test statistic: %.3f' % test_statistic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "H(t) &= 1 - 2 \\sum_{j=1}^\\infty (-1)^{j-1} e^{-2j^2t^2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "let's create an approximation for $H^{-1}(1 - \\alpha)$ -- so that we may find the p-value for $\\alpha$ such that $ \\sqrt{\\frac{n_1 n_2}{n_1 + n_2}} D > H^{-1}(1 - \\alpha) $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def H(t, xtol=1e-8):\n",
    "    assert t != 0, \"t must be non-zero\"\n",
    "    t2 = t * t\n",
    "    j_max = int(np.ceil(np.sqrt(- np.log(xtol / 2) / (2 * t2) )))\n",
    "    jj = np.arange(1, j_max+1)\n",
    "\n",
    "    return 1 + 2 * sum((-1)**(jj) * np.exp(-2 * (jj**2) * t2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9d328c83KxBCAiZhCVtACIsIlQiWFsVWFKmW2tZWrba2KnVrH++qrb2fVu/7tpu1tta6ILdata1Sta2lPAjihisKKCBLAiEsCVsSELJvM7/nj0RNY0Imycycmcn1fr3yypw5ZzKXw8nlyW/OnJ855xARkegX53UAEREJDhW6iEiMUKGLiMQIFbqISIxQoYuIxAgVuohIjFChi7RiZr80sxuOs/77ZvarcGYSCZQKXXoNM9ttZme1ue9yM3u95XYm8E3gwZblOWZW0ubHLAYuNbOscGQW6QoVusjHLgeWO+dqO9rAOVcHPEdz8YtEFBW6yMfOBVYDmFkKzcU9zMyqWr6GtWz3CvAFbyKKdEyFLvKxKUABgHOumuaC3++c69/ytb9lu23AVI8yinQowesAImH2rJk1tVpOAt5tuZ0OVAbwMyqBtGAHE+kpHaFLb/Ml51z6h1/Ata3WfQCkBvAzUoFjIUkn0gMqdJGPbQLGt1ru6FKkE4GNoY8j0jUqdJGPLQfOaLV8CDjBzNoOr5xB8xumIhFFY+giH3sc2GBmfZ1ztc65fDN7Eigys3hgEnAEmA9M9zKoSHtME1yIfMzMfgGUOufu7mD994ARzrkfhjeZSOdU6CIiMUJj6CIiMUKFLiISI1ToIiIxwrOzXDIyMtzo0aO9enoRkai0fv36cudcZnvrPCv00aNHs27dOq+eXkQkKpnZno7WachFRCRGqNBFRGKECl1EJEao0EVEYkSnhW5mj5hZqZlt7mC9mdk9ZlZoZpvM7JTgxxQRkc4EcoT+KDDvOOvPBca1fC0EHuh5LBER6apOC9059yrNV5jryALgcddsDZBuZkODFVBERAITjPPQs4HiVsslLfcdCMLPFhGJeM456hr91DQ0UdPgo7qhiep6H7Utt2vaLE8fNZDZ49r9bFCPBKPQrZ372r2Eo5ktpHlYhpEjRwbhqUVEeq7J5+dwdQOlFfWUVdVRWlFPeVU9lfVNzSVc72su5QYftS3l/PFyc0l35cK118wZG7GFXgKMaLU8HNjf3obOucXAYoC8vDxdt1dEQqqmoYnSinpKK+spq6yntLKu5fvH95VV1nG4uqHdQk5KiCMlKZ5+SQn0S4qnX3ICKUnxDEtPol9SPCnJzetSkuLpm5TQZjmelOTmx6W0enzfxHji49o7Du65YBT6UuB6M1sCzASOOec03CIiIXe0poFtByrJP1jB3iM1zSVdUU9ZVT2lFXVUN/g+8ZiEOCMzNZnM1GSy0/swbUQamal9yGq5Lys1mawBfcjon0RyQrwH/1Xd12mht0zBNQfIMLMS4DYgEcA5t4jmeRjnA4VADfDtUIUVkd6p0ednV3k12w5UkH+wkvwDFWw7UMnBirqPtklJiidrQB8yU5OZPGwAc3IzyUrt06qkk8nsn8zAfknEhegI2WudFrpz7uJO1jvguqAlEpFerayynvyDFeQfqGRby/fC0ioafH4AEuONsZn9+fTYE5gwJJUJQwcwcUgqmanJmMVmUQdKk0SLiCcamvzsKK0kv2XIJP9gJdsOVFBe1fDRNoMHJDNhyABmj89g4pABTBiaypiM/iQl6EPu7VGhi0jYHK1p4OWCUl7YWsrq7WVU1TcBkJwQx/jBqZyZm/XREXfukFRO6J/sceLookIXkZDac7iaVVsP8cK2Q6zd/QE+vyMzNZnzpw5l1tgMJg4dwOgT+pEQr6PunlKhi0hQ+f2ODSVHeWHrIVZtPcSO0ioAcgencs0ZYzlr0mBOzk6L2TcmvaRCF5Eeq23w8UZhOS9sO8QL20opr6onPs6YmTOIi2eM5KyJgxl5Qj+vY8Y8FbqIdEtZZT0v5R9i1dZSXi8so67RT2pyAmfkZjJ30mDmjM8irV+i1zF7FRW6iASs0efnnxv288Tbe3iv+CjOQXZ6X76eN4KzJg1mZs4JOgPFQyp0EelUXaOPp9YV8+DqIvYdrWX84P7c8PnxzJ00mIlDU3v9+d+RQoUuIh2qrGvkz2v28vDrRZRXNTB91EB+9qWTmJObqRKPQCp0EfmEI9UN/PGNXTz25m4q6po4fXwm180Zy4ycQSryCKZCF5GPHDxWx/++VsQTb++lttHHvMlDuPbMsZw8PN3raBIAFbqIsLu8mgdf3ckz60vwO1gwbRjXnDGWcYNTvY4mXaBCF+nF8g9WcP/LO1m2aT8J8XF8/dQRfPf0sYwYpHPGo5EKXaQXenfvB9z/ciEvbCslJSmeq2aP4YrP5pA1oI/X0aQHVOgivcjBY3X8+O+beLmgjPR+ifzHWeP51qxRpPdL8jqaBIEKXaSX+OeGffz02c00+hy3nDuBy04bRUqyKiCW6F9TJMYdrWngJ89uZtmmA0wbkc7vvj6NnIwUr2NJCKjQRWLY6u1l/PCZjRyuauDGueO5Zs5YXaY2hqnQRWJQbYOPXz63jcff2sOJWf156JunMmV4mtexJMRU6CIxZkPxUX7w1w0UlVfznc/k8MN5ufRJjK7Z66V7VOgiMaLR5+cPLxVy38uFDE5N5okrZzLrxAyvY0kYqdBFYkBhaRU/eGoDm0qOccGnsvmvL04mra+uRd7bqNBFopjf73jsrd386rl8+ibFc/83TmH+lKFexxKPqNBFotSBY7Xc/PQmXi8sZ05uJr/+ysn6pGcvp0IXiUJLN+7nJ/94n0af4+cXnMQlM0bqsraiQheJNve/UsivVxRwysh0fvu1aYzWh4SkhQpdJEo457hjRQGLVu/ki1OHcdfXppKoDwlJKyp0kSjg9ztuXbqZP6/ZyyUzR3L7gpOIj9MQi/w7FbpIhGv0+bn56Y08u2E/3z1jDLfMm6DxcmmXCl0kgtU1+vjek++xaushbj4nl2vnjFWZS4dU6CIRqrq+iYV/WscbhYf57y9O5luzRnsdSSJcQO+omNk8Mysws0Izu6Wd9Wlm9i8z22hmW8zs28GPKtJ7HKtp5NKH3+atnYe568KpKnMJSKdH6GYWD9wHzAVKgLVmttQ5t7XVZtcBW51z55tZJlBgZn9xzjWEJLVIDCurrOeyh9+mqKya+78xnXknDfE6kkSJQI7QZwCFzrmiloJeAixos40DUq15cK8/cARoCmpSkV5g39FavvbgW+w5XMPDl+epzKVLAhlDzwaKWy2XADPbbHMvsBTYD6QCX3fO+dv+IDNbCCwEGDlyZHfyisSsorIqLn3obSrrm/jTFTPIGz3I60gSZQI5Qm/vLXXXZvkcYAMwDJgG3GtmAz7xIOcWO+fynHN5mZmZXQ4rEqu27q/gaw++RX2TnyULT1OZS7cEUuglwIhWy8NpPhJv7dvA312zQmAXMCE4EUVi2/o9H3DR4rdIjI/jqas/zeRhmllIuieQQl8LjDOzHDNLAi6ieXiltb3A5wHMbDCQCxQFM6hILHp9RzmXPfw2g1KSePrqTzM2s7/XkSSKdTqG7pxrMrPrgZVAPPCIc26LmV3dsn4RcDvwqJm9T/MQzY+cc+UhzC0S9dbvOcJ3HlvLmIwUHr9iBlmpuvSt9ExAHyxyzi0Hlre5b1Gr2/uBs4MbTSR27Ttay3f/tJ6haX148qrTGJiS5HUkiQH6pKhImNU0NHHlY+uob/SzZGGeylyCRoUuEkZ+v+PGpzZScLCChy8/lROzUr2OJDFEF1MWCaO7X9zBc5sP8p/zJ3JmbpbXcSTGqNBFwuRfG/dzz4s7uHD6cK74bI7XcSQGqdBFwuD9kmPc9PRG8kYN5GcXnKRL4EpIqNBFQqy0oo6rHl9HRv9kFl02neSEeK8jSYzSm6IiIVTX6OOqP62noq6RZ66eRUb/ZK8jSQxToYuEiHOOW/62iY3FR1l06XQmDfvE5Y1EgkpDLiIh8sDqnTy7YT83nT1el8GVsFChi4TAqq2HuHNlAedPHcZ1Z57odRzpJVToIkGWf7CCG5a8x5TsNO786sk6o0XCRoUuEkSHq+q58rF1pCQnsPiyPPok6owWCR+9KSoSJA1Nfq75y7uUVdbz1+9+miFpunqihJcKXSQInHPctnQz7+w6wu8vmsa0EeleR5JeSEMuIkHwzw37efKdYq47cywLpmV7HUd6KRW6SA8dPFbHrf/czPRRA/nB3Fyv40gvpkIX6QHnHD/62yYafY67LpxKfJzOaBHvqNBFemDJ2mJWby/jlnMnMDojxes40sup0EW6qfhIDT9btpVZY0/gstNGeR1HRIUu0h1+v+PmZzZiZvz6qycTp6EWiQAqdJFueOyt3awpOsKt501i+MB+XscRAVToIl22s6yKXz2Xz+cmZHFh3nCv44h8RIUu0gVNPj83Pb2RPonx/OrLU3SdFoko+qSoSBcsfq2I9/Ye5Z6LP0XWAH20XyKLjtBFApR/sILfrdrOF6YM5fyTh3odR+QTVOgiAWho8vODv24krW8it39JkzxLZNKQi0gA7n1pB1sPVPDgZdMZlJLkdRyRdukIXaQTG4uPct8rO/nyKdmcM1lTyUnkUqGLHEddo48bn95IZv9kbjt/stdxRI5LQy4ix/HbVdspLK3i8e/MIK1votdxRI5LR+giHVi7+wj/+1oR35g5ktPHZ3odR6RTARW6mc0zswIzKzSzWzrYZo6ZbTCzLWa2OrgxRcKrpqGJm57eyPCBffnP+RO9jiMSkE6HXMwsHrgPmAuUAGvNbKlzbmurbdKB+4F5zrm9ZpYVqsAi4XDvS4XsOVzDkoWnkZKskUmJDoEcoc8ACp1zRc65BmAJsKDNNpcAf3fO7QVwzpUGN6ZI+Owur+ah13bx5VOyOW3MCV7HEQlYIIWeDRS3Wi5pua+18cBAM3vFzNab2Tfb+0FmttDM1pnZurKysu4lFgmx25dtJTHeuGXeBK+jiHRJIIXe3kfiXJvlBGA68AXgHOCnZjb+Ew9ybrFzLs85l5eZqTeZJPK8nF/Ki/mlfP/z43StFok6gQwOlgAjWi0PB/a3s025c64aqDazV4GpwPagpBQJg4YmP/+zbCtjMlL49mdyvI4j0mWBHKGvBcaZWY6ZJQEXAUvbbPNPYLaZJZhZP2AmsC24UUVC65E3drGrvJpbz59EUoLO6JXo0+kRunOuycyuB1YC8cAjzrktZnZ1y/pFzrltZrYC2AT4gYecc5tDGVwkmA5V1PGHF3dw1sQs5uTqJC2JTgGdj+WcWw4sb3PfojbLdwJ3Bi+aSPjc8Vw+jT7HT74wyesoIt2mvyul11u/5wh/f28fV52ew+iMFK/jiHSbCl16NZ/fcdvSLQwZ0Idr55zodRyRHlGhS6/21LpiNu+r4MfzJ+gToRL1VOjSax2raeTOlQXMGD2IL04d5nUckR5ToUuv9bsXtnO0poHbvjhJU8pJTFChS6+Uf7CCP63ZwyUzRzJ5WJrXcUSCQoUuvY5zjv9eupXUPgncODfX6zgiQaNCl17nuc0HeavoMDeenctATfgsMUSFLr1KbYOPn/+/bUwcOoBLZoz0Oo5IUKnQpVd5YPVO9h2t5b/On0R8nN4IldiiQpdeo/hIDYtW7+SLU4cxUxNXSAxSoUuv8euVBcSb8eP5mrhCYpMKXXqFzfuO8a+N+7niszkMTevrdRyRkFChS69wx4p80vslsvCMMV5HEQkZFbrEvDcLy3ltRznXzTmRAX0SvY4jEjIqdIlpzjnuWJHPsLQ+XPbpUV7HEQkpFbrEtBWbD7Kx5Bg3zB1Pn8R4r+OIhJQKXWJWk8/Pnc8XMC6rP185ZbjXcURCToUuMevp9SUUlVVz0zm5+hCR9AoqdIlJtQ0+7n5hO6eMTOfsSYO9jiMSFip0iUmPvbWbQxX1/GjeBF3rXHoNFbrEnGM1jdz/ciFn5mbqI/7Sq6jQJeY8sHonlfVN3HyOPuIvvYsKXWLKwWN1/PGNXSyYOoxJwwZ4HUckrFToElN+/+IO/M5x49maiUh6HxW6xIydZVU8ta6Yb8wcxYhB/byOIxJ2KnSJGXc9X0ByQhzXnXmi11FEPKFCl5iwsfgoy98/yJWzx5CZmux1HBFPqNAlJvx6ZT6DUpK4anaO11FEPKNCl6j32o4y3ig8zHVnnkiqLo8rvVhAhW5m88yswMwKzeyW42x3qpn5zOyrwYso0jHnHL9eUUB2el8uPW2k13FEPNVpoZtZPHAfcC4wCbjYzCZ1sN0dwMpghxTpyMotB3l/3zFuOGscyQm6PK70boEcoc8ACp1zRc65BmAJsKCd7b4H/A0oDWI+kQ75/Y7frdrBmIwULvhUttdxRDwXSKFnA8Wtlkta7vuImWUDFwCLghdN5PiWvX+AgkOV3DB3PAnxejtIJJDfgvYuVefaLN8N/Mg55zvuDzJbaGbrzGxdWVlZoBlFPqHJ5+fuVdvJHZzKeVOGeh1HJCIkBLBNCTCi1fJwYH+bbfKAJS2XKc0A5ptZk3Pu2dYbOecWA4sB8vLy2v5PQSRgz27YT1F5NYsunU6cJq8QAQIr9LXAODPLAfYBFwGXtN7AOffRyb9m9iiwrG2ZiwRLo8/P71/czknZAzhnsiavEPlQp0Muzrkm4Hqaz17ZBjzlnNtiZleb2dWhDijS1tPrSig+UsuNc3M1eYVIK4EcoeOcWw4sb3Nfu2+AOucu73kskfbVNfr4w0s7+NTIdObkZnodRySi6NQAiSpL3tnLgWN13HS2js5F2lKhS9SobfBx78s7mZkziFljNbWcSFsqdIkaf1qzm/Kqem7U0blIu1ToEhWq6pt44JWdzB6XwYycQV7HEYlIKnSJCo++sYsPaho1tZzIcajQJeIdq21k8atFnDUxi2kj0r2OIxKxVOgS8R5+rYiKuib+Y+54r6OIRDQVukS0I9UNPPLGbuZPGcLkYWlexxGJaCp0iWgPvrqT6oYmbjhLR+cinVGhS8QqrazjsTd3s2DqMMYPTvU6jkjEU6FLxHrglZ00+hz/R0fnIgFRoUtE2ne0lr+s2ctXTskmJyPF6zgiUUGFLhHpnhd2AOjoXKQLVOgScXaWVfHMuyVcetoostP7eh1HJGqo0CXi/G7VdpIT4rj2zLFeRxGJKip0iShb9h9j2aYDXPHZHDL6J3sdRySqqNAlovxmZQFpfRO5cvYYr6OIRB0VukSMtbuP8HJBGVefMZa0volexxGJOip0iQjOOe5cUUBmajLfmjXK6zgiUUmFLhHh1R3lvLP7CN/73In0SwpoqlsRaUOFLp5zznHnynyGD+zLRaeO9DqOSNRSoYvnVmw+yOZ9Fdxw1niSErRLinSXfnvEUz6/4zfPF3BiVn8u+FS213FEopoKXTz1j/f2sbOsmpvOHk98nCZ+FukJFbp4pr7Jx+9WbWdKdhrnTB7idRyRqKdCF8/8dW0x+47WcvM5uZjp6Fykp1To4omahibuebGQmTmDmD0uw+s4IjFBhS6e+OMbuymvqtfRuUgQqdAl7Eor67j/5ULmThpM3uhBXscRiRkqdAm7u1Zup8Hn5//On+h1FJGYokKXsNq87xhPrS/m8lmjGa2p5USCKqBCN7N5ZlZgZoVmdks7679hZptavt40s6nBjyrRzjnH7cu2MrBfEtd/bpzXcURiTqeFbmbxwH3AucAk4GIzm9Rms13AGc65k4HbgcXBDirRb+WWg7y96wg/mDtel8cVCYFAjtBnAIXOuSLnXAOwBFjQegPn3JvOuQ9aFtcAw4MbU6JdfZOPny/fRu7gVC46dYTXcURiUiCFng0Ut1ouabmvI1cAz7W3wswWmtk6M1tXVlYWeEqJen98YzfFR2r5yXkTSYjXWzcioRDIb1Z7Jwm7djc0O5PmQv9Re+udc4udc3nOubzMzMzAU0pUK6us596XCjlrYhazx+nfXSRUAplJoARo/TfycGB/243M7GTgIeBc59zh4MSTWPDbVQXUNfr4T52mKBJSgRyhrwXGmVmOmSUBFwFLW29gZiOBvwOXOee2Bz+mRKst+4+xZG0x35o1mjGZ/b2OIxLTOj1Cd841mdn1wEogHnjEObfFzK5uWb8IuBU4Abi/5WPcTc65vNDFlmjw4WmK6X0T+b5OUxQJuYAmb3TOLQeWt7lvUavbVwJXBjeaRLvntx5iTdERbl8wmbR+Ok1RJNR0uoGERH2Tj18s38a4rP5cPEPzhIqEgwpdQuKxN3ez53ANPz1vkk5TFAkT/aZJ0BUfqeHuF3bwuQlZnD5epymKhIsKXYLK73fc/MxG4sy4/UsneR1HpFdRoUtQ/fntPawpOsJPvjCR7PS+XscR6VVU6BI0ew/X8Mvl+Zw+PpOv63otImGnQpeg8PsdNz2zkYQ4446vTNG0ciIeUKFLUDz+1m7e2XWEn543iaFpGmoR8YIKXXpsz+Fq7lhRwJzcTC7M05WTRbyiQpceaT6rZRMJ8cYvv6yhFhEvqdClRx57q3mo5VYNtYh4ToUu3barvJo7VuTzuQlZfHW6hlpEvKZCl27x+x0/fGYjifFx/OICDbWIRIKArrYo0tYDq3eydvcH3HXhVIak9fE6joigI3TphhWbD/Cb5ws47+ShfPmU400vKyLhpEKXLtlYfJQb/rqBaSPS+c2FUzXUIhJBVOgSsH1Ha7ny8XVk9E9m8WV59EmM9zqSiLSiMXQJSGVdI1c8upa6Bh9PXDmTzNRkryOJSBsqdOlUk8/P9558jx2lVTz67VMZNzjV60gi0g4NuUinbl+2lVcKyrh9wUnMHqcJK0QilQpdjuvRN3bx2Ft7uGp2DpfM1NygIpFMhS4dWrnlIP+zbCtzJw3mlnMneh1HRDqhMXT5BOccj765m9uXbWVKdhq/v2ga8XE6PVEk0qnQ5d80+vzctnQLT7y9l7mTBnP316fRL0m7iUg00G+qfORoTQPX/uVd3tx5mGvmjOXms3OJ05G5SNRQoQsAO8uquPKxdez7oJa7LpzKV3T1RJGoo0IXXt9RzrV/WU9ifBxPXDWTvNGDvI4kIt2gQu/FSivr+MOLhTzxzl5OzOzPQ9/KY8Sgfl7HEpFuUqH3QhV1jSxeXcTDr++i0efn4hkj+NG8CaT2SfQ6moj0gAq9F6lr9PHnNXu47+VCPqhp5LyTh3LT2bmMzkjxOpqIBIEKPcb5/Y5N+47x4rZD/G19CfuP1TF7XAY/PGcCU4aneR1PRIIooEI3s3nA74F44CHn3K/arLeW9fOBGuBy59y7Qc4qAapt8PF6YTkvbjvEi/mllFXWE2cwI2cQv7lwKrNOzPA6ooiEQKeFbmbxwH3AXKAEWGtmS51zW1ttdi4wruVrJvBAy3cJskafn8q6JipqG6msa+JITQO7y6vZ1fK1+3A1JR/U4vM7UpMTOD03k7MmZjFnfBYDU5K8ji8iIRTIEfoMoNA5VwRgZkuABUDrQl8APO6cc8AaM0s3s6HOuQPBDrx6exk/W7a18w0jkGt927lPrnPN3/3O4Rz4/I5Gn5+mlu+NPj91jf52f3ZKUjw5mSlMyU5jwdRhzBxzAqeOHkRSgi7XI9JbBFLo2UBxq+USPnn03d422cC/FbqZLQQWAowc2b0r9/VPTmDc4P7demwkMKz1Qpt1EGeGWfPt+Lg4EuONhHgjMT6OpPg4+icnkNongQF9E0ntk0h6v0RGndCPzP7Jmg5OpJcLpNDbawnXjW1wzi0GFgPk5eV9Yn0gpo8ayPRR07vzUBGRmBbI3+MlwIhWy8OB/d3YRkREQiiQQl8LjDOzHDNLAi4ClrbZZinwTWt2GnAsFOPnIiLSsU6HXJxzTWZ2PbCS5tMWH3HObTGzq1vWLwKW03zKYiHNpy1+O3SRRUSkPQGdh+6cW05zabe+b1Gr2w64LrjRRESkK3ROm4hIjFChi4jECBW6iEiMUKGLiMQIa/sR9LA9sVkZsKebD88AyoMYJ5SUNTSiJWu05ARlDZVgZx3lnMtsb4Vnhd4TZrbOOZfndY5AKGtoREvWaMkJyhoq4cyqIRcRkRihQhcRiRHRWuiLvQ7QBcoaGtGSNVpygrKGStiyRuUYuoiIfFK0HqGLiEgbKnQRkRgRcYVuZvPMrMDMCs3slnbWm5nd07J+k5mdEuhjw5zzGy35NpnZm2Y2tdW63Wb2vpltMLN1ocwZYNY5ZnasJc8GM7s10Md6kPXmVjk3m5nPzAa1rAvb62pmj5hZqZlt7mB9ROynAWaNpH21s6yRtK92ljX8+6pzLmK+aL48705gDJAEbAQmtdlmPvAczbMknQa8Hehjw5xzFjCw5fa5H+ZsWd4NZETQazoHWNadx4Y7a5vtzwde8uh1PR04BdjcwXrP99MuZI2IfTXArBGxrwaStc22YdlXI+0I/aMJqZ1zDcCHE1K39tGE1M65NUC6mQ0N8LFhy+mce9M590HL4hqaZ3HyQk9el3C+pt15vouBJ0OYp0POuVeBI8fZJBL204CyRtC+Gsjr2pGIe13bCMu+GmmF3tFk04FsE8hjg6Wrz3UFzUdrH3LA82a23ponzg6lQLN+2sw2mtlzZja5i48NloCfz8z6AfOAv7W6O5yva2ciYT/tDi/31UBFwr4asHDuqwFNcBFGPZmQOqCJqoMk4OcyszNp/iX5bKu7P+Oc229mWcAqM8tv+b99KASS9V2arw9RZWbzgWeBcQE+Npi68nznA28451ofIYXzde1MJOynXRIB+2ogImVf7Yqw7auRdoTekwmpwzlRdUDPZWYnAw8BC5xzhz+83zm3v+V7KfAPmv9cDJVOszrnKpxzVS23lwOJZpYRyGPDnbWVi2jzJ2yYX9fORMJ+GrAI2Vc7FUH7aleEb18N5ZsGXf2i+S+GIiCHj9/YmNxmmy/w7282vRPoY8OccyTNc6zOanN/CpDa6vabwDyPX9MhfPwhsxnA3pbXN2yvaVf+DYE0mscuU7x6XVueZzQdv3nn+X7ahawRsa8GmDUi9tVAsnqxr0bUkIvrwYTUHT3Ww5y3AicA95sZQJNrvuLaYOAfLZMvm2kAAACFSURBVPclAE8451aEImcXsn4VuMbMmoBa4CLXvLeF7TXtQlaAC4DnnXPVrR4e1tfVzJ6k+YyLDDMrAW4DElvl9Hw/7ULWiNhXA8waEftqgFkhzPuqPvovIhIjIm0MXUREukmFLiISI1ToIiIxQoUuIhIjVOgiIjFChS4iEiNU6CIiMeL/A6/+AJEelm2ZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's plot some values for H(t)\n",
    "\n",
    "tt = np.logspace(-3, 0.25, 100)\n",
    "ht = [H(t) for t in tt]\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.plot(tt, ht)\n",
    "plt.title('H(t)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "def H_inv(q):\n",
    "    def loss_function(x):\n",
    "        return (H(x) - q)**2\n",
    "    \n",
    "    x0 = 1.0\n",
    "    return minimize(loss_function, x0, method='nelder-mead').x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAfW0lEQVR4nO3de3RV9Z338fcXknDJlVwgEAgJN7kLiiBira2XVmpr23FaL7WOdR6mTjvLzmOnnbZrehm7uto+s3q1rWOtj4+t0860OmrrbWyrVapgAbkHJICQkIQQcoUAycn5Pn+ckxhjICdw7vm81jqLnLN39v5uEz78/O7f3tvcHRERSX2jEl2AiIhEhwJdRCRNKNBFRNKEAl1EJE0o0EVE0oQCXUQkTSjQJS2Y2e1m9vdDrLPAzF4zs3YzW3kO+/qtmR0xs/9zmuUlZrbbzMZGsK1JZlZlZmPOth6RXgp0SVpm9oaZXTngs78xs7UDPvsn4GvAXWZ29xk2+QlgH1Dg7q+cYb8fMbOXzazTzF4YuNzd3w+8C/ismRUMsol/Bv6vu588Qy292zoMPA+sGWpdkaEo0CWlmdmtwB3AZcA7gOvN7FOnWb0QqHL34BCbbQa+B3zzdCu4+/bwl0UD6hkD3Ar8Yujq+zwM/N0w1hcZVEaiCxA5W2b2PuDzwDvdvSb82TuBZ82s0d1/PeBbMoChwhx3/314W3871Kq8/e/QCqDV3Wv71VkJPAhcAKwDdhP6v4SPhVdZD8wws+nufmCo+kRORyN0SVnu/qS7z+8N8/Bnje6+dGCYm1khsAw4GMUSaoArzcz6fbaIUGD39x/ARqAYuJvQCL6PuweAauD8KNYmI5ACXZLdY2bW2vsCfjzcDZjZPwBHgTbg/0Wxts8A3wVa+n1WAHT023c5cBHwL+5+yt1fBH47yLY6wt8rctYU6JLsPujuBb0v4IwzWQbj7j8EJgOlwHX9l5nZvWZ2LPz64jA3/TXgq4R6871agNx+76cALe5+vN9ng7VVcoHWYe5f5C0U6DIiuHsD8Aowf8Dnn3T3nPDrG8Pc7Dzg8QEnWbcCc/q9rwcmmFl2v8/K+2/EzDKAWcCWYe5f5C0U6DKSnAKyhlrJzEaH55BnAKPMbKyZZQ6yakZ4m/29ChSYWRlA+CTnBuBrZpZlZpcC7x/wPcuBN3RCVM6VAl1GkiCR/c7fApwAfkJoKuQJ4Kf9VzCz3u28ZdaMu3cRmtHysX4f30Ro9ksz8BXgoQH7uxm4N5IDEDkT0wMuZKQws28AS4EPuHv3OW5rOaEWToG7dwxYVgK8BCx19xODfO9XgVnu/jEzmwj8KbzukBciiZyJRugyktwPjAPqzOzis92ImT0G/IbQzJWOgcvd/Yi7zx0szAdZt9Hd5ynMJRp0YZGMGO6+D7g8Ctv54LlXIxJ9armIiKQJtVxERNJEwlouxcXFXlFRkajdi4ikpI0bNza5e8lgyxIW6BUVFWzYsCFRuxcRSUlmdtrrFdRyERFJEwp0EZE0MWSgm9k0M3s+/JisHWZ25yDrXG5mbWa2Ofz6cmzKFRGR04mkhx4A7nL3TWaWC2w0s+fcfeeA9V5y92ujX6KIiERiyBG6u9e7+6bw1x1AFVAW68JERGR4htVDN7MKQvfCWD/I4pVmtsXMnjazBVGoTUREhiHiaYtmlgM8AnzG3dsHLN4ETHf3Y2a2GngMmD3INtYQfrp5eXn5wMUiInIOIhqhh+8F/QjwsLs/OnC5u7e7+7Hw108BmWZWPMh697n7MndfVlIy6Lx4EZG09v3f7+GlPUdisu1IZrkY8DOgyt2/c5p1SnsflBu+regoQs9wFBGRfu55fg+v7I1NPEbScllF6Ib/28xsc/izLxJ+jJa73wtcD9xhZgFCDwO4wXXXLxGRt3B3unuczNGxuQRoyEB397WADbHOPcA90SpKRCQdBYKhcW5WRmwCXVeKiojESXdP6ImFGaPOOEY+awp0EZE46Q6ERuixarko0EVE4qQrPELPVMtFRCS1BYKhQM8arZaLiEhKU8tFRCRN9LZcMhToIiKprXeWi1ouIiIprjfQ1XIREUlx3T3qoYuIpAWN0EVE0sSbga4euohIStMIXUQkTXRpHrqISHrou1I0Qy0XEZGUppaLiEia6L30X1eKioikuC7NchERSQ9vXvqvEbqISEpTD11EJE30XvqfoZaLiEhq6xuhj9IIXUQkpXX3BMkYZYzSQ6JFRFJbd4/HrH8OCnQRkbjpCgRjNmURFOgiInETCAY1QhcRSQfdAbVcRETSQndPkMwY3ZgLFOgiInHT1aOWi4hIWujuCcbssn9QoIuIxE2gx2N2lSgo0EVE4kYtFxGRNNGtQBcRSQ/dPa4euohIOgiN0NVDFxFJed09HrPHz4ECXUQkbjRtUUQkTSS85WJm08zseTOrMrMdZnbnIOuYmf3AzKrNbKuZXRCbckVEUld3ILazXDIiWCcA3OXum8wsF9hoZs+5+85+61wDzA6/VgA/Cf8pIiJhXYnuobt7vbtvCn/dAVQBZQNWuw54yEPWAQVmNjnq1YqIpLBAMEhWssxyMbMKYCmwfsCiMqCm3/ta3h76mNkaM9tgZhuOHDkyvEpFRFJcrFsuEW/ZzHKAR4DPuHv7wMWDfIu/7QP3+9x9mbsvKykpGV6lIiIprrvHycxIcKCbWSahMH/Y3R8dZJVaYFq/91OBunMvT0QkPbh74u/lYmYG/AyocvfvnGa1J4CPh2e7XAy0uXt9FOsUEUlpgWCoaZE5KnY99EhmuawCbgG2mdnm8GdfBMoB3P1e4ClgNVANdAK3Rb9UEZHU1d0TBIhpy2XIQHf3tQzeI++/jgOfilZRIiLpprsnPEJPhpOiIiJy9npH6EkzbVFERM5OX8tFI3QRkdTWHQi1XHS3RRGRFNfVN0JXy0VEJKUFgr09dI3QRURSWm/LRT10EZEU1xWHeegKdBGROOib5RLDK0UV6CIicRCPK0UV6CIicRDQlaIiIulB0xZFRNKErhQVEUkTCnQRkTTx5jx0tVxERFJat64UFRFJD90BtVxERNJC7wMuMtRyERFJbV06KSoikh5OqeUiIpIeapo7mZI/ltG6l4uISGrb09jBzIk5Md2HAl1EJMaCQWdv43FmT8yN6X4U6CIiMXao9QQnunuYPUkjdBGRlFbdeAyAWWq5iIiktr5AL1Ggi4iktD2NHRTnZDEhOyum+1Ggi4jE2J7GYzFvt4ACXUQkptyd6sZjMZ/hAgp0EZGY2nvkOB0nA8yfkhfzfSnQRURi6M/VTQCsmlkc830p0EVEYmhtdRPTCsdRXjQ+5vtSoIuIxEigJ8i6fUe5dFbsR+egQBcRiZlth9roOBngkji0W0CBLiISM2v3hPrnl8wsisv+FOgiIjHy5LZ6lpYXUJQzJi77U6CLiMRAVX07uxo6+PDSsrjtc8hAN7MHzKzRzLafZvnlZtZmZpvDry9Hv0wRkdTy2GuHyBhlvG/xlLjtMyOCdR4E7gEeOsM6L7n7tVGpSEQkxfUEncc313H5eSUUxvj+Lf0NOUJ39xeB5jjUIiKSFp7beZiG9pNcf+HUuO43Wj30lWa2xcyeNrMFUdqmiEhK+ulL+5hWOI6r5pfGdb/RCPRNwHR3Px/4IfDY6VY0szVmtsHMNhw5ciQKuxYRSS4bD7Sw8UALt6+qjOkDoQdzzoHu7u3ufiz89VNAppkNOove3e9z92XuvqykpORcdy0iknR+/Hw1eWMz+Otl0+K+73MOdDMrNTMLf708vM2j57pdEZFU8+r+Zv6wq5FPXj6T7DGRzDmJriH3aGa/BC4His2sFvgKkAng7vcC1wN3mFkAOAHc4O4es4pFRJKQu/PNp6uYlDeG2y6pTEgNQwa6u984xPJ7CE1rFBEZsR7fXMemg61888OLGJc1OiE16EpREZFz1NrZxd2/28mSaQUJ6Z33in+TR0QkzXz9ySpaT3Tz8w8tivvMlv40QhcROQdPbavnNxtrueOdM+PymLkzUaCLiJylQ60n+MKj2zh/aj53Xjk70eUo0EVEzsbJ7h7u+MVGgkHnezcsJXN04uNUPXQRkWFyd77039vZWtvGfbdcSGVxdqJLAjRCFxEZtm8/u5tHNtXymStnc/WC+N6v5UwU6CIiw/Cztfv5yQt7uWlFOXdekfi+eX8KdBGRCP33a7Xc/budXLOwlLuvW0j4ridJQ4EuIhKB//zLQe76ry2snFHEdz+6JKHzzU9HJ0VFRIZw/0v7+PqTVbxzTgn3fuxCxmYm5tL+oSjQRUROw935znOv88M/VrN6USnf++hSsjKSt7GhQBcRGURnV4B/+vVWntxWz0eWTeUbH1pERhLMNT8TBbqIyACHWk+w5qEN7Kxv5wvXzGXNZTOS7gToYBToIiL9vLC7kbv+awtdgSAP3HoR75o7MdElRUyBLiICdPcE+bf/2c2//2kf503K5Uc3L2XWxNxElzUsCnQRGfGqG49x16+3sKWmlZtXlPMv185P2pksZ6JAF5ERqyfo/PSlfXznudcZnzWaH910Ae9bPDnRZZ01BbqIjEjbD7Xxpce2s6WmlavnT+LrH1rIxNyxiS7rnCjQRWREae3s4t/+ZzcPrz9IUXYW379hCR84f0pKzGIZigJdREaEYND5zw01fPuZXbSd6ObWlRX841VzyB+XmejSokaBLiJpzd15cU8T335mFzvq2rmoYgL/et1C5k1O7OPiYkGBLiJp67WDLXzrmV2s29fM1Anj+N5Hl3DdkvRorwxGgS4iaWdnXTs/+MMentnRQFF2Fl99/3xuXFHOmIzUm4o4HAp0EUkbrx1s4UfPV/P7qkZyxmTwj1fO4fZ3VJIzZmRE3cg4ShFJW+7OK3uP8uMX9rK2uomC8Zn876vmcOvKCvLHp88Jz0go0EUkJZ3s7uGJLXU8sHY/uxo6KMkdw5dWz+OmFeVkj5AR+UAj86hFJGU1dpzkF+sO8vC6Axw93sXc0ly+/VeL+cCSKSl5uX40KdBFJOm5OxsOtPAf6w/yu611BILOFXMn8olVlaycWZS2s1aGS4EuIkmr5XgXj752iF++epDqxmPkjMngxuXl3Laqksri7ESXl3QU6CKSVIJBZ/3+Zn71l4M8vb2BrkCQJdMK+PZfLeZ9iyeP2P54JPRfRkSSQlV9O49tPsRvN9dR13aS3LEZ3HjRNG5YXp6WV3XGggJdRBLmUOsJHt98iMdfq2P34Q4yRhmXzSnh89fM5er5pYzLGtknOYdLgS4icdXa2cWT2+p5/LU6Xn2jGYALp0/g7usWsHrRZIpyxiS4wtSlQBeRmDvScYrndh7mmR0NvLK3ie4eZ2ZJNp+9eg7XLSljWuH4RJeYFhToIhITNc2dPLujgWd3NLDhQAvuML1oPJ9YVcn7z5/Cgil5mm4YZQp0EYkKd2dP4zGe3d7AMzsa2FHXDsC8yXncecVs3ruwlPMm5SrEY2jIQDezB4BrgUZ3XzjIcgO+D6wGOoG/cfdN0S5URJLPye4eXtl3lBd2NfL87iMcbO4EQj3xL66ey3sWlDK9SPPF4yWSEfqDwD3AQ6dZfg0wO/xaAfwk/KeIpKGa5k5e2N3IH3c18vLeo5wKBBmXOZpVs4pYc9kMrpo/iUl5qf1szlQ1ZKC7+4tmVnGGVa4DHnJ3B9aZWYGZTXb3+ijVKCIJdLK7hw1vtPCn10Oj8OrGY0CoH37j8nLeNXciKyoLR/x9VJJBNHroZUBNv/e14c/eFuhmtgZYA1BeXh6FXYtItAWDzo66dtZWN/Hn6iZefaOZrkCQrNGjWDGjMBTi55UwoyQn0aXKANEI9MHOcPhgK7r7fcB9AMuWLRt0HRGJv5rmTtZWN7F2TxMv722ipbMbgLmludxy8XQunVXM8spCXXaf5KLx06kFpvV7PxWoi8J2RSRGapo7WbfvKOv3N7N+/1Fqmk8AMClvDO+eO4l3zC7mkllFTMxVLzyVRCPQnwA+bWa/InQytE39c5Hk4e4cbO5k/b7mvhA/1BoK8AnjM1leWcgnVlXyjtnFzCzJ0bTCFBbJtMVfApcDxWZWC3wFyARw93uBpwhNWawmNG3xtlgVKyJDc3f2Nx0Pjb73HWXdvmYa2k8CUJSdxYoZhay5bAYrZhQyZ2Iuo0YpwNNFJLNcbhxiuQOfilpFIjIspwI97KhrZ9OBFjYeaGHDgRaOdJwCoDhnDCtmFHJxZSEXzyhi1kSNwNOZznCIpJimY6dC4X2whU0HWthS20ZXIAjAtMJxXDKziOWVhayoLGJmSbYCfARRoIsksWAwdDn9xvDoe9PBFvY3HQcgc7SxsCyfj188nWUVE7igfAITdUHPiKZAF0kiRzpOsbW2lS01rbxW08rmmlY6TgaAUP/7gukTuOGiaVw4fQILy/J1MY+8hQJdJEGOnwqw/VAbW2pb2VLTxuaa1r7ZJ6NHGXMm5XLt4iksmz6BC6dPYHrReLVP5IwU6CJxEOgJ8vrhY2yuCY2+t9S28vrhDoLhy+umFY5jaXkBt62q4PxpBSyYksf4LP31lOHRb4xIlHX3BKluPMa2Q23sONTGtkNt7Kxv52R36MRlwfhMzp9awNULSlk6rYDFU/P1lB6JCgW6yDk4Fehhz+FQeG8Pv6oaOvpmnWRnjWbBlHxuXF7OkmkFLJlWQHmhWicSGwp0kQid7O5hV0PHW0berx/uoLsn1DfJHZvBwin53LpyOgvL8llYlk9lUbYu3JG4UaCLDKKzK0BVfTvbatvYXtfO9kNt7Gk8Rk+46V0wPpOFU/K5/dIZLCzLY1FZPtMmjFd4S0Ip0GXEazp2iqr6dqrq29lZ186Ounb2HjnWd8KyKDuLhWX5XDlvEgvL8lhYlk9ZwTi1TSTpKNBlxAj0BNnfdJyd9e1U1XeE/2zvu0weoDRvLAum5HHNosksKstnYVkepXljFd6SEhTokpbaTnSzKxzYveH9+uEOToVPVmaONmZNzOUds4uZPzmPeeFXYXZWgisXOXsKdElpwaBT09IZapfUd/S1TXov0AEozM5i3uTQgxp6g3vWxByyMkYlsHKR6FOgS8ro7Aqwu6GDqnBwV9W3s6uhg2OnQpfGjzKoLM5maXkBN60oZ/7kPOZPyWNi7hi1TGREUKBL0gkGnQPNnexuCLVLdjd0sKuhnQPNnXj4RGXOmAzmTc7lwxeU9Y26z5uUy7gs3dtERi4FuiRUy/EudoUDe1d9B7sOd/B6QwcnunsAMIPKomzmTc7jQ0uncl5pLgum5DF1gmaZiAykQJe4OBXoYW/jcXY1tIfaJg0d7G5o53D7mzNMCrOzmFuayw3LpzGvNI+5k3OZPVGjbpFIKdAlqtyd+raT7BrQLtl35DiB8MTurNGjmDUxh1WziplbmsvccHiX5KjXLXIuFOhy1o6dCvQF9q5weFc1tPfdvxugrGAcc0tzuWr+JM4rzWNeaS4VxdlkjtYME5FoU6DLkHqCzoGjx0O97vp2qsIhXtP85tTAnDEZzC3N5QPnT2Hu5DzmluYyZ1Iu+eMyE1i5yMiiQJe3aO3s6gvuXeFed/+TlL1TAxdPLeAjF04LzS4pzdVJSpEkoEAfoQI9QfY1He+by90b4PVtJ/vWmTA+k3mT87hxeTlzJ+cyrzSP2ZNy9NgzkSSlQB8BOk52h0bb4asod4bDu6vfZfAzS3K4eEZR6CTl5FCvu0QX5IikFAV6GumdYbKzLnznwPDrwNHOvnUmjM9k/pQ8bl05nflT8phbmsfMEl0GL5IOFOgpqvcxZwPDu7Wzu2+diqLxLJiSx19fOJX5U/KYPzmfSXkadYukKwV6CjjR1cPO+ja21baxI9wy2XP4GF09oZbJmIxRzC3N5ZqFpX33LzmvNI+cMfrxiowk+hufZELhHXpCztbatvCTcjre8rCF+VPyuO3SCuZPzmPBlDwqirLJ0LxukRFPgZ5AJ7sHC+83H3NWnJPForJ83rNgEgvL8lk8tUAtExE5LQV6nLg7bxztZHNNC68dbOW1g61U1bf3XQ5flJ3Foqn5XDV/EovK8lk0NV9PyhGRYVGgx0jbiW621ISCe3NNC5trWmkJn7DMzhrN4qkFrLlsBounFrB4aj6T8xXeInJuFOhR4O7Utpzg1f3NvLq/mQ0Hmtl75DgQuv3r7Ik5XD2/lCXlBSwtL2D2xFxG6+nwIhJlCvSz4O7sPXI8HOBHeXV/M3XhKyzzxmawrKKQDy4pY2n5BBZPyydvrO5nIiKxp0CP0KHWE/x5TxMvVTfxyt4mmo51AVCcM4YVlYX8XWUhyysLOW9SLqM0+haRBFCgn0bHyW7+XN3E2uom/lx9lP1NoRZKSe4YLp1VzMqZRVxUUUhlcbZ63yKSFBTo/dQ0d/KHqsP8YVcj6/YdpbvHGZ81motnFHHLxdO5dHYxsyfmKMBFJCmN6EB3d3bWt/P0tgZ+X3WYXQ0dAMwoyea2VZW8e+5ELiifoPuciEhKGJGBXtPcyeObD/H45jr2NB5j9Chj2fQJfGn1PK6YN5EZJTmJLlFEZNgiCnQzey/wfWA0cL+7f3PA8suBx4H94Y8edfd/jWKd5+xkdw9Pbq3nl68eZMOBFgAuqpjA1z+4kNWLJlOYnZXgCkVEzs2QgW5mo4EfAVcBtcBfzOwJd985YNWX3P3aGNR4TurbTvDA2v38emMtrZ3dzCjO5nPvPY8PnD+FqRPGJ7o8EZGoiWSEvhyodvd9AGb2K+A6YGCgJ5X9Tcf59z/t5ZFNtbjDexaWcvOKclbOKNJJTRFJS5EEehlQ0+99LbBikPVWmtkWoA74rLvvGLiCma0B1gCUl5cPv9oItJ3o5rvPvc7P1x0gY5Rx0/Jy/tdlMzQaF5G0F0mgDzac9QHvNwHT3f2Yma0GHgNmv+2b3O8D7gNYtmzZwG2cE3fnNxtr+dYzu2g+3sVNK8q584o5lOSOieZuRESSViSBXgtM6/d+KqFReB93b+/39VNm9mMzK3b3puiUeWbtJ7v53K+38syOBi4oL+DB25azsCw/HrsWEUkakQT6X4DZZlYJHAJuAG7qv4KZlQKH3d3NbDkwCjga7WIH09h+kpvvX8++puN8afU8br+0Upfei8iINGSgu3vAzD4NPEto2uID7r7DzD4ZXn4vcD1wh5kFgBPADe4e1ZbKYI50nOKj963jcPtJfv6J5VwyqzjWuxQRSVoRzUN396eApwZ8dm+/r+8B7oluaUPWxOd+s4W61hM8/LcrWFZRGM/di4gknZS9pv3h9Qd5fvcRvnDNXIW5iAgpHOgPvvwGF06fwMdXViS6FBGRpJCSgd4TdA4e7eSiikKdABURCUvJQK9rPUFXT5DKYl0sJCLSKyUDvfdhE9OLshNciYhI8kjJQH/jaCjQK4sV6CIivVIy0Pc3HWd81mgm6rJ+EZE+KRnobzQdZ3qRnuUpItJfSgb6gaOdOiEqIjJAygV6oCfIweZOKnRCVETkLVIu0A+1niAQdCp0QlRE5C1SLtB7pyxqhC4i8lYpF+g5YzK4ev4kZpQo0EVE+ovobovJZFlFoW7GJSIyiJQboYuIyOAU6CIiaUKBLiKSJhToIiJpQoEuIpImFOgiImlCgS4ikiYU6CIiacLcPTE7NjsCHDjLby8GmqJYTirQMY8MOuaR4VyOebq7lwy2IGGBfi7MbIO7L0t0HfGkYx4ZdMwjQ6yOWS0XEZE0oUAXEUkTqRro9yW6gATQMY8MOuaRISbHnJI9dBERebtUHaGLiMgACnQRkTSR1IFuZu81s91mVm1m/zzIcjOzH4SXbzWzCxJRZzRFcMw3h491q5m9bGbnJ6LOaBrqmPutd5GZ9ZjZ9fGsLxYiOWYzu9zMNpvZDjP7U7xrjLYIfrfzzey3ZrYlfMy3JaLOaDGzB8ys0cy2n2Z59PPL3ZPyBYwG9gIzgCxgCzB/wDqrgacBAy4G1ie67jgc8yXAhPDX14yEY+633h+Bp4DrE113HH7OBcBOoDz8fmKi647DMX8R+Fb46xKgGchKdO3ncMyXARcA20+zPOr5lcwj9OVAtbvvc/cu4FfAdQPWuQ54yEPWAQVmNjnehUbRkMfs7i+7e0v47TpgapxrjLZIfs4A/wA8AjTGs7gYieSYbwIedfeDAO6e6scdyTE7kGtmBuQQCvRAfMuMHnd/kdAxnE7U8yuZA70MqOn3vjb82XDXSSXDPZ7bCf0Ln8qGPGYzKwM+BNwbx7piKZKf8xxggpm9YGYbzezjcasuNiI55nuAeUAdsA24092D8SkvIaKeX8n8kGgb5LOBcywjWSeVRHw8ZvYuQoF+aUwrir1Ijvl7wOfdvSc0eEt5kRxzBnAhcAUwDnjFzNa5++uxLi5GIjnm9wCbgXcDM4HnzOwld2+PdXEJEvX8SuZArwWm9Xs/ldC/3MNdJ5VEdDxmthi4H7jG3Y/GqbZYieSYlwG/Cod5MbDazALu/lh8Soy6SH+3m9z9OHDczF4EzgdSNdAjOebbgG96qMFcbWb7gbnAq/EpMe6inl/J3HL5CzDbzCrNLAu4AXhiwDpPAB8Pny2+GGhz9/p4FxpFQx6zmZUDjwK3pPBorb8hj9ndK929wt0rgN8Af5/CYQ6R/W4/DrzDzDLMbDywAqiKc53RFMkxHyT0fySY2STgPGBfXKuMr6jnV9KO0N09YGafBp4ldIb8AXffYWafDC+/l9CMh9VANdBJ6F/4lBXhMX8ZKAJ+HB6xBjyF71QX4TGnlUiO2d2rzOwZYCsQBO5390Gnv6WCCH/OdwMPmtk2Qu2Iz7t7yt5W18x+CVwOFJtZLfAVIBNil1+69F9EJE0kc8tFRESGQYEuIpImFOgiImlCgS4ikiYU6CIiaUKBLiKSJhToIiJp4v8DzNEECUJ0f8EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's plot some values for H_inv(q)\n",
    "\n",
    "qq = np.linspace(1e-5, 1 - 1e-5, 200)\n",
    "xx = [H_inv(q) for q in qq]\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.plot(qq, xx)\n",
    "plt.title('H^{-1}(q)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test statistic: \t\t0.936\n",
      "H^{-1}(test_statistic): \t1.313\n"
     ]
    }
   ],
   "source": [
    "h_inv_test_statistic = H_inv(test_statistic)\n",
    "\n",
    "print('Test statistic: \\t\\t%.3f' % test_statistic)\n",
    "print('H^{-1}(test_statistic): \\t%.3f' % h_inv_test_statistic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that the inverse of the test statistic is greater than 1, we have that the test from Theorem 16.16 rejects the null hypothesis for any $\\alpha > 0$ -- so it claims the distributions are distinct at any approximate $\\alpha$ level -- so, by Theorem 16.15, they are associated."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}