likelihood-ratio-test-from-first-principles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "690fceab",
   "metadata": {},
   "source": [
    "## Likelihood Ratio Test \n",
    "\n",
    "The Likelihood Ratio Test is nicely described here: \n",
    "https://en.wikipedia.org/wiki/Likelihood-ratio_test \n",
    "\n",
    "Briefly, a model with more parameters (here the \"alternative\") will always fit a given data set \n",
    "at least as well — i.e., have the same or greater log-likelihood — as a \"nested\" model, with fewer parameters (here called \"null\"). Whether the fit is significantly better, and should thus be preferred, is determined by deriving how likely it is that the observed log likelihood differences could have occurred \n",
    "by chance (sampling error) alone, if the null model were true.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6018fefd",
   "metadata": {},
   "source": [
    "## The simulated null distribution of the likelihood-ratio statistic is approximated by the chi-square distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3d5361b2",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZZklEQVR4nO3deVhUZf8G8PvMDLOwgyibIO5buCviUloatmjWm6+piZpZljvpq7y55Erlnpm+mqktllZq/dI0JZc0ldRQKcQNxYVFVEB2mDm/P5DJUdQZmDMHhvtzXXMxc7b5nsGYu+d5znkEURRFEBEREdkJhdwFEBEREVkTww0RERHZFYYbIiIisisMN0RERGRXGG6IiIjIrjDcEBERkV1huCEiIiK7wnBDREREdoXhhoiIiOwKww0RERHZFVnDzf79+9G7d2/4+flBEARs3br1kfvs3bsXbdq0gUajQYMGDbBu3TrJ6yQiIqKqQ9Zwk5OTg5YtW2L58uVmbZ+YmIjnnnsO3bt3R2xsLMaPH4/XX38dO3fulLhSIiIiqiqEyjJxpiAI2LJlC/r27fvAbSZPnoxt27YhLi7OuOyVV15BRkYGduzYYYMqiYiIqLJTyV2AJQ4dOoQePXqYLAsLC8P48eMfuE9BQQEKCgqMrw0GA27evIkaNWpAEASpSiUiIiIrEkURt2/fhp+fHxSKh3c8Valwk5KSAm9vb5Nl3t7eyMrKQl5eHnQ63X37REVFYebMmbYqkYiIiCR0+fJl1K5d+6HbVKlwUx6RkZGIiIgwvs7MzERgYCAuX74MV1dXq7/f8l/Pmbwe9WQDq78HERFRdZOVlYWAgAC4uLg8ctsqFW58fHyQmppqsiw1NRWurq5lttoAgEajgUajuW+5q6urJOFG6+R83/sQERGRdZgzpKRK3ecmNDQU0dHRJst27dqF0NBQmSoiIiKiykbWcJOdnY3Y2FjExsYCKLnUOzY2FklJSQBKupTCw8ON248cORIXLlzAf/7zH5w+fRqffPIJNm3ahAkTJshRPhEREVVCsoabo0ePonXr1mjdujUAICIiAq1bt8b06dMBAMnJycagAwB169bFtm3bsGvXLrRs2RILFy7Ep59+irCwMFnqJyIiosqn0tznxlaysrLg5uaGzMxMScbDLN51xuT1hJ6NrP4eRETVmV6vR1FRkdxlkATUavUDL/O25Pu7Sg0oJiKi6ksURaSkpCAjI0PuUkgiCoUCdevWhVqtrtBxGG6IiKhKKA02tWrVgqOjI2/EamcMBgOuXbuG5ORkBAYGVuj3y3BDRESVnl6vNwabGjVqyF0OSaRmzZq4du0aiouL4eDgUO7jVKlLwYmIqHoqHWPj6OgocyUkpdLuKL1eX6HjMNwQEVGVwa4o+2at3y/DDREREdkVhhsiIiI70K1bN4wfP97s7detWwd3d3fJ6pETBxQTEVGVde+9xaRW3e5dNnfuXGzbtg2xsbFQq9VV5jJ8ttwQERFRmQoLC9GvXz+89dZbcpdiEYYbIiIiiXTr1g1jxozB+PHj4eHhAW9vb6xevRo5OTkYNmwYXFxc0KBBA/z8888m++3btw8dOnSARqOBr68vpkyZguLiYuP6nJwchIeHw9nZGb6+vli4cOF9711QUICJEyfC398fTk5OCAkJwd69ey2qf+bMmZgwYQKCg4PLdf5yYbiR0MW8w1hxYgX+Sv9L7lKIiEgm69evh5eXF2JiYjBmzBi89dZb6NevHzp16oTjx4/j6aefxuDBg5GbmwsAuHr1Kp599lm0b98eJ06cwIoVK7BmzRrMmTPHeMxJkyZh3759+OGHH/DLL79g7969OH78uMn7jh49GocOHcI333yDkydPol+/fujVqxfOnj1r0/OXA8ONhC7k/YZPYj/ByfSTcpdCREQyadmyJaZOnYqGDRsiMjISWq0WXl5eGDFiBBo2bIjp06fjxo0bOHmy5Lvik08+QUBAAD7++GM0adIEffv2xcyZM7Fw4UIYDAZkZ2djzZo1WLBgAZ566ikEBwdj/fr1Ji07SUlJWLt2Lb799lt07doV9evXx8SJE9GlSxesXbtWro/CZjigWEJKQQMAKCgukLkSIiKSS4sWLYzPlUolatSoYdLN4+3tDQBIS0sDAMTHxyM0NNTkni+dO3dGdnY2rly5glu3bqGwsBAhISHG9Z6enmjcuLHx9alTp6DX69GokekA6IKCgmpxh2eGGwmp7oSbPH2ezJUQEZFc7p1GQBAEk2WlIcZgMFjtPbOzs6FUKnHs2DEolUqTdc7OzlZ7n8qK4UZCKqHkNtJsuSEiInM1bdoU33//PURRNAafgwcPwsXFBbVr14anpyccHBxw5MgRBAYGAgBu3bqFM2fO4IknngAAtG7dGnq9Hmlpaejatats5yIXjrmRkPJOuMnX58tcCRERVRVvv/02Ll++jDFjxuD06dP44YcfMGPGDEREREChUMDZ2RnDhw/HpEmT8OuvvyIuLg5Dhw6FQvHPV3qjRo0waNAghIeHY/PmzUhMTERMTAyioqKwbds2s2tJSkpCbGwskpKSoNfrERsbi9jYWGRnZ0tx6lbDlhsJlXZL5Rcz3BARkXn8/f2xfft2TJo0CS1btoSnpyeGDx+OqVOnGreZP38+srOz0bt3b7i4uOCdd95BZmamyXHWrl2LOXPm4J133sHVq1fh5eWFjh074vnnnze7lunTp2P9+vXG161btwYA7NmzB926davYiUpIEEVRlLsIW8rKyoKbmxsyMzPh6upq9ePffbfMuOwfEZO1Hs/Vew7vd33f6u9FRFRd5OfnIzExEXXr1oVWq5W7HJLIw37Plnx/s1tKQmy5ISIisj2GGwkpGW6IiIhsjuFGQioOKCYiIrI5hhsJsVuKiIjI9hhuJGS8FJzhhoiIyGYYbiRkbLlhtxQREZHNMNxISMWWGyIiIptjuJEQ71BMRERkeww3ElJxVnAiIiKbY7iRUGm4KRaLUWQokrkaIiKyZ926dcP48ePN3n7dunVwd3eXrB45cW4pCZV2SwElrTcOaoeHbE1ERBbbE2Xb9+seadv3k9HFixcxe/Zs/Prrr0hJSYGfnx9effVVvPvuu1Cr1Y8+gIwYbiSkhAMECBAhIl+fD2c4y10SERGRWU6fPg2DwYD//e9/aNCgAeLi4jBixAjk5ORgwYIFcpf3UOyWkpAgCNCqSib+yivOk7kaIiKytW7dumHMmDEYP348PDw84O3tjdWrVyMnJwfDhg2Di4sLGjRogJ9//tlkv3379qFDhw7QaDTw9fXFlClTUFxcbFyfk5OD8PBwODs7w9fXFwsXLrzvvQsKCjBx4kT4+/vDyckJISEh2Lt3r9m19+rVC2vXrsXTTz+NevXqoU+fPpg4cSI2b95c7s/DVhhuJKZVloQbDiomIqqe1q9fDy8vL8TExGDMmDF466230K9fP3Tq1AnHjx/H008/jcGDByM3NxcAcPXqVTz77LNo3749Tpw4gRUrVmDNmjWYM2eO8ZiTJk3Cvn378MMPP+CXX37B3r17cfz4cZP3HT16NA4dOoRvvvkGJ0+eRL9+/dCrVy+cPXu23OeSmZkJT0/Pcu9vKww3EittueHl4ERE1VPLli0xdepUNGzYEJGRkdBqtfDy8sKIESPQsGFDTJ8+HTdu3MDJkycBAJ988gkCAgLw8ccfo0mTJujbty9mzpyJhQsXwmAwIDs7G2vWrMGCBQvw1FNPITg4GOvXrzdp2UlKSsLatWvx7bffomvXrqhfvz4mTpyILl26YO3ateU6j3PnzmHZsmV48803rfK5SIljbiSmUZZcMcVuKSKi6qlFixbG50qlEjVq1EBwcLBxmbe3NwAgLS0NABAfH4/Q0FAIgmDcpnPnzsjOzsaVK1dw69YtFBYWIiQkxLje09MTjRs3Nr4+deoU9Ho9GjVqZFJLQUEBatSoYfE5XL16Fb169UK/fv0wYsQIi/e3NYYbielUOgBAgZ7dUkRE1ZGDg+mVsoIgmCwrDTEGg8Fq75mdnQ2lUoljx45BqVSarHN2tuzilmvXrqF79+7o1KkTVq1aZbUapcRwIzFjtxSnYCAiIjM0bdoU33//PURRNAafgwcPwsXFBbVr14anpyccHBxw5MgRBAYGAgBu3bqFM2fO4IknngAAtG7dGnq9HmlpaejatWu5a7l69Sq6d++Otm3bYu3atVAoqsZolqpRZRXGbikiIrLE22+/jcuXL2PMmDE4ffo0fvjhB8yYMQMRERFQKBRwdnbG8OHDMWnSJPz666+Ii4vD0KFDTYJHo0aNMGjQIISHh2Pz5s1ITExETEwMoqKisG3bNrPquHr1Krp164bAwEAsWLAA169fR0pKClJSUqQ6dathy43ESltu2C1FRETm8Pf3x/bt2zFp0iS0bNkSnp6eGD58OKZOnWrcZv78+cjOzkbv3r3h4uKCd955B5mZmSbHWbt2LebMmYN33nkHV69ehZeXFzp27Ijnn3/erDp27dqFc+fO4dy5c6hdu7bJOlEUK36iEhLEyl6hlWVlZcHNzQ2ZmZlwdXW1+vEX7zpj8jpZ/Sl+vvgzJrefjFebvWr19yMiqg7y8/ORmJiIunXrQqvVyl0OSeRhv2dLvr/ZLSUxjaqkW4qXghMREdkGw43ESm/ixwHFREREtsFwI7HSS8EZboiIiGyD4UZi7JYiIiKyLYYbibFbioiIyLYYbiTGuaWIiIhsi+FGYpwVnIiIyLYYbiRW2nKTp+cdiomIiGyB4UZiHHNDRERkWww3EuP0C0REVJaLFy9CEATExsY+cJt169bB3d3dZjXZC4YbiXFWcCIiKq/+/fvjzJkzj96QTHDiTImVdktxVnAiIrKUTqeDTqeTuwyLiaIIvV4PlUqemMGWG4mxW4qIqHozGAz48MMP0aBBA2g0GgQGBmLu3LnG9RcuXED37t3h6OiIli1b4tChQ8Z15nRL7d27Fx06dICTkxPc3d3RuXNnXLp0ybj+/fffh7e3N1xcXDB8+HBMmTIFrVq1Mq7v1q0bxo8fb3LMvn37YujQocbXX3zxBdq1awcXFxf4+Phg4MCBSEtLM6lBEAT8/PPPaNu2LTQaDQ4cOACDwYCoqCjUrVsXOp0OLVu2xHfffWfZB1gObLmRGLuliIisTxRF2VrEdSodBEEwe/vIyEisXr0aixcvRpcuXZCcnIzTp08b17/77rtYsGABGjZsiHfffRcDBgzAuXPnzGr1KC4uRt++fTFixAh8/fXXKCwsRExMjLG+TZs24b333sPy5cvRpUsXfPHFF/joo49Qr149i865qKgIs2fPRuPGjZGWloaIiAgMHToU27dvN9luypQpWLBgAerVqwcPDw9ERUXhyy+/xMqVK9GwYUPs378fr776KmrWrIknnnjCohoswXAjMY3yn+kXRFG06D8IIiIqW15xHkI2hMjy3kcGHoGjg6NZ296+fRtLly7Fxx9/jCFDhgAA6tevjy5duuDixYsAgIkTJ+K5554DAMycORPNmzfHuXPn0KRJk0cePysrC5mZmXj++edRv359AEDTpk2N65csWYLhw4dj+PDhAIA5c+Zg9+7dyM+37H+4X3vtNePzevXq4aOPPkL79u2RnZ0NZ2dn47pZs2ahZ8+eAICCggLMmzcPu3fvRmhoqHHfAwcO4H//+5+k4YbdUhIrnTgTYNcUEVF1Ex8fj4KCAjz11FMP3KZFixbG576+vgBg0uVTKikpCc7OzsbHvHnz4OnpiaFDhyIsLAy9e/fG0qVLkZycbPL+ISGmIbA0aFji2LFj6N27NwIDA+Hi4mIMJklJSSbbtWvXzvj83LlzyM3NRc+ePU3q/vzzz3H+/HmLa7AEW24kVtpyA5R0TZV2UxERUfnpVDocGXhEtvc2e1szBgM7ODgYn5e27hsMhvu28/PzM7ls3NPTEwCwdu1ajB07Fjt27MDGjRsxdepU7Nq1Cx07djSrRoVCAVEUTZYVFRUZn+fk5CAsLAxhYWH46quvULNmTSQlJSEsLAyFhYUm+zk5ORmfZ2dnAwC2bdsGf39/k+00Gg2kxHAjMZVCBZVChWJDMeeXIiKyEkEQzO4aklPDhg2h0+kQHR2N119/vULHUqlUaNCgQZnrWrdujdatWyMyMhKhoaHYsGEDOnbsiKZNm+LIkSMIDw83bnv48GGTfWvWrGnS2qPX6xEXF4fu3bsDAE6fPo0bN27g/fffR0BAAADg6NGjj6y3WbNm0Gg0SEpKkrQLqiwMNzagU+pw23Cbg4qJiKoZrVaLyZMn4z//+Q/UajU6d+6M69ev46+//npoV5W5EhMTsWrVKvTp0wd+fn5ISEjA2bNnjWFm3LhxGDp0KNq1a4fOnTvjq6++wl9//WUyoPjJJ59EREQEtm3bhvr162PRokXIyMgwrg8MDIRarcayZcswcuRIxMXFYfbs2Y+szcXFBRMnTsSECRNgMBjQpUsXZGZm4uDBg3B1dTWOQZICw40NaFQa3C66zZYbIqJqaNq0aVCpVJg+fTquXbsGX19fjBw50irHdnR0xOnTp7F+/XrcuHEDvr6+GDVqFN58800AJTcBPH/+PP7zn/8gPz8f//rXv/DWW29h586dxmO89tprOHHiBMLDw6FSqTBhwgRjqw1Q0rKzbt06/Pe//8VHH32ENm3aYMGCBejTp88j65s9ezZq1qyJqKgoXLhwAe7u7mjTpg3++9//WuX8H0QQ7+1os3NZWVlwc3NDZmYmXF1drX78xbtM7yQ5oWcjPPP9M7iSfQVfPPMFWtVqZfX3JCKyd/n5+UhMTETdunWh1XLsYkW899572Lp160OnfZDLw37Plnx/82opGzDe64YtN0RERJJjuLEBzgxORERkOww3NsCWGyIiqizee++9StklZU2yh5vly5cjKCgIWq0WISEhiImJeej2S5YsQePGjaHT6RAQEIAJEyZYfKdFW+MUDERERLYja7jZuHEjIiIiMGPGDBw/fhwtW7ZEWFhYmXdmBIANGzZgypQpmDFjBuLj47FmzRps3LhR8lHXFVXaLVVQzDsUExFVRDW7BqbasdbvV9Zws2jRIowYMQLDhg1Ds2bNsHLlSjg6OuKzzz4rc/vff/8dnTt3xsCBAxEUFISnn34aAwYMeGRrj9zYLUVEVDGld/HNzc2VuRKSUukdj5VKZYWOI9t9bgoLC3Hs2DFERkYalykUCvTo0cNkuve7derUCV9++SViYmLQoUMHXLhwAdu3b8fgwYMf+D4FBQUoKPinxSQrK8t6J2Gm0ikY5JrBloioqlMqlXB3dze27Ds6OnIiYjtjMBhw/fp1ODo6mjUj+sPIFm7S09Oh1+vh7e1tstzb29tkKvi7DRw4EOnp6ejSpQtEUURxcTFGjhz50G6pqKgozJw506q1W6p0HhJOnElEVH4+Pj4Ayp5UkuyDQqFAYGBghYNrlbpD8d69ezFv3jx88sknCAkJwblz5zBu3DjMnj0b06ZNK3OfyMhIREREGF9nZWUZ58awFQ4oJiKqOEEQ4Ovri1q1aplM7Ej2Q61WQ6Go+IgZ2cKNl5cXlEolUlNTTZanpqYa0/m9pk2bhsGDBxsnHwsODkZOTg7eeOMNvPvuu2V+IBqNRvLZRx+F3VJERNajVCorPCaD7JtsA4rVajXatm2L6Oho4zKDwYDo6GiEhoaWuU9ubu59Aab0H3hlHkHPbikiIiLbkbVbKiIiAkOGDEG7du3QoUMHLFmyBDk5ORg2bBgAIDw8HP7+/oiKigIA9O7dG4sWLULr1q2N3VLTpk1D7969K3WK5x2KiYiIbEfWcNO/f39cv34d06dPR0pKClq1aoUdO3YYBxknJSWZtNRMnToVgiBg6tSpuHr1KmrWrInevXtj7ty5cp2CWTSqO91SenZLERERSU32AcWjR4/G6NGjy1y3d+9ek9cqlQozZszAjBkzbFCZ9ZQOKOZN/IiIiKQn+/QL1YFOWTLmht1SRERE0mO4sYHSbineoZiIiEh6DDc2wAHFREREtsNwYwOcW4qIiMh2GG5sgLOCExER2Q7DjQ2UttzwUnAiIiLpMdzYQGnLTbGhGMWGYpmrISIism8MNzZQ2nIDcAoGIiIiqTHc2EDpxJkAJ88kIiKSGsONDQiCwMvBiYiIbIThxkaMUzCwW4qIiEhSDDc2YrzXDVtuiIiIJMVwYyOl3VIcc0NERCQthhsbYbcUERGRbTDc2EjpFVPsliIiIpIWw42NcH4pIiIi22C4sRGdUgeALTdERERSY7ixEY3qTrcUW26IiIgkxXBjI7yJHxERkW0w3NgIx9wQERHZBsONjbDlhoiIyDYYbmyEdygmIiKyDYYbG2G3FBERkW0w3NgIu6WIiIhsg+HGRthyQ0REZBsMNzbC6ReIiIhsg+HGRnQq3qGYiIjIFhhubISzghMREdkGw42NlHZL5RXnyVwJERGRfWO4sZHSbim23BAREUmL4cZGeCk4ERGRbTDc2AhnBSciIrINhhsbuftqKVEUZa6GiIjIfjHc2Ehpt5QIEYWGQpmrISIisl8MNzZS2i0FcNwNERGRlBhubMRB4QCVoALAcENERCQlhhsb4qBiIiIi6ankLsDeLd51xvjcoHcAwJYbIiIiKbHlxoaUghoAW26IiIikxHBjQ6rScMOWGyIiIslYHG727NkjRR3VgkooGXPDKRiIiIikY3G46dWrF+rXr485c+bg8uXLUtRkt0rDDSfPJCIiko7F4ebq1asYPXo0vvvuO9SrVw9hYWHYtGkTCgt5Y7pHUbJbioiISHIWhxsvLy9MmDABsbGxOHLkCBo1aoS3334bfn5+GDt2LE6cOCFFnXaB3VJERETSq9CA4jZt2iAyMhKjR49GdnY2PvvsM7Rt2xZdu3bFX3/9Za0a7Qa7pYiIiKRXrnBTVFSE7777Ds8++yzq1KmDnTt34uOPP0ZqairOnTuHOnXqoF+/ftautcpjtxQREZH0LL6J35gxY/D1119DFEUMHjwYH374IR577DHjeicnJyxYsAB+fn5WLdQelF4Kzm4pIiIi6Vgcbv7++28sW7YML730EjQaTZnbeHl58ZLxMpS23LBbioiISDoWd0vNmDED/fr1uy/YFBcXY//+/QAAlUqFJ554wjoV2pHSMTe8QzEREZF0LA433bt3x82bN+9bnpmZie7du1ulKHtlvFqqmN1SREREUrE43IiiCEEQ7lt+48YNODk5WaUoe8W5pYiIiKRn9pibl156CQAgCAKGDh1q0i2l1+tx8uRJdOrUyfoV2hFeCk5ERCQ9s8ONm5sbgJKWGxcXF+h0OuM6tVqNjh07YsSIEdav0I7wJn5ERETSMzvcrF27FgAQFBSEiRMnsguqHDgrOBERkfQsvhR8xowZUtRRLXDMDRERkfTMCjdt2rRBdHQ0PDw80Lp16zIHFJc6fvy41YqzN8ZLwdlyQ0REJBmzws0LL7xgHEDct29fKeuxa+yWIiIikp5Z4eburih2S5Ufb+JHREQkPYvvc3P58mVcuXLF+DomJgbjx4/HqlWrrFqYPVKWXgpelAdRFGWuhoiIyD5ZHG4GDhxonDcqJSUFPXr0QExMDN59913MmjXL6gXaE63gAgAoFouRU5QjczVERET2yeJwExcXhw4dOgAANm3ahODgYPz+++/46quvsG7dOmvXZ1dUCo2xa+pW/i2ZqyEiIrJPFoeboqIi4+Di3bt3o0+fPgCAJk2aIDk52brV2SGtwhUAcLPg/vm5iIiIqOIsDjfNmzfHypUr8dtvv2HXrl3o1asXAODatWuoUaOG1Qu0N6Xhhi03RERE0rA43HzwwQf43//+h27dumHAgAFo2bIlAODHH380dldZYvny5QgKCoJWq0VISAhiYmIeun1GRgZGjRoFX19faDQaNGrUCNu3b7f4feWiYbghIiKSlMV3KO7WrRvS09ORlZUFDw8P4/I33ngDjo6OFh1r48aNiIiIwMqVKxESEoIlS5YgLCwMCQkJqFWr1n3bFxYWomfPnqhVqxa+++47+Pv749KlS3B3d7f0NGRj7JbKZ7cUERGRFCwONwCgVCpNgg1QMueUpRYtWoQRI0Zg2LBhAICVK1di27Zt+OyzzzBlypT7tv/ss89w8+ZN/P7773BwcCj3+8qpNNzsTriArJQzAIAJPRvJWRIREZFdsbhbKjU1FYMHD4afnx9UKhWUSqXJw1yFhYU4duwYevTo8U8xCgV69OiBQ4cOlbnPjz/+iNDQUIwaNQre3t547LHHMG/ePOj1+ge+T0FBAbKyskwecioNN/kGeesgIiKyVxa33AwdOhRJSUmYNm0afH19HzrP1MOkp6dDr9fD29vbZLm3tzdOnz5d5j4XLlzAr7/+ikGDBmH79u04d+4c3n77bRQVFT3wzslRUVGYOXNmuWqUgo7hhoiISFIWh5sDBw7gt99+Q6tWrSQo5+EMBgNq1aqFVatWQalUom3btrh69Srmz5//wHATGRmJiIgI4+usrCwEBATYquT7aJVuABhuiIiIpGJxuAkICLDK1AFeXl5QKpVITU01WZ6amgofH58y9/H19YWDg4NJ91fTpk2RkpKCwsJCqNXq+/bRaDTG+/JUBhpFyV2KGW6IiIikYfGYmyVLlmDKlCm4ePFihd5YrVajbdu2iI6ONi4zGAyIjo5GaGhomft07twZ586dg8FgMC47c+YMfH19yww2lRHH3BAREUnL4nDTv39/7N27F/Xr14eLiws8PT1NHpaIiIjA6tWrsX79esTHx+Ott95CTk6O8eqp8PBwREZGGrd/6623cPPmTYwbNw5nzpzBtm3bMG/ePIwaNcrS05CNVlHSLVUs5qNYLJC5GiIiIvtjcbfUkiVLrPbm/fv3x/Xr1zF9+nSkpKSgVatW2LFjh3GQcVJSEhSKf/JXQEAAdu7ciQkTJqBFixbw9/fHuHHjMHnyZKvVJDW14AgFVDCgGPmG23BWVp4uMyIiInsgiNYYQFOFZGVlwc3NDZmZmXB1dbX68RfvOvPIbb5OeR15hlvo4zUfXup6vM8NERHRI1jy/W1xtxQAnD9/HlOnTsWAAQOQlpYGAPj555/x119/ledw1c4/424yZa6EiIjI/lgcbvbt24fg4GAcOXIEmzdvRnZ2NgDgxIkTD7wcm0xxUDEREZF0LA43U6ZMwZw5c7Br1y6TK5SefPJJHD582KrF2SuGGyIiIulYHG5OnTqFF1988b7ltWrVQnp6ulWKsndaJcMNERGRVCwON+7u7khOTr5v+Z9//gl/f3+rFGXv2HJDREQkHYvDzSuvvILJkycjJSUFgiDAYDDg4MGDmDhxIsLDw6Wo0e4w3BAREUnH4nAzb948NGnSBAEBAcjOzkazZs3w+OOPo1OnTpg6daoUNdqd0hv5MdwQERFZn8U38VOr1Vi9ejWmTZuGuLg4ZGdno3Xr1mjYsKEU9dklben8UnqGGyIiImuzONyUCgwMRGBgoDVrqTbYLUVERCQds8JNRESE2QdctGhRuYupLkq7pQrFbBjEYpmrISIisi9mhZs///zT5PXx48dRXFyMxo0bAyiZmVupVKJt27bWr9AOaRTOAAQAIgoM2XKXQ0REZFfMCjd79uwxPl+0aBFcXFywfv16eHh4AABu3bqFYcOGoWvXrtJUaWcUghIahTMKDLfZNUVERGRlFl8ttXDhQkRFRRmDDQB4eHhgzpw5WLhwoVWLs2ecX4qIiEgaFoebrKwsXL9+/b7l169fx+3bt61SVHVQGm7y2HJDRERkVRaHmxdffBHDhg3D5s2bceXKFVy5cgXff/89hg8fjpdeekmKGu0Sr5giIiKShsWXgq9cuRITJ07EwIEDUVRUVHIQlQrDhw/H/PnzrV6gvTKGG97rhoiIyKosDjeOjo745JNPMH/+fJw/fx4AUL9+fTg5OVm9OHtWGm4K2HJDRERkVeW+iZ+TkxNatGhhzVqqFXZLERERScPiMTdkHZxfioiISBoMNzJhyw0REZE0GG5kolXemTyT4YaIiMiqLA43OTk5UtRR7dzdLWUQDTJXQ0REZD8sDjfe3t547bXXcODAASnqqTZKu6VEGHC7kDc/JCIishaLw82XX36Jmzdv4sknn0SjRo3w/vvv49q1a1LUZteUggMcBB0A4Fb+LZmrISIish8Wh5u+ffti69atuHr1KkaOHIkNGzagTp06eP7557F582YUFxdLUaddKm29uVXAcENERGQt5R5QXLNmTURERODkyZNYtGgRdu/ejZdffhl+fn6YPn06cnNzrVmnXSoNNzfzb8pcCRERkf0o9038UlNTsX79eqxbtw6XLl3Cyy+/jOHDh+PKlSv44IMPcPjwYfzyyy/WrNXuGFtu2C1FRERkNRaHm82bN2Pt2rXYuXMnmjVrhrfffhuvvvoq3N3djdt06tQJTZs2tWaddonhhoiIyPosDjfDhg3DK6+8goMHD6J9+/ZlbuPn54d33323wsXZO3ZLERERWZ/F4SY5ORmOjo4P3Uan02HGjBnlLqq60CpLwk1GQYa8hRAREdkRiwcUu7i4IC0t7b7lN27cgFKptEpR1QW7pYiIiKzP4nAjimKZywsKCqBWqytcUHVSepdidksRERFZj9ndUh999BEAQBAEfPrpp3B2djau0+v12L9/P5o0aWL9Cu2YRlEyvxTvc0NERGQ9ZoebxYsXAyhpuVm5cqVJF5RarUZQUBBWrlxp/Qrt2N3dUqIoQhAEmSsiIiKq+swON4mJiQCA7t27Y/PmzfDw8JCsqOpCd6dbqkBfgLziPDg6PHygNhERET2axVdL7dmzR4o6qiWVoIUSDtCjCLcKbjHcEBERWYFZ4SYiIgKzZ8+Gk5MTIiIiHrrtokWLrFJYdSAIArQKV+QYbuBW/i34O/vLXRIREVGVZ1a4+fPPP1FUVGR8/iAcM2I5rbIk3PCKKSIiIuswK9zc3RXFbinr0vBeN0RERFZV7lnByTp4Iz8iIiLrMqvl5qWXXjL7gJs3by53MdVR6RVT6XnpMldCRERkH8wKN25ublLXUW05K2sCAK7lXJO5EiIiIvtgVrhZu3at1HVUWy4qbwDAldtXZK6EiIjIPnDMjcxclCXh5vLtyw+ct4uIiIjMZ1bLTZs2bRAdHQ0PDw+0bt36oZd8Hz9+3GrFVQel4Sa7KBtZhVlw07ALkIiIqCLMCjcvvPACNBoNAKBv375S1lPtqBQaeOm8kJ6Xjiu3rzDcEBERVZBZ4WbGjBllPifrqO1cG+l56bicfRnNvZrLXQ4REVGVZvHcUqWOHj2K+Ph4AECzZs3Qtm1bqxVV3dR2qY3Y67G4evuq3KUQERFVeRaHmytXrmDAgAE4ePAg3N3dAQAZGRno1KkTvvnmG9SuXdvaNdq92i4ln9mVbF4xRUREVFEWXy31+uuvo6ioCPHx8bh58yZu3ryJ+Ph4GAwGvP7661LUaPdqO98JN7wcnIiIqMIsbrnZt28ffv/9dzRu3Ni4rHHjxli2bBm6du1q1eKqC2PLDcMNERFRhVncchMQEGCcIfxuer0efn5+VimquiltuUnOSUaxoVjmaoiIiKo2i8PN/PnzMWbMGBw9etS47OjRoxg3bhwWLFhg1eKqi5qONaFWqKEX9UjJSZG7HCIioirNrG4pDw8Pkxv35eTkICQkBCpVye7FxcVQqVR47bXXeB+cclAICvi7+CMxMxFXsq8Yu6mIiIjIcmaFmyVLlkhcBtV2rl0Sbm5fAXzlroaIiKjqMivcDBkyROo6qj0OKiYiIrKOct/EDwDy8/NRWFhosszV1bVCBVVX/s7+AHivGyIiooqyeEBxTk4ORo8ejVq1asHJyQkeHh4mDyofttwQERFZh8Xh5j//+Q9+/fVXrFixAhqNBp9++ilmzpwJPz8/fP7551LUWC0Yb+THlhsiIqIKsbhb6v/+7//w+eefo1u3bhg2bBi6du2KBg0aoE6dOvjqq68waNAgKeq0e6UtN5kFmcgqzIKrmt17RERE5WFxy83NmzdRr149ACXja27evAkA6NKlC/bv32/d6qoRJwcneGo9AYATaBIREVWAxeGmXr16SExMBAA0adIEmzZtAlDSolM6kSaVD7umiIiIKs7icDNs2DCcOHECADBlyhQsX74cWq0WEyZMwKRJk6xeYHXi73LniikOKiYiIio3i8PNhAkTMHbsWABAjx49EB8fjw0bNuDPP//EuHHjylXE8uXLERQUBK1Wi5CQEMTExJi13zfffANBEOzmrsicHZyIiKjiKnSfGwAICgpCUFBQufffuHEjIiIisHLlSoSEhGDJkiUICwtDQkICatWq9cD9Ll68iIkTJ9rVTOQBLgEAgKvZHHNDRERUXha33ABAdHQ0nn/+edSvXx/169fH888/j927d5ergEWLFmHEiBEYNmwYmjVrhpUrV8LR0RGfffbZA/fR6/UYNGgQZs6caRzcbA+M97rhmBsiIqJyszjcfPLJJ+jVqxdcXFwwbtw4jBs3Dq6urnj22WexfPlyi45VWFiIY8eOoUePHv8UpFCgR48eOHTo0AP3mzVrFmrVqoXhw4c/8j0KCgqQlZVl8qisSrulrmZfhd6gl7kaIiKiqsnibql58+Zh8eLFGD16tHHZ2LFj0blzZ8ybNw+jRo0y+1jp6enQ6/Xw9vY2We7t7Y3Tp0+Xuc+BAwewZs0axMbGmvUeUVFRmDlzptk1yamWYy2oFCoUG4qRlpsGX2fOoElERGQpi1tuMjIy0KtXr/uWP/3008jMzLRKUQ9y+/ZtDB48GKtXr4aXl5dZ+0RGRiIzM9P4uHz5sqQ1VoRSoeQcU0RERBVkcctNnz59sGXLlvsu+/7hhx/w/PPPW3QsLy8vKJVKpKammixPTU2Fj4/PfdufP38eFy9eRO/evY3LDAYDAEClUiEhIQH169c32Uej0UCj0VhUl5xqO9fGpaxLuHL7Ctr7tJe7HCIioirHrHDz0UcfGZ83a9YMc+fOxd69exEaGgoAOHz4MA4ePIh33nnHojdXq9Vo27YtoqOjjZdzGwwGREdHm3R7lWrSpAlOnTplsmzq1Km4ffs2li5dioCAAIvevzIqHVR8+XblbWEiIiKqzMwKN4sXLzZ57eHhgb///ht///23cZm7uzs+++wzTJ061aICIiIiMGTIELRr1w4dOnTAkiVLkJOTg2HDhgEAwsPD4e/vj6ioKGi1Wjz22GMm+5feFfne5VUV71JMRERUMWaFm9LpFqTQv39/XL9+HdOnT0dKSgpatWqFHTt2GAcZJyUlQaEo1xXrVVJpyw3nlyIiIiqfCt3ETxRFAIAgCBUqYvTo0WV2QwHA3r17H7rvunXrKvTelQ3vdUNERFQx5WoS+fzzzxEcHAydTgedTocWLVrgiy++sHZt1VLp1VI3828iuzBb5mqIiIiqHovDzaJFi/DWW2/h2WefxaZNm7Bp0yb06tULI0eOvG9sDlnORe2CWrqSaSfOZZyTuRoiIqKqx+JuqWXLlmHFihUIDw83LuvTpw+aN2+O9957DxMmTLBqgdVRY8/GSLuahvib8WhVq5Xc5RAREVUpFrfcJCcno1OnTvct79SpE5KTk61SVHXXxLMJAOD0zbLv0kxEREQPZnG4adCgATZt2nTf8o0bN6Jhw4ZWKaq6a1qjKQAg/ka8zJUQERFVPRZ3S82cORP9+/fH/v370blzZwDAwYMHER0dXWboIcuVttycyziHIkMRHBQOMldERERUdVjccvOvf/0LMTEx8PLywtatW7F161Z4eXkhJiYGL774ohQ1Vjv+zv5wdnBGkaEIFzIuyF0OERFRlWJRy01RURHefPNNTJs2DV9++aVUNVV7CkEBZyEQ2fgbS3/bi4aOJfcRmtCzkbyFERERVQEWtdw4ODjg+++/l6oWuounQ10AwI0i6e4OTUREZI8s7pbq27cvtm7dKkEpdLcad8LNTYYbIiIii1g8oLhhw4aYNWsWDh48iLZt28LJyclk/dixY61WXHVWw9hycxGiaIAgVJ/5tYiIiCrC4nCzZs0auLu749ixYzh27JjJOkEQGG6sxF1VGwqoUCTm4rY+Da4qH7lLIiIiqhIsDjdSzhBO/1AIKng4BOJG0QXcLEpkuCEiIjJThfo6RFE0zgxO1leDg4qJiIgsVq5ws2bNGjz22GPQarXQarV47LHH8Omnn1q7tmqPV0wRERFZzuJuqenTp2PRokUYM2YMQkNDAQCHDh3ChAkTkJSUhFmzZlm9yOqKLTdERESWszjcrFixAqtXr8aAAQOMy/r06YMWLVpgzJgxDDdW5KmqA0BAnuEW8vQZcpdDRERUJVjcLVVUVIR27drdt7xt27YoLi62SlFUwkGhg6vSFwBbb4iIiMxlccvN4MGDsWLFCixatMhk+apVqzBo0CCrFVadLN515oHrajjURZb+GsMNERGRmSwON0DJgOJffvkFHTt2BAAcOXIESUlJCA8PR0REhHG7ewMQWa6GQ10k5h/knYqJiIjMZHG4iYuLQ5s2bQAA58+fBwB4eXnBy8sLcXFxxu0EQbBSidUbBxUTERFZxuJws2fPHinqoAcovRw8S5+CnKIcODk4PWIPIiKi6q1c3VJkOzqlGxwVnsg13MTsndHw1jQxrpvQs5GMlREREVVOnI2xCjB2TRVfkLkSIiKiyo/hpgow3qm4kOGGiIjoURhuqoBa6sYAgOTCv2SuhIiIqPJjuKkCfNRNIUCBbH0abhenyV0OERFRpcZwUwU4KHTwcmgAAEgujHvE1kRERNUbw00V4at5DACQXMBwQ0RE9DAMN1XEP+HmFERRlLkaIiKiyovhporwdmgCBVTINdxElj5Z7nKIiIgqLYabKkKl0KCWuuSmfeyaIiIiejCGmyrER81xN0RERI/CcFOF+N0Zd5NS+BfH3RARET0Aw00VUlPdCEqokWfIQEbxFbnLISIiqpQYbqoQpeBw192KT8lcDRERUeXEcFPF+GmCAXDcDRER0YMw3FQxpfe7SSn4CwbRIHM1RERElQ/DTRXj5VAfKkGLAjEbZ26dkbscIiKiSofhpopRCCr4qJsCAGKSY2SuhoiIqPJhuKmCSrumYlIYboiIiO7FcFMF+apLBhUfSz2GYkOxzNUQERFVLgw3VZCnQxDUghOyi7JxKp2XhBMREd2N4aYKUghK1Na2AQDsvrRb5mqIiIgqF4abKipIGwoA2HVpF6diICIiugvDTRVVW9sKOpUOyTnJ+OvGX3KXQ0REVGkw3FRRKkGDx2s/DgD45dIvMldDRERUeTDcVGE96/QEAOy6yK4pIiKiUgw3VVhX/67QKrW4kn0FCbcS5C6HiIioUmC4qcIcHRzRxb8LAOCXi+yaIiIiAhhuqrTFu86gKKvkbsUb/96GRb+w9YaIiIjhpooL0LaFAipk6a8ho/iy3OUQERHJjuGmilMrHOGvaQUASMw7JG8xRERElQDDjR0I0pXc0O9S/mGZKyEiIpIfw40dCNS2gwIq3CpOwoXMC3KXQ0REJCuGGzugUTjDT1MyUzjnmiIiouqO4cZOlHZN/XThJ97Qj4iIqjWGGzsRpA2FStAiMTMRR1OPyl0OERGRbBhu7IRa4Yj6uq4AgG8TvpW5GiIiIvkw3NiRJk5PAwB2Je3CjbwbMldDREQkD4YbO1LDoR6CvYJRbCjG1nNb5S6HiIhIFgw3dubfjf8NAPj2zLcwiAaZqyEiIrI9hhs7c/Z8fagFJ1zNvopJP3HsDRERVT+VItwsX74cQUFB0Gq1CAkJQUxMzAO3Xb16Nbp27QoPDw94eHigR48eD92+ulEpNGjg2A0AcDp3p7zFEBERyUD2cLNx40ZERERgxowZOH78OFq2bImwsDCkpaWVuf3evXsxYMAA7NmzB4cOHUJAQACefvppXL161caVV15NHEsGFl/OP4aUnBSZqyEiIrIt2cPNokWLMGLECAwbNgzNmjXDypUr4ejoiM8++6zM7b/66iu8/fbbaNWqFZo0aYJPP/0UBoMB0dHRNq688nJ3qA0fdXOIMGDz2c1yl0NERGRTsoabwsJCHDt2DD169DAuUygU6NGjBw4dMm+G69zcXBQVFcHT07PM9QUFBcjKyjJ5VAell4V/f+Z7FBmKZK6GiIjIdmQNN+np6dDr9fD29jZZ7u3tjZQU87pTJk+eDD8/P5OAdLeoqCi4ubkZHwEBARWuuyqoow2BVuGKtLw0zjdFRETViuzdUhXx/vvv45tvvsGWLVug1WrL3CYyMhKZmZnGx+XLl21cpTyUggOaOj0DAFh1chUvCyciompD1nDj5eUFpVKJ1NRUk+Wpqanw8fF56L4LFizA+++/j19++QUtWrR44HYajQaurq4mj+qimdNzcHFwwbmMc2y9ISKiakPWcKNWq9G2bVuTwcClg4NDQ0MfuN+HH36I2bNnY8eOHWjXrp0tSq2SNAonDGo2CACw8uRKtt4QEVG1IHu3VEREBFavXo3169cjPj4eb731FnJycjBs2DAAQHh4OCIjI43bf/DBB5g2bRo+++wzBAUFISUlBSkpKcjOzpbrFCq1V5u+CicHJ5y9dRa/Jv0qdzlERESSkz3c9O/fHwsWLMD06dPRqlUrxMbGYseOHcZBxklJSUhOTjZuv2LFChQWFuLll1+Gr6+v8bFgwQK5TqFSc9O4YWCTgQCAlSdWQhRFmSsiIiKSliBWs2+7rKwsuLm5ITMzU5LxN4t3nbH6MStiQs9GyMjPQNj3YcgtzsXS7kvxZOCTcpdFRERkEUu+v2VvuSHpuWvdMbApW2+IiKh6UMldAEmrtCUpX98ZKuFLxN+Mx/4r+/FEwBMyV0ZERCQNttxUE1qlK5o69QIALI9dDr1BL3NFRERE0mC4qUYec+oDB8ER8Tfjsfkc55wiIiL7xHBTjeiUbmjj0h8A8NHxj5BZkClzRURERNbHcFPNNHXqhQbuDZBRkIFlfy6TuxwiIiKrY7ipZhSCCpEdSm6K+O2Zb3H65mmZKyIiIrIuhptq6GCcO+pqO8EgGjB653ReGk5ERHaF4aaaau82BCpBg9TCeGxL3CZ3OURERFbDcFNNOSu90NL5XwCAhUcXIruQc3MREZF9YLipxh5z7gNXpQ/S89Kx6NgiucshIiKyCoabakwpOKCT+0gAJYOL913eJ3NFREREFcdwU835aYIxuNlgAMD036fjRt4NmSsiIiKqGIYbwrg249DAvQFu5t/Ee4fe49VTRERUpTHcED759RKaO7wFBVTYe3kv3vphhdwlERERlRvDDQEAajgEoa3rAADAkay1SMpKkrkiIiKi8mG4IaPmTr3ho26OYjEfkb9FolBfKHdJREREFmO4ISOFoMTj7qOhFhxxMv0k5h6Zy/E3RERU5TDckAlnVS1084iAQlBg89nN2HB6g9wlERERWYThhu5TW9sabV1eBQB8EPMhJv30ncwVERERmY/hhsr0mFMf1Nc9DhEG7Lm5EJdvX5a7JCIiIrMw3FCZBEFAZ/eR8HJogAIxG2N/HYucohy5yyIiInokldwFUOWlEjR4yvM/+PH6f3Au4xxe+u4N9KgRCZWgBgBM6NlI5gqJiIjux5YbeignZQ308JwClaDFtcKT2HtrEQxisdxlERERPRDDDT1STXVD9PCcAiUckJT/B37LWA5RNMhdFhERUZnYLUVm8dMEo7vnO4i+OR/n8/ZDJWgxQVwAQRDkLo2Iqoo9Uaavu0fKUwfZPYYbMlugtj2e8BiLvbeWICH3Fyw8uhDvtHuHAYdIalKFgnuPa+6xy7tfVVDWud3LXs7VjjHckEXq6bqgyJCPg5krsP7v9cjX5yOyQySUCqXcpRERVS32HBJlxnBDFmvs1AMi9DiUuRobEzYioyADUV2i4KB0kLs0IqoszGkBMWe/6vZlX93P30oYbqhcmjiF4cVWDRH5WyR2XtyJ24W3sbjbYjg6OMpdGlHlYE//V17eoEIkE4YbKrdeQb3gqnbF+D3j8fu13zHilxFY/tRyuGvd5S6NqPqpTv/HX9nGxdhTkLUTDDdUbot3nQHghafcp2PXjbk4mX4SA7YNwNInl6KRB2/wR3asKnyZVffWlqrwO7Ilc/893PsZVdHPkeGGKqyWuhGe85qLXTfn4Ur2Fby6/VXM7jwbYUFhcpdGVLlU0S8KshJrBc7q1EpXTgw3ZBXuDrXRp+YH2HNrMa4VnMDEfRPx+bHf8fmL03klFVUttg4g/KJ6OClboKw16Lkyqgo1SojhhqxGo3DB057v4ljWVziV8wNOZm/GqF9TMbfzXNTQ1ZC7PCJpVfMvE7qHPf97qAKBnOGGrEohKNHeLRyeDnVxIPMTHLx6EP/68V+Y1XkWHq/9uNzlUXVTBf4I38eevxTvZetzrYqfbVWsuRJguCFJ1HfsCg+HAOy7tRQ38pMwKnoUmjr2Qnu3cEx6Olju8qgyMSeAWCukcMwLUbXAcEOS8XQIQu+aH+Bo1pf4O2cb4nN3ILkwDs+mL0Bzr+Zyl0dEZWFLQfVmJ79/hhuSlEpQo6PbawjQtMH+jI+RUXwFA7cPxMAmAzG69Wg4OTjJXSJVNub8cbVmC0x534+IKi2GG7IJf20rvFhzEQ5nrcGFvAP4Mv5L/HLpF/y3w3/xZOCTnHyTKo4BRF78/KuvStjdy3BDNqNVuqKbxwQ01HXH75mrkZabgvF7x6Nb7W6Y2H4i6rjWkbvE6knKQbdVcUAvEVV5DDdkc/7aVnhRswgnbn+PU9k/YO+VvThw9QD6Ne6HkS1HwlPrKXeJ9qOyBRf+37118HOku/Hfw30UchdA1ZNK0KCt60D0rbkQj9d+HMViMb4+/TWe2/wcPj31KfKL8+UukYiIqii23JCs3B1qw90wDo41nsQfWetxoygRS48vxVfxX2Fo86Ho16gfZxqvKvh/j0RUSQiiKIpyF2FLWVlZcHNzQ2ZmJlxdXa1+/JLJJKk8RNGA83m/4fjtr5Gtvw4A8NR6IrxZOF5p8gqvrHoUhgsiqiwkGF9nyfc3W26o0hAEBRo4PoG6uk44l7sPJ7M342Z+KpYcX4K1f63Fvxv9G/0b94e3k7fcpdpeJbwagYiosmK4oUpHKTigsVMPNHTsjvN5v+HE7e+RWXANq0+txqenPkNdXSiaOz2Pec89J3epRERUCTHcUKWlEJRo6NgN9XVdkZR/FH/n/ISUwr9xIe8ALuQdQOJPX+Jfjf6FXkG94Kx2lrtc6+KN5YiIyo3hhio9haBEkC4EQboQ3Ci6gL+yt+FC3gHE3YhD3KE4fPjHhwgLCsNLDV9Cq5qtKvcNARlIiIgkx3BDVUoNh3p43GMM2ruG41zePpzJ3Y3M4qvYem4rtp7bitrOtfFM3WfwbN1n0cCjgW2L47gYIqJKgVdLWRmvlrItURSRVpiAM7m7kZj/O4rFAuM6D1UgXg1+AU8FPoV6bvUq1qLDFhciIvPJfLUUw42VMdzIp8iQj8sFR3Eh7wCu5P8JA4qN61yVfqij64BAbXvMfeY5KBXKf3ZkcCEisi5eCk5kHQ4KLerpuqCergsKDNm4lHcYGbe+xWnhBrL013AqeytOZW/FgU0fopNvJ3T274xOfp1QU+7CiYjIqhhuyC5pFM5o5NQDHW9cQD6KESekI1ZIwynhOjILMvHzxZ/x88WfAQCNHNzRQeuN9hpvtNXWgptCLXP1RERUEQw3ZPe0UKGd6IN2og/0MCARmfhLkY44IR1JyMKZogycKcrAl7cTIABo4uCBNpqaaHXn4aPi9A9ERFUJww1VWR2TVlm8jxIKNIAHGhg88AIa4jYKYfAuxB8FqYjJT8PF4izEF91CfNEtfJVdMn7KR+mIlhovBKtroJnaE83UnnBSOFj7dIiIyEoYbqhKKE+QMYcL1ECqGj3gjB6ojwzk44xwCzlu+YgtuI6Eogyk6HORkpuEnblJAABBBHzghLZONdFY7VHycHCHp1IrSY1ERGQZhhuSVFmh5HDgG+XazxbcoUUH0RfIALqjLvJRjItCJi4gE5eETFwUsnBLyEcycvBTbg5+yr1o3LemUof6Dm5o4OBm/FnPwQ2uHMNDRGRTDDckO7mCjDm0UKGJWANNUAO4c9OETBTgkpCJy7iNy8JtXBFuI03IxXV9Hq7r83A4P8XkGJ4KLYIcXFBX5YogB1cEqJwRoHJBbZUzHBX8T5CIyNr4l5VsrjKHGXO4QYMWYi20QC1j4MlHMa4hG9eEO487zzOEAtw05ONmQT6OF1y/71iuohr1NK7wUznBX+UMP5UT/JRO8FE5wkfpxPBDRFQO/MtJ5VbVQ4o1aaFCPbijnuhuDDwAkIdipCEHKUIOUoVcpCAH14VcXEcucoViZAmFiC1MR2xhepnHdVOo4aN0hLfSEbVUjqil1KGmUmf86aXUwUOhgUpQ2OZEiYiqAIYbKhODi3XooEIduKGO6GYSegAgB0W4jlykC3m4gTzcEPKQfufnLeQjX9Aj01CITEMhEooyHvgeAgAPhRY1lFp4KjTwvOunh1ILD4UGbgo1PBQauCtLnjsIygcej4ioqmO4sXMMKZWXExzgBDcEiW4lC+4JP7kowi3k46aQjwzkIwMFyBAK7vzMR56yCDcNBTBALOn6MuSb/d6Oggrud0KPq0INF+NPB7gq1HBWOMBFoYaz4ABnxZ2H4ADHO8/VUFTu2deJqFqrFOFm+fLlmD9/PlJSUtCyZUssW7YMHTp0eOD23377LaZNm4aLFy+iYcOG+OCDD/Dss8/asOLKgcHFvjnCAY5wgL/o8s/CuwOQHjBARDYKkYVCZAoFxue3hZKfOShEtlCEbBQiG0XIRRFEAcgVi5GrL8Y1fU65alNBAUeFCk6CCo4KBzgJKugEFXQKFRzvPHdUqKAV7rxWqKAVlNAKpT//ea4RlNDcWa++s47hiYgqQvZws3HjRkRERGDlypUICQnBkiVLEBYWhoSEBNSqVeu+7X///XcMGDAAUVFReP7557Fhwwb07dsXx48fx2OPPSbDGViOoYSsRQEBrtDAFRrUflAIuosBIvJQhJw7j2yhJPDkoRg5KEKeUIxcFCEXxchDMfKEkp/5dx4Fgh4AUAwDsgwlAQp6ac5NDQU0dwJPyU8FHAQlNIICakEJB9z5KShK1qFkvYOgKHlAAdVdzx2Ekteq0ucQSl7fWaYSBKiggPKunw5QQCkooIRgslxp/HnnISiguLNN6TKGMyL5yD4reEhICNq3b4+PP/4YAGAwGBAQEIAxY8ZgypQp923fv39/5OTk4KeffjIu69ixI1q1aoWVK1c+8v2knhX80JqJVj8mUWVhgPhP0IHeGHjyoUfBnWXGh1CMQhhQeOd1YelDKFlWBL1xffGdn6Id5QEBgBKCMfQojM8VJsuUECAAUAilz0vWCRCgQEmAVdy1v0K4s/2d16X7KlDy4Snu7Cfccwzhrv1Ktyvdt6TeO/vdWY4764S7lv3zE8bj3bvOeP7CXdui9D3ufQ5jCLz3GGUd7+79cM+xcM+y+7YVTJffu/29x7x3u7uVdax7tzXnOPdu96Bt7zu2GcHZnP+U7n1v8/d7NPc2QxHqF2rGluarMrOCFxYW4tixY4iM/GdqdIVCgR49euDQoUNl7nPo0CFERESYLAsLC8PWrVvL3L6goAAFBQXG15mZmQBKPiQp5OQVPHojoipOAwU0UMAV1puGQoQIPcQ7YUdEMQzG0FMMA/QQUSQYUAw9iiBCD4PJdvo72xQLhjvrSl6XPv55bYDhzk89ROgF8b7lBgB6iHdel/y8+1G67aPCWLHVPh2iqiU4PwPNn25u1WOWfm+b0yYja7hJT0+HXq+Ht7e3yXJvb2+cPn26zH1SUlLK3D4lJaXM7aOiojBz5sz7lgcEBJSzaiIiInqYeMRjEzZJcuzbt2/Dzc3todvIPuZGapGRkSYtPQaDATdv3kSNGjWs3ieelZWFgIAAXL58WZIur8qG52vfeL72rbqdL1D9ztnezlcURdy+fRt+fn6P3FbWcOPl5QWlUonU1FST5ampqfDx8SlzHx8fH4u212g00Gg0Jsvc3d3LX7QZXF1d7eIfkrl4vvaN52vfqtv5AtXvnO3pfB/VYlNK1tuaqtVqtG3bFtHR0cZlBoMB0dHRCA0teyBSaGioyfYAsGvXrgduT0RERNWL7N1SERERGDJkCNq1a4cOHTpgyZIlyMnJwbBhwwAA4eHh8Pf3R1RUFABg3LhxeOKJJ7Bw4UI899xz+Oabb3D06FGsWsXLq4mIiKgShJv+/fvj+vXrmD59OlJSUtCqVSvs2LHDOGg4KSkJCsU/DUydOnXChg0bMHXqVPz3v/9Fw4YNsXXr1kpxjxuNRoMZM2bc1w1mr3i+9o3na9+q2/kC1e+cq9v53k32+9wQERERWROnEiYiIiK7wnBDREREdoXhhoiIiOwKww0RERHZFYYbK1m+fDmCgoKg1WoREhKCmJgYuUuSTFRUFNq3bw8XFxfUqlULffv2RUJCgtxl2cT7778PQRAwfvx4uUuR1NWrV/Hqq6+iRo0a0Ol0CA4OxtGjR+UuSxJ6vR7Tpk1D3bp1odPpUL9+fcyePdus+Wuqgv3796N3797w8/ODIAj3zcMniiKmT58OX19f6HQ69OjRA2fPnpWnWCt42PkWFRVh8uTJCA4OhpOTE/z8/BAeHo5r167JV3AFPer3e7eRI0dCEAQsWbLEZvXJheHGCjZu3IiIiAjMmDEDx48fR8uWLREWFoa0tDS5S5PEvn37MGrUKBw+fBi7du1CUVERnn76aeTk5MhdmqT++OMP/O9//0OLFi3kLkVSt27dQufOneHg4ICff/4Zf//9NxYuXAgPDw+5S5PEBx98gBUrVuDjjz9GfHw8PvjgA3z44YdYtmyZ3KVZRU5ODlq2bInly5eXuf7DDz/ERx99hJUrV+LIkSNwcnJCWFgY8vPzbVypdTzsfHNzc3H8+HFMmzYNx48fx+bNm5GQkIA+ffrIUKl1POr3W2rLli04fPiwWVMX2AWRKqxDhw7iqFGjjK/1er3o5+cnRkVFyViV7aSlpYkAxH379sldimRu374tNmzYUNy1a5f4xBNPiOPGjZO7JMlMnjxZ7NKli9xl2Mxzzz0nvvbaaybLXnrpJXHQoEEyVSQdAOKWLVuMrw0Gg+jj4yPOnz/fuCwjI0PUaDTi119/LUOF1nXv+ZYlJiZGBCBeunTJNkVJ6EHne+XKFdHf31+Mi4sT69SpIy5evNjmtdkaW24qqLCwEMeOHUOPHj2MyxQKBXr06IFDhw7JWJntZGZmAgA8PT1lrkQ6o0aNwnPPPWfye7ZXP/74I9q1a4d+/fqhVq1aaN26NVavXi13WZLp1KkToqOjcebMGQDAiRMncODAATzzzDMyVya9xMREpKSkmPy7dnNzQ0hISLX6+yUIguRzDsrFYDBg8ODBmDRpEpo3by53OTYj+x2Kq7r09HTo9XrjHZVLeXt74/Tp0zJVZTsGgwHjx49H586dK8VdoqXwzTff4Pjx4/jjjz/kLsUmLly4gBUrViAiIgL//e9/8ccff2Ds2LFQq9UYMmSI3OVZ3ZQpU5CVlYUmTZpAqVRCr9dj7ty5GDRokNylSS4lJQUAyvz7VbrOnuXn52Py5MkYMGCA3Uwsea8PPvgAKpUKY8eOlbsUm2K4oQoZNWoU4uLicODAAblLkcTly5cxbtw47Nq1C1qtVu5ybMJgMKBdu3aYN28eAKB169aIi4vDypUr7TLcbNq0CV999RU2bNiA5s2bIzY2FuPHj4efn59dni+VKCoqwr///W+IoogVK1bIXY4kjh07hqVLl+L48eMQBEHucmyK3VIV5OXlBaVSidTUVJPlqamp8PHxkakq2xg9ejR++ukn7NmzB7Vr15a7HEkcO3YMaWlpaNOmDVQqFVQqFfbt24ePPvoIKpUKer1e7hKtztfXF82aNTNZ1rRpUyQlJclUkbQmTZqEKVOm4JVXXkFwcDAGDx6MCRMmGCfrtWelf6Oq29+v0mBz6dIl7Nq1y25bbX777TekpaUhMDDQ+Pfr0qVLeOeddxAUFCR3eZJiuKkgtVqNtm3bIjo62rjMYDAgOjoaoaGhMlYmHVEUMXr0aGzZsgW//vor6tatK3dJknnqqadw6tQpxMbGGh/t2rXDoEGDEBsbC6VSKXeJVte5c+f7Lu0/c+YM6tSpI1NF0srNzTWZnBcAlEolDAaDTBXZTt26deHj42Py9ysrKwtHjhyx279fpcHm7Nmz2L17N2rUqCF3SZIZPHgwTp48afL3y8/PD5MmTcLOnTvlLk9S7JaygoiICAwZMgTt2rVDhw4dsGTJEuTk5GDYsGFylyaJUaNGYcOGDfjhhx/g4uJi7Jt3c3ODTqeTuTrrcnFxuW8skZOTE2rUqGG3Y4wmTJiATp06Yd68efj3v/+NmJgYrFq1CqtWrZK7NEn07t0bc+fORWBgIJo3b44///wTixYtwmuvvSZ3aVaRnZ2Nc+fOGV8nJiYiNjYWnp6eCAwMxPjx4zFnzhw0bNgQdevWxbRp0+Dn54e+ffvKV3QFPOx8fX198fLLL+P48eP46aefoNfrjX+/PD09oVar5Sq73B71+703vDk4OMDHxweNGze2dam2JfflWvZi2bJlYmBgoKhWq8UOHTqIhw8flrskyQAo87F27Vq5S7MJe78UXBRF8f/+7//Exx57TNRoNGKTJk3EVatWyV2SZLKyssRx48aJgYGBolarFevVqye+++67YkFBgdylWcWePXvK/O91yJAhoiiWXA4+bdo00dvbW9RoNOJTTz0lJiQkyFt0BTzsfBMTEx/492vPnj1yl14uj/r93qu6XAouiKKd3IaTiIiICBxzQ0RERHaG4YaIiIjsCsMNERER2RWGGyIiIrIrDDdERERkVxhuiIiIyK4w3BAREZFdYbghkkC3bt0wfvx44+ugoCAsWbLE+FoQBGzdurXcx3/Y8S5evAhBEBAbG1vu45fXunXr4O7ubtP3vPezprLt3bsXgiAgIyOj3MfgZ01VBcMNkQ388ccfeOONNyQ7fnJyMp555hnJjl8ZPOjLefPmzZg9e3aFjl3RsFmWewOoLY9RVgjp1KkTkpOT4ebm9sj9pfysiWyBc0sR2UDNmjUlPX5VnsG5sLCwQnP6eHp6WrEa+6VWqyv874SfNVUVbLkhsoFH/R/4jBkz4Ovri5MnTwIADhw4gK5du0Kn0yEgIABjx45FTk7OA/cvq+XhwoUL6N69OxwdHdGyZUscOnTIZP3333+P5s2bQ6PRICgoCAsXLjRZf+vWLYSHh8PDwwOOjo545plncPbsWZNt1q1bh8DAQDg6OuLFF1/EjRs3HvlZDB06FH379sXcuXPh5+dnnMDviy++QLt27eDi4gIfHx8MHDgQaWlpAEq62rp37w4A8PDwgCAIGDp0KID7WynMqftuQUFBAIAXX3wRgiAYXwPADz/8gDZt2kCr1aJevXqYOXMmiouLAQCiKOK9995DYGAgNBoN/Pz8MHbsWGNNly5dwoQJEyAIAgRBKPO9y3OMGzduYMCAAfD394ejoyOCg4Px9ddfm3y++/btw9KlS437Xbx48b7WmEuXLqF3797w8PCAk5MTmjdvju3bt1v0WRcUFGDy5MkICAiARqNBgwYNsGbNmgd+1kQ2I+/UVkT26d7JNe+drA6AuGXLFtFgMIijR48Wg4KCxLNnz4qiKIrnzp0TnZycxMWLF4tnzpwRDx48KLZu3VocOnToI48niqJxcsAmTZqIP/30k5iQkCC+/PLLYp06dcSioiJRFEXx6NGjokKhEGfNmiUmJCSIa9euFXU6ncnkp3369BGbNm0q7t+/X4yNjRXDwsLEBg0aiIWFhaIoiuLhw4dFhUIhfvDBB2JCQoK4dOlS0d3dXXRzc3voZzNkyBDR2dlZHDx4sBgXFyfGxcWJoiiKa9asEbdv3y6eP39ePHTokBgaGio+88wzoiiKYnFxsfj999+LAMSEhAQxOTlZzMjIKPOzflTd90pLSzNO/JqcnCympaWJoiiK+/fvF11dXcV169aJ58+fF3/55RcxKChIfO+990RRFMVvv/1WdHV1Fbdv3y5eunRJPHLkiHGC0Rs3boi1a9cWZ82aJSYnJ4vJycllvnd5jnHlyhVx/vz54p9//imeP39e/Oijj0SlUikeOXJEFEVRzMjIEENDQ8URI0YY9ysuLjZOsHjr1i1RFEXxueeeE3v27CmePHlSPH/+vPh///d/4r59+yz6rP/973+LAQEB4ubNm8Xz58+Lu3fvFr/55puH/v6JbIHhhkgC5oSbb7/9Vhw4cKDYtGlT8cqVK8Z1w4cPF9944w2T4/3222+iQqEQ8/LyHni8e8PNp59+alz/119/iQDE+Ph4URRFceDAgWLPnj1N3mPSpElis2bNRFEUxTNnzogAxIMHDxrXp6enizqdTty0aZMoiqI4YMAA8dlnnzU5Rv/+/c0KN97e3o+cdfuPP/4QAYi3b98WRVG878u51N2ftTl1l+Xuz6/UU089Jc6bN89k2RdffCH6+vqKoiiKCxcuFBs1avTA0GTO7MvWOIYolgSVd955x/i6rJnr7/38goODjUHtXuZ81gkJCSIAcdeuXY+sj8jW2C1FJJMJEybgyJEj2L9/P/z9/Y3LT5w4gXXr1sHZ2dn4CAsLg8FgQGJiotnHb9GihfG5r68vABi7eeLj49G5c2eT7Tt37oyzZ89Cr9cjPj4eKpUKISEhxvU1atRA48aNER8fbzzG3esBIDQ01Pg8KSnJ5BzmzZtnXBccHHzfOJtjx46hd+/eCAwMhIuLC5544gnjccxlTt3mOnHiBGbNmmVyDiNGjEBycjJyc3PRr18/5OXloV69ehgxYgS2bNli7LIyV3mOodfrMXv2bAQHB8PT0xPOzs7YuXOnRZ8TAIwdOxZz5sxB586dMWPGDGOXqLliY2OhVCqNvyeiyoThhkgmPXv2xNWrV7Fz506T5dnZ2XjzzTcRGxtrfJw4cQJnz55F/fr1zT6+g4OD8XnpeA2DwWCd4s3g5+dncg4jR440rnNycjLZNicnB2FhYXB1dcVXX32FP/74A1u2bAFQMuBYDtnZ2Zg5c6bJOZw6dQpnz56FVqtFQEAAEhIS8Mknn0Cn0+Htt9/G448/jqKiIrPfozzHmD9/PpYuXYrJkydjz549iI2NRVhYmMWf0+uvv44LFy5g8ODBOHXqFNq1a4dly5aZvb9Op7Po/YhsiVdLEcmkT58+6N27NwYOHAilUolXXnkFANCmTRv8/fffaNCggWTv3bRpUxw8eNBk2cGDB9GoUSMolUo0bdoUxcXFOHLkCDp16gSgZCBrQkICmjVrZjzGkSNHTI5x+PBh43OVSmX2OZw+fRo3btzA+++/j4CAAADA0aNHTbYpbenR6/UPPa9H1V0WBweH+47bpk0bJCQkPPQcdDodevfujd69e2PUqFFo0qQJTp06hTZt2kCtVj+01vIe4+DBg3jhhRfw6quvAigJrGfOnDE5P3PfOyAgACNHjsTIkSMRGRmJ1atXY8yYMWZ91sHBwTAYDNi3bx969OjxyPcisiW23BDJ6MUXX8QXX3yBYcOG4bvvvgMATJ48Gb///jtGjx6N2NhYnD17Fj/88ANGjx5ttfd95513EB0djdmzZ+PMmTNYv349Pv74Y0ycOBEA0LBhQ7zwwgsYMWIEDhw4gBMnTuDVV1+Fv78/XnjhBQAl3Ro7duzAggULcPbsWXz88cfYsWNHueoJDAyEWq3GsmXLcOHCBfz444/33U+lTp06EAQBP/30E65fv47s7Oz7jmNO3WUJCgpCdHQ0UlJScOvWLQDA9OnT8fnnn2PmzJn466+/EB8fj2+++QZTp04FUHKl2Jo1axAXF4cLFy7gyy+/hE6nQ506dYzH3L9/P65evYr09PQy37c8x2jYsCF27dqF33//HfHx8XjzzTeRmpp63/kcOXIEFy9eRHp6epktduPHj8fOnTuRmJiI48ePY8+ePWjatKnZn3VQUBCGDBmC1157DVu3bkViYiL27t2LTZs2PfBzJrIZuQf9ENkjc6+WKrVx40ZRq9WK33//vSiKohgTEyP27NlTdHZ2Fp2cnMQWLVqIc+fONet4pQOK//zzT+P6W7duiQDEPXv2GJd99913YrNmzUQHBwcxMDBQnD9/vsk53Lx5Uxw8eLDo5uYm6nQ6MSwsTDxz5ozJNmvWrBFr164t6nQ6sXfv3uKCBQvMGlD8wgsv3Ld8w4YNYlBQkKjRaMTQ0FDxxx9/vO88Zs2aJfr4+IiCIIhDhgwRRfH+z9qcuu/1448/ig0aNBBVKpVYp04d4/IdO3aInTp1EnU6nejq6ip26NDBeDXTli1bxJCQENHV1VV0cnISO3bsKO7evdu476FDh8QWLVqIGo1GfNCf2vIc48aNG+ILL7wgOjs7i7Vq1RKnTp0qhoeHm3ymCQkJYseOHUWdTicCEBMTE+8bJDx69Gixfv36okajEWvWrCkOHjxYTE9Pt+izzsvLEydMmCD6+vqKarVabNCggfjZZ5899LMmsgVBFEVRtmRFREREZGXsliIiIiK7wnBDREREdoXhhoiIiOwKww0RERHZFYYbIiIisisMN0RERGRXGG6IiIjIrjDcEBERkV1huCEiIiK7wnBDREREdoXhhoiIiOwKww0RERHZlf8H86XQ0Bv7omwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np \n",
    "\n",
    "# a single experiment consists of tossing a p1-coin n times, \n",
    "# and then tossing a p2-coin n times \n",
    "def toss_coins(n, p1, p2, number_experiments):\n",
    "    number_heads_1 = np.random.binomial(n, p1, number_experiments) \n",
    "    number_heads_2 = np.random.binomial(n, p2, number_experiments) \n",
    "    return number_heads_1, number_heads_2\n",
    "\n",
    "def compute_maximum_likelihood_estimates(number_heads, number_heads_1, number_heads_2, n):\n",
    "    p_null = number_heads / (2*n) \n",
    "    p1_alternative = number_heads_1 / n\n",
    "    p2_alternative = number_heads_2 / n\n",
    "    return p_null, p1_alternative, p2_alternative\n",
    "\n",
    "def compute_maximum_log_likelihood(number_heads, number_tosses, p_hat):\n",
    "    number_tails = number_tosses - number_heads\n",
    "    q_hat = 1 - p_hat\n",
    "    return number_heads*np.log(p_hat) + number_tails*np.log(q_hat)\n",
    "\n",
    "def compute_test_statistic(n, p1, p2, number_experiments): \n",
    "    number_heads_1, number_heads_2 = toss_coins(n, p1, p2, number_experiments)\n",
    "    number_heads = number_heads_1 + number_heads_2 \n",
    "\n",
    "    p_null, p1_alternative, p2_alternative = compute_maximum_likelihood_estimates(\n",
    "        number_heads, number_heads_1, number_heads_2, n\n",
    "    )\n",
    "    \n",
    "    maximum_log_likelihood_null = compute_maximum_log_likelihood(\n",
    "        number_heads=number_heads, \n",
    "        number_tosses=2*n, \n",
    "        p_hat=p_null\n",
    "    )\n",
    "    \n",
    "    maximum_log_likelihood1_alternative = compute_maximum_log_likelihood(\n",
    "        number_heads=number_heads_1, \n",
    "        number_tosses=n, \n",
    "        p_hat=p1_alternative\n",
    "    )\n",
    "    maximum_log_likelihood2_alternative = compute_maximum_log_likelihood(\n",
    "        number_heads=number_heads_2, \n",
    "        number_tosses=n, \n",
    "        p_hat=p2_alternative\n",
    "    )\n",
    "    maximum_log_likelihood_alternative = maximum_log_likelihood1_alternative + maximum_log_likelihood2_alternative\n",
    "    \n",
    "    return 2 * (maximum_log_likelihood_alternative - maximum_log_likelihood_null)\n",
    "    \n",
    "import matplotlib.pyplot as plt \n",
    "from scipy.stats import chi2\n",
    "\n",
    "def plot_test_statistic_distributions(n=500, p=0.25, p1=0.2, p2=0.3, number_experiments=100000):\n",
    "    test_statistic_under_model1 = compute_test_statistic(n, p, p, number_experiments)\n",
    "    test_statistic_under_model2 = compute_test_statistic(n, p1, p2, number_experiments)   \n",
    "    \n",
    "    bins = np.linspace(0, 15, 100)\n",
    "    plt.hist(test_statistic_under_model1, bins, alpha=0.5, label='model 1', density=True)\n",
    "    plt.hist(test_statistic_under_model2, bins, alpha=0.5, label='model 2', density=True)\n",
    "\n",
    "    degrees_of_freedom = 1 # there are two parameters in general model, and one parameter in nested model\n",
    "    plt.plot(bins, chi2.pdf(bins, degrees_of_freedom), label='chi-square')\n",
    "\n",
    "    plt.legend(loc='upper right')\n",
    "    plt.xlabel('likelihood-ratio test statistic')\n",
    "    plt.ylabel('probability density')\n",
    "    plt.ylim([0, 1])\n",
    "    \n",
    "plot_test_statistic_distributions()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7215e35",
   "metadata": {},
   "source": [
    "## Resources \n",
    "\n",
    "* Wilks, S. S. (1938). “The Large Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses.” Annals of Mathematical Statistics, 9: 60–62. URL http://projecteuclid.org/euclid.aoms/1177732360\n",
    "* Wilks' Theorem: https://en.wikipedia.org/wiki/Wilks%27_theorem\n",
    "* Likelihood ratio test controls for overfitting: https://stats.stackexchange.com/a/213481\n",
    "* Similar analysis to this notebook: https://towardsdatascience.com/the-likelihood-ratio-test-463455b34de9\n",
    "* Coin-tossing math: https://en.wikipedia.org/wiki/Wilks%27_theorem#Coin_tossing\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dd44c65",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}