chap02.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from abc import ABC, abstractmethod\n",
    "import random\n",
    "import math\n",
    "\n",
    "\n",
    "# 離散確率分布の抽象クラス（現在は1次元のみ）\n",
    "class DiscreteProb(ABC):\n",
    "    @abstractmethod\n",
    "    def pmf(self, x: int):\n",
    "        # 確率量子関数\n",
    "        pass\n",
    "\n",
    "    @abstractmethod\n",
    "    def __call__(self):\n",
    "        # 確率変数\n",
    "        pass\n",
    "\n",
    "\n",
    "class Bern(DiscreteProb):\n",
    "    def __init__(self, mu: float):\n",
    "        self.mu = mu\n",
    "\n",
    "    def pmf(self, x: int):\n",
    "        if x == 0:\n",
    "            return 1 - self.mu\n",
    "        elif x == 1:\n",
    "            return self.mu\n",
    "        else:\n",
    "            return 0\n",
    "\n",
    "    def __call__(self):\n",
    "        return 1 if random.random() < self.mu else 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# サンプリングによる期待値計算とエントロピー計算\n",
    "\n",
    "\n",
    "class SimpleSampler(object):\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    def __call__(self, p, L: int, f=None) -> float:\n",
    "        # 単純なサンプリングを用いて f(x) の期待値を計算する\n",
    "        # # p: 確率分布\n",
    "        # # f: 関数\n",
    "        # L: サンプル数\n",
    "        if f is None:\n",
    "            f = lambda x: x\n",
    "        return sum([f(p()) for _ in range(L)]) / L\n",
    "\n",
    "\n",
    "class Entropy(object):\n",
    "    def __init__(self, sampler=SimpleSampler()):\n",
    "        self.sampler = sampler\n",
    "\n",
    "    def __call__(self, p, n_samples: int = 10000) -> float:\n",
    "        return -1 * self.sampler(p, n_samples, lambda x: math.log(p.pmf(x)))\n",
    "\n",
    "\n",
    "class KLDiv(object):\n",
    "    def __init__(self, sampler=SimpleSampler()):\n",
    "        self.sampler = sampler\n",
    "\n",
    "    def __call__(self, *, p, q, L: int = 10000) -> float:\n",
    "        # p(X) と q(X) のKL距離\n",
    "        #    KL[q(x) || p(x)]]\n",
    "        # を計算する\n",
    "        return abs(-1 * self.sampler(\n",
    "            q, f=lambda x: log(p.pmf(x) / q.pmf(x)), L=L\n",
    "        ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "近似値：0.64078858\n",
      "理論値：0.63651417\n"
     ]
    }
   ],
   "source": [
    "# ================================================\n",
    "# 2.1\n",
    "# ================================================\n",
    "\n",
    "\n",
    "# 2.1.5 サンプリングによる期待値の近似計算\n",
    "# 例題：単純な確率分布のエントロピー計算\n",
    "\n",
    "from math import log\n",
    "\n",
    "\n",
    "class MyDist(DiscreteProb):\n",
    "    def pmf(self, x: int) -> float:\n",
    "        assert x in [0, 1]\n",
    "        if x == 0:\n",
    "            return 2.0 / 3\n",
    "        else:\n",
    "            return 1.0 / 3\n",
    "\n",
    "    def __call__(self) -> int:\n",
    "        if random.random() < 2.0 / 3:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1\n",
    "\n",
    "\n",
    "p = MyDist()\n",
    "H = Entropy()\n",
    "print(f\"近似値：{H(p):.8f}\")\n",
    "print(f\"理論値：{-(1/3 * log(1/3) + 2/3 * log(2/3)):.8f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0]\n",
      "[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "E[Bern(x|0.5)] = 0.5\n",
      "E[Bern(x|0.9)] = 0.8957\n",
      "H[Bern(x|0.5)] = 0.6931471805600546\n",
      "H[Bern(x|0.9)] = 0.3220068589832419\n"
     ]
    },
    {
     "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>0.1</th>\n",
       "      <th>0.5</th>\n",
       "      <th>0.9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.1</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.512364</td>\n",
       "      <td>1.760416</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.5</th>\n",
       "      <td>0.368372</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.370635</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.9</th>\n",
       "      <td>1.766085</td>\n",
       "      <td>0.509353</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          0.1       0.5       0.9\n",
       "0.1  0.000000  0.512364  1.760416\n",
       "0.5  0.368372  0.000000  0.370635\n",
       "0.9  1.766085  0.509353  0.000000"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG2CAYAAACTTOmSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABbUUlEQVR4nO3deVzUdeI/8NcczAwgh9yIKILgiaCohOZRUqZ+Pba26FKzstbs2Ni2sjbt2E0z82eHm5tptV26ld2mJml5kBSIiiIKKCDCcCk3M8zM5/cHR5EXA8O853g9H495FONnZl7jlPPy83kfMkmSJBARERE5IbnoAERERESisAgRERGR02IRIiIiIqfFIkREREROi0WIiIiInBaLEBERETktFiEiIiJyWixCRERE5LRYhIiIiMhpsQgRERGR07KJIrR27VqEhYVBo9EgPj4eaWlplzx28uTJkMlkF9xmzJhhxcRERETkCIQXoc2bNyM5ORnLli1DRkYGYmJiMHXqVJSVlV30+C1btqCkpKT9lpWVBYVCgZtvvtnKyYmIiMjeyURvuhofH48xY8bgjTfeAACYTCaEhobioYcewpNPPnnFx69ZswZLly5FSUkJ3N3dezouERERORClyBfX6/VIT0/HkiVL2u+Ty+VITExEampqp55jw4YNuPXWWy9bgnQ6HXQ6XfvPJpMJVVVV8PX1hUwm6/obICIiIquRJAm1tbXo06cP5HLLXNQSWoQqKipgNBoRGBjY4f7AwEAcP378io9PS0tDVlYWNmzYcNnjli9fjueee65bWYmIiMg2FBUVoW/fvhZ5LqFFqLs2bNiA6OhojB079rLHLVmyBMnJye0/V1dXo1+/figqKoKnp2dPxyQiIiILqKmpQWhoKDw8PCz2nEKLkJ+fHxQKBbRabYf7tVotgoKCLvvY+vp6bNq0Cc8///wVX0etVkOtVl9wv6enJ4sQkZ0xmiTUNRnQ2GxEY7MRTa3/1BtMMEkSIAEmCZDQMvxRrVRA4yKHWqmAWimHm1oBX3c1FHJeFieyV5Yc1iK0CKlUKsTFxSElJQVz5swB0DJ+JyUlBQ8++OBlH/vJJ59Ap9PhzjvvtEJSIrIGSZKgrdHhVEU9TlfW43RFPc6cb0RlnQ6VdXpU1etxrkEPUzeneMhlgF8vNQI81Qj00CDQS4PQ3m7o59N683WDl6uLZd4UEdk04ZfGkpOTMX/+fIwePRpjx47FmjVrUF9fjwULFgAA5s2bh5CQECxfvrzD4zZs2IA5c+bA19dXRGwi6ia9wYQT2lpkFVfjSHE1ss7W4ERpLRqbjZ16vFIug6uLAmoXBVxVcqgUcshlMshlMshkLX9jlCQJeoMJTc1G6Fr/2dBshEkCymp1KKvVIQs1F31+D40SwV4aBHq23II8Nejb2xUJEb7o78sZqkSOQngRSkpKQnl5OZYuXYrS0lLExsZi27Zt7QOoCwsLLxgZnpOTg71792LHjh0iIhORmXQGI3JKa1sKT2vxySmtRbPxwlM7CrkMfXu7YoCfO8J83RHq4wZ/DzV83VXw7aWCj7sKvd1UcFF0bcaIwWhCVb0e2hodymqbUFarQ8n5RhRWNbTeGlFRp0NtkwG1TXU4oa274Dn6+bhhYpQfJkT6Y1yELzw0PHtEZK+EryMkQk1NDby8vFBdXc0xQkQWZjCacLKsDofPnEdmUTUOnzmPE9qLlx4vVxcMD/HE8BAvDO/jhaF9PNHPx63LJcdSGvQGFJ9rRGlNE0qrW8pSaXUTcrS1yCg4B8Pvrs0p5DKM6OuFhHBfJET4YnR/H7iqFALTEzmunvj+ZhFiESLqFkmScLy0Fj+dKMdPJ8uRXnAOTc2mC47r7eaC4SFeiA7xav9n396udreWV53OgNS8Suw5WY6fTpTjdGVDh193Ucgwsl9vXD80ENOjg9HH21VQUiLHwyJkISxCRN3T1GzE7pxy7MzWYs/JcmhrdB1+3UOtxPAQL4wI9UJMX2+7LT2dceZcA1LzKpGaX4mf8ypxtrqpw6/HhnpjenQQpg0PRqiPm6CURI6BRchCWISIzGc0STiQX4kvMovxXVYpapsM7b/m6qLAVeE+mBjlj6sH+iHCvxfkTjg9XZIkFFY1YNfxMmzNKsUvp6vw+z9hY0K9MXNEMM8UEXURi5CFsAgRdV69zoD1e/Lx0YFClNX+duYnyFODGSOCcc2gAIwO6w2NC8fF/FFZbRO2H9Vi6+ESHDhV2WHa/+j+vTFjRDAmRvkj3M/dIc+WEVkai5CFsAgRXZnRJOGzjDNYtT2nvQB5ubpgenQwZsf2wdgwH6c869NV5bU6fJdVgm8OlSDtdFWHXwv0VCMh3BfjIvwwbqAv+vbmJTSii2ERshAWIaLL259bgX9+m41jJS1r7PTzccNjUwfhhmFBUCnFzuhyBCXVjfj2cAl2ZmuRUXgeekPHweUTo/xx/8RwjIvgxtBEv8ciZCEsQkQd1ekM+OV0FX7Or0RqXiUOn6kG0LKo4MPXRmLeuP5QK3npqyc0NRuRUXiuZcB1XiUyCs+1X0IbGuyJ+yaGY8aIYOFLChDZAhYhC2ERIgJMJgkfphXi0/QzyCquhvEPa+PcGd8PjyRGwcddJTCl8ymsbMDGfaew+Zei9lW2+3hp8MS0wZgV04dniMipsQhZCIsQObuiqgb8/dND+Dn/t7Eq/XzccFW4D65qHasS5KURmJDO1evx4YECvLu/ABV1LWO0JkT64YXZwxHmxy0+yDmxCFkIixA5K0mS8L9fi/DCN9mo0xng6qLA366PwrToYIRwOrdNamo24j8/5mPt7lzoDSaolHI8eM1A3D8pnJcryemwCFkIixA5o+LzjfjH50ewK6ccADAmrDdW3RzDDUTtxKmKejzzRRb25lYAAML93bF48kD8X0wwCxE5DRYhC2ERImdRVtOE77JK8e3hEvxS0LK4n0opx9+vH4S7rx4ABae/2xVJkvDVobN44Zvs9stlfr1UuG1sP9wR35+XM8nhsQhZCIsQObK2L8uPDhQi7Q8rG18V7oPnZw9HVKCHuIDUbdWNzfjwQAHeTy1ASeuWHkq5DNOig/FoYiTC/XsJTkjUM1iELIRFiBxVVnE1ln11FOkF59rvG9XPGzNG9MG04UHc1sHBGIwm7Dimxbv7TyPtVMvAd5VCjoUTB+DBayLhquIlM3IsLEIWwiJEjqaqXo+Xt+dg0y+FkCTATaXA/RMj8OfRfTkI2kkcPVuNl7fnYHfrGLAQb1csnTkU1w8N5JR7chgsQhbCIkSOwmSS8FFaIV7enoPqxmYAwOzYPlgybQjHizghSZKw/agWL3xzDMXnGwEA1wzyx/Ozh3Pne3IILEIWwiJEjqCwsgGPf/bbWkBDgj3x3KxhGDvAR3AyEq1Bb8DaXbl466d8NBul9mUSFoznAHmybyxCFsIiRPbMZJLw/s8FWPHdcTQ2G+HqosDjNwzCvIQwfslRB3nldXhqyxEcaB0/FNPXCytuGoEhwfxzj+wTi5CFsAiRvTpdUY/HPzvcPjA2foAPXv5zDPr58rIHXZzJJGHzr0V4cWs2apsMUMpl+MukCDw8JZIb6JLdYRGyEBYhsjdNzUas+zEP/96dB73BBDeVAkumDcYd8f0h51kg6gRtTROWfpmF7Ue1AICxYT5YNzeOe8mRXWERshAWIbInPxzX4tmvjqGwqgEAcPVAPyy/MZqDX6lLvjtSgsc/PYxanQH9fNyw8a7RGBjAdaXIPrAIWQiLENmDoqoGPPf1MezMbvkbfJCnBs/831BMjw7idGjqlpPaWtz93i8oqmqEh0aJtbePwsQof9GxiK6oJ76/eYGYyAZlFVdj1ht7sTNbC6VchvsnhiPlb5MwY0QwSxB1W2SgB754YDzGhPVGbZMBC979Bf9NPS06FpEQLEJENuZQ0Xncvv5nnGtoRnSIF757ZAKWTB8Cd7VSdDRyIL691Pjg3njcNKovjCYJS788irveScMJba3oaERWxSJEZEPSC87hzrcPoKbJgLj+vfHhwnhEcl8w6iFqpQKrbh6BJ6cNhotCht055bhhzU946vMjKK/ViY5HZBUcI8QxQmQj0k5VYcE7aajXGzF2gA823jUGvXgWiKzkdEU9Vnx3HNuOlgIA3FUKLJocgYUTw6FWcs8ysg0cI0TkoH46UY75G1tK0PiBvnh3AUsQWVeYnzvWzY3D/+5PQExfL9TrjVi14wTu2vhL+/YtRI6IZ4R4RogEyi6pwSs7crAzuwwAMDHKH2/NjYPGhX8DJ3FMJglfHirGM18cRZ3OgEGBHnj37jEI9uIGviQWp89bCIsQiZZfXof/t/Mkvj50FgAglwFJY/ph2cyhLEFkM46ercaCd35BWa0OwV4avLtgLAYFccwaicMiZCEsQiRKU7MRy7dm44MDhTCaWv7X+78RwXj0uihE+PcSnI7oQmfONWD+xjTkldfDQ6PEW3NHIyHCV3QsclIsQhbCIkQinD3fiEUfZuBQ0XkAQOKQACRfNwhD+/C/QbJt5xv0WPjfX/HL6XNQKeRYNmsobh/bj2takdWxCFkIixBZ28/5lVj8YQYq6/XwdnPBmqRYTB4UIDoWUac1NRvx6OZMfJfVMqvshmFBWHFTNLzduFcZWQ9njRHZGUmSsGHvKdzx9gFU1usxNNgTXz94NUsQ2R2NiwJrbx+Fp6cPgYtChm1HSzHt1T04kF8pOhpRt7AIEfUQvcGEv/3vEF745hiMJgl/GhmCzxaN42apZLfkchkWTgzHlkXjMcDPHSXVTbht/c9YvSMHBqNJdDyiLmERIuoBdToD7n73F2w5WAyFXIZlM4di9S0xcFVxRhjZv+i+Xvjmoatxc1xfmCTgtR9yMW9jGqrq9aKjEZmNRYjIwsprdbj1rVTsza2Am0qBjXeNwYLxAziwlByKu1qJl2+Owau3xsJNpcD+vErMfH0vsoqrRUcjMguLEJEFna6ox01v7kdWcQ183VX4eOFVmBTlLzoWUY+ZHRuCLxa3XCorPt+Im97cjy0ZZ0THIuo0FiEiCzlypho3vbkfhVUNCPVxxaeLxiEm1Ft0LKIeFxXogS8Wj8e1gwOgM5iQ/L9DeParo2jmuCGyAyxCRN1kMrXMDLtp3f72mWGfLRqHAX7uoqMRWY2XqwvenjcaD187EADw7v7TuOudNFQ3cJ8ysm0sQkTdUFrdhHkb0/DCN8egN5hw7eAAbL7/KgR4aERHI7I6uVyG5OsH4T9z4+CmUmBfbiX+9O99OFVRLzoa0SWxCBF10TeHz2Lqmp+wN7cCGhc5XpgzHBvmj4aHxkV0NCKhpg4Lwqd/GYc+XhrkV9Rjztp92J9bIToW0UWxCBGZSZIkLNlyBA9+dBDVjc0Y0dcL3z48AXOv6s+ZYUSthvbxxBcPjsfIft6obmzGvI1p+OhAoehYRBewiSK0du1ahIWFQaPRID4+HmlpaZc9/vz581i8eDGCg4OhVqsRFRWFrVu3WiktObsNe0/h47RCyGXAw9cOxGeLxnHDVKKLCPDQ4OOFV2FWTB8YTBKe+vwI1u7KFR2LqAOl6ACbN29GcnIy1q1bh/j4eKxZswZTp05FTk4OAgIu3IZAr9fjuuuuQ0BAAD799FOEhISgoKAA3t7e1g9PTueX01VY/t1xAMCzs4ZhXkKY2EBENk7josCrt8Yi3N8da3aexMvbc+Dt5oI74vuLjkYEwAY2XY2Pj8eYMWPwxhtvAABMJhNCQ0Px0EMP4cknn7zg+HXr1uHll1/G8ePH4eLStbEY3HSVuqKiTocZr+2BtkaHmTF98NqtsbwURmSGV3bk4PUfciGTAW/cNgozRgSLjkR2xuE2XdXr9UhPT0diYmL7fXK5HImJiUhNTb3oY7766iskJCRg8eLFCAwMxPDhw/Hiiy/CaDRe8nV0Oh1qamo63IjMYTRJ+OumTGhrdIjwd8fyG6NZgojMlHxdFO6I7wdJAv66+SD2nCwXHYlIbBGqqKiA0WhEYGBgh/sDAwNRWlp60cfk5+fj008/hdFoxNatW/HMM8/glVdewT//+c9Lvs7y5cvh5eXVfgsNDbXo+yDH92rKSezNrYCriwJv3hmHXmrhV5WJ7I5MJsPzs4djxohgNBsl3P9+OjKLzouORU7OJgZLm8NkMiEgIABvvfUW4uLikJSUhKeffhrr1q275GOWLFmC6urq9ltRUZEVE5O9251Thtd/OAkAWH5jNKICPQQnIrJfCrkM/++WWEyI9EOD3oi73klDTmmt6FjkxIQWIT8/PygUCmi12g73a7VaBAUFXfQxwcHBiIqKgkLx2y7eQ4YMQWlpKfT6i+98rFar4enp2eFG1BmpeZV46KODkCTgjvh+mDMyRHQkIrunUsqx7s44xIR643xDM259K5WbtZIwQouQSqVCXFwcUlJS2u8zmUxISUlBQkLCRR8zfvx45ObmwmT6bQ+bEydOIDg4GCqVqsczk/P4+tBZzN+YhlqdAWMH+OCZ/xsqOhKRw3BXK/HegjGI6euFcw3NuO2tn5FeUCU6Fjkh4ZfGkpOTsX79erz33nvIzs7GokWLUF9fjwULFgAA5s2bhyVLlrQfv2jRIlRVVeGRRx7BiRMn8O233+LFF1/E4sWLRb0FckBv78nHQx8fhN5owg3DgvDfu8dC46K48gOJqNO83VT44N54jB3gg1qdAXe+nYa9J7kCNVmX8BGfSUlJKC8vx9KlS1FaWorY2Fhs27atfQB1YWEh5PLf+lpoaCi2b9+ORx99FCNGjEBISAgeeeQRPPHEE6LeAjkQk0nCi1uz8fbeUwCAeQn9sWzmMCjknCFG1BM8NC54b8FY3P9BOn46UY673/0F/75jFBKHBl75wUQWIHwdIRG4jhBdjMkk4bFPDmHLwWIAwBM3DMZfJoVzmjyRFegMRjz88UFsP6qFUi7DG7ePxA3Duc4QdeRw6wgR2ZJXU05iy8FiKOUyrL4lBosmR7AEEVmJWqnA2ttH4U8jQ2AwSXj440zsz+NlMup5LEJEALYeKcGrKS1T5F+8MRo3juorOBGR81Eq5Fh1cwxuGBYEvdGE+/6bztlk1ONYhMjpHT1bjb/97xAA4J6rB+CW0Vxwk0gUhVyGNbfG4qpwH9TpDLjrnV9QUFkvOhY5MBYhcmrltTosfO9XNDYbMTHKH0umDRYdicjpaVwUeGveaAwJ9kRFnQ5zN6ShrLZJdCxyUCxC5LT0BhMWfZCOs9VNCPdzx+u3jYRSwf8liGyBp8YF7909Bv183FBY1YC7Nv6CmqZm0bHIAfFPfXJKBqMJT245jF8LzsFDo8T6+aPh5eoiOhYR/U6Ahwbv3zMWfr1UOFZSg0UfpENvMF35gURmYBEip1NZp8O8jWnYklEMuQx44/ZRiPDvJToWEV1Ef193vLtgLNxVCuzLrcSSLUfghKu+UA9iESKnklVcjVlv7MP+vEq4qxT49x2jMCnKX3QsIrqM4SFeeOOOUVDIZfgs40z7DE8iS2ARIqfxWfoZ3PTmfhSfb8QAP3d8sXg8F2wjshPXDArAC7OHAwDW7DyJT9PPCE5EjkL4FhtEPU2SWrbNWL+nZduMKYMDsDoplmOCiOzM7fH9cOZcA/69Ow9PfnYYQZ4aXB3pJzoW2TmeESKH99Whs+0l6OEpkVg/jwOjiezVY9cPwqyYPjCYJCz6IB3HS2tERyI7xyJEDq2kuhHPfJEFAHhkSiSSr4uCnBuoEtktuVyGl28egfjWHevv2vgLzp5vFB2L7BiLEDksk0nC3z85jJomA2L6euHBaweKjkREFqBWKvDW3NGIDOiF0pomzN+YhuoGrjFEXcMiRA7r/Z8LsDe3AhoXOVYnxcKFiyUSOQwvNxe8e/dYBHlqcLKsDgv/+yuamo2iY5Ed4jcDOaS88jos/y4bALBk2hCuE0TkgEK8XfHu3WPgoVYi7XQVHt2cCaOJawyReViEyOE0G01I3pyJpmYTJkT6Ye5V/UVHIqIeMjjIE/+ZFweVQo7vskrx/NdHueAimYVFiBzO2l25OHSmGp4aJV7+cwwHRxM5uHERfnjllhgAwHupBXjzxzzBiciesAiRQ9mfV4HXWledfWHOcAR5aQQnIiJrmBnTB/+YMQQAsHJbDhdcpE5jESKHcfZ8Ix766CBMEnDjqBDMjg0RHYmIrOjeCeG4b2I4AOCJzw5jV06Z4ERkD1iEyCHoDEYs+jADlfV6DA32xIt/ihYdiYgEePKGwfjTyBAYTRIe+CADmUXnRUciG8ciRA7h2a+O4VDReXi5uuA/c+OgcVGIjkREAsjlMrx00whMiPRDY7MRd7/7C/LL60THIhvGIkR2b1NaIT5OK4RMBrx220iE+riJjkREAqmUcqy7Mw4j+nqhql6PeRvTUFbTJDoW2SgWIbJrh4rOY+mXRwEAf7suCpOi/AUnIiJb4K5WYuNdYxDm64Yz5xrx2KeHOa2eLopFiOzW+QY9HvgwA3qjCYlDAvHAZG6hQUS/8eulxoa7xkCllOOnE+X45nCJ6Ehkg1iEyC5JkoTHPjmM4vONCPN1w+okrhdERBeK8O+FB69p+UvSc18fQ3Uj9ySjjliEyC69u/80dmZroVLI8cbto+CpcREdiYhs1P2TwhHh746KOh1WbjsuOg7ZGBYhsjtZxdVYvrXlD7Onpg/G8BAvwYmIyJaplYr2JTU+PFCI9IJzghORLWERIrtSpzPgwY9axgVdNzQQ88eFiY5ERHYgPtwXt4zuCwB4assRNBtNghORrWARIrshSRKe/vwITlc2oI+XBi//eQRkMo4LIqLOWTJtCHzcVcjR1uLtPadExyEbwSJEduOT9DP4MvMsFHIZXrttJLzdVKIjEZEd6e2uat+P7NWUEyisbBCciGwBixDZhbzyOixrXS8o+boojA7zEZyIiOzRn0aGYFyEL5qaTXhk80HoDEbRkUgwFiGyeXqDCY9sOojGZiPGRfhi0aQI0ZGIyE7JZDKsuHEEPDVKHCw8j+e+PiY6EgnGIkQ275UdOcgqroG3mwtW3xLL9YKIqFv6+brhtdtGQiYDPjpQiM2/FIqORAKxCJFN23uyAv/5KR8A8NJNIxDkpRGciIgcweRBAfjbdVEAgGe+OMpd6p0YixDZrKp6PZL/lwkAuD2+H6YOCxIbiIgcygOTB+L6oYHQG01Y9EE6Kup0oiORACxCZJMkScITnx1GWa0OEf7ueGbGUNGRiMjByOUyvHJLDCL83VFS3YTFH2ZwfSEnxCJENumjtEJ8f0wLF4UMr946Eq4qhehIROSAPDQu+M/c0eilVuLAqSqs+I5bcDgbFiGyOdklNXjhm5aZHI9P5RYaRNSzBgb0wiu3xAAANuw9hW1Z3KXembAIkU2paWrGog/S0dRswqQof9xz9QDRkYjICUwdFoT7JoYDAP7+yWGcrqgXnIishUWIbIYkSXjsf4dwurIBId6uWJPEqfJEZD1/nzoIo/v3Rq3OgAc+zEBTMxdbdAYsQmQz/vNTPnYc00KlkOPfd4xCb3duoUFE1uOikOP120fCx12FYyU1eO7ro6IjkRWwCJFNSM2rxMptLYMUl80aiphQb7GBiMgpBXu54tVbYyGTAR+nFWFLxhnRkaiH2UQRWrt2LcLCwqDRaBAfH4+0tLRLHvvuu+9CJpN1uGk0XGTPnmlrmvDQxxkwScBNo/ri9rH9REciIic2IdIfD18bCQB4+vMsnNDWCk5EPUl4Edq8eTOSk5OxbNkyZGRkICYmBlOnTkVZWdklH+Pp6YmSkpL2W0FBgRUTkyUZjCY8+FEGKur0GBzkgX/OGQ6ZjOOCiEish6dE4uqBfmhsNuKpLUcgSZLoSNRDhBeh1atXY+HChViwYAGGDh2KdevWwc3NDRs3brzkY2QyGYKCgtpvgYGBVkxMlrR2Vx5+OX0OHmol1t0Zx/WCiMgmKOQyrLo5BhoXOX4tOIftR0tFR6IeIrQI6fV6pKenIzExsf0+uVyOxMREpKamXvJxdXV16N+/P0JDQzF79mwcPXr5AW06nQ41NTUdbiTeoaLzeO2HkwCAF+YMR5ifu+BERES/CfLS4L4JLVPqV3x3HHoDV512REKLUEVFBYxG4wVndAIDA1FaevH2PWjQIGzcuBFffvklPvjgA5hMJowbNw5nzlx6QNvy5cvh5eXVfgsNDbXo+yDzNegNeHRzJowmCf83IhizY/uIjkREdIH7J0XA30ON05UNeP9nDsNwRMIvjZkrISEB8+bNQ2xsLCZNmoQtW7bA398f//nPfy75mCVLlqC6urr9VlRUZMXEdDEvbs1GfkU9gjw1+NecaI4LIiKb5K5Wtu9S/1rKSVQ3NAtORJYmtAj5+flBoVBAq9V2uF+r1SIoqHM7jbu4uGDkyJHIzc295DFqtRqenp4dbiTOruNl+ODnQgDAK7fEwMvNRXAiIqJLu3l0KAYFeqC6sRmvt17OJ8chtAipVCrExcUhJSWl/T6TyYSUlBQkJCR06jmMRiOOHDmC4ODgnopJFlRZp8PfPz0MALh7/ACMH+gnOBER0eUp5DI8NWMIAOC91NMoqOT2G45E+KWx5ORkrF+/Hu+99x6ys7OxaNEi1NfXY8GCBQCAefPmYcmSJe3HP//889ixYwfy8/ORkZGBO++8EwUFBbj33ntFvQXqJEmS8NTnR1BRp0NUYC88fsMg0ZGIiDplUpQ/JkT6odkoYeW2HNFxyIKUogMkJSWhvLwcS5cuRWlpKWJjY7Ft27b2AdSFhYWQy3/ra+fOncPChQtRWlqK3r17Iy4uDvv378fQoUNFvQXqpM8yirH9qBYuChn+X1IsNC6cKk9E9uPpGUMw/dU9+PZICe4uqEJcfx/RkcgCZJITrhJVU1MDLy8vVFdXc7yQlZw514Bpa/agVmfA4zcMwgOTB4qORERktic+PYzNvxYhwt8dXz90NdxUws8nOJWe+P4WfmmMHJ/JJOHvnxxGrc6AuP69cf/ECNGRiIi65IlpgxHoqUZeeT2WfclNWR0BixD1uPdSTyM1vxKuLgq8cnMMFHJOlSci++TjrsKapJGQy4BP0s/gi4PFoiNRN7EIUY/KLavDiu9adpV/asYQrh5NRHYvIcIXD7VvynoEpyo4i8yesQhRjzEYTfjb/zKhM5gwMcofd8ZzV3kicgwPXTsQYwf4oF5vxEMfZ0BnMIqORF3EIkQ95t+783DoTDU8NUqsvGkEV48mIoehVMjx6q2x6O3mgqzimvYz32R/WISoR5zU1uK1lN82VA3y0ghORERkWcFerlh1cwwA4J19p7HzmPYKjyBbxCJEFidJEpZ+eRQGk4TEIYGYFcMNVYnIMU0ZEoi7xw8AADzx2WFU1OkEJyJzsQiRxX1zuASp+ZVQK+VYNnMoL4kRkUN7YtogDA7yQGW9Hk9+dgROuDyfXWMRIouq1xnwr2+zAQAPTB6IUB83wYmIiHqWWqnA6lti4aKQYWe2Fp+knxEdiczAIkQW9doPJ1Fa04R+Pm64f1K46DhERFYxtI8nkq9r2T/x+a+PoaiqQXAi6iwWIbKY3LI6bNhzCgCwbOZQ7iVGRE7lvonhGBPWG3U6A/72ySEYTbxEZg9YhMgiJEnCs1+1DJCeMjgAU4YEio5ERGRVCrkMr9wcC3eVAmmnqrBhb77oSNQJLEJkEd9llWJvbgVUSjmWzRwmOg4RkRD9fN3wzP8NBQCs2n4Cx0trBCeiK2ERom5rajbin98cAwD8ZVIE+vlygDQROa+kMaFIHBIAvdGEpz/P4iwyG8ciRN32wc8FOFvdhBBvVzwwmTvLE5Fzk8lk+NefouHqokB6wTl8l1UqOhJdBosQdUu9zoA3d+cBAB6ZEskB0kREAAI9NbhvYsvM2RXfHedeZDaMRYi65d39p1FZr0eYrxtuHBUiOg4Rkc24f1I4AjzUKKxqwPupBaLj0CWwCFGX1TQ1462fWmZF/DUxCkoF/3MiImrjplLib9dHAQBeSzmJc/V6wYnoYvjNRV22Yc8pVDc2IzKgF2ZyPzEiogv8OS4Ug4M8UNNkwGs/nBQdhy6CRYi65Fy9Hhv2tiye+Oh1UVDIuZ8YEdEfKeQyPD1jCADg/dQCnKqoF5yI/ohFiLrkrT35qNMZMCTYEzcMCxIdh4jIZk2I9MfkQf4wmCSs+C5bdBz6AxYhMlt5rQ7v7jsNAPjbdVGQ82wQEdFlPTV9COQyYPtRLdJOVYmOQ7/DIkRme3N3HhqbjYgJ9caUIQGi4xAR2byoQA8kjekHAFj+XTYXWbQhLEJkltMV9fjgQMs00L9dFwWZjGeDiIg649HESGhc5DhYeB4/HC8THYdasQhRp5lMEh7/7DD0BhMmRPphQqSf6EhERHYjwFOD+ePCAACrdpyAibvT2wQWIeq0938uQNqpKripFHjxT9E8G0REZKa/TIyAh1qJ7JIabM0qER2HwCJEnVRU1YCXth0HADw5bTBCfbixKhGRuXq7q3DPhAEAgNXfn4DBaBKciFiE6IokScKTWw6jQW/E2AE+uDO+v+hIRER2656rB6C3mwvyy+vx+cFi0XGcHosQXdGmX4qwL7cSGhc5Vt40gtPliYi6wUPjgr9MigAArNl5khuyCsYiRJd19nwj/vVtywJgj10/CGF+7oITERHZv3kJYfD3UKP4fCM2/1IkOo5TYxGiS5IkCU99fgR1OgNG9vPGgvEDREciInIIrioFHr52IADg9R9y0ajnWSFRWIToknbnlGN3TjlUCjle/vMI7idGRGRBSWP6oW9vV5TX6vDf1NOi4zgtpTkH19TUmP0Cnp6eZj+GxDOZJLy8PQcAcNf4MAwM8BCciIjIsaiUcjwyJRJ///Qw3vwxD7fF94OnxkV0LKdjVhHy9vY2a+0YmUyGEydOIDw83OxgJNbWrBIcK6lBL7USi1oH9RERkWX9aWQI1v2Yh7zyerz9Uz6Srx8kOpLTMasIAcCnn34KHx+fKx4nSRKmT5/epVAklsFowuodJwAACyeEo7e7SnAiIiLHpFTI8dj1g7Dowwy8vfcU5o0Lg18vtehYTsWsItS/f39MnDgRvr6+nTo+PDwcLi48zWdvPss4g/yKevj8buEvIiLqGTcMD0J0iBeOFFfj37vysHTmUNGRnIpZg6VPnTrV6RIEAFlZWQgNDTU7FInT1GzEqztPAgAemByBXmqzTxoSEZEZZDIZ/j615ZLYBz8XoPh8o+BEzoWzxqiDDw8U4mx1E4K9NLjzKq4gTURkDRMi/RA/wAd6owmvtf5llKyjy3/dP3XqFPbs2YOCggI0NDTA398fI0eOREJCAjQajSUzkpXU6Qz4965cAMDDUyKhcVEITkRE5BxkMhkev2EQbnozFZ9mnMF9k8IR4d9LdCynYHYR+vDDD/Hqq6/i119/RWBgIPr06QNXV1dUVVUhLy8PGo0Gd9xxB5544gn0788zCvbknb2nUFmvxwA/d/w5rq/oOERETiWuvw+mDA5AyvEyrP7+BNbePkp0JKdg1qWxkSNH4rXXXsNdd92FgoIClJSUID09HXv37sWxY8dQU1ODL7/8EiaTCaNHj8Ynn3zSU7nJwqobmvHWT/kAgEevi4KLgldNiYis7bHWsULfHi5BVnG14DTOwaxvuxUrVuDAgQN44IEHLjoIWq1WY/LkyVi3bh2OHz/e6fWD1q5di7CwMGg0GsTHxyMtLa1Tj9u0aRNkMhnmzJljztugi3gv9TRqdQYMCvTA/0UHi45DROSUhgR7YlZMHwDA//v+hOA0zsGsIjR16tQrHmMwGFBXVwdfX1/ExcVd8fjNmzcjOTkZy5YtQ0ZGBmJiYjB16lSUlZVd9nGnT5/GY489hgkTJnQ6P11cvc6AjftOAQAWXzuQu8sTEQn06HVRkMuAlONlPCtkBV2+/lFYWHjR25dffokhQ4Z0+nlWr16NhQsXYsGCBRg6dCjWrVsHNzc3bNy48ZKPMRqNuOOOO/Dcc89x1WoL+DitEOcbmhHm64YZPBtERCTUAD/39rNCr6VwBllP6/KssbCwsEtutzFy5MhOPYder0d6ejqWLFnSfp9cLkdiYiJSU1Mv+bjnn38eAQEBuOeee7Bnzx7zglMHOoMR6/e0jA36y6QIbqxKRGQDHrx2IL48dBY7jmmRXVKDIcHct7OndLkIHTx4sMPPRqMReXl5eP755/Hwww936jkqKipgNBoRGBjY4f7AwEAcP378oo/Zu3cvNmzYgMzMzE5n1el00Ol07T93ZfNYR/VZejG0NToEe2lw4yjOFCMisgUDAzwwIzoY3xwuwRs/5GLtHZxB1lO6XIRiYmIuuG/UqFHo27cvbrzxRsybN69bwS6mtrYWc+fOxfr16+Hn59fpxy1fvhzPPfecxfPYO4PRhHU/5gFo2VNMpeRMMSIiW/HQtZH45nAJtmaV4IS2FlGBHqIjOSSLf/P5+Ph0OPtyOX5+flAoFNBqtR3u12q1CAoKuuD4vLw8nD59GjNnzoRSqYRSqcR///tffPXVV1AqlcjLy7vo6yxZsgTV1dXtt6KiIvPfmAP69kgJCqsa4OOuwq1juRUKEZEtGRTkgWnDgyBJwBs/5IqO47C6XIRqamouuJ08eRJPPfUUHnnkkQ73X4pKpUJcXBxSUlLa7zOZTEhJSUFCQsIFxw8ePBhHjhxBZmZm+23WrFm45pprkJmZecl9zdRqNTw9PTvcnJ3JJOHfu1qK493jw+Cm4p5iRES25sFrBwIAvj58FrlldYLTOKYuf/t5e3tfdLC0JEn4/PPP8fzzz0OSJMhkMhiNxks+T3JyMubPn4/Ro0dj7NixWLNmDerr67FgwQIAwLx58xASEoLly5dDo9Fg+PDhF+QAcMH9dHk7s7XI0dbCQ63E3IQw0XGIiOgihvXxQuKQQOzM1uLfu3KxOilWdCSH0+UitGvXLosESEpKQnl5OZYuXYrS0lLExsZi27Zt7QOoCwsLIZdz7IolSZKEtbtbzgbNTegPL1cXwYmIiOhSHp4yEDuztfgisxgPT4lEmJ+76EgORSZJkiQ6hLXV1NTAy8sL1dXVTnmZ7NvDJVj8UQbUSjn2PXkt/HqpRUciIqLLWPBOGnbllCNpdChe+vMI0XGE6Ynvb55qcTJ1OgOe/+YogJZ1g1iCiIhsX9tYoc8zi1Fe27kJSdQ5Zl0aGzBgwCUXUbycv/71r51eW4h61prvT0Bbo0N/XzcsmhwhOg4REXVCXH8fxIZ6I7PoPD74uQCPXhclOpLDMKsIvfvuu116kbCwsC49jizreGkN3tl/GgDw7Kxh0LgoxAYiIqJOu+fqAXjo44P44OcCLJocwT/DLcSsIjRp0qSeykE9zGSS8I/Ps2A0SZg2PAjXDAoQHYmIiMwwbXgQQrxdUXy+EV9lnsUtY7j+myV0eYxQeXn5JX/tyJEjXX1a6iGfZZzBrwXn4KZS4Jn/Gyo6DhERmUmpkGP+uP4AgLf35sMJ5zr1iC4XoejoaHz77bcX3L9q1SqMHTu2W6HIss436LH8u5a92x6ZEok+3q6CExERUVckjekHN5UCJ7R12JtbITqOQ+hyEUpOTsZNN92ERYsWobGxEcXFxZgyZQpWrlyJjz76yJIZqZtWbs9BVb0ekQG9cPfVA0THISKiLvJydcEto1suiW3Ye0pwGsfQ5SL0+OOPIzU1FXv27MGIESMwYsQIqNVqHD58GH/6058smZG64YS2Fh+nFQIA/jlnOFwUXDGBiMieLRgfBpkM2J1TjtyyWtFx7F63vhUHDhyI4cOH4/Tp06ipqUFSUtJFN0slcdb9mAdJAqYOC0R8uK/oOERE1E39fd1x3ZCW3Rc27D0tNowD6HIR2rdvH0aMGIGTJ0/i8OHDePPNN/HQQw8hKSkJ586ds2RG6qK2mQUAsGjyQMFpiIjIUu5pHeawJeMMqur1gtPYty4XoWuvvRZJSUn4+eefMWTIENx77704ePAgCgsLER0dbcmM1EVv78mHwSQhIdwXsaHeouMQEZGFjB3gg+gQL+gMJryfWiA6jl3rchHasWMHVqxYAReX3zbsjIiIwL59+3D//fdbJBx13bl6PTalFQEAV5AmInIwMpkM905oOSv01k95KKtpEpzIfnW5CF1qcUW5XI5nnnmmy4HIMt5LPY3GZiOG9fHEhEg/0XGIiMjCZo7og9hQb9Trje1LpJD5zCpCmzZt6vSxRUVF2Ldvn9mBqPsa9Aa817qVxl8mRXRpfzgiIrJtcrkMz88eBpkM+PxgMdJOVYmOZJfMKkJvvvkmhgwZgpUrVyI7O/uCX6+ursbWrVtx++23Y9SoUaisrLRYUOq8zb8U4VxDM/r5uGHacM7iIyJyVCP6euPWMf0AAEu/zILBaBKcyP6YVYR+/PFHvPTSS/j+++8xfPhweHp6IjIyEtHR0ejbty98fX1x9913o1+/fsjKysKsWbN6KjddQrPRhLf3tCyydd/EcCi5bhARkUP7+9RB8HJ1wfHSWnzUum4cdZ5Zm64CwKxZszBr1ixUVFRg7969KCgoQGNjI/z8/DBy5EiMHDkScjm/fEX5+tBZFJ9vhF8vFf4c11d0HCIi6mE+7io8NnUQnvkiC6u252BGdDB8e6lFx7IbZhehNn5+fpgzZ44Fo1B3SZKE//yYDwBYMH4ANC4KwYmIiMgabh/bDx8fKMSxkhq8vD0HK24aITqS3ej2qRu9Xo8zZ86gsLCww42sL+1UFXK0tXBXKXDnVf1FxyEiIitRtA6cBoDNvxYhs+i82EB2pMtF6OTJk5gwYQJcXV3Rv39/DBgwAAMGDEBYWBgGDODGniJ8eahlFenp0cHwcnW5wtFERORIRof54MZRIZAkYNX2HNFx7EaXL43dddddUCqV+OabbxAcHMwp2oLpDSZsPVICAJgdGyI4DRERiZB8XRQ+P1iMvbkVOF1RjzA/d9GRbF6Xi1BmZibS09MxePBgS+ahLvrpRDnONzTD30ONhAhurkpE5Iz69nbDpCh/7M4px8e/FGLJtCGiI9m8Ll8aGzp0KCoqKiyZhbqh7bLYzBF9oJDz7BwRkbO6fWzLukKf/noGegPXFbqSLhehl156CY8//jh2796NyspK1NTUdLiR9dTrDPj+WCkAYHZsH8FpiIhIpGsHByDIU4PKej22Hy0VHcfmdfnSWGJiIgBgypQpHe6XJAkymQxGo7F7yajTvj+mRVOzCWG+bhjR10t0HCIiEkipkOOWMaF4LeUkPjpQiJkx/Avy5XS5CO3atcuSOagbvswsBtAySJqD1omIKGlMKN744SRS8yuRX16HcP9eoiPZrC4VoebmZjz//PNYt24dIiMjLZ2JzFBZp8NPJ1vGas3iZTEiIgIQ4u2KyYMC8MPxMnycVoinZwwVHclmdWmMkIuLCw4fPmzpLNQFW4+UwGiSEB3ihQg2fiIiatU+aDr9DHQGDle5lC4Plr7zzjuxYcMGS2ahLvgys2W2GAdJExHR700e5I9gLw3ONTRjWxYHTV9Kl8cIGQwGbNy4ETt37kRcXBzc3Tsu2rR69epuh6PLK6pqwK8F5yCTgYPhiIioA6VCjqQxoVizs2XQNBfbvbguF6GsrCyMGjUKAHDixIkOv8YBu9bx9eGWs0EJ4b4I9NQITkNERLYmqXX22IFTVcgtq8PAAA6h+CPOGrNjXx7kZTEiIrq0YC9XXDs4EDuztfjoQCGWzuSg6T/q9u7zubm52L59OxobGwG0rCNEPe9g4TnkaGuhUshxw/Bg0XGIiMhG3XFVy6Dp//1ahJqmZsFpbE+Xi1BlZSWmTJmCqKgoTJ8+HSUlLRt+3nPPPfjb3/5msYB0ca+lnATQMmWeO80TEdGlTIr0R2RAL9TpDPj4QKHoODany0Xo0UcfhYuLCwoLC+Hm5tZ+f1JSErZt22aRcHRxh8+cx66ccshlwIPXDBQdh4iIbJhcLsPCieEAgHf2neb+Y3/Q5SK0Y8cOvPTSS+jbt2+H+yMjI1FQUNDtYHRpr6XkAgDmxIYgzM/9CkcTEZGzmx3bBwEeapTWNOGr1k26qUWXi1B9fX2HM0FtqqqqoFaruxWKLi2ruBo7s7WQy4DF1/JsEBERXZlaqcCC8QMAAOt/yud43t/pchGaMGEC/vvf/7b/LJPJYDKZsHLlSlxzzTUWCUcXev2HlrFBM2P6cCVpIiLqtNvj+8FdpUCOthY/nigXHcdmdHn6/MqVKzFlyhT8+uuv0Ov1ePzxx3H06FFUVVVh3759lsxIrbJLarD9qBYyjg0iIiIzebm64Nax/bBh7ym89VM+Jg8KEB3JJnT5jNDw4cNx4sQJXH311Zg9ezbq6+tx44034uDBg4iIiLBkRmr1xg8tY4OmRwcjMtBDcBoiIrI3d189AAq5DPvzKpFVXC06jk3o8hkhAPDy8sLTTz9tqSx0GSe0tdia1bJEwcPXRgpOQ0RE9ijE2xUzRwTji8yz+M9P+Xj9tpGiIwnX7QUVgZaB0xs3bsTatWtx8uRJSzwl/cHrP+RCkoBpw4MwKIhng4iIqGvum9hy1WbrkRIUVTUITiOe2UWosLAQkyZNgoeHB6677joUFhZi1KhRuPfee/HQQw8hNjYWP/30k1nPuXbtWoSFhUGj0SA+Ph5paWmXPHbLli0YPXo0vL294e7ujtjYWLz//vvmvg27UljZgG9a9xV7kDPFiIioG4b28cSESD8YTRI27jslOo5wZhehxx57DHq9HuvWrYObmxumTp2KyMhIlJSUQKvVYtq0aXj22Wc7/XybN29GcnIyli1bhoyMDMTExGDq1KkoKyu76PE+Pj54+umnkZqaisOHD2PBggVYsGABtm/fbu5bsRtfZhZDkoDxA30xrI+X6DhERGTn7rm6ZSr9Z+ln0NRsFJxGLJlk5mICQUFB+OqrrzB27FhUVVXBz88P+/btQ0JCAgDg0KFDmDJlCioqKjr1fPHx8RgzZgzeeOMNAIDJZEJoaCgeeughPPnkk516jlGjRmHGjBl44YUXOnV8TU0NvLy8UF1dDU9Pz049RhRJkpC4+kfkldfj5T+PwM2jQ0VHIiIiO2cySZiwcheKzzfi1VtjMTs2RHSkTumJ72+zzwiVlZWhf//+AFrOzri5uSEwMLD914OCgnDu3LlOPZder0d6ejoSExN/CySXIzExEampqVd8vCRJSElJQU5ODiZOnHjJ43Q6HWpqajrc7MWxkhrklddDpZRj6vAg0XGIiMgByOUy/DmuZWeIzb8UCU4jVpcGS8tksov+u7kqKipgNBo7FCkACAwMRGlp6SUfV11djV69ekGlUmHGjBl4/fXXcd11113y+OXLl8PLy6v9FhpqP2dVvsxsGRuUOCQAnhpurkpERJZx8+i+kMmA/XmVKKx03kHTXZo+v3Tp0vbtNfR6Pf71r3/By6tl7EpDQ8//Znp4eCAzMxN1dXVISUlBcnIywsPDMXny5Isev2TJEiQnJ7f/XFNTYxdlyGSS8FVrEZoVYx+nLYmIyD707e2Gqwf6Yc/JCvzv1yI8NnWQ6EhCmF2EJk6ciJycnPafx40bh/z8/AuO6Qw/Pz8oFApotdoO92u1WgQFXfoykFwux8CBLbOnYmNjkZ2djeXLl1+yCKnVarvc/yztdBVKa5rgoVFi8iB/0XGIiMjBJI0JxZ6TFfg0/QwevS4KCnnXr/LYK7OL0O7duy324iqVCnFxcUhJScGcOXMAtAyWTklJwYMPPtjp5zGZTNDpdBbLZSu+zCwGAEwfHgyNi0JwGiIicjTXDQ1EbzcXlNY04acT5bhmsPNtu2GRBRW7Izk5GevXr8d7772H7OxsLFq0CPX19ViwYAEAYN68eViyZEn78cuXL8f333+P/Px8ZGdn45VXXsH777+PO++8U9Rb6BE6gxFbj7SMk5od20dwGiIickRqpQJzRrYMvXDWQdNmFaHk5GTU19d3+vglS5agqqrqssckJSVh1apVWLp0KWJjY5GZmYlt27a1D6AuLCxESUlJ+/H19fV44IEHMGzYMIwfPx6fffYZPvjgA9x7773mvBWb99OJClQ3NiPAQ434cF/RcYiIyEEljWkZM7szW4uKOse7unIlZq0jpFAoUFpaCn//zo1X8fT0RGZmJsLDw7scsCfYwzpCiz/KwLeHS3Dv1QPwj/8bKjoOERE5sNlv7MWhM9V4evoQLJxoW9/Zv9cT399mjRGSJAlRUVGdnjJvztkj+k2dzoCdx1oGkNvLIldERGS/bhkTikNnqrH51yLcO2FAt5bGsTdmFaF33nnH7Bf44xpBdGU7jpZCZzAh3M8dw0Ns84wVERE5jpkxffDCN8eQW1aHjMLziOvfW3QkqzGrCM2fP7/TqzLb6iUne/BF69pBs2NDnKqVExGRGJ4aF0yPDsaWjGL875ciFqHL8fb27tSXs9Ho3Ju4dVV5rQ77clv2aZvF2WJERGQlSaNDsSWjGN8eKcGzs4bBVeUcy7aYXYR27drV/u+SJGH69Ol4++23ERLCsSyWsPmXQhhNEkb288YAP3fRcYiIyEmMCfNBqI8riqoaseNYqdOMUTW7CE2aNKnDzwqFAldddZXNzQyzRwajCR8eKAQAzL2qv+A0RETkTORyGW4c2RevppzEp+lnnKYICV9QkX6TcrwMJdVN8HFXYXp0sOg4RETkZG4c1VJ+9uVWoLS6SXAa62ARsiHvpxYAAG4ZHcotNYiIyOr6+7pjTFhvmCTgi9ZtnhydRYoQZzZ1X155HfbmVkAmA+6I7yc6DhEROambRvUFAHyWfgZmrLlst8weI3TjjTd2+LmpqQl/+ctf4O7ecWDvli1bupfMybSdDZoyOAChPm6C0xARkbOaPiIYy746ipNldThSXI0Rfb1FR+pRZhchLy+vDj872manIjToDfgs/QwAYG5CmNgwRETk1Dw1Lpg6LAhfHTqLz9LPsAj9UVdWl6bL++LgWdTqDAjzdcOEgX6i4xARkZO7cVQIvjp0Fl8dOounZwyFSum4Q4od953ZCUmS8N/U0wCAO6/qD7mc462IiEisCZH+CPBQ41xDM3bllImO06NYhAT7teAcjpfWQuMix81xoaLjEBERQSGX4U8jW6bStw3dcFQsQoK1DZKeHRMCLzcXwWmIiIha3BTXMntsV04Zqur1gtP0HBYhgSrqdPguqwQAMDeBK0kTEZHtiAr0QHSIF5qNEr4+dFZ0nB7DIiTQ9qOlaDZKiA7xwvAQrys/gIiIyIraVpr+LMNxL4+xCAm0/agWADAtOkhwEiIiogvNiukDhVyGw2eqkVdeJzpOj2AREqS6sRn7cysAADcMYxEiIiLb49tLjUlR/gCALw465pYbLEKC/HBcC4NJQlRgL4T79xIdh4iI6KLaZo99frDYIbfcYBESZFtWKQCeDSIiItuWOCQQvdRKnDnXiPSCc6LjWByLkAANegN+PFEOAJg6nEWIiIhsl6tKgRtav6u2OODlMRYhAX46UY6mZhNCfVwxNNhTdBwiIqLLurH18ti3h0ugMxgFp7EsFiEB2i6LTR0aBJmMW2oQEZFtiw/3RZCnBtWNzdh1vFx0HItiEbIyvcGElOMt+7bcwMtiRERkBxRyGWbH9gHgeLPHWISsLDW/ErVNBvh7qDGqX2/RcYiIiDrlT62LK/5wvAzVDc2C01gOi5CVtV0Wu35oIHeaJyIiuzE4yBODgzygN5rw7ZES0XEshkXIiowmCd8fa502z8tiRERkZ9rWFHKky2MsQlaUXnAOFXV6eGqUuCrcV3QcIiIis8yODYFMBqSdrkJRVYPoOBbBImRFbZfFEocEwkXB33oiIrIvQV4ajIto+Yv8l5mOcVaI38ZWIkkSth9tnTbPy2JERGSn5sS2XB775rBjjBNiEbKSo2drUHy+Ea4uCkyM9Bcdh4iIqEumDAmETAYcL61FaXWT6DjdxiJkJTtazwZNjPKDq0ohOA0REVHX+LirEBvqDQD48USZ2DAWwCJkJTuOaQEA1w/lZTEiIrJvk6MCAMAhVplmEbKCwsoGHC+thUIuw7WDA0THISIi6pZrBrcM8dibWwG9wSQ4TfewCFnBjta1g8aG+aC3u0pwGiIiou4Z3scLfr1UqNMZkF5wTnScbmERsoK2y2LXDQ0UnISIiKj75HIZJka1nBXanWPf44RYhHpYVb0ev56uAsAiREREjuOaQa3jhFiE6HJSsrUwScDQYE+E+riJjkNERGQREyL9IJcBJ7R1KD7fKDpOl7EI9TBeFiMiIkfk7abCqH69Adj35TEWoR7UqDdiz8mWqYXXD2MRIiIix3LNYPufRs8i1IP2nCxHU7MJId6uGBrsKToOERGRRU1qHTC9P68COoNRcJqusYkitHbtWoSFhUGj0SA+Ph5paWmXPHb9+vWYMGECevfujd69eyMxMfGyx4v0+8tiMplMcBoiIiLLGtbHEwEeajTojfjllH1OoxdehDZv3ozk5GQsW7YMGRkZiImJwdSpU1FWdvHrjbt378Ztt92GXbt2ITU1FaGhobj++utRXGxbu+AajCakZLeuJs3LYkRE5IBkMhkmD7LvafTCi9Dq1auxcOFCLFiwAEOHDsW6devg5uaGjRs3XvT4Dz/8EA888ABiY2MxePBgvP322zCZTEhJSbFy8stLLziHcw3N8HJ1wdgwH9FxiIiIesRkO59GL7QI6fV6pKenIzExsf0+uVyOxMREpKamduo5Ghoa0NzcDB+fS5cNnU6HmpqaDree1nZZbMrgACgVwvsmERFRj7g60g8KuQx55fUoqmoQHcdsQr+hKyoqYDQaERjY8dJRYGAgSktLO/UcTzzxBPr06dOhTP3R8uXL4eXl1X4LDQ3tVu4rkSSpfVsNXhYjIiJH5qlxwej+9juN3q5PVaxYsQKbNm3C559/Do1Gc8njlixZgurq6vZbUVFRj+Y6VlKDoqpGqJXy9iXIiYiIHFXb5bGU4yxCZvHz84NCoYBWq+1wv1arRVBQ0GUfu2rVKqxYsQI7duzAiBEjLnusWq2Gp6dnh1tP+uJgy8DtyYP84aZS9uhrERERiZY4pKUI7c+tRG1Ts+A05hFahFQqFeLi4joMdG4b+JyQkHDJx61cuRIvvPACtm3bhtGjR1sjaqcZjCZ8kXkWAHDjqL6C0xAREfW8yEAPRPi7Q2804Qc7Oysk/NJYcnIy1q9fj/feew/Z2dlYtGgR6uvrsWDBAgDAvHnzsGTJkvbjX3rpJTzzzDPYuHEjwsLCUFpaitLSUtTV1Yl6Cx3sy6tEea0O3m4u7RvSERERObobhrdcyfnuSOfG+NoK4ddtkpKSUF5ejqVLl6K0tBSxsbHYtm1b+wDqwsJCyOW/9bU333wTer0ef/7znzs8z7Jly/Dss89aM/pFfZ5xBgAwK6YPVErhPZOIiMgqpg0Pxtpdedh9ogwNeoPdDA2RSZIkiQ5hbTU1NfDy8kJ1dbVFxwvV6QwY/c/v0dRswucPjMPI1s3oiIiIHJ0kSZiwchfOnGvEm3eMwrToYIu/Rk98f/OUhQV9d6QETc0mhPu5IzbUW3QcIiIiq5HJZJjWdnksy34uj7EIWdDnrbPFbhwVwr3FiIjI6dwwvOUs0A/Hy+xmE1YWIQs5e74RqfmVAIDZsSGC0xAREVnfyFBvBHqqUaczYO/JCtFxOoVFyEK+yCyGJAHxA3wQ6uMmOg4REZHVyeUy3DDMvi6PsQhZgCRJ2JLRclnsJq4dRERETqzt8tj3x7RoNpoEp7kyFiELyCquQW5ZHdRKOaZFX35FbCIiIkc2doAPfN1VqG5sxs+tQ0ZsGYuQBXzWunbQ9cOC4KFxEZyGiIhIHIVc1r7huD1cHmMR6qZmowlfH2rbUoODpImIiNouj+04WgqjybaXK2QR6qa0U1WorNfDr5cKEwb6iY5DREQkXEK4Lzw1SlTU6fHr6SrRcS6LRaib9uW2TA+cGOUPpYK/nURERCqlHIlD7ePyGL+5u2lfXstAsHERPBtERETUZtrvZo/Z8m5eLELdUNPUjCNnzgMAxkX4ig1DRERkQ8YP9IVKKUfx+UacLKsTHeeSWIS64UB+FUwSMMDPHX28XUXHISIishluKiUSwltOEvxwvExwmktjEeqGtvFBPBtERER0oWsHBwAAdrEIOabU1vFB4zlbjIiI6ALXDGopQr8WnEN1Y7PgNBfHItRF5bU65GhrAQBXhfOMEBER0R/183VDhL87jCYJe06Wi45zUSxCXdS20/zQYE/4uKsEpyEiIrJNbZfHbHWcEItQF+3n+CAiIqIruqa1CP2YUw6TDa4yzSLURfs5PoiIiOiKxoT5wEOtRGW9Hodal5yxJSxCXVBU1YDCqgYo5TKMGeAjOg4REZHNclHIMSGq5aSBLc4eYxHqgv15LZfFYkK90UutFJyGiIjItrXNHvshh0XIIbRfFuP4ICIioiua3FqEsoprUFbTJDhNRyxCZpIkqb0IJXB/MSIioivy91Ajpq8XAGCXjZ0VYhEyU25ZHcprdVAr5RjV31t0HCIiIrtwjY1Oo2cRMlPbthpjwnygVioEpyEiIrIPbeOE9p6sgM5gFJzmNyxCZmq7LDZuIMcHERERdVZ0iBf8eqlRrzfil1PnRMdpxyJkBqNJws+tK0qP4/ggIiKiTpPLZZg8yB+AbV0eYxEyw9Gz1ahpMsBDo8TwPp6i4xAREdmVtu02dtvQgGkWITO0nQ2KH+ADpYK/dUREROa4OtIPSrkM+RX1KKxsEB0HAIuQWQ7kVwEA4gdwfBAREZG5PDUuGNW/NwBg9wnbOCvEItRJRpOEtNOtRSic22oQERF1RdvsMVvZboNFqJOyS2pQ22SAh1qJocEcH0RERNQVbQOmU/Mr0dQsfho9i1AntY0PGh3Wm+ODiIiIumhwkAeCPDVoajbhwKkq0XFYhDqr7cOKD+f4ICIioq6SyX6bRm8Ls8dYhDrBZJKQ1laEBnB8EBERUXf8VoTKBSdhEeqU46W1qG5shrtKgeEhXqLjEBER2bXxA1um0Z+qqEdBZb3QLCxCnXDgVMv4oLgwH7hwfBAREVG3eGhcENc2jV7wWSF+q3fCb+sH8bIYERGRJVxjI6tMswhdgel36wddxYHSREREFtE2Tmh/nthp9CxCV3CyrA5V9Xq4uigwoi/HBxEREVnCoMCWafQ6g6l9iRoRWISuoH18UP/eHB9ERERkITKZDNcMFj97jN/sV8DxQURERD1jUlTLOKEfTzh5EVq7di3CwsKg0WgQHx+PtLS0Sx579OhR3HTTTQgLC4NMJsOaNWt6LJckSe1nhK6K4PggIiIiSxo/0Ld9Gv3pCjHT6IUXoc2bNyM5ORnLli1DRkYGYmJiMHXqVJSVXXwUeUNDA8LDw7FixQoEBQX1aLa88jpU1OmhVso5PoiIiMjCPDQuGBPWcsVF1Owx4UVo9erVWLhwIRYsWIChQ4di3bp1cHNzw8aNGy96/JgxY/Dyyy/j1ltvhVqt7tFsP7deFhvVrzfUSkWPvhYREZEzal9lWtDlMaFFSK/XIz09HYmJie33yeVyJCYmIjU11WKvo9PpUFNT0+HWGb/tL8bxQURERD1h/EA/AMDBwvOQJMnqry+0CFVUVMBoNCIwMLDD/YGBgSgtLbXY6yxfvhxeXl7tt9DQ0Cs+RpKk9ul8XD+IiIioZwwK8oBKKUd1YzMKKhus/vrCL41Zw5IlS1BdXd1+KyoquuJjTlXUo7xWB5VSjthQ754PSURE5IRcFHIMDfYEABw6c97qry+0CPn5+UGhUECr1Xa4X6vVWnQgtFqthqenZ4fblbSND4oN9YbGheODiIiIekrbhKTDZ6qt/tpCi5BKpUJcXBxSUlLa7zOZTEhJSUFCQoLAZMD+vAoAwDhOmyciIupRI/p6AwCOCChCSqu/4h8kJydj/vz5GD16NMaOHYs1a9agvr4eCxYsAADMmzcPISEhWL58OYCWAdbHjh1r//fi4mJkZmaiV69eGDhwoEUytYwPajkjlMDxQURERD0qpvWMUNbZahhNEhRymdVeW3gRSkpKQnl5OZYuXYrS0lLExsZi27Zt7QOoCwsLIZf/duLq7NmzGDlyZPvPq1atwqpVqzBp0iTs3r3bIplyy+pQUaeDWilHbD9vizwnERERXVy4fy+4qxSo1xuRW1aHQUEeVntt4UUIAB588EE8+OCDF/21P5absLCwHp9el9o6W2x0GNcPIiIi6mkKuQzDQryQdqoKh86ct2oRcopZY+ZKzWspQrwsRkREZB0x7QOmz1v1dVmE/sBk+m39oAQOlCYiIrIKUQOmWYT+IEdbi3MNzXBTKdo/FCIiIupZMa3fudkltdAbTFZ7XRahP2i7LDY6zAcuCv72EBERWUOojyu83VygN5pwvLRzW2FZAr/p/6BtoDTHBxEREVmPTCZDdEjLOKFDVrw8xiL0O0aThAMcH0RERCRETPs4ofNWe00Wod/JLqlBTZMBvdRKDO9z5W04iIiIyHJEbLXBIvQ7beODxg7wgZLjg4iIiKyqbZLSCW0tGvQGq7wmv+1/h+ODiIiIxAny0iDAQw2TBBw9a50B0yxCrQxGE9JOte4vxvFBREREQrSdFTpUdN4qr8ci1CrrbA3qdAZ4apQYEszxQURERCK0rTB9pNg644RYhFq1jQ+KD/e16q63RERE9JtoKw+YZhFqxfFBRERE4rVdGjtVUY/qxuYefz0WIQDNRhN+Pc3xQURERKL5uKsQ6uMKwDr7jrEIoeU6ZIPeiN5uLhgU6CE6DhERkVNrHzBthYUVWYQA5JTWAgCi+3pDzvFBREREQo1o3WqDZ4SsJLesDgAw0L+X4CRERETUNmD6aAmLkFW0F6EAFiEiIiLR2oapFFU19vgK0yxC+K0IRQayCBEREYnm20sNv14qAMBJbV2PvpbTF6F6nQHF5xsB8NIYERGRrYgMaDkrdEJb26Ov4/RFKK+8pWn69VKht7tKcBoiIiICgEFBLEJW0XZZLIJng4iIiGxGVGBbEeKlsR7FgdJERES2J6p13C7PCPWwk20DpVmEiIiIbEZk6xmhkuqmHt1qw+mLUF77GSGuKE1ERGQrvFxdEOSpAQDklvXcWSGnLkI6gxGnK+sBcOo8ERGRrYlqHTCdU9pz44ScuggVVjbAJAEeaiUCPNSi4xAREdHvRAX0/Dghpy5C+eUtZ4MiAnpBJuMeY0RERLYkygpT6J26CLWtIcSB0kRERLbHGlPonboI5Zdz6jwREZGtajtRUVGnQ1W9vkdew7mLUAUHShMREdkqd7USfXu7Aui5y2NOXYROVTYAAAb6c+o8ERGRLRoU2LPjhJy6CDUbTFAr5QhpbZtERERkWyJZhHpWhH8vKOScMUZERGSLBgW1TqHvobWEnL4IcaA0ERGR7WqfOVZWC0mSLP78Tl+EOHWeiIjIdkX494JcBpxvaEZFrc7iz+/0RYhnhIiIiGyXxkWBMF93AEBu60LIlsQixCJERERk09qWuTnZA5uvOnURUspl6N/aMomIiMg2tU2hzyuz/IBppy5CoT6uUCmd+reAiIjI5rVNoT/JImRZEf68LEZERGTrBrVuvtq2R6glOXURCvdjESIiIrJ1Yb7uUMplqNcZLf7cNlGE1q5di7CwMGg0GsTHxyMtLe2yx3/yyScYPHgwNBoNoqOjsXXr1i69brg/xwcRERHZOpVS3mPf2cKL0ObNm5GcnIxly5YhIyMDMTExmDp1KsrKyi56/P79+3HbbbfhnnvuwcGDBzFnzhzMmTMHWVlZZr82ixAREZF9aBsnZGkyqSeWaTRDfHw8xowZgzfeeAMAYDKZEBoaioceeghPPvnkBccnJSWhvr4e33zzTft9V111FWJjY7Fu3bpOvWZNTQ28vLxQWl6FQL/elnkjRERE1GNeSzmJVd9komjNLaiuroanp6dFnldpkWfpIr1ej/T0dCxZsqT9PrlcjsTERKSmpl70MampqUhOTu5w39SpU/HFF19c8nV0Oh10ut9Wo6yurgYANDfVo6ZG0Y13QERERNYQ4g6YdA0AYNGtNoQWoYqKChiNRgQGBna4PzAwEMePH7/oY0pLSy96fGlp6SVfZ/ny5XjuuecuuD80NLQLqYmIiEikyspKeHl5WeS5hBYha1myZEmHs0jnz59H//79UVhYaLHfSOqampoahIaGoqioyGKnOalr+FnYDn4WtoWfh+2orq5Gv3794OPjY7HnFFqE/Pz8oFAooNVqO9yv1WoRFBR00ccEBQWZdTwAqNVqqNXqC+738vLif9Q2wtPTk5+FjeBnYTv4WdgWfh62Qy633FwvobPGVCoV4uLikJKS0n6fyWRCSkoKEhISLvqYhISEDscDwPfff3/J44mIiIguRfilseTkZMyfPx+jR4/G2LFjsWbNGtTX12PBggUAgHnz5iEkJATLly8HADzyyCOYNGkSXnnlFcyYMQObNm3Cr7/+irfeekvk2yAiIiI7JLwIJSUloby8HEuXLkVpaSliY2Oxbdu29gHRhYWFHU6BjRs3Dh999BH+8Y9/4KmnnkJkZCS++OILDB8+vNOvqVarsWzZsoteLiPr4mdhO/hZ2A5+FraFn4ft6InPQvg6QkRERESiCF9ZmoiIiEgUFiEiIiJyWixCRERE5LRYhIiIiMhpOWwRWrt2LcLCwqDRaBAfH4+0tLTLHv/JJ59g8ODB0Gg0iI6OxtatW62U1PGZ81msX78eEyZMQO/evdG7d28kJiZe8bOjzjP3/4s2mzZtgkwmw5w5c3o2oBMx97M4f/48Fi9ejODgYKjVakRFRfHPKQsx97NYs2YNBg0aBFdXV4SGhuLRRx9FU1OTldI6rp9++gkzZ85Enz59IJPJLruHaJvdu3dj1KhRUKvVGDhwIN59913zX1hyQJs2bZJUKpW0ceNG6ejRo9LChQslb29vSavVXvT4ffv2SQqFQlq5cqV07Ngx6R//+Ifk4uIiHTlyxMrJHY+5n8Xtt98urV27Vjp48KCUnZ0t3XXXXZKXl5d05swZKyd3POZ+Fm1OnTolhYSESBMmTJBmz55tnbAOztzPQqfTSaNHj5amT58u7d27Vzp16pS0e/duKTMz08rJHY+5n8WHH34oqdVq6cMPP5ROnTolbd++XQoODpYeffRRKyd3PFu3bpWefvppacuWLRIA6fPPP7/s8fn5+ZKbm5uUnJwsHTt2THr99dclhUIhbdu2zazXdcgiNHbsWGnx4sXtPxuNRqlPnz7S8uXLL3r8LbfcIs2YMaPDffHx8dL999/fozmdgbmfxR8ZDAbJw8NDeu+993oqotPoymdhMBikcePGSW+//bY0f/58FiELMfezePPNN6Xw8HBJr9dbK6LTMPezWLx4sXTttdd2uC85OVkaP358j+Z0Np0pQo8//rg0bNiwDvclJSVJU6dONeu1HO7SmF6vR3p6OhITE9vvk8vlSExMRGpq6kUfk5qa2uF4AJg6deolj6fO6cpn8UcNDQ1obm626AZ7zqirn8Xzzz+PgIAA3HPPPdaI6RS68ll89dVXSEhIwOLFixEYGIjhw4fjxRdfhNFotFZsh9SVz2LcuHFIT09vv3yWn5+PrVu3Yvr06VbJTL+x1He38JWlLa2iogJGo7F9Zeo2gYGBOH78+EUfU1paetHjS0tLeyynM+jKZ/FHTzzxBPr06XPBf+xknq58Fnv37sWGDRuQmZlphYTOoyufRX5+Pn744Qfccccd2Lp1K3Jzc/HAAw+gubkZy5Yts0Zsh9SVz+L2229HRUUFrr76akiSBIPBgL/85S946qmnrBGZfudS3901NTVobGyEq6trp57H4c4IkeNYsWIFNm3ahM8//xwajUZ0HKdSW1uLuXPnYv369fDz8xMdx+mZTCYEBATgrbfeQlxcHJKSkvD0009j3bp1oqM5nd27d+PFF1/Ev//9b2RkZGDLli349ttv8cILL4iORl3kcGeE/Pz8oFAooNVqO9yv1WoRFBR00ccEBQWZdTx1Tlc+izarVq3CihUrsHPnTowYMaInYzoFcz+LvLw8nD59GjNnzmy/z2QyAQCUSiVycnIQERHRs6EdVFf+vwgODoaLiwsUCkX7fUOGDEFpaSn0ej1UKlWPZnZUXfksnnnmGcydOxf33nsvACA6Ohr19fW477778PTTT3fYG5N61qW+uz09PTt9NghwwDNCKpUKcXFxSElJab/PZDIhJSUFCQkJF31MQkJCh+MB4Pvvv7/k8dQ5XfksAGDlypV44YUXsG3bNowePdoaUR2euZ/F4MGDceTIEWRmZrbfZs2ahWuuuQaZmZkIDQ21ZnyH0pX/L8aPH4/c3Nz2MgoAJ06cQHBwMEtQN3Tls2hoaLig7LQVVIlbd1qVxb67zRvHbR82bdokqdVq6d1335WOHTsm3XfffZK3t7dUWloqSZIkzZ07V3ryySfbj9+3b5+kVCqlVatWSdnZ2dKyZcs4fd5CzP0sVqxYIalUKunTTz+VSkpK2m+1tbWi3oLDMPez+CPOGrMccz+LwsJCycPDQ3rwwQelnJwc6ZtvvpECAgKkf/7zn6LegsMw97NYtmyZ5OHhIX388cdSfn6+tGPHDikiIkK65ZZbRL0Fh1FbWysdPHhQOnjwoARAWr16tXTw4EGpoKBAkiRJevLJJ6W5c+e2H982ff7vf/+7lJ2dLa1du5bT53/v9ddfl/r16yepVCpp7Nix0s8//9z+a5MmTZLmz5/f4fj//e9/UlRUlKRSqaRhw4ZJ3377rZUTOy5zPov+/ftLAC64LVu2zPrBHZC5/1/8HouQZZn7Wezfv1+Kj4+X1Gq1FB4eLv3rX/+SDAaDlVM7JnM+i+bmZunZZ5+VIiIiJI1GI4WGhkoPPPCAdO7cOesHdzC7du266J//bb//8+fPlyZNmnTBY2JjYyWVSiWFh4dL77zzjtmvK5MknssjIiIi5+RwY4SIiIiIOotFiIiIiJwWixARERE5LRYhIiIiclosQkREROS0WISIiIjIabEIERERkdNiESIiIiKnxSJERERETotFiIiIiJwWixAROYSwsDCsWbOmw32xsbF49tlnheQhIvvAIkREREROi0WIiIiInBaLEBERETktFiEiclhGo1F0BCKycSxCROQwtFpt+783NzejqKhIYBoisgcsQkTkMDZu3IidO3fi5MmTePTRR1FdXY28vLwOBYmI6PdYhIjIYcycORMPP/wwoqOjUVVVhX/+85/YsmULdu7cKToaEdkomSRJkugQRETdFRYWhr/+9a/461//KjoKEdkRnhEiIiIip8UiRERERE6Ll8aIiIjIafGMEBERETktFiEiIiJyWixCRERE5LRYhIiIiMhpsQgRERGR02IRIiIiIqfFIkREREROi0WIiIiInBaLEBERETmt/w874A+M0qupswAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ================================================\n",
    "# 2.2.1 ベルヌーイ分布\n",
    "# ================================================\n",
    "\n",
    "b1 = Bern(0.5)\n",
    "print([b1() for _ in range(20)])\n",
    "\n",
    "b2 = Bern(0.9)\n",
    "print([b2() for _ in range(20)])\n",
    "\n",
    "sampler = SimpleSampler()\n",
    "print(f\"E[Bern(x|0.5)] = {sampler(b1, 10000)}\")\n",
    "print(f\"E[Bern(x|0.9)] = {sampler(b2, 10000)}\")\n",
    "\n",
    "H = Entropy()\n",
    "print(f\"H[Bern(x|0.5)] = {H(b1)}\")\n",
    "print(f\"H[Bern(x|0.9)] = {H(b2)}\")\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "x = np.linspace(0, 1, 100)\n",
    "y = [H(Bern(mu), 100000) for mu in x]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, y)\n",
    "ax.set_ylim(0, 0.7)\n",
    "ax.set_xlim(0, 1.0)\n",
    "ax.set_xlabel(\"μ\")\n",
    "ax.set_ylabel(\"H[Bern(x|μ)]\")\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "kl = KLDiv()\n",
    "\n",
    "x = [0.1, 0.5, 0.9]\n",
    "data = pd.DataFrame(\n",
    "    [[kl(p=Bern(mu2), q=Bern(mu1), L=100000) for mu1 in x] for mu2 in x],\n",
    "    index=x,\n",
    "    columns=x\n",
    ")\n",
    "\n",
    "data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 0.35)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGiCAYAAADEJZ3cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlIElEQVR4nO3df0zU9+HH8Rdcy+FP1FLvwNEeqC3rKrCC3Jfvt12befPwazrZXINmCZQ0LrHS79yt7WSrYKffnFrnmJNJ1m9ctVtbZrK6bDV03a34zTKUDGuabtao0YDVO8QFqPgVGu7z/aPrdVdB/aDtvcHnI/mk8Ln358374/Xqsx8+B0mWZVkCAAAwWHKiFwAAAHA1BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAw3qiCpaGhQR6PR6mpqfJ6vWpraxtx7G9+8xsVFRVp2rRpmjRpkgoKCvTiiy/GjXn00UeVlJQUt5WWlo5maQAAYBy6xe4BTU1NCgQCamxslNfrVX19vfx+v44ePaqZM2deNn7GjBn6wQ9+oNzcXKWkpOj3v/+9qqqqNHPmTPn9/ti40tJS/eIXv4h97nQ6R3lKAABgvEmy+8sPvV6v5s+fr+3bt0uSotGosrKy9MQTT2jNmjXXNMd9992nxYsXa/369ZI+vMLS09OjvXv32ls9AAC4Kdi6wjI4OKj29nbV1NTE9iUnJ8vn86m1tfWqx1uWpT/96U86evSoNm3aFPdYS0uLZs6cqenTp+vLX/6yNmzYoNtuu23YeQYGBjQwMBD7PBqN6h//+Iduu+02JSUl2TklAACQIJZl6f3331dmZqaSk698l4qtYOnu7tbQ0JBcLlfcfpfLpXfffXfE43p7ezVr1iwNDAzI4XDoZz/7mb7yla/EHi8tLdXXv/51ZWdn68SJE/r+97+vRYsWqbW1VQ6H47L5gsGgnn32WTtLBwAAhurs7NTnPve5K46xfQ/LaEyZMkWHDx/WhQsXFAqFFAgElJOTo4ceekiStGzZstjYefPmKS8vT7Nnz1ZLS4sWLFhw2Xw1NTUKBAKxz3t7e3XHHXeos7NTU6dO/dTPBwAAXL++vj5lZWVpypQpVx1rK1jS09PlcDgUiUTi9kciEbnd7hGPS05O1pw5cyRJBQUFOnLkiILBYCxYPiknJ0fp6ek6fvz4sMHidDqHvSl36tSpBAsAAGPMtdzOYettzSkpKSosLFQoFIrti0ajCoVCKikpueZ5otFo3D0on3T69GmdP39eGRkZdpYHAADGKdvfEgoEAqqsrFRRUZGKi4tVX1+v/v5+VVVVSZIqKio0a9YsBYNBSR/eb1JUVKTZs2drYGBA+/bt04svvqgdO3ZIki5cuKBnn31WS5culdvt1okTJ/T0009rzpw5cW97BgAANy/bwVJeXq5z586ptrZW4XBYBQUFam5ujt2I29HREXenb39/vx5//HGdPn1aEyZMUG5urn75y1+qvLxckuRwOPT2229r165d6unpUWZmphYuXKj169fzs1gAAICkUfwcFhP19fUpLS1Nvb293MMCAMAYYefvb36XEAAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMN6pgaWhokMfjUWpqqrxer9ra2kYc+5vf/EZFRUWaNm2aJk2apIKCAr344otxYyzLUm1trTIyMjRhwgT5fD4dO3ZsNEsDAADjkO1gaWpqUiAQUF1dnQ4dOqT8/Hz5/X51dXUNO37GjBn6wQ9+oNbWVr399tuqqqpSVVWVXn/99diYzZs3a9u2bWpsbNTBgwc1adIk+f1+Xbp0afRnBgAAxo0ky7IsOwd4vV7Nnz9f27dvlyRFo1FlZWXpiSee0Jo1a65pjvvuu0+LFy/W+vXrZVmWMjMz9d3vfldPPvmkJKm3t1cul0svvPCCli1bdtnxAwMDGhgYiH3e19enrKws9fb2aurUqXZOBwAAJEhfX5/S0tKu6e9vW1dYBgcH1d7eLp/P9/EEycny+XxqbW296vGWZSkUCuno0aP60pe+JEk6efKkwuFw3JxpaWnyer0jzhkMBpWWlhbbsrKy7JwGAAAYY2wFS3d3t4aGhuRyueL2u1wuhcPhEY/r7e3V5MmTlZKSosWLF+unP/2pvvKVr0hS7Dg7c9bU1Ki3tze2dXZ22jkNAAAwxtzyWXyRKVOm6PDhw7pw4YJCoZACgYBycnL00EMPjWo+p9Mpp9N5YxcJAACMZStY0tPT5XA4FIlE4vZHIhG53e4Rj0tOTtacOXMkSQUFBTpy5IiCwaAeeuih2HGRSEQZGRlxcxYUFNhZHgAAGKdsfUsoJSVFhYWFCoVCsX3RaFShUEglJSXXPE80Go3dNJudnS232x03Z19fnw4ePGhrTgAAMH7Z/pZQIBBQZWWlioqKVFxcrPr6evX396uqqkqSVFFRoVmzZikYDEr68AbZoqIizZ49WwMDA9q3b59efPFF7dixQ5KUlJSk1atXa8OGDZo7d66ys7O1du1aZWZmqqys7MadKQAAGLNsB0t5ebnOnTun2tpahcNhFRQUqLm5OXbTbEdHh5KTP75w09/fr8cff1ynT5/WhAkTlJubq1/+8pcqLy+PjXn66afV39+vb33rW+rp6dH999+v5uZmpaam3oBTBAAAY53tn8NiIjvv4wYAAGb41H4OCwAAQCIQLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeKMKloaGBnk8HqWmpsrr9aqtrW3Esc8//7weeOABTZ8+XdOnT5fP57ts/KOPPqqkpKS4rbS0dDRLAwAA45DtYGlqalIgEFBdXZ0OHTqk/Px8+f1+dXV1DTu+paVFy5cv15tvvqnW1lZlZWVp4cKFeu+99+LGlZaW6uzZs7Ht5ZdfHt0ZAQCAcSfJsizLzgFer1fz58/X9u3bJUnRaFRZWVl64okntGbNmqsePzQ0pOnTp2v79u2qqKiQ9OEVlp6eHu3du/ea1jAwMKCBgYHY5319fcrKylJvb6+mTp1q53QAAECC9PX1KS0t7Zr+/rZ1hWVwcFDt7e3y+XwfT5CcLJ/Pp9bW1mua4+LFi/rggw80Y8aMuP0tLS2aOXOm7r77bq1cuVLnz58fcY5gMKi0tLTYlpWVZec0AADAGGMrWLq7uzU0NCSXyxW33+VyKRwOX9Mc3/ve95SZmRkXPaWlpdq9e7dCoZA2bdqk/fv3a9GiRRoaGhp2jpqaGvX29sa2zs5OO6cBAADGmFs+yy+2ceNGvfLKK2ppaVFqamps/7Jly2Ifz5s3T3l5eZo9e7ZaWlq0YMGCy+ZxOp1yOp2fyZoBAEDi2brCkp6eLofDoUgkErc/EonI7XZf8dgtW7Zo48aN+sMf/qC8vLwrjs3JyVF6erqOHz9uZ3kAAGCcsnWFJSUlRYWFhQqFQiorK5P04U23oVBI1dXVIx63efNm/fd//7def/11FRUVXfXrnD59WufPn1dGRoad5QGSJM+a127IPKc2Lr4h8wAArp/ttzUHAgE9//zz2rVrl44cOaKVK1eqv79fVVVVkqSKigrV1NTExm/atElr167Vzp075fF4FA6HFQ6HdeHCBUnShQsX9NRTT+nAgQM6deqUQqGQlixZojlz5sjv99+g0wQAAGOZ7XtYysvLde7cOdXW1iocDqugoEDNzc2xG3E7OjqUnPxxB+3YsUODg4P6xje+ETdPXV2d1q1bJ4fDobffflu7du1ST0+PMjMztXDhQq1fv577VAAAgKRR/BwWE9l5HzfGP74lBABjw6f2c1gAAAASgWABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGC8WxK9ANzcPGteu+45Tm1cfANWcnVjaa0AMN6M6gpLQ0ODPB6PUlNT5fV61dbWNuLY559/Xg888ICmT5+u6dOny+fzXTbesizV1tYqIyNDEyZMkM/n07Fjx0azNAAAMA7ZDpampiYFAgHV1dXp0KFDys/Pl9/vV1dX17DjW1patHz5cr355ptqbW1VVlaWFi5cqPfeey82ZvPmzdq2bZsaGxt18OBBTZo0SX6/X5cuXRr9mQEAgHHDdrBs3bpVK1asUFVVle655x41NjZq4sSJ2rlz57Djf/WrX+nxxx9XQUGBcnNz9T//8z+KRqMKhUKSPry6Ul9fr2eeeUZLlixRXl6edu/erTNnzmjv3r3XdXIAAGB8sBUsg4ODam9vl8/n+3iC5GT5fD61trZe0xwXL17UBx98oBkzZkiSTp48qXA4HDdnWlqavF7viHMODAyor68vbgMAAOOXrWDp7u7W0NCQXC5X3H6Xy6VwOHxNc3zve99TZmZmLFA+Os7OnMFgUGlpabEtKyvLzmkAAIAx5jN9W/PGjRv1yiuv6NVXX1Vqauqo56mpqVFvb29s6+zsvIGrBAAAprH1tub09HQ5HA5FIpG4/ZFIRG63+4rHbtmyRRs3btQf//hH5eXlxfZ/dFwkElFGRkbcnAUFBcPO5XQ65XQ67SwdAACMYbausKSkpKiwsDB2w6yk2A20JSUlIx63efNmrV+/Xs3NzSoqKop7LDs7W263O27Ovr4+HTx48IpzAgCAm4ftHxwXCARUWVmpoqIiFRcXq76+Xv39/aqqqpIkVVRUaNasWQoGg5KkTZs2qba2Vi+99JI8Hk/svpTJkydr8uTJSkpK0urVq7VhwwbNnTtX2dnZWrt2rTIzM1VWVnbjzhQAAIxZtoOlvLxc586dU21trcLhsAoKCtTc3By7abajo0PJyR9fuNmxY4cGBwf1jW98I26euro6rVu3TpL09NNPq7+/X9/61rfU09Oj+++/X83Nzdd1nwsAABg/RvWj+aurq1VdXT3sYy0tLXGfnzp16qrzJSUl6Yc//KF++MMfjmY5AABgnOOXHwIAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeKMKloaGBnk8HqWmpsrr9aqtrW3EsX/729+0dOlSeTweJSUlqb6+/rIx69atU1JSUtyWm5s7mqUBAIBxyHawNDU1KRAIqK6uTocOHVJ+fr78fr+6urqGHX/x4kXl5ORo48aNcrvdI877hS98QWfPno1tf/7zn+0uDQAAjFO2g2Xr1q1asWKFqqqqdM8996ixsVETJ07Uzp07hx0/f/58Pffcc1q2bJmcTueI895yyy1yu92xLT093e7SAADAOGUrWAYHB9Xe3i6fz/fxBMnJ8vl8am1tva6FHDt2TJmZmcrJydE3v/lNdXR0jDh2YGBAfX19cRsAABi/bAVLd3e3hoaG5HK54va7XC6Fw+FRL8Lr9eqFF15Qc3OzduzYoZMnT+qBBx7Q+++/P+z4YDCotLS02JaVlTXqrw0AAMxnxLuEFi1apEceeUR5eXny+/3at2+fenp69Otf/3rY8TU1Nert7Y1tnZ2dn/GKAQDAZ+kWO4PT09PlcDgUiUTi9kcikSveUGvXtGnTdNddd+n48ePDPu50Oq94PwwAABhfbF1hSUlJUWFhoUKhUGxfNBpVKBRSSUnJDVvUhQsXdOLECWVkZNywOQEAwNhl6wqLJAUCAVVWVqqoqEjFxcWqr69Xf3+/qqqqJEkVFRWaNWuWgsGgpA9v1P373/8e+/i9997T4cOHNXnyZM2ZM0eS9OSTT+rhhx/WnXfeqTNnzqiurk4Oh0PLly+/UecJAADGMNvBUl5ernPnzqm2tlbhcFgFBQVqbm6O3Yjb0dGh5OSPL9ycOXNGX/ziF2Ofb9myRVu2bNGDDz6olpYWSdLp06e1fPlynT9/Xrfffrvuv/9+HThwQLfffvt1nh4AABgPbAeLJFVXV6u6unrYxz6KkI94PB5ZlnXF+V555ZXRLAMAANwkjHiXEAAAwJUQLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjHdLoheAscGz5rXrnuPUxsU3YCXjD3+2AHB1XGEBAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGG1WwNDQ0yOPxKDU1VV6vV21tbSOO/dvf/qalS5fK4/EoKSlJ9fX11z0nAAC4udgOlqamJgUCAdXV1enQoUPKz8+X3+9XV1fXsOMvXryonJwcbdy4UW63+4bMCQAAbi62g2Xr1q1asWKFqqqqdM8996ixsVETJ07Uzp07hx0/f/58Pffcc1q2bJmcTucNmRMAANxcbAXL4OCg2tvb5fP5Pp4gOVk+n0+tra2jWsBo5hwYGFBfX1/cBgAAxi9bwdLd3a2hoSG5XK64/S6XS+FweFQLGM2cwWBQaWlpsS0rK2tUXxsAAIwNY/JdQjU1Nert7Y1tnZ2diV4SAAD4FN1iZ3B6erocDocikUjc/kgkMuINtZ/GnE6nc8T7YQAAwPhj6wpLSkqKCgsLFQqFYvui0ahCoZBKSkpGtYBPY04AADC+2LrCIkmBQECVlZUqKipScXGx6uvr1d/fr6qqKklSRUWFZs2apWAwKOnDm2r//ve/xz5+7733dPjwYU2ePFlz5sy5pjkBAMDNzXawlJeX69y5c6qtrVU4HFZBQYGam5tjN812dHQoOfnjCzdnzpzRF7/4xdjnW7Zs0ZYtW/Tggw+qpaXlmuYEAAA3N9vBIknV1dWqrq4e9rGPIuQjHo9HlmVd15wAAODmNibfJQQAAG4uBAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMN6ogqWhoUEej0epqanyer1qa2u74vg9e/YoNzdXqampmjdvnvbt2xf3+KOPPqqkpKS4rbS0dDRLAwAA45DtYGlqalIgEFBdXZ0OHTqk/Px8+f1+dXV1DTv+L3/5i5YvX67HHntMb731lsrKylRWVqZ33nknblxpaanOnj0b215++eXRnREAABh3bAfL1q1btWLFClVVVemee+5RY2OjJk6cqJ07dw47/ic/+YlKS0v11FNP6fOf/7zWr1+v++67T9u3b48b53Q65Xa7Y9v06dNHXMPAwID6+vriNgAAMH7ZCpbBwUG1t7fL5/N9PEFysnw+n1pbW4c9prW1NW68JPn9/svGt7S0aObMmbr77ru1cuVKnT9/fsR1BINBpaWlxbasrCw7pwEAAMYYW8HS3d2toaEhuVyuuP0ul0vhcHjYY8Lh8FXHl5aWavfu3QqFQtq0aZP279+vRYsWaWhoaNg5a2pq1NvbG9s6OzvtnAYAABhjbkn0AiRp2bJlsY/nzZunvLw8zZ49Wy0tLVqwYMFl451Op5xO52e5RAAAkEC2rrCkp6fL4XAoEonE7Y9EInK73cMe43a7bY2XpJycHKWnp+v48eN2lgcAAMYpW8GSkpKiwsJChUKh2L5oNKpQKKSSkpJhjykpKYkbL0lvvPHGiOMl6fTp0zp//rwyMjLsLA8AAIxTtt8lFAgE9Pzzz2vXrl06cuSIVq5cqf7+flVVVUmSKioqVFNTExv/7W9/W83NzfrRj36kd999V+vWrdNf//pXVVdXS5IuXLigp556SgcOHNCpU6cUCoW0ZMkSzZkzR36//wadJgAAGMts38NSXl6uc+fOqba2VuFwWAUFBWpubo7dWNvR0aHk5I876N///d/10ksv6ZlnntH3v/99zZ07V3v37tW9994rSXI4HHr77be1a9cu9fT0KDMzUwsXLtT69eu5TwUAAEga5U231dXVsSskn9TS0nLZvkceeUSPPPLIsOMnTJig119/fTTLwDA8a167IfOc2rj4hsyDxODfAwDjDb9LCAAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPFuSfQCbmaeNa9d9xynNi6+ASsBrg3/zgJIFK6wAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4owqWhoYGeTwepaamyuv1qq2t7Yrj9+zZo9zcXKWmpmrevHnat29f3OOWZam2tlYZGRmaMGGCfD6fjh07NpqlAQCAcch2sDQ1NSkQCKiurk6HDh1Sfn6+/H6/urq6hh3/l7/8RcuXL9djjz2mt956S2VlZSorK9M777wTG7N582Zt27ZNjY2NOnjwoCZNmiS/369Lly6N/swAAMC4YfsHx23dulUrVqxQVVWVJKmxsVGvvfaadu7cqTVr1lw2/ic/+YlKS0v11FNPSZLWr1+vN954Q9u3b1djY6Msy1J9fb2eeeYZLVmyRJK0e/duuVwu7d27V8uWLbtszoGBAQ0MDMQ+7+3tlST19fXZPZ1rcm/d69c9xzvP+i/bFx24eN3zfvKcb8Scn9a8wz0//BnwZyB9eq8xAGb76L8HlmVdfbBlw8DAgOVwOKxXX301bn9FRYX11a9+ddhjsrKyrB//+Mdx+2pra628vDzLsizrxIkTliTrrbfeihvzpS99yfqv//qvYeesq6uzJLGxsbGxsbGNg62zs/OqDWLrCkt3d7eGhobkcrni9rtcLr377rvDHhMOh4cdHw6HY49/tG+kMZ9UU1OjQCAQ+zwajeof//iHbrvtNiUlJdk5pRuir69PWVlZ6uzs1NSpUz/zrw97eL7GFp6vsYfnbGxJ5PNlWZbef/99ZWZmXnXsmPxdQk6nU06nM27ftGnTErOYfzF16lRenGMIz9fYwvM19vCcjS2Jer7S0tKuaZytm27T09PlcDgUiUTi9kciEbnd7mGPcbvdVxz/0T/tzAkAAG4utoIlJSVFhYWFCoVCsX3RaFShUEglJSXDHlNSUhI3XpLeeOON2Pjs7Gy53e64MX19fTp48OCIcwIAgJuL7W8JBQIBVVZWqqioSMXFxaqvr1d/f3/sXUMVFRWaNWuWgsGgJOnb3/62HnzwQf3oRz/S4sWL9corr+ivf/2rfv7zn0uSkpKStHr1am3YsEFz585Vdna21q5dq8zMTJWVld24M/0UOZ1O1dXVXfZtKpiJ52ts4fkae3jOxpax8nwlWda1vJco3vbt2/Xcc88pHA6roKBA27Ztk9frlSQ99NBD8ng8euGFF2Lj9+zZo2eeeUanTp3S3LlztXnzZv3nf/5n7HHLslRXV6ef//zn6unp0f3336+f/exnuuuuu67/DAEAwJg3qmABAAD4LPG7hAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYrlNDQ4M8Ho9SU1Pl9XrV1taW6CVhBOvWrVNSUlLclpubm+hl4Z/+93//Vw8//LAyMzOVlJSkvXv3xj1uWZZqa2uVkZGhCRMmyOfz6dixY4lZLCRd/Tl79NFHL3vNlZaWJmaxUDAY1Pz58zVlyhTNnDlTZWVlOnr0aNyYS5cuadWqVbrttts0efJkLV269LIf7JooBMt1aGpqUiAQUF1dnQ4dOqT8/Hz5/X51dXUlemkYwRe+8AWdPXs2tv35z39O9JLwT/39/crPz1dDQ8Owj2/evFnbtm1TY2OjDh48qEmTJsnv9+vSpUuf8Urxkas9Z5JUWloa95p7+eWXP8MV4l/t379fq1at0oEDB/TGG2/ogw8+0MKFC9Xf3x8b853vfEe/+93vtGfPHu3fv19nzpzR17/+9QSu+l9c9dcjYkTFxcXWqlWrYp8PDQ1ZmZmZVjAYTOCqMJK6ujorPz8/0cvANZAU91vho9Go5Xa7reeeey62r6enx3I6ndbLL7+cgBXikz75nFmWZVVWVlpLlixJyHpwdV1dXZYka//+/ZZlffiauvXWW609e/bExhw5csSSZLW2tiZqmTFcYRmlwcFBtbe3y+fzxfYlJyfL5/OptbU1gSvDlRw7dkyZmZnKycnRN7/5TXV0dCR6SbgGJ0+eVDgcjnu9paWlyev18nozXEtLi2bOnKm7775bK1eu1Pnz5xO9JPxTb2+vJGnGjBmSpPb2dn3wwQdxr7Pc3FzdcccdRrzOCJZR6u7u1tDQkFwuV9x+l8ulcDicoFXhSrxer1544QU1Nzdrx44dOnnypB544AG9//77iV4aruKj1xSvt7GltLRUu3fvVigU0qZNm7R//34tWrRIQ0NDiV7aTS8ajWr16tX6j//4D917772SPnydpaSkaNq0aXFjTXmd2f5dQsBYtWjRotjHeXl58nq9uvPOO/XrX/9ajz32WAJXBoxPy5Yti308b9485eXlafbs2WppadGCBQsSuDKsWrVK77zzzpi6j48rLKOUnp4uh8Nx2d3TkUhEbrc7QauCHdOmTdNdd92l48ePJ3opuIqPXlO83sa2nJwcpaen85pLsOrqav3+97/Xm2++qc997nOx/W63W4ODg+rp6Ykbb8rrjGAZpZSUFBUWFioUCsX2RaNRhUIhlZSUJHBluFYXLlzQiRMnlJGRkeil4Cqys7PldrvjXm99fX06ePAgr7cx5PTp0zp//jyvuQSxLEvV1dV69dVX9ac//UnZ2dlxjxcWFurWW2+Ne50dPXpUHR0dRrzO+JbQdQgEAqqsrFRRUZGKi4tVX1+v/v5+VVVVJXppGMaTTz6phx9+WHfeeafOnDmjuro6ORwOLV++PNFLgz4MyH/9P++TJ0/q8OHDmjFjhu644w6tXr1aGzZs0Ny5c5Wdna21a9cqMzNTZWVliVv0Te5Kz9mMGTP07LPPaunSpXK73Tpx4oSefvppzZkzR36/P4GrvnmtWrVKL730kn77299qypQpsftS0tLSNGHCBKWlpemxxx5TIBDQjBkzNHXqVD3xxBMqKSnRv/3bvyV49eJtzdfrpz/9qXXHHXdYKSkpVnFxsXXgwIFELwkjKC8vtzIyMqyUlBRr1qxZVnl5uXX8+PFELwv/9Oabb1qSLtsqKysty/rwrc1r1661XC6X5XQ6rQULFlhHjx5N7KJvcld6zi5evGgtXLjQuv32261bb73VuvPOO60VK1ZY4XA40cu+aQ33XEmyfvGLX8TG/N///Z/1+OOPW9OnT7cmTpxofe1rX7POnj2buEX/iyTLsqzPPpMAAACuHfewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMN7/A9SYzzrWFwupAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ================================================\n",
    "# 2.2.2 二項分布\n",
    "# ================================================\n",
    "\n",
    "class Bin(DiscreteProb):\n",
    "    def __init__(self, M: int, mu: float):\n",
    "        self.M = M\n",
    "        self.mu = mu\n",
    "\n",
    "    def pmf(self, m: int) -> float:\n",
    "        return math.comb(self.M, m) * self.mu**m * (1 - self.mu)**(self.M - m)\n",
    "\n",
    "    def __call__(self) -> int:\n",
    "        return sum([1 if random.random() < self.mu else 0 for _ in range(self.M)])\n",
    "\n",
    "\n",
    "for (M, mu) in [(10, 0.5), (10, 0.2), (10, 0.9), (100, 0.5)]:\n",
    "    b = Bin(M=M, mu=mu)\n",
    "    x = range(0, 21)\n",
    "    y = [b.pmf(m) for m in x]\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.bar(x, y)\n",
    "    ax.set_ylim(0, 0.35)\n",
    "    ax.title(f\"M={M}, μ={mu}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.11.7 (main, Dec  4 2023, 18:10:11) [Clang 15.0.0 (clang-1500.1.0.2.5)]\n",
      "Requirement already satisfied: matplotlib in /opt/homebrew/lib/python3.11/site-packages (3.8.1)\n",
      "Requirement already satisfied: pandas in /opt/homebrew/lib/python3.11/site-packages (2.2.3)\n",
      "Requirement already satisfied: numpy in /opt/homebrew/lib/python3.11/site-packages (1.26.1)\n",
      "Requirement already satisfied: contourpy>=1.0.1 in /opt/homebrew/lib/python3.11/site-packages (from matplotlib) (1.1.1)\n",
      "Requirement already satisfied: cycler>=0.10 in /opt/homebrew/lib/python3.11/site-packages (from matplotlib) (0.12.1)\n",
      "Requirement already satisfied: fonttools>=4.22.0 in /opt/homebrew/lib/python3.11/site-packages (from matplotlib) (4.43.1)\n",
      "Requirement already satisfied: kiwisolver>=1.3.1 in /opt/homebrew/lib/python3.11/site-packages (from matplotlib) (1.4.5)\n",
      "Requirement already satisfied: packaging>=20.0 in /Users/keisukefukuda/Library/Python/3.11/lib/python/site-packages (from matplotlib) (23.2)\n",
      "Requirement already satisfied: pillow>=8 in /opt/homebrew/lib/python3.11/site-packages (from matplotlib) (10.1.0)\n",
      "Requirement already satisfied: pyparsing>=2.3.1 in /opt/homebrew/lib/python3.11/site-packages (from matplotlib) (3.1.1)\n",
      "Requirement already satisfied: python-dateutil>=2.7 in /Users/keisukefukuda/Library/Python/3.11/lib/python/site-packages (from matplotlib) (2.8.2)\n",
      "Requirement already satisfied: pytz>=2020.1 in /opt/homebrew/lib/python3.11/site-packages (from pandas) (2024.2)\n",
      "Requirement already satisfied: tzdata>=2022.7 in /opt/homebrew/lib/python3.11/site-packages (from pandas) (2024.2)\n",
      "Requirement already satisfied: six>=1.5 in /opt/homebrew/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.3.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m24.3.1\u001b[0m\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpython3.11 -m pip install --upgrade pip\u001b[0m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=['/opt/homebrew/opt/python@3.11/bin/python3.11', '-m', 'pip', 'install', 'matplotlib', 'pandas', 'numpy'], returncode=0)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "math.comb(10, 3)\n",
    "\n",
    "import sys\n",
    "print(sys.version)\n",
    "from subprocess import run\n",
    "run([sys.executable, \"-m\", \"pip\", \"install\", \"matplotlib\", \"pandas\", \"numpy\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.11",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}