SectionD10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c70d3b81",
   "metadata": {},
   "source": [
    "# Entropy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b450f8be",
   "metadata": {},
   "source": [
    "## Definition of statistical entropy\n",
    "\n",
    "The entropy $S$ of an isolated system at thermodynamic equilibrium is the natural logarithm of W, the number of distinct microscopic states available to the system given the macroscopic constraints such as a fixed total energy E."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "206edfd3",
   "metadata": {},
   "source": [
    "$$S=k \\log W$$\n",
    "\n",
    "- S: entropy\n",
    "- W: number of possbile arrangements\n",
    "- k: Boltzmann constant used to allow expression in terms of energy. Set to 1 below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1d1506a",
   "metadata": {},
   "source": [
    "The Boltzmann equation relates the microscopic details of the system via W to its macroscopic state via the entropy S. Named after Ludwig Boltzmann who related probability to entropy in 1877, but $k$ was defined by Max Planck around 1900."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f26b601",
   "metadata": {},
   "source": [
    "$$\n",
    "\\Delta S=\\int \\frac{\\mathrm{d} Q}{T}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0286321",
   "metadata": {},
   "source": [
    "- The Boltzmann constant ($k_{B}$ or $k$) relates the statistical entropy to the classical macroscopic entropy. \n",
    "- It relates average relative kinetic energy of particles with the thermodynamic temperature of the system: $$1.380649×10^{−23} J⋅K^{−1}$$ \n",
    "- It has the same units as the classical entropy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddda42b5",
   "metadata": {},
   "source": [
    "### Measurement of the classical entropy of a system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c2103a31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"./SectionD10/images/MacroscopicEntropy.png\" width=\"900\" height=\"900\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url=\"./SectionD10/images/MacroscopicEntropy.png\", width=900, height=900)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e6ab74d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"./SectionD10/images/StatisticalEntropy.png\" width=\"900\" height=\"900\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url=\"./SectionD10/images/StatisticalEntropy.png\", width=900, height=900)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3227159c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"./SectionD10/images/entropy.png\" width=\"900\" height=\"900\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url=\"./SectionD10/images/entropy.png\", width=900, height=900)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccb255b4",
   "metadata": {},
   "source": [
    "### Entropy quantifies the possible arrangement of objects\n",
    "\n",
    "- molecules\n",
    "- characters in text\n",
    "- pixels\n",
    "- bubbles in sourdough bread\n",
    "- books in a bookcase"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f9bbb63",
   "metadata": {},
   "source": [
    "## Rescaled, dimensionless Entropy == Information Entropy\n",
    "\n",
    "Set $k=1$:\n",
    "$$\n",
    "S^{\\prime}=\\ln W, \\quad \\Delta S^{\\prime}=\\int \\frac{\\mathrm{d} Q}{k T}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8625e3d1",
   "metadata": {},
   "source": [
    "$$S' = H(p)=-\\mathbb{E}[\\log p]=-\\sum_i^n p_i \\log p_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9751a67c",
   "metadata": {},
   "source": [
    "For a proability distribution p with n possible distinct outcomes, each event $i$ has a probability of $p_{i}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d5ada05",
   "metadata": {},
   "source": [
    "### Multiplicity, W\n",
    "\n",
    "\n",
    "$$\n",
    "W=\\frac{N !}{n_{1} ! n_{2} ! \\cdots n_{t} !}\n",
    "$$\n",
    "\n",
    "\n",
    "- $N:$ total number of events\n",
    "- $n_{i}$: number of times distinct event $i$ occurs\n",
    "- $t$: number of distinct events\n",
    "\n",
    "\n",
    "\n",
    "With large $N$, we can apply Stirling's approximation:\n",
    "$$x ! \\approx\\left(\\frac{x}{e}\\right)^x$$\n",
    "\n",
    "\n",
    "Because $p_i=\\frac{n_i}{N}$, we can rewrite the above as\n",
    "\n",
    "$$\n",
    "W=\\frac{1}{p_1^{n_1} p_2^{n_2} \\cdots p_t^{n_t}}\n",
    "$$\n",
    "\n",
    "Take the logarithm (????), \n",
    "\n",
    "$$\n",
    "\\log W=-\\sum_i^n p_i \\log p_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "75edd9bc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:27:35) [Clang 14.0.6 ]\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print(sys.version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "38d6cda8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/blaine/opt/anaconda3/envs/stat-rethink2-pymc_v4/lib/python3.10/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.5\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n",
      "2022-12-24 06:04:54.801900: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n"
     ]
    }
   ],
   "source": [
    "# Basic\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "from patsy import bs, dmatrix\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Exploratory Analysis of Bayesian Models\n",
    "import arviz as az\n",
    "\n",
    "# Probabilistic programming languages\n",
    "import bambi as bmb\n",
    "import pymc as pm\n",
    "import tensorflow_probability as tfp\n",
    "\n",
    "tfd =tfp.distributions\n",
    "\n",
    "# Computational Backend\n",
    "import theano\n",
    "import theano.tensor as tt\n",
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1d0d6181",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLsAAAH/CAYAAAC/998mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJe0lEQVR4nO3df3RV9bkn/uck5AcoCUUgAQ0Qq21pcVkbqgW/3KtWUNprcXRaRkd+LMEuqt5ezO10QFpFp+qd/kDsKFirFr1VL3VKrbdw26ZjtbTQ1lKwHXTZ2w4YNIkgbRMwNoFwvn9Yco0JkJMEzsnO67XWWcuz89n7PDmf7vD0fT5n71Q6nU4HAAAAACRAXrYLAAAAAIC+IuwCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAgAAACAxhF0AAAAAJIawCwAAAIDEEHYBvMW6devi/e9/fxQVFUVlZWV8+ctfjmXLlkUqlcp2aQAAx93ReqMdO3ZEKpWK1atXd9o3lUrFsmXLjm/BABExKNsFAOSK//N//k/MnDkzJk+eHP/yL/8SbW1t8cUvfjFeffXVbJcGAHDc6Y2A/krYBfBXS5cujbKysqipqYni4uKIiLjoooti/Pjx2S0MACAL9EZAf+VrjAAR8frrr8ezzz4bl112WXszFxExdOjQuOSSS7JYGQDA8ac3AvozYRdARPzpT3+KgwcPRnl5eaefdbUNACDJ9EZAfybsAoiId7zjHZFKpaKhoaHTz7raBgCQZN3tjQ6t+mppaekwZs+ePce2QIAjEHYBRMQJJ5wQZ599dqxduzb+8pe/tG/fu3dv/Ou//msWKwMAOP662xuVlZVFcXFx/OY3v+mw/3e/+93jVivA27lAPcBf/Y//8T/i4osvjmnTpsU//uM/RltbW/zP//k/44QTTog//vGP2S4PAOC46k5vlEql4qqrrooHH3ww3vnOd8aZZ54Zv/zlL+PRRx/NcvXAQGZlF8BfTZs2LZ544oloamqKWbNmRXV1dVx++eVx9dVXZ7s0AIDjrru90Ve+8pW46qqr4otf/GLMnDkzNm3aFN/73veyVDVARCqdTqezXQRALlu2bFnccsst4c8lAIDeCMh9VnYBAAAAkBjCLgAAAAASw9cYAQAAAEgMK7sAAAAASAxhFwAAAACJIewCAAAAIDEGZbuA7jh48GDU1dXF0KFDI5VKZbscAIA+lU6nY+/evTFmzJjIy+vZZ5H6JQAg6brbM/WLsKuuri4qKiqyXQYAwDG1c+fOOOWUU3q0r34JABgojtYz9Yuwa+jQoRHx5i9TUlKS5WoAAPpWU1NTVFRUtPc8PaFfAgCSrrs9U78Iuw4txS8pKdG8AQCJ1ZuvH+qXAICB4mg9kwvUAwAAAJAYwi4AAAAAEkPYBTmsufVAjF+8LsYvXhfNrQeyXQ4AAADkvH5xzS4AILe0tbXF/v37s11Gv1JQUBD5+fnZLgMAOE7S6XQcOHAg2trasl1Kv5Gfnx+DBg3q1XVMI4RdAECG9u3bFy+//HKk0+lsl9KvpFKpOOWUU+LEE0/MdikAwDHW2toa9fX10dzcnO1S+p0hQ4bE6NGjo7CwsMfHEHYBAN3W1tYWL7/8cgwZMiRGjhzZ60/dBop0Oh27d++Ol19+OU4//XQrvAAgwQ4ePBjbt2+P/Pz8GDNmTBQWFuqZuiGdTkdra2vs3r07tm/fHqeffnrk5fXs6lvCLgCg2/bv3x/pdDpGjhwZgwcPznY5/crIkSNjx44dsX//fmEXACRYa2trHDx4MCoqKmLIkCHZLqdfGTx4cBQUFMRLL70Ura2tUVxc3KPjuEA9AJAxn05mznsGAANLT1clDXR98b555wEAAABIDGEXAAAAAIkh7AIAAAAgMYRdAMBxN2/evLj00ks7bX/66acjlUrFn//85z57rX/4h3+IqqqqKCoqive///1HHb9jx45IpVJdPh5//PH2cbfddltMmTIlhgwZEsOGDeuzegEADjlePdNzzz0XV1xxRVRUVMTgwYNjwoQJcdddd3V7/3Q6HTNmzIhUKhVPPPFEh5/96U9/itmzZ0dpaWmUlpbG7Nmz+7TX64qwCwBItHQ6HVdffXXMmjWrW+MrKiqivr6+w+OWW26JE044IWbMmNE+rrW1NT7+8Y/Hpz71qWNVOgDAcbF58+YYOXJkfPOb34xt27bF0qVLY8mSJXH33Xd3a/8VK1Yc9mY8V155ZWzdujW+//3vx/e///3YunVrzJ49uy/L72TQMT06AECWffWrX42IiN27d8dvfvObo47Pz8+P8vLyDtu+853vxKxZs+LEE09s33bLLbdERMTq1av7rlgAgCy4+uqrOzw/9dRTY9OmTbF27dq4/vrrj7jvc889F8uXL49nn302Ro8e3eFnL7zwQnz/+9+Pn//853HOOedERMTXv/71mDx5crz44ovx7ne/u29/kb+ysgsA6FdmzJgRJ5544hEffWnz5s2xdevWmD9/fp8eFwDgWOptz9TY2BjDhw8/4pjm5ua44oor4u677+70YWFExKZNm6K0tLQ96IqI+NCHPhSlpaWxcePGnv1i3WBlF/SR5tYD8d6bfhAREc/felEMKXR6ARzJ9773vU5NVltb21H3u//+++ONN944VmV18sADD8SECRNiypQpx+01AQAOyUbPtGnTpvjWt74V69atO+K4G264IaZMmRIzZ87s8ucNDQ0xatSoTttHjRoVDQ0NPaqtO/y/cQAgK84///xYtWpVh22/+MUv4qqrrjrifieffPKxLKuDN954Ix599NH4/Oc/f9xeEwDgrY53z7Rt27aYOXNm3HTTTTFt2rTDjnvyySfjqaeeii1bthzxeF1dyyudTh/2Gl99QdgFAGTFCSecEKeddlqHbS+//PJR95sxY0Zs2LDhiGP27dvXq9oO+d//+39Hc3NzzJkzp0+OBwCQqePZMz3//PNxwQUXxDXXXBOf+9znjrjvU089FX/4wx863ZX68ssvj6lTp8bTTz8d5eXl8eqrr3bad/fu3VFWVnbU36GnhF0AQL9yPL/G+MADD8THPvaxGDly5HF5PQCAvpJpz7Rt27a44IILYu7cuXHbbbcddfzixYtjwYIFHbadccYZceedd8Yll1wSERGTJ0+OxsbG+OUvfxlnn312RLy5Kq2xsfGYXiJC2AUA9CuZLsn//e9/H/v27YuGhoZ44403YuvWrRER8d73vjcKCwvjlVdeiQ9/+MPx8MMPtzdhh/b7yU9+EuvXr+/yuLW1tfHHP/4xamtro62trf24p512Wp9fJB8AIFOZ9Ezbtm2L888/P6ZPnx7V1dXt19PKz89v/9Dv7T1TeXl5lxelHzt2bFRWVkZExIQJE+Liiy+Oa665Jr72ta9FRMQnP/nJ+Lu/+7tjdifGCGEXAJBwCxYsiGeeeab9+VlnnRUREdu3b4/x48fH/v3748UXX4zm5uYO+z344INx8sknx/Tp07s87k033RQPPfRQp+P++Mc/jvPOO6+PfwsAgGPn8ccfj927d8cjjzwSjzzySPv2cePGxY4dOyIiDtszHc0jjzwSn/70p9t7qo997GNx991391ntXUml0+n0MX2FPtDU1BSlpaXR2NgYJSUl2S4HunQs7sboDo9ArvnLX/4S27dvj8rKyiguLs52Of3Kkd67vuh19EsAkBv0S73TFz1T3rEuEgAAAACOF2EXAAAAAIkh7AIAAAAgMYRdAAAAACSGsAsAyFg/uL9NzvGeAcDA4t/+numL903YBQB0W35+fkREtLa2ZrmS/ufQe3boPQQAkqmgoCAiIpqbm7NcSf906H079D72xKC+KgYASL5BgwbFkCFDYvfu3VFQUBB5eT43646DBw/G7t27Y8iQITFokPYLAJIsPz8/hg0bFrt27YqIiCFDhkQqlcpyVbkvnU5Hc3Nz7Nq1K4YNG9arDwh1WwBAt6VSqRg9enRs3749XnrppWyX06/k5eXF2LFjNbsAMACUl5dHRLQHXnTfsGHD2t+/nhJ2AQAZKSwsjNNPP91XGTNUWFhoJRwADBCHPiAcNWpU7N+/P9vl9BsFBQV9cskHYRcAkLG8vLwoLi7OdhkAADktPz/f9TqzoEcfL65cuTIqKyujuLg4qqqqYsOGDYcdu3bt2pg2bVqMHDkySkpKYvLkyfGDH/ygxwUDAAAAwOFkHHatWbMmFi1aFEuXLo0tW7bE1KlTY8aMGVFbW9vl+J/85Ccxbdq0WL9+fWzevDnOP//8uOSSS2LLli29Lh4AAAAA3irjsGv58uUxf/78WLBgQUyYMCFWrFgRFRUVsWrVqi7Hr1ixIj772c/GBz/4wTj99NPj9ttvj9NPPz3+9V//tdfFAwAAAMBbZXTNrtbW1ti8eXMsXry4w/bp06fHxo0bu3WMgwcPxt69e2P48OGHHdPS0hItLS3tz5uamjIpEwAg8fRLAABdy2hl12uvvRZtbW1RVlbWYXtZWVk0NDR06xhf+cpX4vXXX49PfOIThx1zxx13RGlpafujoqIikzIBABJPvwQA0LUeXaA+lUp1eJ5Opztt68pjjz0Wy5YtizVr1sSoUaMOO27JkiXR2NjY/ti5c2dPygQASCz9EgBA1zL6GuOIESMiPz+/0yquXbt2dVrt9XZr1qyJ+fPnx+OPPx4XXnjhEccWFRVFUVFRJqUBAAwo+iUAgK5ltLKrsLAwqqqqoqampsP2mpqamDJlymH3e+yxx2LevHnx6KOPxkc/+tGeVQoAAAAAR5HRyq6IiOrq6pg9e3ZMmjQpJk+eHPfdd1/U1tbGwoULI+LNJfWvvPJKPPzwwxHxZtA1Z86cuOuuu+JDH/pQ+6qwwYMHR2lpaR/+KgAAAAAMdBmHXbNmzYo9e/bErbfeGvX19TFx4sRYv359jBs3LiIi6uvro7a2tn381772tThw4EBcd911cd1117Vvnzt3bqxevbr3vwGQkebWA/Hem34QERHP33pRDCnM+M8AAAAA5Kwe/b/ca6+9Nq699touf/b2AOvpp5/uyUsAAAAAQMZ6dDdGAAAAAMhFwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAgAAACAxhF0AAAAAJIawiwGrufVAjF+8LsYvXhfNrQeyXQ4AAADQB4RdAAAAACSGsAsAAACAxBB2AQAAAJAYwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEiMHoVdK1eujMrKyiguLo6qqqrYsGHDYcfW19fHlVdeGe9+97sjLy8vFi1a1NNaAQAAAOCIMg671qxZE4sWLYqlS5fGli1bYurUqTFjxoyora3tcnxLS0uMHDkyli5dGmeeeWavCwYAAACAw8k47Fq+fHnMnz8/FixYEBMmTIgVK1ZERUVFrFq1qsvx48ePj7vuuivmzJkTpaWlvS4YAAAAAA5nUCaDW1tbY/PmzbF48eIO26dPnx4bN27ss6JaWlqipaWl/XlTU1OfHRsAIAn0SwAAXctoZddrr70WbW1tUVZW1mF7WVlZNDQ09FlRd9xxR5SWlrY/Kioq+uzYAABJoF8CAOhajy5Qn0qlOjxPp9OdtvXGkiVLorGxsf2xc+fOPjs2AEAS6JcAALqW0dcYR4wYEfn5+Z1Wce3atavTaq/eKCoqiqKioj47HgBA0uiXAAC6ltHKrsLCwqiqqoqampoO22tqamLKlCl9WhgAAAAAZCrjrzFWV1fH/fffHw8++GC88MILccMNN0RtbW0sXLgwIt5cUj9nzpwO+2zdujW2bt0a+/bti927d8fWrVvj+eef75vfAACA46K59UCMX7wuxi9eF82tBxwzh443kI/ZH2o8VscEoGsZfY0xImLWrFmxZ8+euPXWW6O+vj4mTpwY69evj3HjxkVERH19fdTW1nbY56yzzmr/782bN8ejjz4a48aNix07dvSuegAAAAB4i4zDroiIa6+9Nq699touf7Z69epO29LpdE9eBgAAAAAy0qO7MQIAAABALhJ2AQAAAJAYwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIu+oXm1gMxfvG6GL94XTS3Hsh2OQAAAECOEnYBvSaMBAAAIFcIuwAAAABIDGEXAAAAAIkh7AIAAAAgMYRdAAAAACSGsAsAAACAxBB2AQAAAJAYwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2EWfa249EOMXr4vxi9dFc+uBbJcDAAAADCDCLgAAAAASQ9gFAAAAQGL0KOxauXJlVFZWRnFxcVRVVcWGDRuOOP6ZZ56JqqqqKC4ujlNPPTXuvffeHhULAAAAAEeScdi1Zs2aWLRoUSxdujS2bNkSU6dOjRkzZkRtbW2X47dv3x4f+chHYurUqbFly5a48cYb49Of/nR8+9vf7nXxAAAAAPBWGYddy5cvj/nz58eCBQtiwoQJsWLFiqioqIhVq1Z1Of7ee++NsWPHxooVK2LChAmxYMGCuPrqq+PLX/5yr4un91xMHgAAAEiSjMKu1tbW2Lx5c0yfPr3D9unTp8fGjRu73GfTpk2dxl900UXxq1/9Kvbv39/lPi0tLdHU1NThAQDAf9AvAQB0LZVOp9PdHVxXVxcnn3xy/OxnP4spU6a0b7/99tvjoYceihdffLHTPu9617ti3rx5ceONN7Zv27hxY5x77rlRV1cXo0eP7rTPsmXL4pZbbum0vbGxMUpKSrpbLgBAv9DU1BSlpaUZ9Tr6JQBgoOluz9SjC9SnUqkOz9PpdKdtRxvf1fZDlixZEo2Nje2PnTt39qRMAIDE0i8BAHRtUCaDR4wYEfn5+dHQ0NBh+65du6KsrKzLfcrLy7scP2jQoDjppJO63KeoqCiKiooyKQ0AYEDRLwEAdC2jlV2FhYVRVVUVNTU1HbbX1NR0+FrjW02ePLnT+B/+8IcxadKkKCgoyLBcAAAAADi8jL/GWF1dHffff388+OCD8cILL8QNN9wQtbW1sXDhwoh4c0n9nDlz2scvXLgwXnrppaiuro4XXnghHnzwwXjggQfiM5/5TN/9FgAAAAAQGX6NMSJi1qxZsWfPnrj11lujvr4+Jk6cGOvXr49x48ZFRER9fX3U1ta2j6+srIz169fHDTfcEPfcc0+MGTMmvvrVr8bll1/e7dc8dI0vdxkCAJLoUI+TwX2DOtEvAQBJ192eKaO7MWbLyy+/HBUVFdkuAwDgmNq5c2eccsopPdpXvwQADBRH65n6Rdh18ODBqKuri6FDhx7xro+91dTUFBUVFbFz50637M5B5ie3mZ/cZn5yl7nJbcdrftLpdOzduzfGjBkTeXk9ulm2fomIMD+5zvzkNvOT28xPbsu1ninjrzFmQ15eXo8/5eyJkpISJ08OMz+5zfzkNvOTu8xNbjse81NaWtqr/fVLvJX5yW3mJ7eZn9xmfnJbrvRMPfvoEAAAAABykLALAAAAgMQQdr1FUVFR3HzzzVFUVJTtUuiC+clt5ie3mZ/cZW5ym/npzHuS28xPbjM/uc385Dbzk9tybX76xQXqAQAAAKA7rOwCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAuiG5ubmbJcAAJAz9EZALhN2AbzNsmXLIpVKxa9//ev4z//5P8c73vGOeOc735ntsgAAsuJwvdF5550X5513Xqfx8+bNi/Hjxx/3OgEOEXYBHMZll10Wp512Wjz++ONx7733ZrscAICs0hsB/cWgbBcAkKvmzp0bt9xyS7bLAADICW/vje68884sVgNweFZ2ARzG5Zdfnu0SAAByht4I6C+EXQCHMXr06GyXAACQM/RGQH8h7AI4jFQqle0SAAByxtt7o+Li4mhpaek07rXXXjteJQF0SdgFAABAxsaPHx+/+93vOgRee/bsiY0bN2axKgBhFwAAAD0we/bs+OMf/xhXXXVV/PCHP4zHHnssLrzwwigpKcl2acAAJ+wCAAAgY+eee2489NBDsW3btpg5c2Z84QtfiCVLlsR5552X7dKAAS6VTqfT2S4CAAAAAPqClV0AAAAAJIawCwAAAIDEEHYBAAAAkBjCLgAAAAASQ9gFAAAAQGIIuwAAAABIjEHZLqA7Dh48GHV1dTF06NBIpVLZLgcAoE+l0+nYu3dvjBkzJvLyevZZpH4JAEi67vZM/SLsqquri4qKimyXAQBwTO3cuTNOOeWUHu2rXwIABoqj9Uz9IuwaOnRoRLz5y5SUlGS5GgCAvtXU1BQVFRXtPU9P6JcAgKTrbs/UL8KuQ0vxS0pKNG8AQGL15uuH+iUAYKA4Ws/kAvUAAAAAJIawCwAAAIDEEHYBAAAAkBj94ppdAEBuOXjwYLS2tma7jH6lsLDwiLfIBgCSp62tLfbv35/tMvqNgoKCyM/P7/VxhF0AQEZaW1tj+/btcfDgwWyX0q/k5eVFZWVlFBYWZrsUAOAYS6fT0dDQEH/+85+zXUq/M2zYsCgvL+/VjXuEXQBAt6XT6aivr4/8/PyoqKiwUqmbDh48GHV1dVFfXx9jx47tVfMGAOS+Q0HXqFGjYsiQIf7t74Z0Oh3Nzc2xa9euiIgYPXp0j48l7AIAuu3AgQPR3NwcY8aMiSFDhmS7nH5l5MiRUVdXFwcOHIiCgoJslwMAHCNtbW3tQddJJ52U7XL6lcGDB0dExK5du2LUqFE9/kqjj2MBgG5ra2uLiPBVvB449J4deg8BgGQ6dI0uHwz2zKH3rTfXOhN2AQAZsxQ/c94zABhY/NvfM33xvgm7AAAAAEgMYRcAAAAAiSHsAgCOu3nz5sWll17aafvTTz8dqVSqz27T/dxzz8UVV1wRFRUVMXjw4JgwYULcdddd3d4/nU7HjBkzIpVKxRNPPNHhZ7fddltMmTIlhgwZEsOGDeuTegEA3up49Ux79uyJiy++OMaMGRNFRUVRUVER119/fTQ1NR1xv/vuuy/OO++8KCkpOWo9LS0t8f73vz9SqVRs3bq1T+o+HGEXAJBYmzdvjpEjR8Y3v/nN2LZtWyxdujSWLFkSd999d7f2X7FixWGvG9Ha2hof//jH41Of+lRflgwAcNzl5eXFzJkz48knn4zf/e53sXr16vjRj34UCxcuPOJ+zc3NcfHFF8eNN9541Nf47Gc/G2PGjOmrko9o0HF5FQCALLj66qs7PD/11FNj06ZNsXbt2rj++uuPuO9zzz0Xy5cvj2effTZGjx7d6ee33HJLRESsXr26z+oFAMiGd7zjHR0+wBs3blxce+218aUvfemI+y1atCgi3lxpdiT/9m//Fj/84Q/j29/+dvzbv/1bb8s9KmEXANCvzJgxIzZs2HDEMfv27TvszxobG2P48OFH3L+5uTmuuOKKuPvuu6O8vLxHdQIAZFNveqa6urpYu3Zt/O3f/m2v63j11VfjmmuuiSeeeCKGDBnS6+N1h7ALAMiK733ve3HiiSd22NbW1nbU/e6///544403evSamzZtim9961uxbt26I4674YYbYsqUKTFz5swevQ4AQF85nj3TFVdcEd/97nfjjTfeiEsuuSTuv//+jPZ/u3Q6HfPmzYuFCxfGpEmTYseOHb06XncJuwCArDj//PNj1apVHbb94he/iKuuuuqI+5188sk9er1t27bFzJkz46abbopp06YddtyTTz4ZTz31VGzZsqVHrwMA0JeOZ8905513xs033xwvvvhi3HjjjVFdXR0rV67M+DiH/K//9b+iqakplixZ0uNj9ISwCwDIihNOOCFOO+20Dttefvnlo+7XkyX5zz//fFxwwQVxzTXXxOc+97kj7vvUU0/FH/7wh053WLz88stj6tSpR70mBQBAXzqePVN5eXmUl5fHe97znjjppJNi6tSp8fnPf77L65d2x1NPPRU///nPo6ioqMP2SZMmxX/9r/81HnrooR4d92iEXQBAv5Lpkvxt27bFBRdcEHPnzo3bbrvtqOMXL14cCxYs6LDtjDPOiDvvvDMuueSSjOsFAMiG3lz6IeLNryBGRLS0tPT4GF/96lfjC1/4Qvvzurq6uOiii2LNmjVxzjnn9Pi4RyPsAgD6lUyW5G/bti3OP//8mD59elRXV0dDQ0NEROTn58fIkSMjIuKVV16JD3/4w/Hwww/H2Wef3f6J5tuNHTs2Kisr25/X1tbGH//4x6itrY22trbYunVrREScdtppna6rAQBwvGXSM61fvz5effXV+OAHPxgnnnhiPP/88/HZz342zj333Bg/fnxEdO6ZIiIaGhqioaEhfv/730dExG9/+9sYOnRojB07NoYPHx5jx47t8DqHeqR3vvOdccopp/TBb9k1YRcAkFiPP/547N69Ox555JF45JFH2rePGzeu/QKp+/fvjxdffDGam5szOvZNN93UYen9WWedFRERP/7xj+O8887rde0AAMfL4MGD4+tf/3rccMMN0dLSEhUVFXHZZZfF4sWL28d01TPde++9ccstt7Q//5u/+ZuIiPjGN74R8+bNO271v10qfWhdWg5ramqK0tLSaGxsjJKSkmyXAwAD1l/+8pfYvn17VFZWRnFxcbbL6VeO9N71Ra+jXwKA3KBf6p2+6JnyjnWRAAAAAHC8CLsAAAAASAxhFwAAAACJIewCAAAAIDGEXQBAxvrB/W1yjvcMAAYW//b3TF+8b8IuAKDb8vPzIyKitbU1y5X0P4fes0PvIQCQTAUFBRER0dzcnOVK+qdD79uh97EnBvVVMQBA8g0aNCiGDBkSu3fvjoKCgsjL87lZdxw8eDB2794dQ4YMiUGDtF8AkGT5+fkxbNiw2LVrV0REDBkyJFKpVJaryn3pdDqam5tj165dMWzYsF59QKjbAgC6LZVKxejRo2P79u3x0ksvZbucfiUvLy/Gjh2r2QWAAaC8vDwioj3wovuGDRvW/v71lLALAMhIYWFhnH766b7KmKHCwkIr4QBggDj0AeGoUaNi//792S6n3ygoKOiTSz4IuwCAjOXl5UVxcXG2ywAAyGn5+fmu15kFPfp4ceXKlVFZWRnFxcVRVVUVGzZsOOzYtWvXxrRp02LkyJFRUlISkydPjh/84Ac9LhgAAAAADifjsGvNmjWxaNGiWLp0aWzZsiWmTp0aM2bMiNra2i7H/+QnP4lp06bF+vXrY/PmzXH++efHJZdcElu2bOl18QAAAADwVql0Op3OZIdzzjknPvCBD8SqVavat02YMCEuvfTSuOOOO7p1jPe9730xa9asuOmmm7o1vqmpKUpLS6OxsTFKSkoyKRcAIOf1Ra+jXwIAkq67/U5G1+xqbW2NzZs3x+LFiztsnz59emzcuLFbxzh48GDs3bs3hg8fftgxLS0t0dLS0v68qakpkzIBABJPvwQA0LWMvsb42muvRVtbW5SVlXXYXlZWFg0NDd06xle+8pV4/fXX4xOf+MRhx9xxxx1RWlra/qioqMikTACAxNMvAQB0rUcXqE+lUh2ep9PpTtu68thjj8WyZctizZo1MWrUqMOOW7JkSTQ2NrY/du7c2ZMyAQASS78EANC1jL7GOGLEiMjPz++0imvXrl2dVnu93Zo1a2L+/Pnx+OOPx4UXXnjEsUVFRVFUVJRJaQAAA4p+CQCgaxmt7CosLIyqqqqoqanpsL2mpiamTJly2P0ee+yxmDdvXjz66KPx0Y9+tGeVAgAAAMBRZLSyKyKiuro6Zs+eHZMmTYrJkyfHfffdF7W1tbFw4cKIeHNJ/SuvvBIPP/xwRLwZdM2ZMyfuuuuu+NCHPtS+Kmzw4MFRWlrah78KAAAAAANdxmHXrFmzYs+ePXHrrbdGfX19TJw4MdavXx/jxo2LiIj6+vqora1tH/+1r30tDhw4ENddd11cd9117dvnzp0bq1ev7v1vAAAAAAB/lUqn0+lsF3E0TU1NUVpaGo2NjVFSUpLtcgAA+lRf9Dr6JQAg6brb7/TobowAAAAAkIuEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAgAAACAxhF0AAAAAJIawCwAAAIDEEHYBAAAAkBjCLgAAAAASQ9gFAAAAQGIIuwAAAABIDGEXAAAAAIkh7AIAAAAgMYRdkMOaWw/E+MXrYvziddHceiDb5QB9zDkOAAB9T9gFAAAAQGIIuwAAAABIDGEXAAAAAIkh7AIAAAAgMYRdAAAAACSGsAsAAACAxBB2QR9pbj0Q4xevi/GL10Vz64FslwMcA85zAADIfcIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAgAAACAxehR2rVy5MiorK6O4uDiqqqpiw4YNhx1bX18fV155Zbz73e+OvLy8WLRoUU9rBQAAAIAjyjjsWrNmTSxatCiWLl0aW7ZsialTp8aMGTOitra2y/EtLS0xcuTIWLp0aZx55pm9LhgAAAAADifjsGv58uUxf/78WLBgQUyYMCFWrFgRFRUVsWrVqi7Hjx8/Pu66666YM2dOlJaW9rpgAAAAADicQZkMbm1tjc2bN8fixYs7bJ8+fXps3Lixz4pqaWmJlpaW9udNTU19dmwAgCTQLwEAdC2jlV2vvfZatLW1RVlZWYftZWVl0dDQ0GdF3XHHHVFaWtr+qKio6LNjAwAkgX4JAKBrPbpAfSqV6vA8nU532tYbS5YsicbGxvbHzp07++zYAABJoF8CAOhaRl9jHDFiROTn53daxbVr165Oq716o6ioKIqKivrseAAASaNfAgDoWkYruwoLC6Oqqipqamo6bK+pqYkpU6b0aWHAsdHceiDGL14X4xevi+bWA9kuB+hjznEAAAa6jFZ2RURUV1fH7NmzY9KkSTF58uS47777ora2NhYuXBgRby6pf+WVV+Lhhx9u32fr1q0REbFv377YvXt3bN26NQoLC+O9731v3/wWAAAAABA9CLtmzZoVe/bsiVtvvTXq6+tj4sSJsX79+hg3blxERNTX10dtbW2Hfc4666z2/968eXM8+uijMW7cuNixY0fvqgcAAACAt8g47IqIuPbaa+Paa6/t8merV6/utC2dTvfkZQAAAAAgIz26GyMAAAAA5CJhFwAAAACJIewCAAAAIDGEXQAAAAAkhrCLAau59UCMX7wuxi9eF82tB7JdDtDHnOMAADAwCbsAAAAASAxhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAgAAACAxhF0AAHRLc+uBGL94XYxfvC6aWw84Zg4dbyAfsz/UOJCP2R9q7C/H7A81DuRjHosa6TlhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIu+gV3toBkc44DAAB9RdgF9JqgApLNOQ4AQH8i7AIAAAAgMYRdAAAAACSGsAsAAACAxBB2AQAAAJAYwi4AAAAAEkPYBQAAAEBi9CjsWrlyZVRWVkZxcXFUVVXFhg0bjjj+mWeeiaqqqiguLo5TTz017r333h4VCwAAAABHknHYtWbNmli0aFEsXbo0tmzZElOnTo0ZM2ZEbW1tl+O3b98eH/nIR2Lq1KmxZcuWuPHGG+PTn/50fPvb3+518eSm5tYDMX7xuhi/eF00tx7IdjlAH3OOAwAAuSzjsGv58uUxf/78WLBgQUyYMCFWrFgRFRUVsWrVqi7H33vvvTF27NhYsWJFTJgwIRYsWBBXX311fPnLX+518QAAAADwVhmFXa2trbF58+aYPn16h+3Tp0+PjRs3drnPpk2bOo2/6KKL4le/+lXs37+/y31aWlqiqampw4NjwwoNSD7nOSSTfgkAoGupdDqd7u7gurq6OPnkk+NnP/tZTJkypX377bffHg899FC8+OKLnfZ517veFfPmzYsbb7yxfdvGjRvj3HPPjbq6uhg9enSnfZYtWxa33HJLp+2NjY1RUlLS3XIBAPqFpqamKC0tzajX0S8BAANNd3umHl2gPpVKdXieTqc7bTva+K62H7JkyZJobGxsf+zcubMnZQIAJJZ+CQCga4MyGTxixIjIz8+PhoaGDtt37doVZWVlXe5TXl7e5fhBgwbFSSed1OU+RUVFUVRUlElpAAADin4JAKBrGa3sKiwsjKqqqqipqemwvaampsPXGt9q8uTJncb/8Ic/jEmTJkVBQUGG5QIAAADA4WX8Ncbq6uq4//7748EHH4wXXnghbrjhhqitrY2FCxdGxJtL6ufMmdM+fuHChfHSSy9FdXV1vPDCC/Hggw/GAw88EJ/5zGf67rcAAAAAgMjwa4wREbNmzYo9e/bErbfeGvX19TFx4sRYv359jBs3LiIi6uvro7a2tn18ZWVlrF+/Pm644Ya45557YsyYMfHVr341Lr/88m6/5qFrfLnLEACQRId6nAzuG9SJfgkASLru9kwZ3Y0xW15++eWoqKjIdhkAAMfUzp0745RTTunRvvolAGCgOFrP1C/CroMHD0ZdXV0MHTr0iHd97K2mpqaoqKiInTt3umV3DjI/uc385Dbzk7vMTW47XvOTTqdj7969MWbMmMjL69HNsvVLRIT5yXXmJ7eZn9xmfnJbrvVMGX+NMRvy8vJ6/ClnT5SUlDh5cpj5yW3mJ7eZn9xlbnLb8Zif0tLSXu2vX+KtzE9uMz+5zfzkNvOT23KlZ+rZR4cAAAAAkIOEXQAAAAAkhrDrLYqKiuLmm2+OoqKibJdCF8xPbjM/uc385C5zk9vMT2fek9xmfnKb+clt5ie3mZ/clmvz0y8uUA8AAAAA3WFlFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAgAAACAxhF0Ab7F79+745Cc/GRUVFVFUVBQjR46Mc889N370ox9luzQAgOPuaL3R+PHjY968eZ32O++88+K88847vsUC/NWgbBcAkEtmz54dv/71r+O2226Ld73rXfHnP/85fv3rX8eePXuyXRoAwHGnNwL6I2EXwFv87Gc/iwULFsQ111zTvm3mzJlZrAgAIHv0RkB/JOwCeIuzzz47Vq9eHSeddFJceOGFUVVVFQUFBdkuCwAgK/RGQH/kml0Ab7FmzZqYO3du3H///TF58uQYPnx4zJkzJxoaGrJdGgDAcac3AvojYRfAW4wYMSJWrFgRO3bsiJdeeinuuOOOWLt2bZcXXgUASLqj9UbFxcXR0tLSab/XXnvtOFcK8B+EXQCHMXbs2Lj++utj2rRp8etf/zrb5QAAZFVXvdH48ePjN7/5TYdxv/vd7+LFF1/MRokAEeGaXQDtGhsb4/zzz48rr7wy3vOe98TQoUPj2Wefje9///tx2WWXZbs8AIDjqju90ezZs+Oqq66Ka6+9Ni6//PJ46aWX4otf/GKMHDkyy9UDA5mwC+CviouL45xzzol//ud/jh07dsT+/ftj7Nix8d//+3+Pz372s9kuDwDguOpOb3TllVdGXV1d3HvvvfGNb3wjJk6cGKtWrYpbbrkly9UDA1kqnU6ns10EAAAAAPQF1+wCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJjULYL6I6DBw9GXV1dDB06NFKpVLbLAQDoU+l0Ovbu3RtjxoyJvLyefRapXwIAkq67PVO/CLvq6uqioqIi22UAABxTO3fujFNOOaVH++qXAICB4mg9U78Iu4YOHRoRb/4yJSUlWa4GAKBvNTU1RUVFRXvP0xP6JQAg6brbM/WLsOvQUvySkhLNGwCQWL35+qF+CQAYKI7WM7lAPQAAAACJIewCAAAAIDGEXQAAAAAkRr+4ZhcAkFsOHjwYra2t2S6jXyksLDziLbIBgORpa2uL/fv3Z7uMfqOgoCDy8/N7fRxhFwCQkdbW1ti+fXscPHgw26X0K3l5eVFZWRmFhYXZLgUAOMbS6XQ0NDTEn//852yX0u8MGzYsysvLe3XjHmEXANBt6XQ66uvrIz8/PyoqKqxU6qaDBw9GXV1d1NfXx9ixY3vVvAEAue9Q0DVq1KgYMmSIf/u7IZ1OR3Nzc+zatSsiIkaPHt3jYwm7AIBuO3DgQDQ3N8eYMWNiyJAh2S6nXxk5cmTU1dXFgQMHoqCgINvlAADHSFtbW3vQddJJJ2W7nH5l8ODBERGxa9euGDVqVI+/0ujjWACg29ra2iIifBWvBw69Z4feQwAgmQ5do8sHgz1z6H3rzbXOhF0AQMYsxc+c9wwABhb/9vdMX7xvwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwBw3M2bNy8uvfTSTtuffvrpSKVS8ec//7lPXue5556LK664IioqKmLw4MExYcKEuOuuu7q9fzqdjhkzZkQqlYonnniiw89uu+22mDJlSgwZMiSGDRvWJ/UCALzV8eqZ9uzZExdffHGMGTMmioqKoqKiIq6//vpoamo64n4tLS3x93//9zFixIg44YQT4mMf+1i8/PLLHcb8+te/jmnTpsWwYcPipJNOik9+8pOxb9++Pqn7cIRdAEBibd68OUaOHBnf/OY3Y9u2bbF06dJYsmRJ3H333d3af8WKFYe9SGpra2t8/OMfj0996lN9WTIAwHGXl5cXM2fOjCeffDJ+97vfxerVq+NHP/pRLFy48Ij7LVq0KL7zne/Ev/zLv8RPf/rT2LdvX/zd3/1d+92n6+rq4sILL4zTTjstfvGLX8T3v//92LZtW8ybN++Y/j6DjunRAQCy6Oqrr+7w/NRTT41NmzbF2rVr4/rrrz/ivs8991wsX748nn322Rg9enSnn99yyy0REbF69eo+qxcAIBve8Y53dPgAb9y4cXHttdfGl770pcPu09jYGA888ED88z//c1x44YUREfHNb34zKioq4kc/+lFcdNFF8b3vfS8KCgrinnvuiby8N9db3XPPPXHWWWfF73//+zjttNOOye9jZRcA0K/MmDEjTjzxxCM+jqSxsTGGDx9+xDHNzc1xxRVXxN133x3l5eV9WT4AwHHRm56prq4u1q5dG3/7t3972DGbN2+O/fv3x/Tp09u3jRkzJiZOnBgbN26MiDe/5lhYWNgedEVEDB48OCIifvrTn/b2VzwsK7sAgKz43ve+16nJOrTk/Ujuv//+eOONN3r0mps2bYpvfetbsW7duiOOu+GGG2LKlCkxc+bMHr0OAEBfOZ490xVXXBHf/e5344033ohLLrkk7r///sOObWhoiMLCwnjHO97RYXtZWVk0NDRERMQFF1wQ1dXV8aUvfSn+4R/+IV5//fW48cYbIyKivr4+o9oyIewCALLi/PPPj1WrVnXY9otf/CKuuuqqI+538skn9+j1tm3bFjNnzoybbroppk2bdthxTz75ZDz11FOxZcuWHr0OAEBfOp4905133hk333xzvPjii3HjjTdGdXV1rFy5MqNjpNPp9muevu9974uHHnooqqurY8mSJZGfnx+f/vSno6ysLPLz8zOur7uEXQBAVpxwwgmdrtPw9rv3dGXGjBmxYcOGI455+x1+nn/++bjgggvimmuuic997nNH3Pepp56KP/zhD53usHj55ZfH1KlT4+mnnz5qjQAAfeV49kzl5eVRXl4e73nPe+Kkk06KqVOnxuc///kur19aXl4era2t8ac//anD6q5du3bFlClT2p9feeWVceWVV8arr74aJ5xwQqRSqVi+fHlUVlYe9XfoKWEXANCvZLokf9u2bXHBBRfE3Llz47bbbjvq+MWLF8eCBQs6bDvjjDPizjvvjEsuuSTjegEAsqE3l36IeHOFVsSb193qSlVVVRQUFERNTU184hOfiIg3v5r4f//v/40vfvGLncaXlZVFRMSDDz4YxcXFR1xp31vCLgCgX8lkSf62bdvi/PPPj+nTp0d1dXX79SPy8/Nj5MiRERHxyiuvxIc//OF4+OGH4+yzz27/RPPtxo4d2+ETyNra2vjjH/8YtbW10dbWFlu3bo2IiNNOO+2oF8kHADjWMumZ1q9fH6+++mp88IMfjBNPPDGef/75+OxnPxvnnntujB8/PiI690ylpaUxf/78+Md//Mc46aSTYvjw4fGZz3wmzjjjjPa7M0ZE3H333TFlypQ48cQTo6amJv7bf/tv8U//9E+dVtH3JWEXAJBYjz/+eOzevTseeeSReOSRR9q3jxs3Lnbs2BEREfv3748XX3wxmpubMzr2TTfdFA899FD787POOisiIn784x/Heeed1+vaAQCOl8GDB8fXv/71uOGGG6KlpSUqKirisssui8WLF7eP6apnuvPOO2PQoEHxiU98It5444348Ic/HKtXr+5wPa5f/vKXcfPNN8e+ffviPe95T3zta1+L2bNnH9PfJ5U+tC4thzU1NUVpaWk0NjZGSUlJtssBgAHrL3/5S2zfvj0qKyujuLg42+X0K0d67/qi19EvAUBu0C/1Tl/0THnHukgAAAAAOF6EXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuACBj/eD+NjnHewYAA8vBgwezXUK/1Bfv26A+qAMAGCAKCgoilUrF7t27Y+TIkZFKpbJdUr+QTqdj9+7dkUqloqCgINvlAADHUGFhYeTl5UVdXV2MHDkyCgsL9UzdkE6no7W1NXbv3h15eXlRWFjY42MJuwCAbsvPz49TTjklXn755dixY0e2y+lXUqlUnHLKKZGfn5/tUgCAYygvLy8qKyujvr4+6urqsl1OvzNkyJAYO3Zs5OX1/MuIwi4AICMnnnhinH766bF///5sl9KvFBQUCLoAYIAoLCyMsWPHxoEDB6KtrS3b5fQb+fn5MWjQoF6vhBN2AQAZy8/PF9wAABzBocsXuITB8ecC9QAAAAAkRo/CrpUrV0ZlZWUUFxdHVVVVbNiw4bBj165dG9OmTYuRI0dGSUlJTJ48OX7wgx/0uGAAAAAAOJyMw641a9bEokWLYunSpbFly5aYOnVqzJgxI2pra7sc/5Of/CSmTZsW69evj82bN8f5558fl1xySWzZsqXXxQMAAADAW6XS6XQ6kx3OOeec+MAHPhCrVq1q3zZhwoS49NJL44477ujWMd73vvfFrFmz4qabburW+KampigtLY3GxsYoKSnJpFwAgJzXF72OfgkASLru9jsZXaC+tbU1Nm/eHIsXL+6wffr06bFx48ZuHePgwYOxd+/eGD58+GHHtLS0REtLS/vzpqamTMoEAEg8/RIAQNcy+hrja6+9Fm1tbVFWVtZhe1lZWTQ0NHTrGF/5ylfi9ddfj0984hOHHXPHHXdEaWlp+6OioiKTMgEAEk+/BADQtR5doD6VSnV4nk6nO23rymOPPRbLli2LNWvWxKhRow47bsmSJdHY2Nj+2LlzZ0/KBABILP0SAEDXMvoa44gRIyI/P7/TKq5du3Z1Wu31dmvWrIn58+fH448/HhdeeOERxxYVFUVRUVEmpQEADCj6JQCArmW0squwsDCqqqqipqamw/aampqYMmXKYfd77LHHYt68efHoo4/GRz/60Z5VCgAAAABHkdHKroiI6urqmD17dkyaNCkmT54c9913X9TW1sbChQsj4s0l9a+88ko8/PDDEfFm0DVnzpy466674kMf+lD7qrDBgwdHaWlpH/4qAAAAAAx0GYdds2bNij179sStt94a9fX1MXHixFi/fn2MGzcuIiLq6+ujtra2ffzXvva1OHDgQFx33XVx3XXXtW+fO3durF69uve/AQAAAAD8VSqdTqezXcTRNDU1RWlpaTQ2NkZJSUm2ywEA6FN90evolwCApOtuv9OjuzECAAAAQC4SdgEAAACQGMIuAAAAABJD2AUAAABAYgi7AAAAAEgMYRcAAAAAiSHsAgAAACAxhF0AAAAAJIawCwAAAIDEEHYBAAAAkBjCLgAAAAASQ9gFAAAAQGIIuwAAAABIDGEXAAAAAIkh7AIAAAAgMYRdAAAAACSGsAsAAACAxBB2QQ5rbj0Q4xevi/GL10Vz64FslwPQZ/x9AwDgWBF2AQAAAJAYwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwAAAACJIeyCPuLOYkBS+fsGAEB/IuwCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABJD2AUAAABAYvQo7Fq5cmVUVlZGcXFxVFVVxYYNGw47tr6+Pq688sp497vfHXl5ebFo0aKe1goAAAAAR5Rx2LVmzZpYtGhRLF26NLZs2RJTp06NGTNmRG1tbZfjW1paYuTIkbF06dI488wze10wAAAAABxOxmHX8uXLY/78+bFgwYKYMGFCrFixIioqKmLVqlVdjh8/fnzcddddMWfOnCgtLe11wQAAAABwOIMyGdza2hqbN2+OxYsXd9g+ffr02LhxY58V1dLSEi0tLe3Pm5qa+uzYAABJoF8CAOhaRiu7XnvttWhra4uysrIO28vKyqKhoaHPirrjjjuitLS0/VFRUdFnxwYASAL9EgBA13p0gfpUKtXheTqd7rStN5YsWRKNjY3tj507d/bZsQEAkkC/BADQtYy+xjhixIjIz8/vtIpr165dnVZ79UZRUVEUFRX12fGA/9DceiDee9MPIiLi+VsviiGFGf0ZAMhZA+3vm34JAKBrGa3sKiwsjKqqqqipqemwvaamJqZMmdKnhQEAAABApjL+yLO6ujpmz54dkyZNismTJ8d9990XtbW1sXDhwoh4c0n9K6+8Eg8//HD7Plu3bo2IiH379sXu3btj69atUVhYGO9973v75rcAAAAAgOhB2DVr1qzYs2dP3HrrrVFfXx8TJ06M9evXx7hx4yIior6+Pmprazvsc9ZZZ7X/9+bNm+PRRx+NcePGxY4dO3pXPQAAAAC8RY8uZnHttdfGtdde2+XPVq9e3WlbOp3uycsAAAAAQEZ6dDdGAAAAAMhFwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwNWc+uBGL94XYxfvC6aWw9kuxyAPuFvGwAAA52wCwAAAIDEEHYBAAAAkBjCLgAAAAASQ9gFAAAAQGIIuwAAAABIDGEXAAAAAIkh7AIAAAAgMYRdAAAAACSGsAsAAACAxBB2AQAAAJAYwi4AALqlufVAjF+8LsYvXhfNrQccM4eON5CP2R9qPFbHpO/0hznvDzUO5GP2hxr70zF7S9gFAAAAQGIIuwAAAABIDGEXAAAAAIkh7AIAAAAgMYRdAAAAACSGsAsAAACAxBB20S/k4q1M+Q/mB3rGuZPbzA8AQP8k7AIAAAAgMYRdAAAAACSGsAsAAACAxBB2AQAAAJAYwi4AAAAAEkPYBQAAAEBiCLsAAAAASIwehV0rV66MysrKKC4ujqqqqtiwYcMRxz/zzDNRVVUVxcXFceqpp8a9997bo2LpH5pbD8T4xeti/OJ10dx6INvlAPQZf98AACD3ZRx2rVmzJhYtWhRLly6NLVu2xNSpU2PGjBlRW1vb5fjt27fHRz7ykZg6dWps2bIlbrzxxvj0pz8d3/72t3tdPAAAAAC8VcZh1/Lly2P+/PmxYMGCmDBhQqxYsSIqKipi1apVXY6/9957Y+zYsbFixYqYMGFCLFiwIK6++ur48pe/3OviAQAAAOCtMgq7WltbY/PmzTF9+vQO26dPnx4bN27scp9NmzZ1Gn/RRRfFr371q9i/f3+X+7S0tERTU1OHB8fm6zO+kgPkAn/fIHP6JQCArqXS6XS6u4Pr6uri5JNPjp/97GcxZcqU9u233357PPTQQ/Hiiy922udd73pXzJs3L2688cb2bRs3boxzzz036urqYvTo0Z32WbZsWdxyyy2dtjc2NkZJSUl3ywUA6BeampqitLQ0o15HvwQADDTd7Zl6dIH6VCrV4Xk6ne607Wjju9p+yJIlS6KxsbH9sXPnzp6UCQCQWPolAICuDcpk8IgRIyI/Pz8aGho6bN+1a1eUlZV1uU95eXmX4wcNGhQnnXRSl/sUFRVFUVFRJqUBAAwo+iUAgK5ltLKrsLAwqqqqoqampsP2mpqaDl9rfKvJkyd3Gv/DH/4wJk2aFAUFBRmWCwAAAACHl9HKroiI6urqmD17dkyaNCkmT54c9913X9TW1sbChQsj4s0l9a+88ko8/PDDERGxcOHCuPvuu6O6ujquueaa2LRpUzzwwAPx2GOPdfs1D33t0YVXAYAkOtTjZHAp1U70SwBA0nW3Z8o47Jo1a1bs2bMnbr311qivr4+JEyfG+vXrY9y4cRERUV9fH7W1te3jKysrY/369XHDDTfEPffcE2PGjImvfvWrcfnll3f7Nffu3RsRERUVFZmWCwDQb+zduzdKS0t7vG+EfgkASL6j9UwZ3Y0xWw4ePBh1dXUxdOjQI14Iv7eampqioqIidu7c6S5GOcj85Dbzk9vMT+4yN7nteM1POp2OvXv3xpgxYyIvr0f3D9IvERHmJ9eZn9xmfnKb+cltudYzZbyyKxvy8vLilFNOOW6vV1JS4uTJYeYnt5mf3GZ+cpe5yW3HY356uqLrEP0Sb2V+cpv5yW3mJ7eZn9yWKz1Tzz46BAAAAIAcJOwCAAAAIDGEXW9RVFQUN998cxQVFWW7FLpgfnKb+clt5id3mZvcZn46857kNvOT28xPbjM/uc385LZcm59+cYF6AAAAAOgOK7sAAAAASAxhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdv3VypUro7KyMoqLi6Oqqio2bNiQ7ZKIiGXLlkUqlerwKC8vz3ZZA9ZPfvKTuOSSS2LMmDGRSqXiiSee6PDzdDody5YtizFjxsTgwYPjvPPOi23btmWn2AHoaPMzb968TufThz70oewUOwDdcccd8cEPfjCGDh0ao0aNiksvvTRefPHFDmOcQ9nRnblx/vwHPVNu0jPlFj1TbtMz5S79Um7rTz2TsCsi1qxZE4sWLYqlS5fGli1bYurUqTFjxoyora3NdmlExPve976or69vf/z2t7/NdkkD1uuvvx5nnnlm3H333V3+/Itf/GIsX7487r777nj22WejvLw8pk2bFnv37j3OlQ5MR5ufiIiLL764w/m0fv3641jhwPbMM8/EddddFz//+c+jpqYmDhw4ENOnT4/XX3+9fYxzKDu6MzcRzp8IPVOu0zPlDj1TbtMz5S79Um7rVz1TmvTZZ5+dXrhwYYdt73nPe9KLFy/OUkUccvPNN6fPPPPMbJdBFyIi/Z3vfKf9+cGDB9Pl5eXpf/qnf2rf9pe//CVdWlqavvfee7NQ4cD29vlJp9PpuXPnpmfOnJmVeuhs165d6YhIP/PMM+l02jmUS94+N+m08+cQPVPu0jPlLj1TbtMz5Tb9Um7L5Z5pwK/sam1tjc2bN8f06dM7bJ8+fXps3LgxS1XxVv/+7/8eY8aMicrKyvgv/+W/xP/7f/8v2yXRhe3bt0dDQ0OHc6moqCj+9m//1rmUQ55++ukYNWpUvOtd74prrrkmdu3ale2SBqzGxsaIiBg+fHhEOIdyydvn5pCBfv7omXKfnql/8Pe+fxjof/NzhX4pt+VyzzTgw67XXnst2traoqysrMP2srKyaGhoyFJVHHLOOefEww8/HD/4wQ/i61//ejQ0NMSUKVNiz5492S6Ntzl0vjiXcteMGTPikUceiaeeeiq+8pWvxLPPPhsXXHBBtLS0ZLu0ASedTkd1dXX8f//f/xcTJ06MCOdQruhqbiKcPxF6plynZ+o//L3Pff7m5wb9Um7L9Z5p0HF9tRyWSqU6PE+n0522cfzNmDGj/b/POOOMmDx5crzzne+Mhx56KKqrq7NYGYfjXMpds2bNav/viRMnxqRJk2LcuHGxbt26uOyyy7JY2cBz/fXXx29+85v46U9/2ulnzqHsOtzcOH/+g/+N5iY9U//jXMpd/ubnBv1Sbsv1nmnAr+waMWJE5Ofnd0qBd+3a1SktJvtOOOGEOOOMM+Lf//3fs10Kb3Pojk/Opf5j9OjRMW7cOOfTcfb3f//38eSTT8aPf/zjOOWUU9q3O4ey73Bz05WBeP7omfoXPVPu8ve+/xmIf/OzTb+U2/pDzzTgw67CwsKoqqqKmpqaDttrampiypQpWaqKw2lpaYkXXnghRo8ene1SeJvKysooLy/vcC61trbGM88841zKUXv27ImdO3c6n46TdDod119/faxduzaeeuqpqKys7PBz51D2HG1uujIQzx89U/+iZ8pd/t73PwPxb3626JdyW3/qmXyNMSKqq6tj9uzZMWnSpJg8eXLcd999UVtbGwsXLsx2aQPeZz7zmbjkkkti7NixsWvXrvjCF74QTU1NMXfu3GyXNiDt27cvfv/737c/3759e2zdujWGDx8eY8eOjUWLFsXtt98ep59+epx++ulx++23x5AhQ+LKK6/MYtUDx5HmZ/jw4bFs2bK4/PLLY/To0bFjx4648cYbY8SIEfGf/tN/ymLVA8d1110Xjz76aHz3u9+NoUOHtn8iWVpaGoMHD45UKuUcypKjzc2+ffucP3+lZ8pdeqbcomfKbXqm3KVfym39qmfKxi0gc9E999yTHjduXLqwsDD9gQ98oMOtM8meWbNmpUePHp0uKChIjxkzJn3ZZZelt23blu2yBqwf//jH6Yjo9Jg7d246nX7zVsA333xzury8PF1UVJT+m7/5m/Rvf/vb7BY9gBxpfpqbm9PTp09Pjxw5Ml1QUJAeO3Zseu7cuena2tpslz1gdDU3EZH+xje+0T7GOZQdR5sb509HeqbcpGfKLXqm3KZnyl36pdzWn3qm1F8LBgAAAIB+b8BfswsAAACA5BB2AQAAAJAYwi4AAAAAEkPYBQAAAEBiCLsAAAAASAxhFwAAAACJIewCAAAAIDGEXQAAAAAkhrALAAAAgMQQdgEAAACQGMIuAAAAABLj/weIHZR0gUW0iwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1200x500 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = range(0, 26)\n",
    "q_pmf =stats.binom(10,0.75).pmf(x)\n",
    "qu_pmf =stats.randint(0,np.max(np.nonzero(q_pmf))+1).pmf(x)\n",
    "r_pmf =(q_pmf+np.roll(q_pmf,12))/2\n",
    "ru_pmf =stats.randint(0,np.max(np.nonzero(r_pmf))+1).pmf(x)\n",
    "s_pmf =(q_pmf+np.roll(q_pmf,15))/2\n",
    "su_pmf =(qu_pmf+np.roll(qu_pmf,15))/2\n",
    "\n",
    "_, ax=plt.subplots(3,2,figsize=(12,5),sharex=True,sharey=True,\n",
    "      constrained_layout=True)\n",
    "ax = np.ravel(ax)\n",
    "\n",
    "zipped = zip([q_pmf, qu_pmf,r_pmf,ru_pmf,s_pmf,su_pmf],\n",
    "         [\"q\", \"qu\", \"r\", \"ru\", \"s\", \"su\"])\n",
    "\n",
    "for idx, (dist,label) in enumerate(zipped):\n",
    "    ax[idx].vlines(x, 0,dist,label=f\"H ={stats.entropy(dist):.2f}\")\n",
    "    ax[idx].set_title(label)\n",
    "    ax[idx].legend(loc=1, handlelength=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7732bc8",
   "metadata": {},
   "source": [
    "26 possible events.\n",
    "\n",
    "- q: $q \\sim \\operatorname{binom}(n=10, p=0.75)$ (has 11 events and  the lowest entropy)\n",
    "- qu: uniform distribution with 11 events (Need to increase n 3 times for $q$ to get a higher entropy than qu).\n",
    "    \n",
    "- r: generated by duplicating $q$ and shifting to the right and then normalizing as the probabilities sum to 1.\n",
    "- ru: uniform distribution with 23 events.\n",
    "\n",
    "- s: like $r$ but with a larger valley\n",
    "- su: uniform distribution with 22 events. H of su is less than H of ru."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "def241fc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([8.67361738e-19, 7.80625564e-17, 3.39572120e-15, 9.50801937e-14,\n",
       "       1.92537392e-12, 3.00358332e-11, 3.75447915e-10, 3.86174998e-09,\n",
       "       3.33075936e-08, 2.44255686e-07, 1.53881082e-06, 8.39351359e-06,\n",
       "       3.98691895e-05, 1.65610480e-04, 6.03295319e-04, 1.93054502e-03,\n",
       "       5.42965787e-03, 1.34144488e-02, 2.90646392e-02, 5.50698426e-02,\n",
       "       9.08652404e-02, 1.29807486e-01, 1.59309188e-01, 1.66235674e-01,\n",
       "       1.45456215e-01, 1.04728475e-01])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_pmf =stats.binom(30,0.75).pmf(x)\n",
    "q_pmf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "33ac5109",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0671606564787637"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = stats.entropy(q_pmf)\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d7e50b1",
   "metadata": {},
   "source": [
    "# Kullback-Leibler Divergence"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cca0147",
   "metadata": {},
   "source": [
    "- When we do not know distribution $p$ or when $p$ is complex, we may want to approximate it with distribution $q$.\n",
    "- How much information do we lose by using $q$ to represent $p$?\n",
    "- This is equivalent to how much extra uncertainty do we introduce when using $q$ to represent $p$.\n",
    "- KL divergence gives us 0 when $p=q$ and positive values otherwise.\n",
    "- KL divergence is the average value of the difference of $\\textrm{log} p$ and $\\textrm{log} q$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f221d1a1",
   "metadata": {},
   "source": [
    "Events appear according to $p$ so compute the expectation (mean value) with respect to $p$.\n",
    "\n",
    "$$\n",
    "\\mathbb{K} \\mathbb{L}(p \\| q)=\\mathbb{E}_p[\\log p-\\log q]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "935a453d",
   "metadata": {},
   "source": [
    "For a discrete distribution, the above equation can be expressed as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eae4378",
   "metadata": {},
   "source": [
    "$$\\mathbb{K} \\mathbb{L}(p \\| q)=\\sum_i^n p_i\\left(\\log p_i-\\log q_i\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c392f823",
   "metadata": {},
   "source": [
    "Using the properties of logarithms, rewrite in this common representation of $\\mathbb{K} \\mathbb{L}(p \\| q)$\n",
    "\n",
    "$$\\mathbb{K} \\mathbb{L}(p \\| q)=\\sum_i^n p_i \\log \\frac{p_i}{q_i}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "027e578d",
   "metadata": {},
   "source": [
    "## Relate the KL divergence, $D_{\\mathrm{KL}}(p \\| q)$,  to information entropy, $H(p, q)$\n",
    "\n",
    "Using the properties of logarithms, rearrange and then expand. Note the negative sign."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5852014a",
   "metadata": {},
   "source": [
    "$$\n",
    "\\mathbb{K} \\mathbb{L}(p \\| q)=-\\sum_i^n p_i\\left(\\log q_i-\\log p_i\\right)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbb{K} \\mathbb{L}(p \\| q)=-\\overbrace{\\sum_i^n p_i \\log q_i}^{H(p, q)}-\\overbrace{\\left(-\\sum_i^n p_i \\log p_i\\right)}^{H(p)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c962d79a",
   "metadata": {},
   "source": [
    "Note that $H(p, q)$ is the entropy of q evaluated at the values of p. It is the *cross entropy*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac9a0039",
   "metadata": {},
   "source": [
    "$$\n",
    "H(p, q)=-\\mathbb{E}_p[\\log q]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c314ca7",
   "metadata": {},
   "source": [
    "Note that when $p=q$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c375d189",
   "metadata": {},
   "source": [
    "$$\\mathbb{K} \\mathbb{L}(p \\| q)=\\mathbb{K} \\mathbb{L}(p \\| p)=\\sum_i p_i\\left(\\log \\left(p_i\\right)-\\log \\left(p_i\\right)\\right)=0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc72ca7",
   "metadata": {},
   "source": [
    "Now replace $\\mathbb{K} \\mathbb{L}(p \\| q)$ with $D_{\\mathrm{KL}}(p \\| q)$ and rearrange. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4904f26b",
   "metadata": {},
   "source": [
    "$$\n",
    "H(p, q)=H(p)+D_{\\mathrm{KL}}(p \\| q)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7da71b42",
   "metadata": {},
   "source": [
    "So the KL divergence can be interpreted as the extra entropy when using $q$ to represent $p$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "067172ee",
   "metadata": {},
   "source": [
    "##  Compute $D_{\\mathrm{KL}}(p \\| q)$\n",
    "\n",
    "- using all pairwise combinations.\n",
    "- white represents infinite values  $D_{\\mathrm{KL}}(p \\| q) = \\infty$\n",
    "- black represents  $D_{\\mathrm{KL}}(p \\| q) = 0$\n",
    "- Note that the plot is not symmetric because usually $D_{\\mathrm{KL}}(p \\| q) \\neq D_{\\mathrm{KL}}(q \\| p)$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "611005a7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAGdCAYAAAAotLvzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlYUlEQVR4nO3df3RU5b3v8c9OAmMgP2iA/AKMHEEuSi5dV7yCWgwquabr6BJarsV1ilTaoxVsMRItYq+htxIXNkqvHF0t1x9gddUer1DXQYEUDPgLllLxFy5FG6+RJlelSSaJOkmYff9IGZnya2b2ntnPnnm/1nrWYnb27P3dpfjN93me/TyWbdu2AACAMbK8DgAAAEQjOQMAYBiSMwAAhiE5AwBgGJIzAACGITkDAGAYkjMAAIYhOQMAYJgcrwP4R+FwWH/961+Vn58vy7K8DgcAECfbttXd3a3y8nJlZXlXA667t0Q5OZZ+8JN2z2JIlGXaCmGffPKJxo0b53UYAACHWltbNXbsWE/u/cknn+isiacrK0v64MO/qrS01JM4EmVccu7q6tKIESN0kb6tHA3xOpyU+mPXBq9DAADHgsGgxo0bp87OThUWFnoSw4/+pVDBnrBCIVvjxuTo/v/d6UkciTKuW/tIV3aOhijHyqzkXFBQ4HUIAOAar4YmDxw4oMee6ta+7aerr9/W+TWtuuWOj3TGGWd4Ek8imBAGAEgr/+O2/6J/+U6+zjpzqKb8p4C++895qv/Zf/Y6rLiQnAEAaeONN97Qpud69fPaosixO5eN1JN/7NG7777rYWTxITkDANLGiroLdf33CzVuzNfDov9UMUQLry7QHXX/1cPI4kNyBgCkhZdeekk7X/5Cy3/yjWN+dsfNRdqyo1evvfaaB5HFj+QMAPA927Z1+7LZWvqv39DoUcfOdS4rydHi60ZoRd3FHkQXP5IzAMD3mpqa9PZ7IdXeMOKE59Td+A3t3vuVdu7cmbrAEkRyBgD4mm3bur3uSt22pEiFBdknPG9kUbaW/fgbur3uchm2xMcxSM4AAF97+umn9df/N6AbF556wZOf/GiEDvylX5s3b05BZIkjOQMAfGtgYEB3/Gy+7ri5SMOGnTql5edl6faffkMrbvuOwuFwCiJMDMkZAOBbG/5trPr6pevmx75M6L9+v1AdnWH94Q9/SGJkzpCcAQC+FAqF9ItfHdKdy4o0dGjsS4WedlqWfn5LkX5++wL19/cnMcLEkZwBAL70m1+NU35elubPyY/7u9f+9wJlZUkP/9qbXbNOheQMAPCd3t5e3bXmb/rFbSOVnR3/Bhs5OZZW3jpS//Pev+mrr75KQoTOkJwBAL6z5pcVGn/6EF3534YnfI3v/nOeRo/M1r/dPc7FyNxBcgYA+M7/+Y8eLVlU6GhbyqwsSzf9cISe+o8eFyNzB8kZAOA7ds5k5Q13nsLyhluyc6a6EJG7jl2AFAAAH7BlKyxn7yrbMnOlMJIzAMCXwnZYh21nyTls6DKeJGcAgC+FJYUdVr6mVs6MOQMAYBgqZwCALw2OODvs1ja0ciY5AwB8KSzpsMMxY1PHnOnWBgDAMFTOAABfCstO2wlhJGcAgC+FZeuww+Rq6o7OSenW7u3t1YIFC5SXl6eysjI1NjaqqqpKS5cuTcbtAAAZyNbX1XOizdTKOSnJua6uTs8//7w2btyobdu2qbm5WXv37j3uuaFQSMFgMKoBAGCihoYGnXfeecrPz1dxcbGuuuoqvffee1HnLFy4UJZlRbXp06fHdR/Xk3NPT48eeugh/epXv9Ls2bNVWVmp9evX6/Dhw8c9v6GhQYWFhZE2bpx5u4MAAMwTtm0ddtjina29c+dOLV68WLt371ZTU5MGBgZUXV2t3t7eqPMuv/xytbW1Rdqzzz4b131cH3P+8MMP1dfXpxkzZkSOFRUVadKkScc9f/ny5aqtrY18DgaDJGgAwCkNdms7E+/3t2zZEvX5kUceUXFxsfbu3auZM2dGjgcCAZWWliYcl+vJ2Y7zt5BAIKBAIOB2GACANHfYpQlhtm0fM6Qaa27q6uqSNFiEHq25uVnFxcUaMWKELr74Yt11110qLi6OOS7Xu7UnTJigIUOGaPfu3ZFjHR0dev/9992+FQAAjh08eDBqeLWwsFANDQ2n/J5t26qtrdVFF12kKVOmRI7X1NTo8ccf144dO9TY2KhXX31Vl1xyiUKhUMwxuV455+XladGiRaqrq9PIkSNVUlKiFStWKCuL9U4AAO4J29Jhh5Otw5LGjBmj/fv3Rx2PpWpesmSJ3nzzTb344otRx6+++urIn6dMmaJp06apoqJCmzdv1ty5c2OKKynvOd9zzz3q6enRlVdeqfz8fN1yyy2R0h8AADe4MeZsS7IsSwUFBXF976abbtIzzzyjXbt2aezYsSc9t6ysTBUVFTpw4EDM109KOZuXl6fHHntMvb29am9vV11dXTJuAwBAStm2rSVLlujpp5/Wjh07NH78+FN+59ChQ2ptbVVZWVnM96GvGQDgS2FZOuywxVt5L168WL/73e/0xBNPKD8/X+3t7Wpvb9eXX34pafB14mXLlumVV17RRx99pObmZl1xxRUaNWqU5syZE/N9WL4TAOBLtj047uxEvN9/8MEHJUlVVVVRxx955BEtXLhQ2dnZeuutt7RhwwZ1dnaqrKxMs2bN0pNPPqn8/PyY75Oy5Nzc3JyqWwEAMsCR6teJcJzfP9Xrwrm5udq6dauTkCTRrQ0AgHHo1gYA+FJYSnnlnCokZwCAL9myFLadJdc4F7VMGZIzAMCXwh6MOacKY84AABiGyhkA4EuDlbOzGtPUypnkDADwpbDtwpizocmZbm0AAAxD5QwA8CVepQIAwDBhO0uHbYdjzg67xZOF5AwA8CVblsJpOiGMMWcAAAxD5QwA8KV0XoSE5AwA8KWwLMacAQAwSViW48qX95wBAEBMqJwBAL7E8p0AABjGlfecDe1ANjMqAAAyGJUzAMCX3FiExNQJYSRnAIAvhW1Lhx2+CsWrVHHa+P5bKsjPrF732VnzvA4h5bLPPsvrEDxhf/SJ1yF4wjpjrNchpJz9l4+9DiHlBuy+lNwnnSeEZVb2AwDAB4ytnAEAOJmwshRmhTAAAMzhTre2mR3IJGcAgC+5MiHMpVjcZuavDAAAZDAqZwCAL7nznrOZNSrJGQDgS4Pd2uk5IczMXxkAAMhgVM4AAF9K5/2cSc4AAF8Ky/muVE5neycLyRkA4Ess3wkAAFKGyhkA4Eu2bTmebc2YMwAALkrnbm2SMwDAl8K25XjjC9vh95PFzKgAAMhgVM4AAF8a7NZOz40vSM4AAF+ybTf2czazA9nMqAAAyGBUzgAAXwpLdGsDAGCSsJx3a7NlJAAALmLLSAAAkDJUzgAAX7LZMhIAALOkc7c2yRkA4EthOd/4wtS1tRlzBgDAMFTOAABfSuddqeJ+qt7eXi1YsEB5eXkqKytTY2OjqqqqtHTpUkmSZVnatGlT1HdGjBihRx991IVwAQAYdGQ/ZyfNNnTMOe7kXFdXp+eff14bN27Utm3b1NzcrL179yYcQCgUUjAYjGoAAGSyuLq1e3p69NBDD2nDhg2aPXu2JGn9+vUaO3ZswgE0NDRo5cqVCX8fAJCZwrIUdtitPfgqle1OQC6K66k+/PBD9fX1acaMGZFjRUVFmjRpUsIBLF++XF1dXZHW2tqa8LUAAJlj8FUqZy0tXqWy7VP/dmFZ1jHn9ff3n/D8QCCgQCAQTxgAAAwuQsKrVNKECRM0ZMgQ7d69O3Kso6ND77//fuTz6NGj1dbWFvl84MABffHFFy6ECgBAZogrOefl5WnRokWqq6vT9u3b9fbbb2vhwoXKyvr6MpdcconWrl2rP//5z3rttdd0ww03aMiQIa4HDgDIbIOztbMctXhnazc0NOi8885Tfn6+iouLddVVV+m99977h7hs1dfXq7y8XLm5uaqqqtI777wT133iHkm/5557NHPmTF155ZW67LLLdNFFF+ncc8+N/LyxsVHjxo3TzJkzdc0112jZsmUaNmxYvLcBAOCkDsty3OLt1t65c6cWL16s3bt3q6mpSQMDA6qurlZvb2/knNWrV+vee+/V2rVr9eqrr6q0tFSzZ89Wd3d3zPeJexGSvLw8PfbYY3rssccixzZv3hz5c3l5ubZu3Rr1nc7OznhvAwDASdm287WxY5hKFWXLli1Rnx955BEVFxdr7969mjlzpmzb1po1a7RixQrNnTtX0uBbTSUlJXriiSd0/fXXx3Qflu8EAGQ027aPWW8jFArF9N2uri5Jg28uSVJLS4va29tVXV0dOScQCOjiiy/Wyy+/HHNMJGcAgC/ZDsebw3aWwsrSwYMHVVhYGNUaGhpiuL+t2tpaXXTRRZoyZYokqb29XZJUUlISdW5JSUnkZ7FwZW3t5uZmNy4DAEDMwnL+KpQtacyYMdq/f3/U8Vhe8V2yZInefPNNvfjii8f8zLKi47Jt+5hjJ8PGFwCAjGZZlgoKCuL6zk033aRnnnlGu3btilols7S0VNJgBV1WVhY5/umnnx5TTZ8M3doAAF8KK/UrhNm2rSVLlujpp5/Wjh07NH78+Kifjx8/XqWlpWpqaooc6+vr086dO3XBBRfEfB8qZwCALx0Zc3Z2jfiS8+LFi/XEE0/oj3/8o/Lz8yPjyIWFhcrNzZVlWVq6dKlWrVqliRMnauLEiVq1apWGDRuma665Jub7kJwBAL4UduFVqnjHrB988EFJUlVVVdTxRx55RAsXLpQk3Xrrrfryyy914403qqOjQ+eff762bdum/Pz8mO9DcgYAIEax7jFRX1+v+vr6hO9DcgYA+JKdwApfx7uGiUjOAABfcmNXqnhXCEsVkjMAwJfCf9/4wtk1zHxpycyoAADIYFTOAABfshN4T/mYa7gUi9tIzgAAX3Jn+U4zJ4TRrQ0AgGGonAEAvuTGbO2wbRlZpZKcAQC+5MqYs6GDziRnAIAvhV1IzmFZynYpHjeZWM0DAJDRqJwBAL5ky/nGF/HuSpUqJGcAgC+5s7a2mUjOAABfcmXM2dDKmTFnAAAMY2zlPON//VDZgdO8DiOl3g7f7HUIKTc7a57XIXgia9gwr0PwxJa37/I6BKRAMBhUYeEfkn4f23ZhzNnQFcKMTc4AAJyMG93apr7nTLc2AACGoXIGAPiSW8t3mojkDADwJdt24z1lkjMAAK5x4z3nsEuxuI0xZwAADEPlDADwJXdma9OtDQCAq5wmV5IzAAAucqVydikWtzHmDACAYaicAQD+ZFt0awMAYJKwC4uQmJqc6dYGAMAwVM4AAH+ynW9cYeqEMJIzAMCXwpLjFcJIzgAAuCmNJ4Qx5gwAgGGonAEAvpTOs7VJzgAAX7KZEAYAgGHc2M/Z0OzMmDMAAIahcgYA+JKdxrO1Sc4AAF+y3ZgQ5vA96WShWxsAAMNQOQMAfMmV2dqGTggjOQMAfMl2YbY2yRkAAFc5nxAmxpwBAEAsUpqc+/r6Unk7AEAas11qJkpqt3ZVVZWmTJmioUOHasOGDTrnnHO0c+fOZN4SAJAheM/ZgfXr1+vHP/6xXnrpJdnHGXkPhUIKhUKRz8FgMNkhAQDSgRulr6Glc9KT84QJE7R69eoT/ryhoUErV65MdhgAAPhG0secp02bdtKfL1++XF1dXZHW2tqa7JAAAGniSNd2os1USa+chw8fftKfBwIBBQKBZIcBAEgz6bwICa9SAQBgGBYhAQD4ErO1AQAwjW0NNkcyMDk3Nzcn8/IAgAxmy4UxY8acAQBALOjWBgD4E4uQAABgGDcmhLkUitvo1gYA+JcHO1/s2rVLV1xxhcrLy2VZljZt2hT184ULF8qyrKg2ffr0uO5BcgYAIA69vb2aOnWq1q5de8JzLr/8crW1tUXas88+G9c96NYGAPjS4AphqX/PuaamRjU1NSc9JxAIqLS0NNGwqJwBAD7lUpe2bdsKBoNR7ejdEhPR3Nys4uJinXXWWfrRj36kTz/9NK7vk5wBABnt4MGDKiwsjGoNDQ0JX6+mpkaPP/64duzYocbGRr366qu65JJL4kr4dGsDAHzKkvMVviyNGTNG+/fvjzrqZEOmq6++OvLnKVOmaNq0aaqoqNDmzZs1d+7cmK5BcgYA+JcL7zlblqWCggJXwjmesrIyVVRU6MCBAzF/h+QMAPAnnyxCcujQIbW2tqqsrCzm75CcAQCIQ09Pjz744IPI55aWFu3bt09FRUUqKipSfX29vvOd76isrEwfffSRbr/9do0aNUpz5syJ+R4kZwCAP7myK1X8XnvtNc2aNSvyuba2VpJ07bXX6sEHH9Rbb72lDRs2qLOzU2VlZZo1a5aefPJJ5efnx3wPkjMAwJe82pWqqqpK9kluvHXrVgcBDSI5AwD8ySdjzongPWcAAAxD5QwA8CdbLow5p37MOhYkZwCAL1n2YHOEbm0AABALKmcAgH+laeVMcgYA+JMb7zl78J50LEjOAAB/cuNVKkMx5gwAgGGMrZx7JgwoK3fA6zCQZE3hf/c6BE/MzprndQhAemDMGQAAg6TxCmEkZwCAP3m08UUqMOYMAIBhqJwBAL5kyfkKYY5XGEsSkjMAwJ/SeMyZbm0AAAxDcgYAwDB0awMA/MmFXalMnetNcgYA+Fearq1NtzYAAIahcgYA+FMaz9YmOQMA/CmNd6UiOQMAfMlyYUKYqcmdMWcAAAxD5QwA8K80rZxJzgAAf0rjMWe6tQEAMAyVMwDAl9yYEMauVAAAuMpK2xXCSM4AAH9izBkAAKQKlTMAwJcYcwYAwES85wwAgEHcWL7TUIw5AwBgGCpnAIA/sWUkAACGSePk7Fq3dl9fn1uXAgAgoyWcnKuqqrRkyRLV1tZq1KhRmjhxoizL0r59+yLndHZ2yrIsNTc3uxAqAABfs/T161QJN68f4gQcVc7r169XTk6OXnrpJW3dujWha4RCIQWDwagGAEAmc5ScJ0yYoNWrV2vSpEk67bTTErpGQ0ODCgsLI23cuHFOQgIAZArbpWYgR8l52rRpjgNYvny5urq6Iq21tdXxNQEA8DNHs7WHDx8e+XNW1mCet+2vfw3p7+8/5TUCgYACgYCTMAAAGSidl+90bbb26NGjJUltbW2RY0dPDgMAwHVp2KUtufiec25urqZPn667775bZ5xxhj7//HPdcccdbl0eAIBovOccm4cfflj9/f2aNm2afvrTn+qXv/ylm5cHACAjJFw5H+/d5cmTJ+uVV16JOnb0GDQAAK5xY+MLQ1MUy3cCAPzL0OTqFLtSAQBgGCpnAIAvpfOrVCRnAIA/pfFsbZIzAMC/DE2uTjHmDACAYaicAQC+xJgzAACmYcwZAADDGL4+thOMOQMAEIddu3bpiiuuUHl5uSzL0qZNm6J+btu26uvrVV5ertzcXFVVVemdd96J6x4kZwCAbx0Zd060JVJ59/b2aurUqVq7du1xf7569Wrde++9Wrt2rV599VWVlpZq9uzZ6u7ujvkedGsDAPzJozHnmpoa1dTUHP9ytq01a9ZoxYoVmjt3riRp/fr1Kikp0RNPPKHrr78+pntQOQMAMppt2woGg1EtFAoldK2Wlha1t7eruro6ciwQCOjiiy/Wyy+/HPN1SM4AAF9y2qVt2ZIl6eDBgyosLIxqDQ0NCcXU3t4uSSopKYk6XlJSEvlZLOjWBgD4lwvd2mPGjNH+/fujDgcCAUeXtSwr+ja2fcyxkyE5AwD8yaUxZ8uyVFBQ4EZEKi0tlTRYQZeVlUWOf/rpp8dU0ydDtzYAAC4ZP368SktL1dTUFDnW19ennTt36oILLoj5OlTOAABfsv7enF4jXj09Pfrggw8in1taWrRv3z4VFRXp9NNP19KlS7Vq1SpNnDhREydO1KpVqzRs2DBdc801Md+D5AwA8CePVgh77bXXNGvWrMjn2tpaSdK1116rRx99VLfeequ+/PJL3Xjjjero6ND555+vbdu2KT8/P+Z7kJwBAL7leOOKBL5fVVUl2z7xFy3LUn19verr6xMOizFnAAAMQ+UMAPAndqUCAMAw7EoFAABSxdjKufCdHGUPNTY8wJGm8L97HYInZmfN8zqElMs67TSvQ0i5AbsvJfex5HxCmOMJZUlC9gMA+BNjzgAAmCWyJ7OTa7gTiusYcwYAwDBUzgAAf6JbGwAAs7gxIYzkDACAm3jPGQAApAqVMwDAnxhzBgDAPOm6CAnd2gAAGIbKGQDgT2k8IYzkDADwJcu2ZdnOsrPT7ycLyRkA4F9pOiGMMWcAAAxD5QwA8CU3Nr4wFckZAOBPvOcMAIBh3Ngy0tDkzJgzAACGoXIGAPiXoZWvUyRnAIAvuTEhjG5tAAAQEypnAIA/MVsbAACzWDK3W9opkjMAwJ9sSU7XxjZ0bW3GnAEAMAyVMwDAn9J4tjbJGQDgX4YmV6dIzgAAX7LCg83pNUyUtDHnp556SpWVlcrNzdXIkSN12WWXqbe3N1m3AwAgbSSlcm5ra9P8+fO1evVqzZkzR93d3XrhhRdkH2dWXCgUUigUinwOBoPJCAkAkI7o1o5dW1ubBgYGNHfuXFVUVEiSKisrj3tuQ0ODVq5cmYwwAABpzJX9nA1N7knp1p46daouvfRSVVZWat68eVq3bp06OjqOe+7y5cvV1dUVaa2trckICQAA30hKcs7OzlZTU5Oee+45nX322br//vs1adIktbS0HHNuIBBQQUFBVAMA4JRs24Xm9UMcX9ImhFmWpQsvvFArV67U66+/rqFDh2rjxo3Juh0AINPYX3dtJ9y8foYTSMqY8549e7R9+3ZVV1eruLhYe/bs0WeffabJkycn43YAgEzlxpizgRk6Kcm5oKBAu3bt0po1axQMBlVRUaHGxkbV1NQk43YAAKSVpCTnyZMna8uWLcm4NAAAktyZrW1lUuUMAEDSHZnU5ewiroTiNpIzAMCXXNnP2czczJaRAACYhsoZAOBPttyZrW0gkjMAwJdcmxBmILq1AQAwDJUzAMCfbElhh6Wv49neyUFyBgD4E2POAACYx/GYszthuI4xZwAADEPlDADwJzdWCKNbGwAA97jxKhXJGQAAN7kxIcxQjDkDAGAYkjMAwJcs2bJs5y0e9fX1siwrqpWWlrr+bHRrAwD8Kfz35kQC3eLnnHOO/vSnP0U+Z2dnOwziWCRnAADikJOTk5Rq+Wh0awMAfMuNbm3bthUMBqNaKBQ64T0PHDig8vJyjR8/Xt/73vf0l7/8xfXnIjkDAPzJdqcdPHhQhYWFUa2hoeG4tzz//PO1YcMGbd26VevWrVN7e7suuOACHTp0yNVHo1sbAOBPbixCImnMmDHav39/1LFAIHDcc2tqaiJ/rqys1IwZM3TmmWdq/fr1qq2tdRzLESRnAEBGsyxLBQUFCX13+PDhqqys1IEDB1yNiW5tAIAvWfp6lbBEm9NFTEKhkN59912VlZW58UgRJGcAgD8d6dZ22uKwbNky7dy5Uy0tLdqzZ4+++93vKhgM6tprr3X10Yzt1u7+p7CyTnP6AhsAk+SMr/A6hJQbaPm/XoeQcmG7PyX3scKDzdE14qycP/nkE82fP1+ff/65Ro8erenTp2v37t2qqHD3/9vGJmcAAEzz+9//PiX3ITkDAHzKjS0jzdw5g+QMAPAndqUCAACpQuUMAPAnW3HvKvWP4p0QliokZwCATzHmDACAWTzaMjIVGHMGAMAwVM4AAF+yZDsec6ZbGwAAN9kyNrk6RXIGAPiTG1tGGprcGXMGAMAwVM4AAH+ylbaztUnOAABfsmznE8IcTyhLErq1AQAwDJUzAMCf3JitbWbhTHIGAPiVC7O1Dc3OJGcAgD/xKhUAAEgVKmcAgD/xKhUAAGax2M8ZAADDMOYMAABShcoZAOBPti2F07NyJjkDAPzJlW5td0JxG93aAAAYhsoZAOBfrBAGAIBBmK0d7amnnlJlZaVyc3M1cuRIXXbZZert7VVVVZWWLl0ade5VV12lhQsXuhAqAABHCdvOm5m5Of7Kua2tTfPnz9fq1as1Z84cdXd364UXXpCd4G8foVBIoVAo8jkYDCZ0HQAA0kVCyXlgYEBz585VRUWFJKmysjLhABoaGrRy5cqEvw8AyFS2ZDtcv9Pp95Mk7m7tqVOn6tJLL1VlZaXmzZundevWqaOjI+EAli9frq6urkhrbW1N+FoAgAxyZMzZUfP6IY4v7uScnZ2tpqYmPffcczr77LN1//33a9KkSWppaVFWVtYx3dv9/f0nvV4gEFBBQUFUAwDglGy5MOZsZnZOaEKYZVm68MILtXLlSr3++usaOnSoNm7cqNGjR6utrS1y3uHDh/X222+7FiwAAJkg7jHnPXv2aPv27aqurlZxcbH27Nmjzz77TJMnT9bw4cNVW1urzZs368wzz9R9992nzs7OJIQNAMh4afwqVdzJuaCgQLt27dKaNWsUDAZVUVGhxsZG1dTUqL+/X2+88YYWLFignJwc3XzzzZo1a1Yy4gYAZDo3krOh4k7OkydP1pYtW477syFDhuiBBx7QAw884DgwAAAyFSuEAQD8i25tAAAMEg4PNidIzgAAuCiNJ4SxZSQAAIahcgYA+FeaVs4kZwCAPx1Z5csJM3MzyRkA4E+2bNkON65w+v1kYcwZAADDUDkDAPyJbm0AAEzDq1QAACBFqJwBAP4UtlkhDAAAs6RvtzbJGQDgS3bYlu2wcrYNTc6MOQMAYBgqZwCAP6XxxhckZwCAP9luvOdMcgYAwD22LTldftPQ5MyYMwAAhqFyBgD4km3bstO0W5vKGQDgT0e6tR21xJLzAw88oPHjx+u0007TueeeqxdeeMHVRyM5AwAQhyeffFJLly7VihUr9Prrr+tb3/qWampq9PHHH7t2D5IzAMCf/t6t7aglUDnfe++9WrRokX74wx9q8uTJWrNmjcaNG6cHH3zQtUczbsz5yP9Q4a++8jiS1AsGg16HACTVQDjkdQgpN2D3ex1Cyg1o8JmTvfpW0P6bsmzL0TW61aFcu+CY//4GAgEFAoFjzu/r69PevXv1s5/9LOp4dXW1Xn75ZUexHM245Nzd3S1J+mTlLz2OJPUKl9/hdQgA4Jru7m4VFhYm5drf/va3talvk0L6m6PrDFOWcnNzj4nzzjvvVH19/THnf/755zp8+LBKSkqijpeUlKi9vd1RLEczLjmXl5ertbVV+fn5sixnvxHFKxgMaty4cWptbVVBQUFK7+2lTHzuTHxmiefOpOf28plt21Z3d7fKy8uTdo+77rpLd911lyvXCoVCCoWie3WOVzUf7R/zk23bruYs45JzVlaWxo4d62kMBQUFGfMP+GiZ+NyZ+MwSz51JvHrmZFXMyXCiLuzjGTVqlLKzs4+pkj/99NNjqmknmBAGAECMhg4dqnPPPVdNTU1Rx5uamnTBBRe4dh/jKmcAAExWW1ur73//+5o2bZpmzJih3/72t/r44491ww03uHYPkvNRAoGA7rzzzpi7N9JFJj53Jj6zxHNn0nNn4jOnytVXX61Dhw7pF7/4hdra2jRlyhQ9++yzqqiocO0elm3qTtMAAGQoxpwBADAMyRkAAMOQnAEAMAzJGQAAw5CcAQAwDMkZ0OBi9ukuE54RSBckZ0m9vb1asGCB8vLyVFZWpsbGRlVVVWnp0qVeh+aqUz2nZVnatGlT1HdGjBihRx99NOWxJltVVZWWLFmi2tpajRo1SrNnz/Y6JNf94zNOnDhRlmVp3759kXM6OztlWZaam5s9izPZnnrqKVVWVio3N1cjR47UZZddpt7eXq/DctWJnvF4/x276qqrtHDhQk/iROxIzpLq6ur0/PPPa+PGjdq2bZuam5u1d+9er8NyXaY8Z6zWr1+vnJwcvfTSS/rNb37jdThJcfQzbt261etwUq6trU3z58/Xddddp3fffVfNzc2aO3du0rcyTKVMeMZMlPErhPX09Oihhx7Shg0bItXT+vXrPd98w22Z8pzxmDBhglavXu11GEl19DN+9NFH3gbjgba2Ng0MDGju3LmR1ZsqKys9jspdmfCMmSjjK+cPP/xQfX19mjFjRuRYUVGRJk2a5GFU7suU54zHtGnTvA4h6TLhGU9m6tSpuvTSS1VZWal58+Zp3bp16ujo8DosV2XCM2aijE/OmdL1E8tzWpZ1zHn9/f3JCslzw4cP9zqEpDv6GbOyBv+5H/13nM5/v5KUnZ2tpqYmPffcczr77LN1//33a9KkSWppafE6NNec7BmzsrIy6t90Osn45DxhwgQNGTJEu3fvjhzr6OjQ+++/72FU7ovlOUePHq22trbI5wMHDuiLL75IaZxIntGjR0tS1N/x0ZPD0pVlWbrwwgu1cuVKvf766xo6dKg2btzodViuOtEz/uO/6cOHD+vtt9/2MFLEKuPHnPPy8rRo0SLV1dVp5MiRKikp0YoVKyJVRrqI5TkvueQSrV27VtOnT1c4HNZtt92mIUOGeBg13JSbm6vp06fr7rvv1hlnnKHPP/9cd9xxh9dhJdWePXu0fft2VVdXq7i4WHv27NFnn32myZMnex2aa072jMOHD1dtba02b96sM888U/fdd586Ozu9DhkxyPjkLEn33HOPenp6dOWVVyo/P1+33HKLurq6vA7Ldad6zsbGRv3gBz/QzJkzVV5erl//+tcZPZs7HT388MO67rrrNG3aNE2aNEmrV69WdXW112ElTUFBgXbt2qU1a9YoGAyqoqJCjY2Nqqmp8To015zsGfv7+/XGG29owYIFysnJ0c0336xZs2Z5HTJiwJaRJ1BVVaVvfvObWrNmjdehJFWmPCcA+El69d0CAJAGSM4AABiGbm0AAAxD5QwAgGFIzgAAGIbkDACAYUjOAAAYhuQMAIBhSM4AABiG5AwAgGFIzgAAGIbkDACAYf4/Z7pCbLNhqgcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dists = [q_pmf, qu_pmf, r_pmf, ru_pmf, s_pmf, su_pmf]\n",
    "names = [\"q\", \"qu\", \"r\", \"ru\", \"s\", \"su\"]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "KL_matrix = np.zeros((6, 6))\n",
    "for i, dist_i in enumerate(dists):\n",
    "    for j, dist_j in enumerate(dists):\n",
    "        KL_matrix[i, j] = stats.entropy(dist_i, dist_j)\n",
    "\n",
    "ax.set_xticks(np.arange(len(names)))\n",
    "ax.set_yticks(np.arange(len(names)))\n",
    "ax.set_xticklabels(names)\n",
    "ax.set_yticklabels(names)\n",
    "plt.set_cmap(\"viridis\")\n",
    "cmap = plt.cm.get_cmap()\n",
    "cmap.set_bad('w', 0.3)\n",
    "im = ax.imshow(KL_matrix)\n",
    "fig.colorbar(im, extend=\"max\");\n",
    "\n",
    "plt.savefig(\"./SectionD10/images/chp11_KL_heatmap.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "777f98b4",
   "metadata": {},
   "source": [
    "- The above is the KL divergence for all pairwise combinations of the distributions q, qu,r,ru,s,and su.\n",
    "\n",
    "- The dark blue corresponds to KL divergence of 0.\n",
    "\n",
    "- The white corresponds to KL divergence of 1.\n",
    "\n",
    "- Note that the matrix is not symmetric.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77156d4e",
   "metadata": {},
   "source": [
    "### Note that the computation of $D_{\\mathrm{KL}}(p \\| q)$ uses the following conventions:\n",
    "\n",
    "$$\n",
    "0 \\log \\frac{0}{0}=0, \\quad 0 \\log \\frac{0}{q(x)}=0, \\quad p(x) \\log \\frac{p(x)}{0}=\\infty\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fffab56",
   "metadata": {},
   "source": [
    "### Ranking of models by KL-divergence is equivalent to ranking on the log-score\n",
    "\n",
    "\n",
    "#### Section 2.5 model comparison\n",
    "\n",
    "- do not want models that are so simple that miss important information\n",
    "- do not want models that are so complex that they fit noise\n",
    "- The sweet spot is hard to find\n",
    "  - no single criterion defines an optimal solution \n",
    "  - the optimal solution may not exist\n",
    "  - must choose from a limited set of models evaluated over a finite dataset\n",
    "\n",
    "\n",
    "- One approach is to compute the generalization error (i.e., out-of-sample predictive accuracy). Compute by cross validation.\n",
    "\n",
    "- logarithmic scoring rules have nice [theoretical properties](https://doi.org/10.1198/016214506000001437)\n",
    "\n",
    "\n",
    "$$\\mathrm{ELPD}=\\sum_{i=1}^n \\int p_t\\left(\\tilde{y}_i\\right) \\log p\\left(\\tilde{y}_i \\mid y_i\\right) d \\tilde{y}_i$$\n",
    "\n",
    "- ELPD: expected log *pointwise* predictive density\n",
    "- \"expected\" because we integrate over all possible datasets that can be generated by the process\n",
    "- \"pointwise\" because the computations are performed per observation $y_i$ over $n$ observations.\n",
    "\n",
    "- $p_t\\left(\\tilde{y}_i\\right)$: is the distribution of the true data-generating process for $\\tilde{y}_i$ (the future data given the current posterior which in turns depends on the model (prior and likelihood) and the current data)\n",
    "- $p\\left(\\tilde{y}_i \\mid y_i\\right)$: posterior predictive distribution $p(\\tilde{{y}} \\mid {y})=\\int_{{\\Theta}} p(\\tilde{{y}} \\mid {\\theta}) p({\\theta} \\mid {y}) d {\\theta}$\n",
    "\n",
    "Use the *deviance* instead:\n",
    "\n",
    "\n",
    "$$\\sum_{i=1}^n \\log \\int p\\left(y_i \\mid \\boldsymbol{\\theta}\\right) p(\\boldsymbol{\\theta} \\mid y) d \\boldsymbol{\\theta}$$\n",
    "\n",
    "- Note that this quantity is proportional to the quadratic mean error when the likelihood is Gaussian.\n",
    "- Note that we are using all of the data to fit the model, so we will over estimate ELPD."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a4cc6e7",
   "metadata": {},
   "source": [
    "##### Leave-One-Out Cross Validation\n",
    "\n",
    "\n",
    "$$\\mathrm{ELPD}_{\\text {LOO-CV }}=\\sum_{i=1}^n \\log \\int p\\left(y_i \\mid \\boldsymbol{\\theta}\\right) p\\left(\\boldsymbol{\\theta} \\mid y_{-i}\\right) d \\boldsymbol{\\theta}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "748bf9d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [σ]\n",
      "/Users/blaine/opt/anaconda3/envs/stat-rethink2-pymc_v4/lib/python3.10/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.5\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n",
      "/Users/blaine/opt/anaconda3/envs/stat-rethink2-pymc_v4/lib/python3.10/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.5\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n",
      "/Users/blaine/opt/anaconda3/envs/stat-rethink2-pymc_v4/lib/python3.10/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.5\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n",
      "/Users/blaine/opt/anaconda3/envs/stat-rethink2-pymc_v4/lib/python3.10/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.5\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:01&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 17 seconds.\n",
      "Sampling: [y]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [4000/4000 00:00&lt;00:00]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "KeyError",
     "evalue": "'y'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[49], line 22\u001b[0m\n\u001b[1;32m     20\u001b[0m     y \u001b[38;5;241m=\u001b[39m pm\u001b[38;5;241m.\u001b[39mNormal(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124my\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;241m0\u001b[39m, σ, observed\u001b[38;5;241m=\u001b[39my_obs)\n\u001b[1;32m     21\u001b[0m     idataB \u001b[38;5;241m=\u001b[39m pm\u001b[38;5;241m.\u001b[39msample(return_inferencedata\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[0;32m---> 22\u001b[0m     idataB\u001b[38;5;241m.\u001b[39madd_groups({\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mposterior_predictive\u001b[39m\u001b[38;5;124m\"\u001b[39m: {\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124my\u001b[39m\u001b[38;5;124m\"\u001b[39m:\u001b[43mpm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msample_posterior_predictive\u001b[49m\u001b[43m(\u001b[49m\u001b[43midataB\u001b[49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43my\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m[\u001b[38;5;28;01mNone\u001b[39;00m,:]}})\n\u001b[1;32m     23\u001b[0m     idatas_cmp[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmB\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m idataB\n\u001b[1;32m     25\u001b[0m \u001b[38;5;66;03m# Generate data from Normal likelihood model\u001b[39;00m\n\u001b[1;32m     26\u001b[0m \u001b[38;5;66;03m# with random mean and random standard deviation\u001b[39;00m\n",
      "File \u001b[0;32m~/opt/anaconda3/envs/stat-rethink2-pymc_v4/lib/python3.10/site-packages/arviz/data/inference_data.py:259\u001b[0m, in \u001b[0;36mInferenceData.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m    257\u001b[0m \u001b[38;5;124;03m\"\"\"Get item by key.\"\"\"\u001b[39;00m\n\u001b[1;32m    258\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m key \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_groups_all:\n\u001b[0;32m--> 259\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m(key)\n\u001b[1;32m    260\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mgetattr\u001b[39m(\u001b[38;5;28mself\u001b[39m, key)\n",
      "\u001b[0;31mKeyError\u001b[0m: 'y'"
     ]
    }
   ],
   "source": [
    "np.random.seed(90210)\n",
    "\n",
    "y_obs =  np.random.normal(0, 1, size=100)\n",
    "idatas_cmp = {}\n",
    "\n",
    "# # Generate data from Skewnormal likelihood model\n",
    "# # with fixed mean and skewness and random standard deviation\n",
    "# with pm.Model() as mA:\n",
    "#     σ = pm.HalfNormal(\"σ\", 1)\n",
    "#     y = pm.SkewNormal(\"y\", 0, σ, alpha=1, observed=y_obs)\n",
    "#     idataA = pm.sample(return_inferencedata=True)\n",
    "# # add_groups modifies an existing az. InferenceData   \n",
    "#     idataA.add_groups({\"posterior_predictive\": {\"y\":pm.sample_posterior_predictive(idataA)[\"y\"][None,:]}})\n",
    "#     idatas_cmp[\"mA\"] = idataA\n",
    "\n",
    "# Generate data from Normal likelihood model\n",
    "# with fixed mean and random standard deviation\n",
    "with pm.Model() as mB:\n",
    "    σ = pm.HalfNormal(\"σ\", 1)\n",
    "    y = pm.Normal(\"y\", 0, σ, observed=y_obs)\n",
    "    idataB = pm.sample(return_inferencedata=True)\n",
    "    idataB.add_groups({\"posterior_predictive\": {\"y\":pm.sample_posterior_predictive(idataB)[\"y\"][None,:]}})\n",
    "    idatas_cmp[\"mB\"] = idataB\n",
    "\n",
    "# Generate data from Normal likelihood model\n",
    "# with random mean and random standard deviation\n",
    "with pm.Model() as mC:\n",
    "    μ = pm.Normal(\"μ\", 0, 1)\n",
    "    σ = pm.HalfNormal(\"σ\", 1)\n",
    "    y = pm.Normal(\"y\", μ, σ, observed=y_obs)\n",
    "    idataC = pm.sample(return_inferencedata=True)\n",
    "    idataC.add_groups({\"posterior_predictive\": {\"y\":pm.sample_posterior_predictive(idataC)[\"y\"][None,:]}})\n",
    "    idatas_cmp[\"mC\"] = idataC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a074b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmp = az.compare(idatas_cmp)\n",
    "cmp.round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "b05437dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"./SectionD10/images/table.png\" width=\"900\" height=\"900\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url=\"./SectionD10/images/table.png\", width=900, height=900)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7aa4744e",
   "metadata": {},
   "source": [
    "- rank: 0 is the model with the best predictive accuracy\n",
    "- loo: the LELPD value\n",
    "- p_loo: values of the penalization term. Roughly the effective number of parameters.\n",
    "- d_loo: difference n LOO for the top-ranked model and the LOO value for each model\n",
    "- weight: probability of each model given the data\n",
    "- se: the standard error of the ELPD\n",
    "- dse: the standard errors of the difference between two values of ELPD. Usually 0 for the top-ranked model\n",
    "- warning: If True, the LOO approximation may not be reliable\n",
    "- loo_scale: log, deviance, and negative-log are the options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "883eceea",
   "metadata": {},
   "outputs": [],
   "source": [
    "az.plot_compare(cmp, figsize=(9, 3))\n",
    "plt.savefig(\"img/chp02/compare_dummy.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "6667f1d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"./SectionD10/images/compareModelsbyLOO.png\" width=\"900\" height=\"900\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url=\"./SectionD10/images/compareModelsbyLOO.png\", width=900, height=900)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc65c3e7",
   "metadata": {},
   "source": [
    "#### Pointwise ELPD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "d2c81ab1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"./SectionD10/images/pointwiseELPD.png\" width=\"900\" height=\"900\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url=\"./SectionD10/images/pointwiseELPD.png\", width=900, height=900)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e46ad9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "az.plot_elpd(idatas_cmp, figsize=(10, 5), plot_kwargs={\"marker\":\".\"}, threshold=2);\n",
    "plt.savefig(\"img/chp02/elpd_dummy.png\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3df09d87",
   "metadata": {},
   "source": [
    "#### Back to section 11.3\n",
    "\n",
    "- Assume that the true model is $M_0$\n",
    "- Assume $k$ model posteriors $\\left\\{q_{M_1}, q_{M_2}, \\cdots q_{M_k}\\right\\}$\n",
    "- Note that $p_{M_0}$ is the same for all comparison and cancels out (????).\n",
    "\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbb{K} \\mathbb{L}\\left(p_{M_0} \\| q_{M_1}\\right)= & \\mathbb{E}\\left[\\log p_{M_0}\\right]-\\mathbb{E}\\left[\\log q_{M_1}\\right] \\\\\n",
    "\\mathbb{K} \\mathbb{L}\\left(p_{M_0} \\| q_{M_2}\\right)= & \\mathbb{E}\\left[\\log p_{M_0}\\right]-\\mathbb{E}\\left[\\log q_{M_2}\\right] \\\\\n",
    "& \\cdots \\\\\n",
    "\\mathbb{K} \\mathbb{L}\\left(p_{M_0} \\| q_{M_k}\\right)= & \\mathbb{E}\\left[\\log p_{M_0}\\right]-\\mathbb{E}\\left[\\log q_{M_k}\\right]\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54ed5319",
   "metadata": {},
   "source": [
    "$\\mathrm{ELPD}=\\sum_{i=1}^n \\int p_t\\left(\\tilde{y}_i\\right) \\log p\\left(\\tilde{y}_i \\mid y_i\\right) d \\tilde{y}_i$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89d774fb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac5f4487",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faa7aa7e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5442e51e",
   "metadata": {},
   "source": [
    "## Measure of uncertainty needs the following properties:\n",
    "\n",
    "1. continuous\n",
    "2. increase with the number of possible events\n",
    "3. additive"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "stat-rethink-pymc4",
   "language": "python",
   "name": "stat-rethink-pymc4"
  },
  "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.8"
  },
  "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
}