{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d65ab641",
   "metadata": {},
   "source": [
    "# 10 - Grande Entropia e Modelos Lineares Generalizados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2cea1405",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import networkx as nx\n",
    "# from causalgraphicalmodels import CausalGraphicalModel\n",
    "\n",
    "import arviz as az\n",
    "# ArviZ ships with style sheets!\n",
    "# https://python.arviz.org/en/stable/examples/styles.html#example-styles\n",
    "az.style.use(\"arviz-darkgrid\")\n",
    "\n",
    "import xarray as xr\n",
    "\n",
    "import stan\n",
    "import nest_asyncio\n",
    "\n",
    "plt.style.use('default')\n",
    "plt.rcParams['axes.facecolor'] = 'lightgray'\n",
    "\n",
    "# To DAG's\n",
    "import daft\n",
    "from causalgraphicalmodels import CausalGraphicalModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "449b2877",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add fonts to matplotlib to run xkcd\n",
    "\n",
    "from matplotlib import font_manager\n",
    "\n",
    "font_dirs = [\"fonts/\"]  # The path to the custom font file.\n",
    "font_files = font_manager.findSystemFonts(fontpaths=font_dirs)\n",
    "\n",
    "for font_file in font_files:\n",
    "    font_manager.fontManager.addfont(font_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2774f695",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To make plots like drawing \n",
    "# plt.xkcd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6f4405d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To running the stan in jupyter notebook\n",
    "nest_asyncio.apply()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c807993a",
   "metadata": {},
   "source": [
    "### R Code 10.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b1b1e4ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "      <th>D</th>\n",
       "      <th>E</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>8</td>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   A  B  C  D  E\n",
       "0  0  0  0  1  2\n",
       "1  0  1  2  2  2\n",
       "2  1  8  6  4  2\n",
       "3  0  1  2  2  2\n",
       "4  0  0  0  1  2"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "buckets = {\n",
    "    \"A\": (0, 0, 1, 0, 0),\n",
    "    \"B\": (0, 1, 8, 1, 0),\n",
    "    \"C\": (0, 2, 6, 2, 0),\n",
    "    \"D\": (1, 2, 4, 2, 1),\n",
    "    \"E\": (2, 2, 2, 2, 2),\n",
    "}\n",
    "\n",
    "df = pd.DataFrame.from_dict(buckets)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8e7607d",
   "metadata": {},
   "source": [
    "### R Code 10.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "19b6ae3d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "      <th>D</th>\n",
       "      <th>E</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.6</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     A    B    C    D    E\n",
       "0  0.0  0.0  0.0  0.1  0.2\n",
       "1  0.0  0.1  0.2  0.2  0.2\n",
       "2  1.0  0.8  0.6  0.4  0.2\n",
       "3  0.0  0.1  0.2  0.2  0.2\n",
       "4  0.0  0.0  0.0  0.1  0.2"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Normalize\n",
    "\n",
    "p_norm = df / df.sum(axis=0)\n",
    "p_norm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4f9da3c",
   "metadata": {},
   "source": [
    "### R Code 10.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "74fdb4ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy(bucket):\n",
    "    uncertainty = []\n",
    "    \n",
    "    for q in bucket:\n",
    "        if q == 0:\n",
    "            uncertainty.append(q)\n",
    "        else:\n",
    "            uncertainty.append(q * np.log(q))\n",
    "    \n",
    "    return (-1) * np.sum(uncertainty)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5ca8a893",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "      <th>D</th>\n",
       "      <th>E</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.639032</td>\n",
       "      <td>0.950271</td>\n",
       "      <td>1.470808</td>\n",
       "      <td>1.609438</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     A         B         C         D         E\n",
       "0 -0.0  0.639032  0.950271  1.470808  1.609438"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = [entropy(p_norm[key]) for key in p_norm.keys()]\n",
    "\n",
    "df_H = pd.DataFrame(H).T\n",
    "df_H.columns = p_norm.keys()\n",
    "df_H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d185249",
   "metadata": {},
   "source": [
    "### R Code 10.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3cfacebb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 0.44998097, 0.7138867 , 1.05400644, 1.16386767])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ways = (1, 90, 1260, 37800, 113400)\n",
    "\n",
    "logwayspp = np.log(ways) / 10\n",
    "\n",
    "logwayspp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2088a2f6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAGDCAYAAAASzPzoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAABMZUlEQVR4nO3deXhV5b328e8vAzKjMkhAFPTgFBFRcJZJUasiUAW1aq3jUavWUlv11FpPtcdaezqrFa1VX60KqASn4pAozkWcWrFVBBQUZxQMAULyvH9kkxORISg7O8P3c11c7r3Ws/e+d7oK3KxnPStSSkiSJEmS1Njl5TqAJEmSJEn1YYGVJEmSJDUJFlhJkiRJUpNggZUkSZIkNQkWWEmSJElSk2CBlSRJkiQ1CRZYSZIaiYg4LiIeynWOuiKid0SkiCjIdRZJkiywkqQmLyLmRURFRHxe59cf6/naxyLi1GxnrI+U0m0ppYO+ymsj4tKIqKzz/V+LiCM3dsYNyGPxlSRtdBZYSVJzMTKl1L7Or7M3xps2sQJ256rvD5wH3BoRW+Q4kyRJG40FVpLUrEXEdyLiyYj4VUQsioi5EfGNzL6fA/sDf6x71jZz5vC7EfEG8EZm22kRMTsiPomIqRHRo85npIg4NyLmRMRHEXFVRORFRKvM+H51xnaLiKUR0XVtWVd73zMi4o2I+DQiro6IqM/3TilNA5YA267pveu8/39kHreJiP+NiLci4rPMz6zNGjIemTnjvXPmO14YEW9GxMcRMTEiNs8MnZ7576eZn+3eEfEfEfF45v0/iog76/NdJElaxQIrSWoJ9gT+DXQBfgn8OSIipfRj4Ang7DWctR2ded1OETEcuAIYBxQBbwF3rPYZY4CBwG7AKODklNKKzLjj64w7Fng0pfRhPbMfDgwCdsl8/sHre0HUOAxoBcyq5+f8Ctgd2AfYHPgRUL3a+54EXAkcmFL6J3AONT+nIUAPYBFwdWb44Mx/N838bJ8BLgMeAjYDtgT+UM9skiQBFlhJUvMxJXOWctWv0+rseyuldH1KqQq4mZoSur6ptVeklD5JKVUAxwE3ppReSCktBy4C9o6I3nXGX5kZ/zbwW2qKKpnPO7bOmdMTgP+3Ad/rFymlTzPvWwbsuo6x4yLiU+BzYCrwPymlT9f3ARGRB5wMfC+l9E5KqSql9HTmu65yHvBDYGhKaXZm2xnAj1NKCzJjLwWOWse060pga6BHSmlZSunJtYyTJGmNLLCSpOZidEpp0zq/rq+z771VD1JKSzMP26/n/ebXedyDmrOuq97jc+BjoOdaxr+VeQ0ppeeApcDQiNgB+A9qymV9vVfn8dL15J6Y+e7tqJk6/O2I+M96fEYXoDXw5jrG/BC4OqW0oM62rYF7Vv2jAfAaUMXa/3HgR0AAf4+IVyPi5HpkkySplgVWktTSpXpsf5easgZARLQDOgPv1BnTq87jrTKvWeVmaqYRnwBMTikt+zqB6yOlNA94EBiZ2VQOtF21PyK61xn+EbCMzPWya3EQcPFqKxvPB76x2j8ctE4pvcMafq4ppfdSSqellHoA/wlcs+oaXEmS6sMCK0lq6d4HtlnPmNuBkyJi14jYBPgf4LlMSVzlhxGxWUT0Ar4H1F2g6FZqrpE9HrhloyVfh4jYEjgEeDWz6WWgOPMdWlMz3ReAlFI1cCPw64joERH5mUWXNqnzlq9m3u/qiDgis+1PwM8jYuvMZ3aNiFGZfR9Scw1t7c82IsZmckHN9bKJ1a6zlSRpXSywkqTm4t744n1g76nn635HzXWbiyLi92sakFJ6BPgJcBewkJozlcesNqwEmAm8BNwP/LnO6+cDL1BT2J6o/1faYEev+v7ADOAp4L8zGV4HfgY8Qs3Kyqtff3o+8I/M6z6hZrGmL/w9IaX0MjWLSl2fWcn5d9RMh34oIpYAz1Kz8NWqqdo/B57KTDHei5rFqJ7L5JtKzTW3czbuj0CS1JxFSmubOSVJkuojIhLQt87iRmsacyPwbkrp4oZLJklS89KUbs4uSVKTlFmt+JvAgBxHkSSpSXMKsSRJWRQRlwH/BK5KKc3NdR5JkpoypxBLkiRJkpoEz8BKkiRJkpoEC6wkSZIkqUlocos4bbbZZqlHjx65jiFJkiRJyoJZs2Z9lFLquqZ9Ta7A9ujRgzvvvHP9AyVJkiRJTU6/fv3eWts+pxBLkiRJkpoEC6wkSZIkqUmwwEqSJEmSmgQLrCRJkiSpSbDASpIkSZKaBAusJEmSJKlJsMBKkiRJkpqEJncfWEmSJElS/S1cUsk9sxZTNq+cZZWJ1oXBsN7tGLNTR4o6FOY63gaxwEqSJElSM/X8OxVcMf1DVlYnqlLNtorKxLTZn/PonHIuGtyVgT3b5DbkBnAKsSRJkiQ1QwuXVHLF9A9ZXvV/5XWVqgTLqxJXTP+QhUsqcxPwK7DASpIkSVIzdM+sxaysTuscs7I6MeW1xQ2U6OvLWoGNiBsj4oOI+Oc6xgyNiJci4tWIeDxbWSRJkiSppSmbV/6lM6+rq0pQNre8YQJtBNk8A3sTcMjadkbEpsA1wBEppWJgbBazSJIkSVKLsqxyPe01o6Ke4xqDrBXYlNJ04JN1DPkWcHdK6e3M+A+ylUWSJEmSWpI33ngDVi6v19g2hZHlNBtPLq+B3Q7YLCIei4iZEfHttQ2MiNMj4vmIeH7RokUNGFGSJEmSGr/PP/+c++67jxdffBGAwsJCKt94ClLVOl+XHzCsT7uGiLhR5PI2OgXA7sABQBvgmYh4NqX0+uoDU0oTgAkAxcXFTef8tiRJkiRlyeeff85jjz3GtGnTePrpp1mxYgVHHnkkAwYMoHfv3tx4wXGcc/97LF/HhbAFecHoHTs2YOqvJ5cFdgHwcUqpHCiPiOlAf+BLBVaSJEmSBJWVlRQWFgJwwgknMHv2bLp168a4ceM46KCD6N+/f+3YHh1bcdHgrl+6DyzUnHktyAsuGtyVog6FDf01vrJcFtgS4I8RUQC0AvYEfpPDPJIkSZLU6NQ90/rqq68ybdo0CgsLOe+88+jYsSP9+/cnL2/NV4cO7NmGPx5exJTXFlM2t5yKykSbwmBYn3aM3rFjkyqvkMUCGxG3A0OBLhGxAPgpUAiQUvpTSum1iPgb8ApQDdyQUlrrLXckSZIkqSX55z//yXXXXVc7Pbhbt24cfPDBVFRUUFhYyJAhQ+r1PkUdCjlzj86cuUfnLCfOvqwV2JTSsfUYcxVwVbYySJIkSVJTsepM6/bbb0/fvn1Zvnw5s2bNYty4cRx88MHssssuaz3T2lLkcgqxJEmSJLVoq0rrQw89xFNPPcWKFSs49dRT+d73vseAAQN4+OGHW3xprcsCK0mSJEkNqKqqivz8fKqrqzniiCP48MMPaxdiWnWmFbC4roEFVpIkSZKyrO6Z1vnz53P33XeTl5fH+eefT48ePZweXE8WWEmSJEnKkpkzZ3LzzTfXTg/u1q0bBx10ECtWrGCTTTbh0EMPzXXEJsUCK0mSJEkbyaozrbvvvjtFRUV89NFHvPrqq1+4T6tnWr86C6wkSZIkfQ1rWojp/PPP58QTT+SAAw5gxIgRltaNxAIrSZIkSRsopUREUFFRwQEHHMDSpUtrF2JadaYVoKDAyrUx+dOUJEmSpHooLy/nscceY9q0aSxbtowJEybQpk0bfvCDH9C3b1+nBzcAC6wkSZIkrcMzzzzD7bff/oWFmA4++GCqq6vJy8tj3LhxuY7YYlhgJUmSJKmOVde07r///nTq1Ik5c+a4EFMjYYGVJEmS1OKtaSGm//mf/2HkyJGMHTuWY4891tLaCFhgJUmSJLVIqxZi+uijjzj44INrpwevvhBTq1atcpxUq1hgJUmSJLUYdc+0tmvXjiuuuIIuXbpw1llnsdtuuzk9uJGzwEqSJElqUhYuqeSeWYspm1fOsspE68JgWO92jNmpI0UdCtf4munTpzN58uQvLMR0xBFH1O4/5ZRTGiq+vgYLrCRJkqQm4/l3Krhi+oesrE5UpZptFZWJabM/59E55Vw0uCsDe7ahvLycxx9/nBEjRlBYWMhLL73kQkzNgAVWkiRJUpOwcEklV0z/kOWrmmsdVQmqqhKXl71H0cs38vfSB1ixYgWdOnVi33335bTTTuPss8+2tDZxFlhJkiRJTcI9sxazsvrL5bWuyqoq5rbq86WFmNq0adMQEZVlFlhJkiRJTULZvHLWcPL1CyK/kI67HMgFR2/VMKHUoDx/LkmSJKlJqKhcT3vdwHFqeiywkiRJkhqtlP6vjOZVrajXa9oURrbiKMcssJIkSZIanUWLFnHzzTczcuRIXn31VQD227IVeaz77Gp+wLA+7RoionLAa2AlSZIkNQopJV544QUmTZrEQw89RGVlJQMGDKCyshKAE/bowXP3LVzjKsSrFOQFo3fs2FCR1cAssJIkSZJyqrq6mry8PCoqKjj77LOJCI466ijGjh1L3759a8cVdSjkosFdv3QfWKg581qQF1w0uCtFHQpz8C3UECywkiRJkhpcSol//OMfTJo0iddff5077riDtm3bcu2117LddtvRtm3bNb5uYM82/PHwIqa8tpiyueVUVCbaFAbD+rRj9I4dLa/NnAVWkiRJUoP5/PPPuf/++5k0aRL//ve/adOmDYcddhgVFRW0bduWXXfddb3vUdShkDP36MyZe3TOfmA1KhZYSZIkSVlXVVVFfn4+06dP5/LLL2eHHXbgJz/5CYceeijt27fPdTw1ERZYSZIkSVmxdOlSpk2bxqRJkxg+fDinnnoqBx54ILfddhv9+vUjwtvdaMNYYCVJkiRtVG+88QaTJk3i3nvv5fPPP2fbbbdliy22AKBVq1bssssuOU6opiprBTYibgQOBz5IKe28jnGDgGeAY1JKk7OVR5IkSVL2rFy5koKCmnrxm9/8hmeffZaDDjqIcePGMWDAAM+2aqPI5hnYm4A/AresbUBE5ANXAg9lMYckSZKkLJk7dy6TJk3i/vvv5/bbb6dHjx5ceOGFdOjQgc022yzX8dTMZK3AppSmR0Tv9Qw7B7gLGJStHJIkSZI2rsrKSh599FEmTpzIjBkzKCgoYPjw4VRWVgKw1VZb5TihmqucXQMbET2BMcAw1lNgI+J04HSAoqKi7IeTJEmS9CWVlZUUFhby2WefcdFFF7HFFlvwve99j9GjR9OlS5dcx1MLkMtFnH4LXJBSql7ffPiU0gRgAkBxcXHKfjRJkiRJUFNaH3/8cSZOnEhVVRV//vOf6dKlC7fffjvbbbcdeXl5uY6oFiSXBXYgcEemvHYBDo2IlSmlKTnMJEmSJAlYuHAhkydP5p577uHDDz9kiy224KijjqK6upq8vDx22GGHXEdUC5SzAptS6rPqcUTcBNxneZUkSZJyp6qqiurqagoLCyktLeX6669n//3355JLLmG//farXWVYypVs3kbndmAo0CUiFgA/BQoBUkp/ytbnSpIkSdowH3zwAXfffTd33XUXZ511FmPGjGHUqFEMGzaMHj165DqeVCubqxAfuwFjv5OtHJIkSZK+LKXE008/zaRJk3jssceoqqpi7733pmfPngC0b9+e9u3b5zil9EXOAZAkSZJakIqKCtq0aUNE8Ic//IGFCxfy7W9/m6OOOsrb36jRs8BKkiRJzVxKiRkzZjBx4kSefvppHnzwQTp16sSvfvUrunXrRqtWrXIdUaoXC6wkSZLUTC1evJh77rmHyZMnM2/ePDp27MioUaNYuXIlAFtuuWWOE0obxgIrSZIkNSMpJZYuXUq7du346KOP+NWvfsWuu+7Kz3/+cw466CBat26d64jSV2aBlSRJkpqBxYsXc9999zFp0iS22WYb/vd//5dtttmG+++/32tb1WxYYCVJkqQmbNasWdxxxx08+OCDLFu2jOLiYgYPHly73/Kq5sQCK0mSJDUx5eXltG7dmvz8fB555BH+9re/cdhhhzF27FiKi4tzHU/KGgusJEmS1ES89tprTJo0ifvvv58rr7ySoUOHcuKJJ3LyySd7z1a1CBZYSZIkqRGrrKysvbb1H//4B5tssgmHHHJI7QrCnTp1ynFCqeFYYCVJkqRG6NNPP2XTTTclLy+P6667jtatW3PhhRdy+OGHW1rVYllgJUmSpEZi+fLlPPTQQ0yaNIm33nqLRx55hMLCQm655Ra6du1KROQ6opRTFlhJkiQpx959911uu+02SkpK+Oyzz9h666056aSTWLlyJYWFhXTr1i3XEaVGwQIrSZIk5UBlZSVLly6lU6dOvPvuu/z1r39l+PDhjBs3jkGDBpGXl5friFKjY4GVJEmSGtD8+fOZPHkyU6ZM4eCDD+a//uu/2H333XnkkUfo3LlzruNJjZoFVpIkSWoA06dP57bbbuPpp58mLy+PIUOGMHz4cAAiwvIq1YMFVpIkScqSDz74oHbxpdLSUt58803OOussxowZQ/fu3XMdT2pyIqWU6wwbpLi4ON155525jiFJkiStUVVVFU899RQTJ07kiSee4NZbb6Vfv34sXryYtm3bUlDgOSRpXfr16zczpTRwTfv8f48kSZK0EZSXl/P//t//4+6772bhwoV06dKFU045hS222AKAjh075jih1PRZYCVJkqSvqLq6mg8++IDu3buTn5/Prbfeyo477sgPf/hDhg4dSmFhYa4jSs2KBVaSJEnaQB9//DFTpkxh8uTJ5Ofnc++999K6dWseeOABz7RKWWSBlSRJkupp1qxZ/OUvf+GRRx5h5cqVDBw4kLFjx1JdXU1+fr7lVcoyC6wkSZK0Dp9++ikRQadOnZg/fz5PP/00xxxzDGPHjmWbbbbJdTypRbHASpIkSatJKfHSSy8xceJEHnroIU477TTOOOMMhg8fzpAhQ2jdunWuI0otkgVWkiRJqmPixIncfvvtzJ49m3bt2jFmzBgOPPBAAAoLC12YScohC6wkSZJatJQSc+fOrZ0O/MQTT7DJJptw6aWX8o1vfIO2bdvmOKGkVSywkiRJapHKy8t54IEHmDRpEq+99hoPPPAAvXr14sorr7S0So2UBVaSJEktygcffMB1113Hfffdx9KlS9luu+24+OKL2XzzzQEsr1IjlrUCGxE3AocDH6SUdl7D/uOAC4AAlgBnppRezlYeSZIktVwVFRV89NFH9OrVi4KCAh588EEOPPBAxo0bxy677EJE5DqipHrI5hnYm4A/AresZf9cYEhKaVFEfAOYAOyZxTySJElqYWbPns2kSZO499576du3LzfffDObb745paWlriQsNUFZK7AppekR0Xsd+5+u8/RZYMtsZZEkSVLL8uSTT3L99dfzwgsvUFhYyIgRIxg7dmztfsur1DQ1lmtgTwEezHUISZIkNV3z5s2jW7dutG3blvnz5/PRRx8xfvx4Ro0aVXt9q6SmLecFNiKGUVNg91vHmNOB0wGKiooaKJkkSZIau8rKSh599FEmT57Mc889xyWXXMLYsWM56qijOProo8nLy8t1REkbUU4LbETsAtwAfCOl9PHaxqWUJlBzjSzFxcWpgeJJkiSpkVq5ciVXX301d999N5988gk9evTg3HPPZdiwYQAUFhbmOKGkbMhZgY2IrYC7gRNSSq/nKockSZKahpUrV/Lvf/+b4uJiCgoKmDlzJrvssgvjxo1jn332IT8/P9cRJWVZNm+jczswFOgSEQuAnwKFACmlPwGXAJ2BazLLlq9MKQ3MVh5JkiQ1TQsXLuSuu+7innvu4dNPP+XRRx9l00035c9//rNnWqUWJpurEB+7nv2nAqdm6/MlSZLUtM2ePZvf/va3PPHEE6SU2HfffRk7dizt27cHnCYstUQ5X8RJkiRJWuXDDz+koqKCrbbaioKCAmbNmsUpp5zCkUceSc+ePXMdT1KOWWAlSZKUU9XV1Tz77LNMnjyZsrIyhgwZwm9/+1t69+7Nww8/7LWtkmpZYCVJkpQzd911F3/+85+ZP38+m266KccffzxHHXVU7X7Lq6S6LLCSJElqMCklZs6cya677kpBQQHvv/8+3bp147vf/S4jRoygVatWuY4oqRGzwEqSJCnrPvvsM0pKSpg8eTJz587lt7/9LQcccABnnHEGZ511Vq7jSWoiLLCSJEnKmiVLlnDFFVcwbdo0VqxYwS677MJll13GPvvsA0BeXl6OE0pqSiywkiRJ2qiWLFnC7NmzGTBgAO3atWP27NmMGTOGsWPHsv322+c6nqQmzAIrSZKkry2lxKuvvsqkSZN48MEHadWqFaWlpbRq1Yo777yTiMh1REnNgAVWkiRJX8uMGTO46qqreO2112jTpg2HHnooY8eOrV2QyfIqaWOxwEqSJGmD/etf/6J9+/ZsueWWFBYWsnLlSn784x9z2GGH0aFDh1zHk9RMWWAlSZJULxUVFUybNo1JkybxyiuvcPTRR3PxxRfTv39/7rrrLs+0Sso6C6wkSZLW6/e//z133HEHS5YsoU+fPlxwwQWMHDkScIqwpIZjgZUkSdKXrFixgieeeILhw4cTEVRUVLDffvsxbtw4dt99d0urpJywwEqSJKnWW2+9xeTJk5kyZQqffvopt956K/379+dHP/qRpVVSzllgJUmSxMKFC/nJT37Cc889R0FBAcOGDWPs2LH069cPcJqwpMbBAitJktRCLViwgPfee4+BAwey+eabs3jxYs455xzGjBlD165dcx1Pkr7EAitJktSCrFy5kunTpzNx4kSefvpptt56a6ZOncomm2zCxIkTcx1PktbJAitJktRC3H///fz617/mgw8+oFu3bpxxxhl885vfdHqwpCbDAitJktRMVVVV8dRTT7HDDjvQrVs32rRpQ9++ffnxj3/M4MGDKSjwr4KSmhZ/15IkSWpmPvroI+6++27uuusu3n33Xc455xxOP/10hg8fzvDhw3MdT5K+MgusJElSM1FdXc2FF17Iww8/zMqVK9lzzz0ZP368pVVSs2GBlSRJasI++eQTnn32WQ499FDy8vJo27Ytxx13HEcddRS9e/fOdTxJ2qgssJIkSU1MSonnn3+eSZMm8cgjj1BZWcluu+1G9+7dufTSS3MdT5KyxgIrSZLUhMyaNYsLL7yQuXPn0qFDB8aNG8fYsWPp3r17rqNJUtZZYCVJkhqxlBIvv/wyVVVV7L777vTs2ZPNN9+ck08+mYMPPpg2bdrkOqIkNRgLrCRJUiO0ZMkS7rvvPiZNmsQbb7zBnnvuyQ033ECnTp246aabch1PknLCAitJktTI3HDDDUyYMIGKigp23HFHLrnkEg499NBcx5KknLPASpIk5djSpUt58MEHOeigg+jQoQNdu3blG9/4BuPGjaO4uDjX8SSp0chagY2IG4HDgQ9SSjuvYX8AvwMOBZYC30kpvZCtPJIkSY3Nv//9byZNmsR9991HeXk5rVq1YuTIkYwaNYpRo0blOp4kNTrZPAN7E/BH4Ja17P8G0Dfza0/g2sx/JUmSmrWlS5dy2mmn8corr9CqVSsOOeQQxo4dS//+/XMdTZIatawV2JTS9IjovY4ho4BbUkoJeDYiNo2IopTSwmxlkiRJypU5c+bwj3/8g1GjRtG2bVt69+7NIYccwhFHHEGnTp1yHU+SmoRcXgPbE5hf5/mCzDYLrCRJahZWrFjBww8/zKRJk5g5cyZt27ZlxIgRtG3blp///Oe5jidJTU6TWMQpIk4HTgcoKirKcRpJkqT1e+KJJ/jxj3/MokWL6NWrF+PHj689+ypJ+mpyWWDfAXrVeb5lZtuXpJQmABMAiouLU/ajSZIkbZjKykrKysrYYost6N+/P71792b33Xdn7Nix7LXXXuTl5eU6oiQ1ebkssFOBsyPiDmoWb/rM618lSVJT884773DXXXdx99138/HHHzNmzBj69+9Pr169+M1vfpPreJLUrGTzNjq3A0OBLhGxAPgpUAiQUvoT8AA1t9CZTc1tdE7KVhZJkqRsuPTSS7n77ruJCAYPHszYsWPZd999cx1LkpqtbK5CfOx69ifgu9n6fEmSpPpauKSSe2YtpmxeOcsqE60Lg2G92zFmp44UdSisHffee+9x3333ceKJJ1JYWMgOO+zAf/7nf3LkkUfSvXv3HH4DSWoZmsQiTpIkSdny/DsVXDH9Q1ZWJ6oyK21UVCamzf6cR+eUc8F+m7N83otMnDiR6dOnk1Kif//+DBo0iGOOOSa34SWphbHASpKkFmvhkkqumP4hy6u+vEZkVYKqqsR/P/ou7/75Z3TMW8HJJ5/MkUceyZZbbpmDtJIkC6wkSWqx7pm1mJXV677BQeQXcMi5v+Rno3emsLBwnWMlSdnleu6SJKnFKptXzhpOvn5R5PNmdWfLqyQ1AhZYSZLU4ixfvpzq6mqWVdbv9vIV9RwnScouC6wkSWoRPvvsM+677z7Gjx/P4MGDefnll2ldGPV6bZt6jpMkZZfXwEqSpGbtvffe4yc/+QnPP/88K1eupEuXLhx66KF07NiRYb3bMW325+ucRpwfMKxPu4YLLElaKwusJElqNlJKvPHGG5SVlbHZZpsxbtw4OnfuzJIlS/j2t7/N8OHD6devH3l5NZPQxnSr5NE55VSto8EW5AWjd+zYUF9BkrQOFlhJktTkvfTSSzz88MOUlpayYMECAA4//HDGjRtHYWEhd9xxxxpfV9ShkIsGd/3SfWCh5sxrQV5w0eCuFHVwASdJagwssJIkqclZtmwZL774InvvvTcAt9xyC4899hh77rknJ598MkOHDqVr1671eq+BPdvwx8OLmPLaYsrmllNRmWhTGAzr047RO3a0vEpSIxIpNa1V9YqLi9Odd96Z6xiSJKmBffrpp0yfPp3S0lKefvppKioqeOCBB+jVqxcLFy6kY8eOtGvntaqS1NT169dvZkpp4Jr2eQZWkiQ1WiklIoInnniCc845h6qqKrp168YRRxzB8OHD6d69OwBFRUU5TipJaggWWEmS1GiklPj3v/9NaWkpZWVljB49muOOO45+/fpx8sknc8ABB7DTTjsR4W1tJKklqleBjYj8lFJVtsNIkqSWKaXEVVddxaOPPsq7775LRDBgwIDa61g33XRTzj333BynlCTlWn3PwL4REXcBf0kpzcpmIEmS1PwtXbqUZ555hrfffpuTTjqJiOCNN95gu+2244wzzmDw4MF07tw51zElSY1MfQtsf+AY4IaIyANuBO5IKS3OWjJJktSsfPLJJzz++OOUlpbyzDPPsHz5cjp37sxxxx1Hq1atmDBhglODJUnrVK8Cm1JaAlwPXB8RQ4C/Ar+JiMnAZSml2VnMKEmSmqi3336bLl260LZtW+6++25+97vfUVRUxFFHHcXw4cMZMGAAhYU1t6mxvEqS1qfe18AChwEnAb2B/wVuA/YHHgC2y1I+SZLUhKSUmDVrFo8++ihlZWXMnj2bX/3qVxx88MGMGjWK/fbbj+23396yKkn6Sup9DSxQBlyVUnq6zvbJETF448eSJElNzccff8zRRx/N+++/T15eHrvvvjsXXHABAwYMAKBr1661izJJkvRV1LfA7pJS+nxNO1JKLgkoSVILU15ezpNPPklpaSmbbbYZF154IZtvvjn7778/u+66K4MHD2azzTbLdUxJUjNT3wLbLSJuB/YGqoFngO+nlOZkLZkkSWp0/va3v1FSUsJzzz1HZWUlm222GSNHjgRqrmH96U9/muOEkqTmrL4F9q/A1cCYzPNjgNuBPbMRSpIkNQ5z585l+vTpHH/88eTn5/Pyyy8zb948jj32WIYPH86uu+5Kfn5+rmNKklqISCmtf1DEKymlXVbb9nJKqX/Wkq1FcXFxuvPOOxv6YyVJahGqq6v55z//SWlpKaWlpcydOxeAO+64g+LiYpYvX06rVq1chEmSlDX9+vWbmVIauKZ99T0D+2BEXAjcASTgaOCBiNgcIKX0yUZJKkmSGtyKFStYtmwZHTt2ZMaMGZx66qkUFBQwcOBAjj32WIYNG0b37t0B2GSTTXKcVpLUktW3wI7L/Pc/V9t+DDWFdpuNlkiSJGXd4sWLaxdhevLJJzn66KP5/ve/z2677cYvfvEL9ttvPzp16pTrmJIkfUG9CmxKqU+2g0iSpIYxfvx4ysrKWLlyJZ07d+aQQw5h//33B6CwsJDDDjssxwklSVqzehXYiCgEzgRW3fP1MeC6lFJllnJJkqSvKaXEm2++SVlZGW+88Qa//OUvAejZsycnnHACw4cPZ5dddiEvLy/HSSVJqp/6TiG+FigErsk8PyGz7dRshJIkSV/d7NmzKSkpobS0lLfffhuAfv36UV5eTrt27fjBD36Q44SSJH019S2wg1Zbcbg0Il5e34si4hDgd0A+cENK6Rer7d8KuBnYNDPmwpTSA/XMJEmSgGXLlvHss8+y00470a1bN1577TVuvfVW9txzT0488USGDh1Kt27dch1TkqSvrb4Ftioitk0pvQkQEdsAVet6QUTkU3Pv2BHAAmBGRExNKc2qM+xiYGJK6dqI2Al4AOi9gd9BkqQW57PPPmP69OmUlpby1FNPUVFRwQUXXMDxxx/PgQceyLBhw2jfvn2uY0qStFHVt8CeD5RFxBwggK2Bk9bzmj2A2SmlOQARcQcwCqhbYBPQMfO4E/BuPfNIktTirFixglatWvH5558zbNgwKisr6datG0cccQTDhw9n0KBBALRp0ybHSSVJyo71FtjMmdT+QF9g+8zmf6eUlq/npT2B+XWeLwD2XG3MpcBDEXEO0A44sB6ZJUlqEVJKvP7665SWllJWVkanTp24/vrrad++PRdeeCE77rgjxcXFLsIkSWox1ltgU0pVEXFsSuk3wCsb+fOPBW5KKf1vROwN/L+I2DmlVF13UEScDpwOUFRUtJEjSJLU+Pz1r3/llltu4Z133iEi2HXXXRk8eHDt/nHjxq3j1ZIkNU/1nUL8VET8EbgTKF+1MaX0wjpe8w7Qq87zLTPb6joFOCTzXs9ERGugC/BB3UEppQnABIDi4uJUz8ySJDUJFRUVPP3005SVlXHBBRfQoUMHUkpsu+22nHbaaQwZMoQuXbrkOqYkSTlX3wK7a+a/P6uzLQHD1/GaGUDfiOhDTXE9BvjWamPeBg4AboqIHYHWwIf1zCRJUpO1ZMkSHnnkEcrKynjmmWdYtmwZHTp0YOzYsfTv35/jjjuO4447LtcxJUlqVOpbYE9ZtRjTKpmViNcqpbQyIs4GplFzi5wbU0qvRsTPgOdTSlOBHwDXR8T3qSnE30kpeYZVktQszZ8/n6qqKnr37s0HH3zAJZdcQvfu3fnmN7/J8OHD2W233SgsLMx1TEmSGq2oT1+MiBdSSruttm1mSmn3rCVbi+Li4nTnnXc29MdKkrTBUkrMmjWrdhGmN954g8MPP5wrrriidoGm7bbbjojIdVRJkhqNfv36zUwpDVzTvnWegY2IHYBioFNEfLPOro7UTPeVJEl1pJRqC+lJJ53EzJkzycvLY7fdduOHP/whw4fXXH0TEWy//fbreitJkrSa9U0h3h44HNgUGFln+xLgtCxlkiSpSVm6dClPPvkkZWVlvPLKK5SUlFBQUMCoUaMYPXo0Q4YMYbPNNst1TEmSmrx1FtiUUglQEhF7p5SeaaBMkiQ1CS+99BLXX389zz77LCtWrGDTTTdl8ODBlJeX06lTJ8aMGZPriJIkNSv1XcRpdkT8F9C77mtSSidnI5QkSY3RvHnzKC0tZd9992X77bdn2bJlvPnmmxx99NEMGzaMAQMGUFBQ3z9aJUnShqrvn7IlwBPAI0BV9uJIktR4VFdX8+qrr1JaWkppaSlz5tQsyJ+fn8/222/PnnvuyYMPPugiTJIkNZD6Fti2KaULsppEkqRGoLKykoULF7LVVltRVVXF6aefTkVFBQMHDmTcuHEMGzaMHj16AFhcJUlqYPUtsPdFxKEppQeymkaSpBxYsmQJTz75JKWlpTz55JN069aNkpISCgsL+cMf/kDfvn3p1KlTrmNKktTi1bfAfg+4KCJWAJVAACml1DFrySRJagATJkzg2muvZeXKlWy++eYcfPDBDB8+vPZ2OAMHrvE2dJIkKQfqW2A7AccBfVJKP4uIrYCi7MWSJGnjSikxd+7c2utZr7zySnr16sUOO+zACSecwLBhw9hll13Iz8/PdVRJkrQW9S2wVwPVwHDgZ9TcB/YuYFCWckmStFF89NFH3HLLLZSVlTFv3jwAdt55ZxYtWkSvXr0YPHgwgwcPzm1ISZJUL/UtsHumlHaLiBcBUkqLIqJVFnNJkvSVLF++nOeee45WrVqx1157UVBQwO23385uu+3G8ccfz9ChQ9liiy1yHVOSJH0F9S2wlRGRDySAiOhKzRlZSZJy7rPPPmP69OmUlZXx5JNPUlFRwf77789ee+3FpptuyvTp02nTpk2uY0qSpK+pvgX298A9QLeI+DlwFHBx1lJJkrQeixYtYrPNNgPgvPPO4/nnn6dr166MHDmS4cOHM2jQ/13lYnmVJKl5qFeBTSndFhEzgQOoWYF4dErptawmkySpjpQSb7zxRu0iTG+88QaPP/44HTt25Oyzz6awsJCdd96ZvLy8XEeVJElZUt8zsKSU/gX8K4tZJElao+eee45LL72UBQsWEBH079+fc889t3b/7rvvnsN0kiSpodS7wEqS1BAqKip45plnKC0tZcSIEQwZMoRu3brRp08fTjnlFIYOHUqXLl1yHVOSJOWABVaSlHNVVVXcd999lJaW8vTTT7Ns2TI6dOjAzjvvDECfPn245pprcpxSkiTlmgVWkpQTCxYs4K233mLfffclLy+P6667juXLlzN69GiGDx/OwIEDKSwszHVMSZLUiFhgJUkNIqXEv/71r9pFmF5//XU6derEY489RkFBATfddBNdu3YlInIdVZIkNVIWWElSvSxcUsk9sxZTNq+cZZWJ1oXBsN7tGLNTR4o6rPlMaWVlJfn5+eTl5XH11Vdz3XXXkZeXx6677sr555/P8OHDKSio+aOoW7duDfl1JElSE2SBlSSt1/PvVHDF9A9ZWZ2oSjXbKioT02Z/zqNzyrlocFcG9qy51+rSpUt56qmnKCsr4/HHH+f3v/89u+++OyNGjKBnz54MGTKEzTffPIffRpIkNVUWWEnSOi1cUskV0z9k+armWkdVgqqqxBXTP+Rn+2zCtVddxjPPPMOKFSvo1KkTQ4cOpX379gBsv/32bL/99g0dX5IkNSMWWEnSOt0zazErq79cXutaWZ14bGEe77//PuPGjWP48OEMGDCgdnqwJEnSxuDfLCRJ61Q2r5w1nHz9gqoEj79VwaRJkxomlCRJapHych1AktS4LatcT3vNqKjnOEmSpK/KAitJ+pKqqiqeeuopfvSjH1G1fGm9XtOm0NvfSJKk7HIKsSSp1qJFi7j55pu59957+eCDD+jUqRPb9xnBR623W+c04vyAYX3aNVxQSZLUInkGVpJauMWLFzNnzhwA8vLy+Otf/8qOO+7Ir3/9a0pLS7n8+KEU5K377GpBXjB6x44NEVeSJLVgWT0DGxGHAL8D8oEbUkq/WMOYccClQAJeTil9K5uZJEk1U4Sfe+45pkyZQmlpKdtvvz233XYbnTp1oqysjHbt/u9salEruGhw1y/dBxZqzrwW5AUXDe5KUYfCHHwTSZLUkmStwEZEPnA1MAJYAMyIiKkppVl1xvQFLgL2TSktiohu2cojSaoxefJkrr322topwt/85jcZNWpU7f665XWVgT3b8MfDi5jy2mLK5pZTUZloUxgM69OO0Tt2tLxKkqQGkc0zsHsAs1NKcwAi4g5gFDCrzpjTgKtTSosAUkofZDGPJLVIixcvZtq0aRx00EF06tSJ/Px8dtxxRy688EKGDBlCq1at6vU+RR0KOXOPzpy5R+csJ5YkSVqzbBbYnsD8Os8XAHuuNmY7gIh4ipppxpemlP62+htFxOnA6QBFRUVZCStJzcnqU4SXL1/OJptswhFHHMGYMWMYM2ZMriNKkiRtsFyvQlwA9AWGAlsC0yOiX0rp07qDUkoTgAkAxcXF3mhQktZhyZIljBkzhvfff5+OHTsyZswYRo8ezU477ZTraJIkSV9LNgvsO0CvOs+3zGyrawHwXEqpEpgbEa9TU2hnZDGXJDUrq6YIv//++5x99tl06NCBQw89lJ133pmhQ4fWe4qwJElSY5fNAjsD6BsRfagprscAq68wPAU4FvhLRHShZkrxnCxmkqRmYU1ThHfccUfOOOMMCgoKGD9+fK4jSpIkbXRZK7AppZURcTYwjZrrW29MKb0aET8Dnk8pTc3sOygiZgFVwA9TSh9nK5MkNXUpJSKCv/zlL/zud7/70hThiHXfr1WSJKkpy+o1sCmlB4AHVtt2SZ3HCRif+SVJWoNVU4RLSko4+eSTGT58OIcddhhbbbWVU4QlSVKLkutFnCRJa1BdXc2zzz77hSnC2267LTX/7lezIrurskuSpJbGAitJjchnn31Gp06dAPjZz37GkiVLGD16NGPGjHGKsCRJavEssJKUY0uWLOFvf/sbJSUlvP322zz66KMUFhZyzTXX0LNnTzbZZJNcR5QkSWoULLCSlCOvv/46N9xwwxemCJ900klUVlZSWFjINttsk+uIkiRJjYoFVpIa0Lx582jVqhU9evTg888/58knn2T06NGMHj2a4uJipwhLkiStgwVWkrJsyZIltasIv/TSSxx99NFcfPHFDBgwgLKyMqcIS5Ik1ZMFVpKy6LLLLqOkpITly5ezzTbbMH78eA477DAAIsLyKkmStAEssJK0Ec2bN4/HHnuME088kYigXbt2jB49mlGjRrHzzjs7RViSJOlrsMBK0te0+hThvLw8hgwZQp8+fRg/fnyu40mSJDUbFlhJ+hpeeuklTj311Nopwt///vc5/PDD6datW66jSZIkNTsWWEnaAG+99RYlJSUUFRUxduxYdthhB4466igOO+wwpwhLkiRlmQVWktZjyZIlPPTQQ5SUlPDiiy+Sl5fH2LFjAWjdujUXXnhhjhNKkiS1DBZYSVqDlFLt2dT/+q//4rHHHnOKsCRJUo5ZYCWpjlVThO+//35uuukmioqKOOOMMzj99NOdIixJkpRjFlhJLd7SpUt58MEHvzBFeJ999qG8vByA4uLiHCeUJEkSWGAltVDV1dUsWrSIzp07s3TpUi677DK23nprpwhLkiQ1YhZYSS3K22+/zZQpU7j33nvZeuutueGGG+jSpQtTpkxh6623doqwJElSI2aBldQilJaWcvPNN/PCCy/UThEePXp07f7evXvnLJskSZLqxwIrqVmqrq7m73//O7vssgtt27ZlwYIFLFq0yCnCkiRJTZgFVlKz8vbbb1NSUsLUqVN57733uPzyyxk1ahTHHnssJ5xwglOEJUmSmjALrKRmoby8nLPOOqt2ivDee+/ND37wA4YNGwZAYWFhjhNKkiTp67LASmqSqqurmTFjBm+99Rbjxo2jXbt2dO3alfPOO4/DDz+cLbbYItcRJUmStJFZYCU1KaumCN97770sXLiQbt268c1vfpOCggJ+9atf5TqeJEmSssgCK6nJuO222/jFL35RO0V4/PjxDBs2jIICfyuTJElqCfxbn6RGadUqwiUlJYwZM4Y99tiDvffe2ynCkiRJLZgFVlKjMn/+fKZMmVI7RbhDhw7sscceAGyzzTZss802OU4oSZKkXLHASsq5qqoq8vPzqa6u5jvf+Q4fffRR7RThoUOH0rp161xHlCRJUiNggZWUE3WnCP/jH/+gpKSE/Px8rrjiCrbeemunCEuSJOlLslpgI+IQ4HdAPnBDSukXaxl3JDAZGJRSej6bmSTl1nvvvcekSZO+MEX4kEMOYenSpV+YLixJkiStLmsFNiLygauBEcACYEZETE0pzVptXAfge8Bz2coiKbfKy8uprKxk00035c033+T6669nn3324fvf/z7Dhg1zirAkSZLqJZtnYPcAZqeU5gBExB3AKGDWauMuA64EfpjFLJIaWHV1NTNmzKCkpIRHHnmEY445hvHjx7PXXnvx0EMP0b1791xHlCRJUhOTzQLbE5hf5/kCYM+6AyJiN6BXSun+iFhrgY2I04HTAYqKirIQVdLGdMMNNzBx4sTaKcKHH344Bx98MAD5+fmWV0mSJH0lOVvEKSLygF8D31nf2JTSBGACQHFxccpuMkkbqry8nKeffpoRI0YAMGfOHLbZZhunCEuSJGmjymaBfQfoVef5lpltq3QAdgYeiwiA7sDUiDjChZykxm/1KcIVFRVMmTKFbbfdlssvv5y8vLxcR5QkSVIzk80COwPoGxF9qCmuxwDfWrUzpfQZ0GXV84h4DDjf8io1frNmzeK8885j4cKFtG/fnsMOO4xRo0axzTbbAFheJUmSlBVZK7AppZURcTYwjZrb6NyYUno1In4GPJ9Smpqtz5a0cZWXl/PQQw/RoUMHDjzwQLbeemu22247zjvvPIYPH+4UYUmSJDWIrF4Dm1J6AHhgtW2XrGXs0GxmkbRhqquref7555kyZUrtFOEDDzyQAw88kHbt2vHHP/4x1xElSZLUwuRsESdJjduFF17Igw8++IUpwv379891LEmSJLVgFlhJtVOE77vvPq688kq6dOnC6NGjGTp0qFOEJUmS1GhYYKUWatUU4ZKSEh5++GEqKiro3bs3CxcupEuXLuyzzz65jihJkiR9gQVWamEqKyspLCzkvffe45RTTvnSFOHMba0kSZKkRscCK7UAq6YIl5SU0KFDB/7whz/Qo0cPrrvuOgYMGECbNm1yHVGSJElaLwus1Iy98sor3HnnnV+YIjxkyJDa/U4TliRJUlNigZWamfnz59OtWzc22WQTnn32WUpLS50iLEmSpGbBAis1A0uXLq2dIvz8889z1VVXccghh3Dcccfx7W9/21WEJUmS1CxYYKUmbOnSpfz85z+vnSK89dZbc+655zJgwAAA2rVrl+OEkiRJ0sZjgZWamAULFvDmm28yZMgQ2rRpw5w5czj00EMZPXq0U4QlSZLUrFlgpSZg9SnCnTp1oqysjMLCQv76179aWiVJktQiWGClRm7q1KlcfvnlX5giPHLkSAoLCwEsr5IkSWoxLLBSI7NgwQLuvfdeBg8eTHFxMdtuu61ThCVJkiQssFKjsHTpUh5++GFKSkqYMWMGEUHbtm0pLi6u/SVJkiS1dBZYKceqq6sZPXo0Cxcu/MIU4e7du+c6miRJktSoWGClBrZgwQKmTp3KCy+8wPXXX09eXh7nnXcePXr0cIqwJEmStA4WWKkBrGmK8J577slnn33GpptuyqGHHprriJIkSVKjZ4GVsqS6uprKyko22WQTnnnmGS6++GK22morzjnnHEaOHElRUVGuI0qSJElNigVW2shWrSJcUlLCqFGjOPPMMxk8eDA333wzAwYMcIqwJEmS9BVZYKWN5L777uPuu+/+whThHXfcEYDCwkJ22223HCeUJEmSmjYLrPQVVVdX869//YuddtoJgIcffpj333/fKcKSJElSllhgpQ1Ud4rwO++8wwMPPECvXr247LLL6NChg1OEJUmSpCyxwEr1NGfOHC6//PLaKcJ77LEH3/3ud+nSpQsAHTt2zHFCSZIkqXmzwEprUV1dzQsvvADAwIED2Xzzzfnkk084++yzGTlyJD169MhxQkmSJKllscBKq3nnnXeYOnUqU6dOZcGCBey9994MHDiQTTfdlClTpuQ6niRJktRiWWClOv7nf/6H22+/vXaK8FlnncUBBxyQ61iSJEmSsMCqBUspMXPmTO69915+8IMf0LFjR3bbbTc6d+7sFGFJkiSpEbLAqsVZfYpw27ZtOfzwwxk0aBCHHHJIruNJkiRJWousFtiIOAT4HZAP3JBS+sVq+8cDpwIrgQ+Bk1NKb2Uzk1q2hQsX8o1vfIOUEnvuuWftFOG2bdvmOpokSZKk9chagY2IfOBqYASwAJgREVNTSrPqDHsRGJhSWhoRZwK/BI7OVia1LKumCJeUlJCXl8d///d/U1RUxCWXXMI+++zjFGFJkiSpicnmGdg9gNkppTkAEXEHMAqoLbAppbI6458Fjs9iHrUQ7777LiUlJV+YInzEEUeQUiIiOOqoo3IdUZIkSdJXkM0C2xOYX+f5AmDPdYw/BXhwTTsi4nTgdICioqKNlU/NyNKlS2nVqhUFBQVMnjyZ66+/3inCkiRJUjPTKBZxiojjgYHAkDXtTylNACYAFBcXpwaMpkYspcQLL7xASUkJ06ZN46qrrmLw4MF861vf4qijjnKKsCRJktTMZLPAvgP0qvN8y8y2L4iIA4EfA0NSSsuzmEfNxLJly7jpppsoKSmpnSJ88MEH156d79KlS44TSpIkScqGbBbYGUDfiOhDTXE9BvhW3QERMQC4DjgkpfRBFrOoiVu6dClz586luLiYVq1aMWXKFLbcckunCEuSJEktSNYKbEppZUScDUyj5jY6N6aUXo2InwHPp5SmAlcB7YFJEQHwdkrpiGxlUtOy+hTh1q1b88gjj1BYWMjdd99taZUkSZJamKxeA5tSegB4YLVtl9R5fGA2P19NV1lZGVdddRXz58+vnSI8atQoCgpqDlnLqyRJktTyNIpFnKSlS5fy6KOPsvPOO9OnTx/at29Pjx49OPPMM50iLEmSJAmwwCqHUkq8+OKLtVOEy8vLOeOMM/jud7/LoEGDGDRoUK4jSpIkSWpELLDKiZQSRx99NK+99toXpgjvtttuuY4mSZIkqZGywKpBVFRU8MgjjzBz5kx++tOfEhEcdthhHH/88Rx44IFOEZYkSZK0XhZYZc2apghvueWWLFq0iM0335wTTzwx1xElSZIkNSEWWG10KSUigoceeojzzz+fNm3acPDBBzN69Gh22203MrdMkiRJkqQNYoHVRrFqinBJSQn7778/J554IoMHD+bnP/+5U4QlSZIkbRQWWH0tL730Evfcc0/tFOGePXvSrl07ANq0acMRRxyR44SSJEmSmgsLrDbYokWL2GyzzQC49tprefHFF7+winBeXl6OE0qSJElqjiywqpe6U4RnzpzJtGnT6NatGz/5yU/YfPPNnSIsSZIkKesssFqnd999l+uuu+4LU4RPP/10CgsLAdhyyy1znFCSJElSS2GB1ZcsXLiQzz//nL59+5Kfn8/DDz/MiBEjnCIsSZIkKacssAK+OEX473//O3vttRcTJkxgiy224LHHHqNVq1a5jihJkiSphbPAigkTJnDjjTfWThE+88wzGTlyZO1+y6skSZKkxsAC2wItXLiQ+++/n+OOO442bdrQsWNHDjzwQEaNGsXuu+/uFGFJkiRJjZIFtoWoqKjg0UcfpaSkhOeee46UEjvssAP77bcfxxxzTK7jSZIkSdJ6WWBbgPfee4/Ro0d/aYqwKwhLkiRJakossM3QwoULuffee1mxYgVnn302W2yxBcceeyz77LOPU4QlSZIkNVkW2GZiTVOEBw8eTEqJiOB73/teriNKkiRJ0tdigW3CUkoARATXXHMNN910Ez179uSMM85g5MiR9OrVK8cJJUmSJGnjscA2Qe+99x5Tp05l6tSp/PSnP2XQoEGMGzeOwYMHO0VYkiRJUrNlgW0iKisreeihhygpKeHZZ58lpcSgQYNqy2qvXr084ypJkiSpWbPANmIpJT744AO22GILAH75y1/Spk0bpwhLkiRJapEssI3Qe++9x7333ktJSQnLly9n2rRpFBYWcuutt9KzZ0+nCEuSJElqkSywjcjMmTO57rrraqcIDxw4kFGjRlFVVUVeXp5nXCVJkiS1aBbYHEop8fLLL9O9e3e6d+9OeXk5b7/9tlOEJUmSJGkNLLAbycIlldwzazFl88pZVploXRgM692OMTt1pKhD4RfG1p0i/NZbb3Haaadx7rnnst9++/HAAw84RViSJEmS1sACuxE8/04FV0z/kJXViaqaW7NSUZmYNvtzHp1TzkWDuzKwZxtSSpx77rk8/vjjtVOETz31VEaMGAFgcZUkSZKkdchqgY2IQ4DfAfnADSmlX6y2fxPgFmB34GPg6JTSvGxm2tgWLqnkiukfsnxVc62jKkFVVeKy0nf50+itKOpQyFZbbeUUYUmSJEn6CrJWYCMiH7gaGAEsAGZExNSU0qw6w04BFqWU/iMijgGuBI7OVqZsuGfWYlZWf7m81rWyGu548QO+P7gnP/zhDxsomSRJkiQ1L9mcs7oHMDulNCeltAK4Axi12phRwM2Zx5OBAyIisphpoyubV84aTr5+QeQX8MzCqoYJJEmSJEnNVDYLbE9gfp3nCzLb1jgmpbQS+AzovPobRcTpEfF8RDy/aNGiLMX9apZVrqe9ZlTUc5wkSZIkac2axKpBKaUJKaWBKaWBm222Wa7jfEHrwvqdMG5Tz3GSJEmSpDXLZoF9B6i7StGWmW1rHBMRBUAnahZzajKG9W5H/nq6aX7AsD7tGiaQJEmSJDVT2SywM4C+EdEnIloBxwBTVxszFTgx8/gooDSl1KTm2o7ZqSMFeetusAV5wegdOzZQIkmSJElqnrJWYDPXtJ4NTANeAyamlF6NiJ9FxBGZYX8GOkfEbGA8cGG28mRLUYdCLhrclU3y40tnYvMDNskPLhrclaIOhbkJKEmSJEnNRFbvA5tSegB4YLVtl9R5vAwYm80MDWFgzzb88fAipry2mLK55VRUJtoUBsP6tGP0jh0tr5IkSZK0EWS1wLYkRR0KOXOPzpy5x5cWUZYkSZIkbQRNYhViSZIkSZIssJIkSZKkJsECK0mSJElqEiywkiRJkqQmwQIrSZIkSWoSLLCSJEmSpCbBAitJkiRJahIssJIkSZKkJiFSSrnOsEEi4kPgrVznWI8uwEe5DqEmwWNFG8LjRfXlsaIN4fGi+vJY0Yb4OsfL1imlrmva0eQKbFMQEc+nlAbmOocaP48VbQiPF9WXx4o2hMeL6stjRRsiW8eLU4glSZIkSU2CBVaSJEmS1CRYYLNjQq4DqMnwWNGG8HhRfXmsaEN4vKi+PFa0IbJyvHgNrCRJkiSpSfAMrCRJkiSpSbDAfg0RcUhE/DsiZkfEhWvYv0lE3JnZ/1xE9M5BTDUC9ThWxkfErIh4JSIejYitc5FTube+Y6XOuCMjIkWEq0G2YPU5XiJiXOb3l1cj4q8NnVGNQz3+HNoqIsoi4sXMn0WH5iKnci8iboyIDyLin2vZHxHx+8yx9EpE7NbQGdV41ON4OS5znPwjIp6OiP5f9zMtsF9RROQDVwPfAHYCjo2InVYbdgqwKKX0H8BvgCsbNqUag3oeKy8CA1NKuwCTgV82bEo1BvU8VoiIDsD3gOcaNqEak/ocLxHRF7gI2DelVAyc19A5lXv1/L3lYmBiSmkAcAxwTcOmVCNyE3DIOvZ/A+ib+XU6cG0DZFLjdRPrPl7mAkNSSv2Ay9gI18VaYL+6PYDZKaU5KaUVwB3AqNXGjAJuzjyeDBwQEdGAGdU4rPdYSSmVpZSWZp4+C2zZwBnVONTn9xWo+QPgSmBZQ4ZTo1Of4+U04OqU0iKAlNIHDZxRjUN9jpUEdMw87gS824D51IiklKYDn6xjyCjgllTjWWDTiChqmHRqbNZ3vKSUnl71ZxAb6e+4Ftivricwv87zBZltaxyTUloJfAZ0bpB0akzqc6zUdQrwYFYTqbFa77GSmarVK6V0f0MGU6NUn99btgO2i4inIuLZiFjXv5Kr+arPsXIpcHxELAAeAM5pmGhqgjb07zXSKhvl77gFGyGIpI0kIo4HBgJDcp1FjU9E5AG/Br6T4yhqOgqomeY3lJp/9Z4eEf1SSp/mMpQapWOBm1JK/xsRewP/LyJ2TilV5zqYpKYvIoZRU2D3+7rv5RnYr+4doFed51tmtq1xTEQUUDMl5+MGSafGpD7HChFxIPBj4IiU0vIGyqbGZX3HSgdgZ+CxiJgH7AVMdSGnFqs+v7csAKamlCpTSnOB16kptGpZ6nOsnAJMBEgpPQO0Bro0SDo1NfX6e420SkTsAtwAjEopfe0uZIH96mYAfSOiT0S0ombBg6mrjZkKnJh5fBRQmrzxbku03mMlIgYA11FTXr1GreVa57GSUvospdQlpdQ7pdSbmmtJjkgpPZ+buMqx+vw5NIWas69ERBdqphTPacCMahzqc6y8DRwAEBE7UlNgP2zQlGoqpgLfzqxGvBfwWUppYa5DqXGKiK2Au4ETUkqvb4z3dArxV5RSWhkRZwPTgHzgxpTSqxHxM+D5lNJU4M/UTMGZTc3FzcfkLrFypZ7HylVAe2BSZp2vt1NKR+QstHKinseKBNT7eJkGHBQRs4Aq4Icb41+/1bTU81j5AXB9RHyfmgWdvuM/urdMEXE7Nf/w1SVzTfRPgUKAlNKfqLlG+lBgNrAUOCk3SdUY1ON4uYSaNYCuyfwdd2VK6WvNHAt/b5IkSZIkNQVOIZYkSZIkNQkWWEmSJElSk2CBlSRJkiQ1CRZYSZIkSVKTYIGVJEmSJDUJFlhJUrMTEZ9/zddPjohtNlaeenze2RFxckN93oaIiJsi4qg1bP9ORPxxLa/5Wj9/SZLWxgIrSVIdEVEM5KeU5jTgx94InLOx3zQivN+7JKlZscBKkpqtqHFVRPwzIv4REUdntudFxDUR8a+IeDgiHqhzlvE4oCQzbmxE/Drz+HsRMSfzeJuIeCrz+JKImJH5jAmZz9w2Il6ok6PvqucR8YuImBURr0TErwBSSkuBeRGxxxq+w00R8aeIeD4iXo+IwzPb8zPfbUbmvf4zs31oRDwREVOBWWt4v88j4jcR8WpEPBoRXTPbt42Iv0XEzMzrd6jzsgNX//yMXhHxWES8ERE/Xcv/Bj+sk/G/1/e/mSRJ62KBlSQ1Z98EdgX6AwcCV0VEUWZ7b2An4ARg7zqv2ReYmXn8BLB/5vH+wMcR0TPzeHpm+x9TSoNSSjsDbYDDU0pvAp9FxK6ZMScBf4mIzsAYoDiltAtweZ3Pfb7OZ62uN7AHcBjwp4hoDZwCfJZSGgQMAk6LiD6Z8bsB30spbbeG92oHPJ9SKgYeB1YVzwnAOSml3YHzgWvW8/lkth0J7AKMjYiBdT8oIg4C+mbG7QrsHhGD1/IdJUlaLwusJKk52w+4PaVUlVJ6n5rCNiizfVJKqTql9B5QVuc1RcCHAJl97SOiA9AL+CswmJqi+URm/LCIeC4i/gEMB4oz228AToqIfODozGs/A5YBf46IbwJL63zuB0CPtXyPiZmsbwBzgB2Ag4BvR8RLwHNAZ2rKIsDfU0pz1/Je1cCdmce3AvtFRHtgH2BS5v2uy/wc1vX5AA+nlD5OKVUAd1Pzc63roMyvF4EXMq/riyRJX5HXxkiS9EUVQOs6z5+m5gzqv6kprSdTc8b2B5kzkdcAA1NK8yPi0jqvvYuas5ulwMyU0scAmWnCBwBHAWdTU3rJvK5iLZnSGp4HNWdMp9XdERFDgfJ6f9ua98oDPk0p7boBn7+u7bVxgCtSStdtQB5JktbKM7CSpObsCeDozPWiXak5e/p34CngyMy1sFsAQ+u85jXgP1Z7j/OpmTL8IjAMWJ5S+oz/K6sfZc5i1q7Wm1JaBkwDrgX+ApAZ0yml9ADwfWqmNq+yHfDPtXyPsZms2wLbUFOmpwFnRkRh5r23i4h29fiZ5NXJ+S3gyZTSYmBuRIzNvFdERN1sa/p8gBERsXlEtAFGU/NzrWsacHLmexMRPSOiWz0ySpK0Rp6BlSQ1Z/dQc7b0ZWrODv4opfReRNxFzVnQWcB8aqa3fpZ5zf3UFNpHMs+foGb68PSUUlVEzAf+BZBS+jQirqemeL4HzFjt82+j5prXhzLPOwAlmTO3AYyvM3Zf4NK1fI+3qSneHYEzUkrLIuIGaq5NfSEigpppz6Pr8TMpB/aIiIupmbZ8dGb7ccC1me2FwB3U/NzW9vlktt0FbAncmlJ6vu4HpZQeiogdgWcy4z8Hjs98riRJGyxSWn22jyRJzV9EtE8pfZ5ZWOnvwL6ZctuGmmti900pVX3NzzifmjOuP1nPuAHA+JTSCWvYdxNwX0pp8tfJUuf9Pk8ptd8Y7yVJUkPzDKwkqaW6LyI2BVoBl2UWbCKlVJG5JUxPas48fiURcQ+wLf93jeu6dAHWWXIlSZJnYCVJkiRJTYSLOEmSJEmSmgQLrCRJkiSpSbDASpIkSZKaBAusJEmSJKlJsMBKkiRJkpoEC6wkSZIkqUn4//sAaHVI0vzGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(16, 6))\n",
    "\n",
    "plt.plot(logwayspp, df_H.T.values, '--', c='black')\n",
    "plt.plot(logwayspp, df_H.T.values, 'o', ms=10)\n",
    "\n",
    "plt.title('Entropy in Buckets')\n",
    "plt.xlabel('log(ways) per pebble')\n",
    "plt.ylabel('entropy')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c75a7a9",
   "metadata": {},
   "source": [
    "### R Code 10.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1c47744a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0, 1.0, 1.0, 1.0]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Build list of canditate distributions\n",
    "p = [\n",
    "    [1/4, 1/4, 1/4, 1/4],\n",
    "    [2/6, 1/6, 1/6, 2/6],\n",
    "    [1/6, 2/6, 2/6, 1/6],\n",
    "    [1/8, 4/8, 2/8, 1/8],\n",
    "]\n",
    "\n",
    "# Compute the expected values of each\n",
    "result = [np.sum(np.multiply(p_i, [0, 1, 1, 2])) for p_i in p]\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33067d14",
   "metadata": {},
   "source": [
    "### R Code 10.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0810f027",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3862943611198906\n",
      "1.3296613488547582\n",
      "1.3296613488547582\n",
      "1.2130075659799042\n"
     ]
    }
   ],
   "source": [
    "# Compute the entropy of each distribution\n",
    "\n",
    "for p_i in p:\n",
    "    print(-np.sum(p_i * np.log(p_i)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "957ee617",
   "metadata": {},
   "source": [
    "### R Code 10.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "93d786e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.09, 0.21, 0.21, 0.49])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 0.7\n",
    "\n",
    "A = [\n",
    "    (1-p)**2,\n",
    "    p*(1-p),\n",
    "    (1-p)*p,\n",
    "    (p)**2,\n",
    "]\n",
    "\n",
    "np.round(A, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c96f28",
   "metadata": {},
   "source": [
    "### R Code 10.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8c4be3b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.221728604109787"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- np.sum(A * np.log(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc6b416a",
   "metadata": {},
   "source": [
    "### R Code 10.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8ae91c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_p(G=1.4):\n",
    "    x = np.random.uniform(0, 1, size=4)\n",
    "    x[3] = 0  # Removing the last random number x4\n",
    "    \n",
    "    x[3] = ( G * np.sum(x) - x[1] - x[2] ) / (2 - G)\n",
    "    \n",
    "    p = x / np.sum(x)\n",
    "    \n",
    "    return [-np.sum(p * np.log(p)), p]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdfdd872",
   "metadata": {},
   "source": [
    "### R Code 10.10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1cb0f5dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAGDCAYAAADd4sgrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbFklEQVR4nO3debSlVXkn4N8rhQOiUopGGUQ72hq6UXBVq8QYjcbEASUxdgKxsB3pZdQYE5OgbSd20t3aSZbRtFNQCSC2GFEM0g5xwhFtC8EJQWmEgBOgpYgDBHz7j3Oqcy2KuqeK2vfcc3mete6653zf/r79nlp73bq/u/fZp7o7AAAAwBg3m3cBAAAAsJYJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQBsU1V9saoeupr7raquqnuMrQgAbhzBGwCWqKqLqupHVXXVkq9XznjtGVX19NE17gpVdbdpaN3yGr9VVa+uqt23tOnuf9PdZ6x0bfPqFwBGEbwB4Poe2917Lvl69q64aVWt2xX32cX26u49kxyU5NAkz5pzPQCw5gjeADCjqnpyVX2sqv6qqjZX1Ver6lHTc/8tyYOTvHLpLPl0VvlZVfWVJF+ZHntGVV1QVd+pqtOqap8lfXRV/W5VXVhVV1TVX1bVzarq5tP2By1pe6eq+mFV3fHGvrbuvizJ+5IcuOT+F1XVL08fv7iq/r6qTqyq70+Xg29Y0vbnpjP+352ee9ySc8dPZ9PfPf23+XhV3bmqXj79dzyvqg65gX7vX1VnTu/7jap6ZVXd/Ma+XgBYSYI3AOyYByQ5P8neSf4iyRuqqrr7PyX5aJJnb2OW/Nem1x1YVQ9L8pIkv5nkLkkuTnLyVn38epINSe6X5PAkT+3ua6btNi5pd2SSD3T35Tf2RU3D/68m+eR2mj1uWsNeSU5LsuWPC7sneWeSf0xypyTPSfKmqrrXkmt/M8mLMvl3uzrJmUk+M31+SpKX3UCf1yV53rTdoUkenuR3dvT1AcA8Cd4AcH3vmM6wbvl6xpJzF3f367r7uiQnZBKef2aZ+72ku7/T3T9K8sQkx3X3Z7r76iQvSHJoVd1tSfv/MW3/T0lenknAzrS/I6uqps+PSvLGG/NCk1xRVd9N8rUkP8gkBN+Qj3X3u6av/Y1J7js9/sAkeyZ5aXdf090fTHL6krqT5NTuPqu7f5zk1CQ/7u4Tp/d6S5JDsg3Taz7Z3dd290VJ/jbJQ3b2xQLAPAjeAHB9v9bdey35et2Sc9/c8qC7fzh9uOcy97tkyeN9Mpnl3nKPq5J8O8m+N9D+4uk16e5PJflhkodW1b2T3COTmefrmS733rJx2oO3U9ve3b1Xkj2SfDzJe7fT9ptLHv8wyS2n71vfJ8kl3f2Trepe+pq+teTxj7bxfJv/hlX1r6vq9Kr6ZlVdmeS/ZzL7DQALQ/AGgF2nZzj+9SQHbHlSVbdOcodMZpy32H/J47tOr9nihEyWmx+V5JTpDPL1O5zsDL5lc7iPLlv4ZDb++CQPrKodDbZfT7J/VS39veKu+enXtLNek+S8JPfs7tsmeWGS2v4lALC6CN4AsOt8K8m/WqbNm5M8paoOrqpbZDKD+6npMuot/rCq1lfV/kmem8lS7C1OyuQ94BuTnLirCp/WclQms9rf3sHLt8zE/1FV7T79DO7H5vrvXd8Zt0lyZZKrprP8z9wF9wSAFSV4A8D1vXOrz/E+dcbrXpHkCdOduv9mWw26+/1J/nOStyX5RpKfTXLEVs3+IclZSc5J8r+TvGHJ9ZdksilZZ7KZ24313aq6KpM/Ghya5HHdfUMz99s03fjtsUkeleSKJK9O8qTuPm8X1Pf8JL+d5PtJXpef/iMEACyE2sH/WwGAgaqqM1lWfcF22hyX5Ovd/aKVqwwA2Fnr5l0AADC76e7nj88N7AIOAKw+lpoDwIKoqj9P8oUkf9ndX513PQDAbCw1BwAAgIHMeAMAAMBAgjcAAAAMtKo2V1u/fn3vs88+8y4DAAAAdsi55557RXffcVvnVlXw3mefffKWt/h4TgAAABbLQQcddPENnbPUHAAAAAYSvAEAAGAgwRsAAAAGErwBAABgIMEbAAAABhK8AQAAYCDBGwAAAAYSvAEAAGAgwRsAAAAGErwBAABgIMEbAAAABhoavKtqr6o6parOq6ovVdWhI/sDAACA1Wbd4Pu/Isl7uvsJVXXzJHsM7g8AAABWlWHBu6pul+QXkzw5Sbr7miTXjOoPAABgtTrspIvnXcJCOn3jAfMuYZcYudT87kkuT/J3VXV2Vb2+qm49sD8AAABYdUYG73VJ7pfkNd19SJIfJDlm60ZVdXRVbaqqTZs3bx5YDgAAAKy8kcH70iSXdvenps9PySSI/5TuPra7N3T3hvXr1w8sBwAAAFbesODd3d9McklV3Wt66OFJzh3VHwAAAKxGo3c1f06SN013NL8wyVMG9wcAAACrytDg3d3nJNkwsg8AAABYzUa+xxsAAABu8gRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgdaNvHlVXZTk+0muS3Jtd28Y2R8AAACsNkOD99QvdfcVK9APAAAArDqWmgMAAMBAo4N3J/nHqjqrqo7eVoOqOrqqNlXVps2bNw8uBwAAAFbW6OD9C919vySPSvKsqvrFrRt097HdvaG7N6xfv35wOQAAALCyhgbv7v7a9PtlSU5Ncv+R/QEAAMBqMyx4V9Wtq+o2Wx4n+ZUkXxjVHwAAAKxGI3c1/5kkp1bVln7+V3e/Z2B/AAAAsOoMC97dfWGS+466PwAAACwCHycGAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAOtm3cBAADA6nDYSRfPu4SFdPrGA+ZdAqucGW8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgYYH76rararOrqrTR/cFAAAAq81KzHg/N8mXVqAfAAAAWHWGBu+q2i/JY5K8fmQ/AAAAsFqNnvF+eZI/SvKTwf0AAADAqjQseFfVYUku6+6zlml3dFVtqqpNmzdvHlUOAAAAzMXIGe8HJXlcVV2U5OQkD6uqk7Zu1N3HdveG7t6wfv36geUAAADAyhsWvLv7Bd29X3ffLckRST7Y3RtH9QcAAACrkc/xBgAAgIHWrUQn3X1GkjNWoi8AAABYTcx4AwAAwECCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQII3AAAADCR4AwAAwECCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQII3AAAADCR4AwAAwECCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQII3AAAADCR4AwAAwECCNwAAAAwkeAMAAMBAMwXvqnp7VT2mqgR1AAAA2AGzBulXJ/ntJF+pqpdW1b0G1gQAAABrxkzBu7vf391PTHK/JBcleX9VfaKqnlJVu48sEAAAABbZzEvHq+oOSZ6c5OlJzk7yikyC+PuGVAYAAABrwLpZGlXVqUnuleSNSR7b3d+YnnpLVW0aVRwAAAAsupmCd5LXdfe7lh6oqlt099XdvWFAXQAAALAmzLrU/L9u49iZu7IQAAAAWIu2O+NdVXdOsm+SW1XVIUlqeuq2SfYYXBsAAAAsvOWWmv9qJhuq7ZfkZUuOfz/JCwfVBAAAAGvGdoN3d5+Q5ISq+o3uftsK1QQAAABrxnJLzTd290lJ7lZVv7/1+e5+2TYuAwAAAKaWW2p+6+n3PUcXAgAAAGvRckvN/3b6/b+sTDkAAACwtsz0cWJV9RdVdduq2r2qPlBVl1fVxtHFAQAAwKKb9XO8f6W7r0xyWJKLktwjyR+OKgoAAADWilmD95Yl6Y9J8tbu/t6gegAAAGBNWW5ztS1Or6rzkvwoyTOr6o5JfjyuLAAAAFgbZprx7u5jkvx8kg3d/c9JfpDk8JGFAQAAwFow64x3ktw7k8/zXnrNibu4HgAAAFhTZgreVfXGJD+b5Jwk100PdwRvAAAA2K5ZZ7w3JDmwu3tkMQAAALDWzLqr+ReS3HlkIQAAALAWzTrjvXeSc6vq/yS5esvB7n7cDV1QVbdM8pEkt5j2c0p3/+mNqBUAAAAWzqzB+8U7ce+rkzysu6+qqt2TfKyq3t3dn9yJewEAAMBCmil4d/eHq+qAJPfs7vdX1R5Jdlvmmk5y1fTp7tMv7xEHAADgJmXWXc2fkeToJLfPZHfzfZO8NsnDl7lutyRnJblHkld196e20ebo6b1zl7vcZUdqBwDgJuawky6edwkL6/SNB8y7BLjJmnVztWcleVCSK5Oku7+S5E7LXdTd13X3wUn2S3L/qvq322hzbHdv6O4N69evn7lwAAAAWASzBu+ru/uaLU+qal12YNl4d383yYeSPHKHqgMAAIAFN2vw/nBVvTDJrarqEUnemuSd27ugqu5YVXtNH98qySOSnHcjagUAAICFM2vwPibJ5Uk+n+Q/JnlXkhctc81dknyoqj6X5NNJ3tfdp+9soQAAALCIZt3V/CdV9Y4k7+juy2e85nNJDrkRtQEAAMDC2+6Md028uKquSHJ+kvOr6vKq+pOVKQ8AAAAW23JLzZ+XyW7m/667b9/dt0/ygCQPqqrnDa8OAAAAFtxywfuoJEd291e3HOjuC5NsTPKkkYUBAADAWrBc8N69u6/Y+uD0fd67jykJAAAA1o7lgvc1O3kOAAAAyPK7mt+3qq7cxvFKcssB9QAAAMCast3g3d27rVQhAAAAsBYtt9QcAAAAuBEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICBBG8AAAAYSPAGAACAgYYF76rav6o+VFXnVtUXq+q5o/oCAACA1WrdwHtfm+QPuvszVXWbJGdV1fu6+9yBfQIAAMCqMmzGu7u/0d2fmT7+fpIvJdl3VH8AAACwGo2c8f7/qupuSQ5J8qmV6A8AYJ4OO+nieZewsE7feMC8SwDY5YZvrlZVeyZ5W5Lf6+4rt3H+6KraVFWbNm/ePLocAAAAWFFDg3dV7Z5J6H5Td799W226+9ju3tDdG9avXz+yHAAAAFhxI3c1ryRvSPKl7n7ZqH4AAABgNRs54/2gJEcleVhVnTP9evTA/gAAAGDVGba5Wnd/LEmNuj8AAAAsguGbqwEAAMBNmeANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADCd4AAAAwkOANAAAAAwneAAAAMJDgDQAAAAMJ3gAAADCQ4A0AAAADrZt3AQDAfBx20sXzLmFhnb7xgHmXAMACGTbjXVXHVdVlVfWFUX0AAADAajdyqfnxSR458P4AAACw6g0L3t39kSTfGXV/AAAAWARz31ytqo6uqk1VtWnz5s3zLgcAAAB2qbkH7+4+trs3dPeG9evXz7scAAAA2KXmHrwBAABgLRO8AQAAYKCRHyf25iRnJrlXVV1aVU8b1RcAAACsVutG3bi7jxx1bwAAAFgUlpoDAADAQII3AAAADCR4AwAAwECCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQII3AAAADCR4AwAAwECCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQII3AAAADCR4AwAAwECCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQII3AAAADCR4AwAAwEDr5l0AAGzPYSddPO8SFtbpGw+YdwkAQMx4AwAAwFCCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQII3AAAADCR4AwAAwECCNwAAAAwkeAMAAMBAgjcAAAAMJHgDAADAQOvmXQDAWnDYSRfPu4SFdPrGA+ZdAgDAcGa8AQAAYCDBGwAAAAYSvAEAAGAgwRsAAAAGErwBAABgIMEbAAAABhK8AQAAYCCf4w03IT5reuf4rGkAAG6MoTPeVfXIqjq/qi6oqmNG9gUAAACr0bAZ76raLcmrkjwiyaVJPl1Vp3X3uaP6ZG0wK7tzzMoCAMDqNHLG+/5JLujuC7v7miQnJzl8YH8AAACw6ox8j/e+SS5Z8vzSJA8Y2N+KMiu7c8zKAgAANzXV3WNuXPWEJI/s7qdPnx+V5AHd/eyt2h2d5Ojp03slOX9IQay0vZNcMe8iWLOML0YyvhjNGGMk44uRjK/tO6C777itEyNnvL+WZP8lz/ebHvsp3X1skmMH1sEcVNWm7t4w7zpYm4wvRjK+GM0YYyTji5GMr5038j3en05yz6q6e1XdPMkRSU4b2B8AAACsOsNmvLv72qp6dpL3JtktyXHd/cVR/QEAAMBqNHKpebr7XUneNbIPVi1vH2Ak44uRjC9GM8YYyfhiJONrJw3bXA0AAAAY+x5vAAAAuMkTvNlpVfXIqjq/qi6oqmO2cf6vq+qc6deXq+q7cyiTBTbDGLtrVX2oqs6uqs9V1aPnUSeLaYbxdUBVfWA6ts6oqv3mUSeLqaqOq6rLquoLN3C+qupvpuPvc1V1v5WukcU1w/i6d1WdWVVXV9XzV7o+FtsM4+uJ059bn6+qT1TVfVe6xkUkeLNTqmq3JK9K8qgkByY5sqoOXNqmu5/X3Qd398FJ/meSt694oSysWcZYkhcl+fvuPiSTT0549cpWyaKacXz9VZITu/s+Sf4syUtWtkoW3PFJHrmd849Kcs/p19FJXrMCNbF2HJ/tj6/vJPndTH6OwY46PtsfX19N8pDuPijJn8f7vmcieLOz7p/kgu6+sLuvSXJyksO30/7IJG9ekcpYK2YZY53kttPHt0vy9RWsj8U2y/g6MMkHp48/tI3zcIO6+yOZhJ8bcngmf9jp7v5kkr2q6i4rUx2Lbrnx1d2Xdfenk/zzylXFWjHD+PpEd2+ePv1kEivCZiB4s7P2TXLJkueXTo9dT1UdkOTu+ZdfYGEWs4yxFyfZWFWXZvIJCs9ZmdJYA2YZX59N8vjp419PcpuqusMK1MZNw8z/jwKsYk9L8u55F7EIBG9WwhFJTunu6+ZdCGvOkUmO7+79kjw6yRurys81dpXnJ3lIVZ2d5CFJvpbEzzEASFJVv5RJ8P7jedeyCIZ+jjdr2teS7L/k+X7TY9tyRJJnDa+ItWaWMfa0TN+D1N1nVtUtk+yd5LIVqZBFtuz46u6vZzrjXVV7JvmN7v7uShXImrcj/48CrCpVdZ8kr0/yqO7+9rzrWQRmhthZn05yz6q6e1XdPJNwfdrWjarq3knWJzlzhetj8c0yxv4pycOTpKp+Lsktk1y+olWyqJYdX1W195IVFC9IctwK18jadlqSJ013N39gku919zfmXRTAcqrqrplsmnxUd3953vUsCjPe7JTuvraqnp3kvUl2S3Jcd3+xqv4syabu3vIL7BFJTu7unletLKYZx9gfJHldVT0vk43WnmysMYsZx9dDk7ykqjrJR2LlDjugqt6cyRjae7oPxZ8m2T1Juvu1mexL8egkFyT5YZKnzKdSFtFy46uq7pxkUyYbkP6kqn4vyYHdfeV8KmaRzPDz60+S3CHJq6sqSa7t7g3zqXZxlN9RAQAAYBxLzQEAAGAgwRsAAAAGErwBAABgIMEbAAAABhK8AQAAYCAfJwYAq1xVXZfk80sOndzdL91O+4cmuaa7PzG4NABgBoI3AKx+P+rug3eg/UOTXJXkesG7qtZ197W7qC4AYAY+xxsAVrmquqq799zG8YuSnJDksUl2T/Lvk/w4ySeTXJfk8iTPSfK06fFDknw8yYlJXptkjyT/N8lTu3tzVZ2R5LNJHpLJH+efmmRTkvOT/Hx3X15VN0vy5SSHdvflg14yAKwp3uMNAKvfrarqnCVfv7Xk3BXdfb8kr0ny/O6+KJNQ/dfdfXB3f3Tabr9MwvPvZxK8/7i775PJEvY/XXK/Paaz67+T5Lju/kmSk5I8cXr+l5N8VugGgNlZag4Aq9/2lpq/ffr9rCSP38493trd11XV7ZLs1d0fnh4/Iclbl7R7c5J090eq6rZVtVeS45L8Q5KXZzIL/nc78yIA4KbKjDcALLarp9+vy/b/oP6DGe+39XvQursvSfKtqnpYkvsnefeOlQgAN22CNwCsPd9Pcpttneju7yXZXFUPnh46KsmHlzT5rSSpql9I8r1p+yR5fSZLzt/a3dcNqRoA1ihLzQFg9btVVZ2z5Pl7uvuY7bR/Z5JTqurwTDZX29p/SPLaqtojyYVJnrLk3I+r6uxMNmt76pLjp2WyxNwycwDYQXY1BwCSJNNdzZ/f3Zu2cW5DJhu2Pfh6FwIA22XGGwDYrqo6Jskz8y87mwMAO8CMNwAAAAxkczUAAAAYSPAGAACAgQRvAAAAGEjwBgAAgIEEbwAAABhI8AYAAICB/h/pZwgbFV+hawAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1224x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "H = pd.DataFrame([ sim_p(1.4) for _ in range(10000)], columns=('entropies', 'distributions'))\n",
    "\n",
    "plt.figure(figsize=(17, 6))\n",
    "\n",
    "plt.hist(H.entropies, density=True, rwidth=0.9)\n",
    "\n",
    "plt.title('Entropy - Binomial')\n",
    "plt.xlabel('Entropy')\n",
    "plt.ylabel('Density')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "914a55af",
   "metadata": {},
   "source": [
    "### R Code 10.11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d11e9d8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# entropies = H.entropies\n",
    "# distributions = H.distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d46899f3",
   "metadata": {},
   "source": [
    "### R Code 10.12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b924401c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2217235826761963"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.entropies.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9386d7e",
   "metadata": {},
   "source": [
    "### R Code 10.13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8fdfa11a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([array([0.08969631, 0.20938934, 0.21121804, 0.48969631])],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.loc[H.entropies == H.entropies.max(), 'distributions'].values\n"
   ]
  }
 ],
 "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}