bayes_cunjugate_prior.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9b55cb09-d779-4618-9c32-f83f4f5a1f4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a044bf21-ffed-4137-9bad-d002d26b5be9",
   "metadata": {},
   "source": [
    "# 正規分布の平均パラメータの事後分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5343ddda-06c6-42ef-b7b7-0f81e0c5357e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc472b6-64c2-4afb-92e6-71f43b05e392",
   "metadata": {},
   "source": [
    "## トイデータ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3737d117-ef11-4975-8183-e26f5ea934dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'x')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEGCAYAAABB6hAxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAS6klEQVR4nO3db6hldb3H8fenmSnlKgnOgYb54+miTyryT+eaXeEihmAp+iADgzKjGOgmGTcI7YGRXLj1pKKMZFBJ+6dhEZMpIWhUXJw8M42mTl3mhhdHBMexxqSyO97vfXCW3XN357jXzG+f2fvs837BwvXnO3t/f7POWfNx7bX2SlUhSZKkY/OacTcgSZK0mhmmJEmSGhimJEmSGhimJEmSGhimJEmSGqwf1xtv3LixZmdnx/X2kiRJve3evfu5qppZatvYwtTs7Czz8/PjentJkqTekvzXctv8mE+SJKmBYUqSJKmBYUqSJKmBYUqSJKmBYUqSJKmBYUqSJKlB7zCVZF2SXya5Z4ltr0tyV5L9SXYlmR1pl5IkSRPqaM5MXQvsW2bbh4HfVdXpwBeBz7c2JkmStBr0ClNJtgCXALcsU3I5cHs3fzfwziRpb0+SJGmy9f0G9C8BnwJOXmb7ZuApgKo6kuQwcCrw3OKiJNuB7QDbtm07hnYlrYTZ63407hZG5snPXTLuFiStMUPPTCW5FHi2qna3vllV7aiquaqam5lZ8vE2kiRJq0qfj/nOBy5L8iRwJ3Bhkm8O1DwNbAVIsh54PXBohH1KkiRNpKFhqqqur6otVTULXAk8UFXvHyjbCXywm7+iq6mRdipJkjSB+l4z9TeS3AjMV9VO4FbgG0n2A8+zELokSZKm3lGFqar6CfCTbv6GRev/DLx3lI1JkiStBn4DuiRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUgPDlCRJUoOhYSrJCUl+keSRJI8n+ewSNVcnOZhkbzd9ZGXalSRJmizre9S8BFxYVS8m2QD8PMl9VfXQQN1dVXXN6FuUJEmaXEPDVFUV8GK3uKGbaiWbkiRJWi16XTOVZF2SvcCzwP1VtWuJsvckeTTJ3Um2jrJJSZKkSdUrTFXVy1V1FrAFODfJWwZKfgjMVtVbgfuB25d6nSTbk8wnmT948GBD25IkSZPhqO7mq6rfAw8CFw+sP1RVL3WLtwBvW+bP76iquaqam5mZOYZ2JUmSJkufu/lmkpzSzZ8IXAT8eqBm06LFy4B9I+xRkiRpYvW5m28TcHuSdSyEr+9W1T1JbgTmq2on8PEklwFHgOeBq1eqYUmSpEnS526+R4Gzl1h/w6L564HrR9uaJEnS5PMb0CVJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoYpiRJkhoMDVNJTkjyiySPJHk8yWeXqHldkruS7E+yK8nsinQrSZI0YfqcmXoJuLCqzgTOAi5Oct5AzYeB31XV6cAXgc+PtEtJkqQJNTRM1YIXu8UN3VQDZZcDt3fzdwPvTJKRdSlJkjSh1vcpSrIO2A2cDny1qnYNlGwGngKoqiNJDgOnAs8NvM52YDvAtm3b2jqXJsDsdT8adwsaME375MnPXTLuFjTF/F0ZnV4XoFfVy1V1FrAFODfJW47lzapqR1XNVdXczMzMsbyEJEnSRDmqu/mq6vfAg8DFA5ueBrYCJFkPvB44NIL+JEmSJlqfu/lmkpzSzZ8IXAT8eqBsJ/DBbv4K4IGqGryuSpIkaer0uWZqE3B7d93Ua4DvVtU9SW4E5qtqJ3Ar8I0k+4HngStXrGNJkqQJMjRMVdWjwNlLrL9h0fyfgfeOtjVJkqTJ5zegS5IkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNTBMSZIkNRgappJsTfJgkieSPJ7k2iVqLkhyOMnebrphZdqVJEmaLOt71BwBPllVe5KcDOxOcn9VPTFQ97OqunT0LUqSJE2uoWemquqZqtrTzf8B2AdsXunGJEmSVoOjumYqySxwNrBric3vSPJIkvuSvHmZP789yXyS+YMHDx59t5IkSROmd5hKchLwPeATVfXCwOY9wGlVdSbwFeAHS71GVe2oqrmqmpuZmTnGliVJkiZHrzCVZAMLQepbVfX9we1V9UJVvdjN3wtsSLJxpJ1KkiRNoD538wW4FdhXVV9YpuYNXR1Jzu1e99AoG5UkSZpEfe7mOx/4APCrJHu7dZ8GtgFU1c3AFcBHkxwB/gRcWVU1+nYlSZImy9AwVVU/BzKk5ibgplE1JUmStFr4DeiSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNhoapJFuTPJjkiSSPJ7l2iZok+XKS/UkeTXLOyrQrSZI0Wdb3qDkCfLKq9iQ5Gdid5P6qemJRzbuAM7rp7cDXuv9KkiRNtaFnpqrqmara083/AdgHbB4ouxy4oxY8BJySZNPIu5UkSZowfc5M/VWSWeBsYNfAps3AU4uWD3Trnhn489uB7QDbtm07ylaPzex1Pzou77PSnvzcJeNuQZKOmcdiTbPeF6AnOQn4HvCJqnrhWN6sqnZU1VxVzc3MzBzLS0iSJE2UXmEqyQYWgtS3qur7S5Q8DWxdtLylWydJkjTV+tzNF+BWYF9VfWGZsp3AVd1dfecBh6vqmWVqJUmSpkafa6bOBz4A/CrJ3m7dp4FtAFV1M3Av8G5gP/BH4EMj71SSJGkCDQ1TVfVzIENqCvjYqJqSJElaLfwGdEmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAaGKUmSpAZDw1SS25I8m+SxZbZfkORwkr3ddMPo25QkSZpM63vUfB24CbjjVWp+VlWXjqQjSZKkVWTomamq+inw/HHoRZIkadUZ1TVT70jySJL7krx5uaIk25PMJ5k/ePDgiN5akiRpfEYRpvYAp1XVmcBXgB8sV1hVO6pqrqrmZmZmRvDWkiRJ49Ucpqrqhap6sZu/F9iQZGNzZ5IkSatAc5hK8oYk6ebP7V7zUOvrSpIkrQZD7+ZL8h3gAmBjkgPAZ4ANAFV1M3AF8NEkR4A/AVdWVa1Yx5IkSRNkaJiqqvcN2X4TC1+dIEmStOb4DeiSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNDFOSJEkNhoapJLcleTbJY8tsT5IvJ9mf5NEk54y+TUmSpMnU58zU14GLX2X7u4Azumk78LX2tiRJklaHoWGqqn4KPP8qJZcDd9SCh4BTkmwaVYOSJEmTbP0IXmMz8NSi5QPdumcGC5NsZ+HsFdu2bRvBW2s1mr3uR+NuQVoV/F2ZPO4TLeW4XoBeVTuqaq6q5mZmZo7nW0uSJK2IUYSpp4Gti5a3dOskSZKm3ijC1E7gqu6uvvOAw1X1Nx/xSZIkTaOh10wl+Q5wAbAxyQHgM8AGgKq6GbgXeDewH/gj8KGValaSJGnSDA1TVfW+IdsL+NjIOpIkSVpF/AZ0SZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBoYpSZKkBr3CVJKLk/wmyf4k1y2x/eokB5Ps7aaPjL5VSZKkybN+WEGSdcBXgYuAA8DDSXZW1RMDpXdV1TUr0KMkSdLE6nNm6lxgf1X9tqr+AtwJXL6ybUmSJK0OfcLUZuCpRcsHunWD3pPk0SR3J9m61Asl2Z5kPsn8wYMHj6FdSZKkyTKqC9B/CMxW1VuB+4Hblyqqqh1VNVdVczMzMyN6a0mSpPHpE6aeBhafadrSrfurqjpUVS91i7cAbxtNe5IkSZOtT5h6GDgjyRuTvBa4Eti5uCDJpkWLlwH7RteiJEnS5Bp6N19VHUlyDfBjYB1wW1U9nuRGYL6qdgIfT3IZcAR4Hrh6BXuWJEmaGEPDFEBV3QvcO7DuhkXz1wPXj7Y1SZKkyec3oEuSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDUwTEmSJDXoFaaSXJzkN0n2J7luie2vS3JXt31XktmRdypJkjSBhoapJOuArwLvAt4EvC/JmwbKPgz8rqpOB74IfH7UjUqSJE2iPmemzgX2V9Vvq+ovwJ3A5QM1lwO3d/N3A+9MktG1KUmSNJnW96jZDDy1aPkA8PblaqrqSJLDwKnAc4uLkmwHtneLLyY5NFizRmzkKMed6TnXd9RjnyJrdexrddzg2Nfi2NfquGGMYz9O/0aettyGPmFqZKpqB7DjleUk81U1dzx7mARrddzg2Nfi2NfquMGxr8Wxr9Vxw9oee5+P+Z4Gti5a3tKtW7ImyXrg9cChUTQoSZI0yfqEqYeBM5K8MclrgSuBnQM1O4EPdvNXAA9UVY2uTUmSpMk09GO+7hqoa4AfA+uA26rq8SQ3AvNVtRO4FfhGkv3A8ywErj52DC+ZSmt13ODY16K1Om5w7GvRWh03rOGxxxNIkiRJx85vQJckSWpgmJIkSWqw4mFqLT+KpsfYr05yMMnebvrIOPoctSS3JXk2yWPLbE+SL3d/L48mOed497hSeoz9giSHF+3zG453jyshydYkDyZ5IsnjSa5domYq93vPsU/rfj8hyS+SPNKN/bNL1EzdMb7nuKfy+P6KJOuS/DLJPUtsm7p9PlRVrdjEwgXr/wn8PfBa4BHgTQM1/wzc3M1fCdy1kj0dr6nn2K8Gbhp3rysw9n8CzgEeW2b7u4H7gADnAbvG3fNxHPsFwD3j7nMFxr0JOKebPxn4jyV+3qdyv/cc+7Tu9wAndfMbgF3AeQM1U3eM7znuqTy+LxrfvwDfXurnehr3+bBppc9MreVH0fQZ+1Sqqp+ycFfnci4H7qgFDwGnJNl0fLpbWT3GPpWq6pmq2tPN/wHYx8KTERabyv3ec+xTqduXL3aLG7pp8K6mqTvG9xz31EqyBbgEuGWZkqnb58OsdJha6lE0gweZ//coGuCVR9Gsdn3GDvCe7iOPu5NsXWL7NOr7dzOt3tF9PHBfkjePu5lR607pn83C/60vNvX7/VXGDlO637uPe/YCzwL3V9Wy+32ajvE9xg3Te3z/EvAp4H+W2T6V+/zVeAH6eP0QmK2qtwL3839JXtNrD3BaVZ0JfAX4wXjbGa0kJwHfAz5RVS+Mu5/jacjYp3a/V9XLVXUWC0/HODfJW8bc0nHRY9xTeXxPcinwbFXtHncvk2Slw9RafhTN0LFX1aGqeqlbvAV423Hqbdz6/FxMpap64ZWPB6rqXmBDko1jbmskkmxgIUx8q6q+v0TJ1O73YWOf5v3+iqr6PfAgcPHApmk9xgPLj3uKj+/nA5cleZKFy1cuTPLNgZqp3udLWekwtZYfRTN07APXi1zGwrUWa8FO4Kru7q7zgMNV9cy4mzoekrzhlWsHkpzLwu/gqj/IdGO6FdhXVV9Ypmwq93ufsU/xfp9Jcko3fyJwEfDrgbKpO8b3Gfe0Ht+r6vqq2lJVsyz8u/ZAVb1/oGzq9vkwQx8n06JW9lE0E63n2D+e5DLgCAtjv3psDY9Qku+wcPfSxiQHgM+wcIEmVXUzcC8Ld3btB/4IfGg8nY5ej7FfAXw0yRHgT8CVU3KQOR/4APCr7joSgE8D22Dq93ufsU/rft8E3J5kHQsB8btVdc8aOMb3GfdUHt+Xswb2+avycTKSJEkNvABdkiSpgWFKkiSpgWFKkiSpgWFKkiSpgWFKkiSpgWFKkiSpgWFKkiSpgWFK0lRI8g/dQ2VPSPJ3SR5fK8+JkzRefmmnpKmR5F+BE4ATgQNV9W9jbknSGmCYkjQ1uudgPgz8GfjHqnp5zC1JWgP8mE/SNDkVOAk4mYUzVJK04jwzJWlqJNkJ3Am8EdhUVdeMuSVJa8D6cTcgSaOQ5Crgv6vq20nWAf+e5MKqemDcvUmabp6ZkiRJauA1U5IkSQ0MU5IkSQ0MU5IkSQ0MU5IkSQ0MU5IkSQ0MU5IkSQ0MU5IkSQ3+F6a0zqwa+ja+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu_set = 2.5\n",
    "var_set = 1.0\n",
    "n = 20\n",
    "sample = stats.norm(mu_set, np.sqrt(var_set)).rvs(n)\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.hist(sample, bins=10);\n",
    "ax.set_xlabel(\"x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e95af448-0128-4157-9a84-77eb13cbf6d1",
   "metadata": {},
   "source": [
    "## 尤度関数\n",
    "\n",
    "$$\n",
    "L\\left( \\mu \\right) = \\prod^n_{i=1} \\mathcal{N}(x_i | \\mu, \\sigma^2) = \\prod^n_{i=1} \\frac{1}{\\sqrt{2\\pi}\\sigma} \\exp \\left[ - \\frac{(x_i - \\mu)^2}{2\\sigma^2} \\right]\\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "63729686-2f6e-455b-ae16-97d39c9c2a7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood_norm(mu, data):\n",
    "    ps = [stats.norm(mu, np.sqrt(var_set)).pdf(x) for x in data]\n",
    "    return np.prod(ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "823cc7a6-97be-4422-8b30-a5674f75e8d6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAETCAYAAAABJsPtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkn0lEQVR4nO3de3TcZ33n8c93RndpZPmii23JV8m3JrbsugEawiUXNqEtKbRQaOFs2XSz5UCX7vaUpduePac923b3tNsDZQttSmk4hQ1tIBRIKCGQpCmEksjxNb7EkuOLZOsS2bpb1/nuH5oxRpVtyZr5/eY3836d4xNpNJ7nkzmJ9NHze57nZ+4uAAAAZEcs7AAAAAD5jLIFAACQRZQtAACALKJsAQAAZBFlCwAAIIsoWwAAAFmU82XLzD5nZr1mdiRDr/ctMxsws8ev8fU/N7ORTIwFAACQ82VL0sOS7s3g6/2JpA/M9wUz2ytpeQbHAgAABS7ny5a7Pyfp4tWPmdnm1AzVPjP7FzPbtojX+66k4bmPm1lcs0XsY0vNDAAAkFYUdoCb9JCkX3f3k2b2OkmflnTnEl/zI5K+7u4XzGzJAQEAAKQIli0zq5L005IevaoUlaa+9i5JfzDPX+ty9393nddcI+ndkt6S0bAAAKDgRa5safbS54C7t879grs/Jumxm3jN3ZKaJbWnClyFmbW7e/NSggIAAOT8mq253H1I0qtm9m5Jslm7lviaT7h7g7tvcPcNksYoWgAAIBNyvmyZ2SOSfiBpq5l1mtkDkn5F0gNmdlDSy5LuX8Tr/YukRyXdlXq9a15eBAAAWCpz97AzAAAA5K2cn9kCAACIskAXyJtZjaTPSrpFkkv6D+7+g2s9f9WqVb5hw4ZgwgEAACzBvn37XnP32rmPB70b8ZOSvuXuv2hmJZIqrvfkDRs2qK2tLZhkAAAAS2BmZ+Z7PLCyZWbLJL1J0q9KkrtPSpoManwAAIAwBLlma6OkPkl/a2b7zeyzZlY590lm9qCZtZlZW19fX4DxAAAAMi/IslUkaY+kz7j7bkmjkj4+90nu/pC773X3vbW1/+ayJwAAQKQEWbY6JXW6+w9Tn39Zs+ULAAAgbwVWtty9W9I5M9uaeuguSUeDGh8AACAMQe9G/A1JX0ztRDwl6YMBjw8AABCoQMuWux+QtDfIMQEAAMIU9MwWABSk0Ylpnekf05n+UZ29OKY7Wmq1Y0112LEABICyBQBZ4u76yCP79cKrF9U3PPFjX/vmkW597cO3h5QMQJAoWwCQJS+dvaQnDl3QXdvqtGf9cm1YWan1Kyv07Ile/em3X9GJ7mFtbUiEHRNAllG2ACBLvnbgvEqLYvrk+3arqvRH325XLyvTJ75zUo+2ndPv/eyOEBMCCEKQ52wBQMGYmknqiUMXdPeO+h8rWpK0sqpUd2+v12P7uzQ5nQwpIYCgULYAIAu+3/6a+kcndf+uNfN+/Zd+qkkXRyf19PGegJMBCBplCwCy4OsHz6u6rEhv3jr/bcfuaFml+upS/UNbZ8DJAASNsgUAGTY+NaMnj3TrvltWq7QoPu9ziuIx/cKeRj17olfdg+MBJwQQJMoWAGTYd4/1anRyRve3zn8JMe09e5uUdOkrLzG7BeQzyhYAZNjXDnSpLlGq121aed3nbVhVqds2rtCjbefk7gGlAxA0yhYAZNDg5Sk9e6JPP7drjeIxu+Hz37O3Saf7x/TCqxcDSAcgDJQtAMigJ490a3ImqXdcYxfiXG+/tUFVpUUslAfyGGULADLoawe7tGFlhXY2LlvQ8ytKivRzu1brm4cvaHh8KsvpAISBsgUAGdI7NK7nO/r1jta1MrvxJcS09+xt0uWpGT1+6EIW0wEIC2ULADLkG4cuyF0LvoSY1tpUo+a6Kn39wPksJQMQJsoWAGTI1w+e10+sqVZzXdWi/p6Z6fbNK3Wwc0AzSXYlAvmGsgUAGXBh8LIOnhvQzy1yVittZ2ONxiZn1N47kuFkAMJG2QKADDh4bkCS9PobnK11LbuaamZfp3MgM4EA5AzKFgBkwKHOQRXFTNsaEjf19zetqlSitEiHKFtA3qFsAUAGHO4a1Jb6hMqK578X4o3EYqZb1i7TwXODGU4GIGyULQBYInfX4a5B3bp2YWdrXcuuphod7x7SxPRMhpIByAWULQBYos5LlzUwNqVbF3iQ6bXsalymqRnXsQvDGUoGIBdQtgBgiQ53zV76W+rM1s70IvnUYnsA+YGyBQBLdLhrUMVx07bVN7c4Pm3NsjKtqiplRyKQZyhbALBER1KL40uLbm5xfJqZaVfjMh3qZJE8kE8CLVtmdtrMDpvZATNrC3JsAMgGd9ehzsEF33j6RnY21qijb4SbUgN5JIyZrbe6e6u77w1hbADIqM5LlzV4eUq3LHG9VtqupmVy/9E6MADRx2VEAFiC9CW/nWtrMvJ6Oxtrfux1AURf0GXLJX3bzPaZ2YPzPcHMHjSzNjNr6+vrCzgeACxOenH8lobF3Xz6WlZUlqhpRTknyQN5JOiy9UZ33yPpPkkfNrM3zX2Cuz/k7nvdfW9tbW3A8QBgcQ53DWhbQ/WSF8dfbVdjDSfJA3kk0LLl7l2pf/ZK+qqk24IcHwAyyd11uHMwY+u10nY11qhr4LJeG5nI6OsCCEdgZcvMKs0skf5Y0tskHQlqfADItLMXxzQ0Pp2xnYhp6dfjUiKQH4Kc2aqX9D0zOyjpBUlPuPu3AhwfADIqUyfHz3XL2mWKmXSAS4lAXigKaiB3PyVpV1DjAUC2He4cVEk8pi31Szs5fq7K0iK11CWY2QLyBEc/AMBNOtw1qG2rEyopyvy30p2pk+TdPeOvDSBYlC0AuAnursNdmV8cn7azqUYXRyfVeelyVl4fQHAoWwBwE870j2l4fFo7s1S2WlOHm3JTaiD6KFsAcBMOpRbHZ2tma2tDQiXxGCfJA3mAsgUAN+FI16BKijK/OD6tpCim7WuqdeDcQFZeH0BwKFsAcBMOdQ5oe0N2Fsen3bq2WsfOD7FIHog4yhYALFIy6Xq5a0i3Zvgw07m21Cc0PDGtniFOkgeijLIFAIt0un9UwxPTGT/MdK7mutmbW5/sHc7qOACyi7IFAIt0vHu2/OxYnd2y1VI3ux7sZM9IVscBkF2ULQBYpPbeEZn9aOYpW1ZVlWh5RbFO9lK2gCijbAHAIp3sHdHamnKVl8SzOo6ZqaUuoXYuIwKRRtkCgEU62TOslizPaqU111fplZ4RdiQCEUbZAoBFmEm6Tr02qpYsna81V0tdlQYvT6lvhB2JQFRRtgBgEc5dHNPkdFLNtcHMbKUXybezSB6ILMoWACxCerF6c31AZas+ffwDZQuIKsoWACxCe7psBbRmqy5RqkRZEWdtARFG2QKARTjZO6z66lJVlxUHMp6ZaUt9grO2gAijbAHAInT0jlxZRxWUlrqqKzNqAKKHsgUAC+TuOtk7EtglxLTmuir1j06qnx2JQCRRtgBggc4PjmtscibwspU+ZoLZLSCaKFsAsEDpshPUgaZp6fFeoWwBkUTZAoAFOtkzuyMwqANN01YvK1NlSVztPexIBKKIsgUAC9TeO6IVlSVaUVkS6Lhmpub6BGdtARFF2QKABWoPYXF82pa6KsoWEFGULQBYgPROxKDXa6W11Fepb3hCA2OToYwP4OYFXrbMLG5m+83s8aDHBoCb1TcyocHLU6HNbF25RyKzW0DkhDGz9VFJx0IYFwBu2o92Iga7OD4tXfK4lAhET6Bly8waJf2MpM8GOS4ALNWVshXQDajnWltTrvLiOLftASIo6JmtT0j6mKTktZ5gZg+aWZuZtfX19QUWDACu52TPiBKlRapLlIYyfixmaq6r4obUQAQFVrbM7Gcl9br7vus9z90fcve97r63trY2oHQAcH3tvSNqrq+SmYWWoaWuipktIIKCnNm6XdI7zOy0pC9JutPMvhDg+ABw08LciZjWUp9Q99C4hsanQs0BYHECK1vu/jvu3ujuGyS9V9LT7v7+oMYHgJs1MDap10YmQtuJmJYue+xIBKKFc7YA4AbC3omYll6c386lRCBSisIY1N2flfRsGGMDwGKlj1sIe2arcXmFSotiLJIHIoaZLQC4gZM9IyovjmttTXmoOeIx0+ZabtsDRA1lCwBuoL1vRJvrKhWLhbcTMa2lnh2JQNRQtgDgBtp7hkNfr5XWUlelroHLGp2YDjsKgAWibAHAdYxMTOv84Hjo67XS0jlO9Y2GnATAQlG2AOA6OnJkcXza5trZHB19XEoEooKyBQDXkSs7EdPWr6xUPGactQVECGULAK6jvXdExXHT+hUVYUeRJJUUxbR+RQUzW0CEULYA4Do6+ka0YWWliuK58+1yc10VM1tAhOTOdw8AyEEdvSM5cwkxbXNtlU73j2p6Jhl2FAALQNkCgGuYnE7qzMWxK4vSc0VzXZWmZlxnL46FHQXAAlC2AOAazvSPaibpOTizVSlJ6uD4ByASKFsAcA3tObYTMW1zKg/rtoBooGwBwDWky8ym1ExSrqguK1ZdopQdiUBEULYA4Bo6+ka0tqZcFSVFYUf5N5rZkQhEBmULAK5h9gbUuXUJMW1zbZU6+kbk7mFHAXADlC0AmEcy6eroHb2yGD3XNNdVaXh8Wn3DE2FHAXADlC0AmMeFoXFdnprJucXxaenjKNpZtwXkPMoWAMzjyk7EHDtjKy1dAjtYtwXkPMoWAMwjXbZydc1WfXWpqkqLOGsLiADKFgDMo6NvRDUVxVpZWRJ2lHmZmTbXVrIjEYgAyhYAzKO9d0TNtVUys7CjXFN6RyKA3EbZAoB5dPSO5Nw9EefaXFelC4PjGpmYDjsKgOugbAHAHJdGJ9U/OpmzOxHT0mXwFLNbQE6jbAHAHOlLc7letpq5RyIQCYsuW2ZWaWbxbIQBgFxwZSdijl9GXL+yQkUxY90WkONuWLbMLGZmv2xmT5hZr6Tjki6Y2VEz+xMza17IQGZWZmYvmNlBM3vZzH5/qeEBIBs6+kZUWhTT2uXlYUe5ruJ4TOtXVjCzBeS4hcxsPSNps6TfkdTg7k3uXifpjZL+VdL/NrP3L+B1JiTd6e67JLVKutfMXn9zsQEge9p7R7SptkrxWO7uREyb3ZHIWVtALlvIrezvdvepuQ+6+0VJX5H0FTMrvtGL+OzdUtO/fhWn/nAHVQA5p71vRLsaa8KOsSCb66r09PFeTc0kVRxnGS6Qi274f2a6aJnZYTP7opn9NzO7z8wazex3r37OjZhZ3MwOSOqV9JS7/3Ce5zxoZm1m1tbX17eofxkAWKrxqRl1Xrqc84vj05prqzSddJ29OBZ2FADXsJhfg94s6a8lXZb0XklHJL19MYO5+4y7t0pqlHSbmd0yz3Mecve97r63trZ2MS8PAEt2qm9U7rm/EzFtMzsSgZy3kMuIkq5cNnw29Udm1iLp925mUHcfMLNnJN2r2dIGADmhvS8aOxHTNtdWShI7EoEctuCZLTPbcvXn7n5S0s5F/P1aM6tJfVwu6R7N7mwEgJzR0TuimEkbV1WGHWVBEmXFqq8uZWYLyGELntmS9FdmtllSl6RDksokHTGzCndfyGKB1ZI+nzqjKybpH9z98UUnBoAsau8bUdOKCpUVR+c4weY6diQCuWwxlxHfKklmtk5S+viGXZIOmFnS3bfd4O8fkrT75qMCQPZF4Z6Ic22urdJjL3XJ3XP6xtlAobph2TIzSx3bIEly97OSzkr6xlXPqc5OPAAIzkzSdeq1Ub1pS7Q25zTXVWlkYlo9QxNqWFYWdhwAcyzoUFMz+43UjNYVZlZiZnea2eclvSs78QAgOJ2XxjQ5nVRzxGa2WuoSkqQTPcMhJwEwn4WUrXslzUh6xMzOp27T86qkk5LeJ+kT7v5wFjMCQCCu3BOxLhqL49O2NaTKVvdQyEkAzOeGlxHdfVzSpyV9OnVS/CpJl919IMvZACBQUbkB9VzLK0tUlyjViW52JAK5aFH3dnD3KXe/kC5aZvb9rKQCgBCc6BlWXaJUNRUlYUdZtK0NCZ3oYWYLyEVLvZHWmoykAIAccKJ7WNtWR3O/z9b6hE72jGgmyS1ngVxzw7JlZp9K3a/wDWaWmPNl/q8GkBemZ5I62TtyZf1T1GxtSGhiOqnT/Zy3BeSahZyzdVjSrZJ+RdItZjaUeuywpGh+VwKAOU73j2pyOhnZsrWtYXZG7pXu4citOQPy3UIWyD909edm1qjZ8rVT0pNZygUAgTp2YfbYhK0RLVvNdVUyk453D+u+W1eHHQfAVRZzux5Jkrt3SuqU9E+ZjwMA4TjRPax4zNRcF81ZofKSuDasrNSJbs7aAnLNUhfIA0BeON49pE2rKlVaFJ17Is61tT6hVzjYFMg5lC0A0Ozlt6juREzb0pDQ6f5RjU/NhB0FwFUoWwAK3vD4lDovXY7s4vi0bQ0JJV062cPhpkAuoWwBKHjpS29b66NdttKL+49z2x4gp1C2ABS846lF5dtWR7tsbVhZqZKiGOu2gBxD2QJQ8I5fGFaitEhra8rDjrIk8Zippa7qSnkEkBsoWwAK3onuYW1tSMjMwo6yZFsbEhz/AOQYyhaAgubuOtY9FNnDTOfa1pBQ7/CELo1Ohh0FQAplC0BBuzA4ruHx6cgf+5C2NXXbnhOs2wJyBmULQEFL79yL+rEPaekdlVxKBHIHZQtAQUsvJs+Xy4j11aVaVl7MInkgh1C2ABS0E93DWltTruqy4rCjZISZaWsDt+0BcgllC0BBO35hOG8uIaZtrU/ole5huXvYUQCIsgWggE1OJ9XRN5I3lxDTtjYkNDwxra6By2FHASDKFoAC1tE3oumk581OxLT0TB2L5IHcEFjZMrMmM3vGzI6a2ctm9tGgxgaA+aTLSL5dRtySLlus2wJyQlGAY01L+i13f8nMEpL2mdlT7n40wAwAcMWx7iEVx00bV1WGHSWjqsuKtWZZGTNbQI4IbGbL3S+4+0upj4clHZO0NqjxAWCuE93Daq5LqDiefysquG0PkDtC+Q5jZhsk7Zb0wzDGBwApP3cipm1tqFZH34imZpJhRwEKXuBly8yqJH1F0m+6+9A8X3/QzNrMrK2vry/oeAAKxODYlLqHxvO2bG1rSGhqxvXqa6NhRwEKXqBly8yKNVu0vujuj833HHd/yN33uvve2traIOMBKCDp2/Tk27EPaTvWzO6wPNw5GHISAEHuRjRJfyPpmLv/WVDjAsB80rez2Z5nxz6kba6tUlVpkfafuxR2FKDgBTmzdbukD0i608wOpP68PcDxAeCK491DqqkoVl2iNOwoWRGPmXY1LdOBcwNhRwEKXmBHP7j79yRZUOMBwPXsPzugnY01mp10z0+tTTX6q38+pcuTMyoviYcdByhY+bffGQBuYHh8Sid6hrVnXU3YUbKqtWm5ppOuI+dZtwWEibIFoOAcPDcod2nPuuVhR8mq1qYaSdKBswOh5gAKHWULQMF56ewlmUmteT6zVZsoVePychbJAyGjbAEoOC+dvaSWuipVlxWHHSXrWptqmNkCQkbZAlBQkknX/rMDeX8JMW33uuU6PziunqHxsKMABYuyBaCgvNo/qsHLUwVTttLrtvYzuwWEhrIFoKC8dGZ2/dKe9TXhBgnIT6ypVnHcOG8LCBFlC0BBeensgKrLirRpVVXYUQJRVhzXjtXVOsAieSA0lC0ABWX/2UtqXbdcsVj+HmY6V2tTjQ51Dmom6WFHAQoSZQtAwSiUw0zn2r1uucYmZ/RKz3DYUYCCRNkCUDAK5TDTuVgkD4SLsgWgYBTKYaZzrV9ZoeUVxazbAkJC2QJQMArpMNOrmdns4absSARCQdkCUBAK7TDTuVqblutk74iGx6fCjgIUHMoWgIJQaIeZzrV7XY3cpUOdg2FHAQoOZQtAQSi0w0zn2pVaJM+lRCB4lC0ABaHQDjOda1l5sTbXVmr/WRbJA0GjbAEoCIV4mOlcrU3LdeDcgNw53BQIEmULQN4r1MNM52pdV6PXRibVeely2FGAgkLZApD3CvUw07l2p9Zt7TvDpUQgSJQtAHmvUA8znWv76mqtqCzRP7/SF3YUoKBQtgDkvbYzhXmY6VzxmOktW2v1zIlebkoNBIiyBSCvjU5M619P9euOltqwo+SEO7fVaWBsil2JQIAoWwDy2nOv9GlyOql7dtSHHSUn3NFSq6KY6bvHe8OOAhQMyhaAvPbUsR7VVBRr7/rCXhyftqy8WD+1YYWePkbZAoJC2QKQt6Znknr6eK/u3Fqnojjf7tLu2l6nEz3D6rw0FnYUoCAE9t3HzD5nZr1mdiSoMQEUtrYzlzQwNsUlxDnu3FYnSXqGS4lAIIL8Ve9hSfcGOB6AAvfU0R6VxGN60xYWx19tU22VNq6qZN0WEJDAypa7PyfpYlDjAShs7q6njvbop5tXqrK0KOw4OefObXV6vqNfY5PTYUcB8l7OLWIwswfNrM3M2vr6OHgPwM052TuisxfHuIR4DXdtq9PkdFLfb+8POwqQ93KubLn7Q+6+19331tYy9Q/g5jx1tEeSdPd2ytZ89m5YoURpkZ4+3hN2FCDv5VzZAoBM+PbRHu1qqlF9dVnYUXJSSdHsWrbvHuuVO6fJA9lE2QKQd3qGxnXw3IDu2V4XdpScdue2OvUOT+jl80NhRwHyWpBHPzwi6QeStppZp5k9ENTYAArLd47NXhq7Z0dDyEly21u21spM+i4HnAJZFeRuxPe5+2p3L3b3Rnf/m6DGBlBYvnO0R+tWVGhLfVXYUXLayqpS7W6qYd0WkGVcRgSQV0YnpvX9jn7ds6NeZhZ2nJx31/Z6HewcVO/weNhRgLxF2QKQV7jx9OKkT5PnXolA9lC2AOSVp45y4+nF2NaQ0KbaSv2/F86yKxHIEsoWgLwxeHlK3z7ao7u21XPj6QUyMz3wxo061DmoF17lJh9ANvDdCEDe+MK/ntHIxLQ+ePuGsKNEyrt2N2p5RbE++71Xw44C5CXKFoC8MD41o7/9/mnd0bJKt6xdFnacSCkviev9r1+v7xzr0auvjYYdB8g7lC0AeeHL+zr12siEPvSWzWFHiaQPvGG9imMxfY7ZLSDjKFsAIm96JqmHnjulXU01esOmlWHHiaS6RJnub12jR/ed08DYZNhxgLxC2QIQef90pFtnL47pQ2/ezNlaS/Brd2zS+FRSX/zh2bCjAHmFsgUg0txdn3m2Q5tqK/U2ztZakq0NCd3RskoPP39aE9MzYccB8gZlC0CkPXfyNR29MKRff9NmxWLMai3Vf7xjk/qGJ/T4wQthRwHyBmULQKT95bMdaqgu0/2714QdJS/c0bJKW+qr9Nf/copDToEMoWwBiKz9Zy/pB6f69Wt3bFRpUTzsOHnBzPRrb9yk493Der6jP+w4QF6gbAGIrL/85w5VlxXpvbetCztKXrl/9xrVJkr1h08cY+0WkAGULQCR9M3DF/Tkyz361ds3qqq0KOw4eaW0KK4/fuetOnphSH/8zeNhxwEij7IFIHI6+kb0sS8f0q6mGn34rRximg1376jXB2/foIefP62njvaEHQeINMoWgEgZnZjWr//dPpUUxfSZX9nDWq0s+vh923TL2mr99pcP6vzA5bDjAJFF2QIQGe6ujz92WB19I/rz9+7WmprysCPltdKiuD71vj2amk7qo1/ar+mZZNiRgEiibAGIjM8/f1rfOHhev/W2rXpjy6qw4xSEjasq9YfvvFUvnr6kP//uybDjAJFE2QIQCfvOXNT/fOKY7t5epw+9mXVaQfr53Wv1iz/ZqE89067vt78WdhwgcihbAHLewXMD+tAXXtLa5eX6P+9p5aT4EPzB/T+hTasq9cDnX9SjbefCjgNECmULQM6aSbr+4pl2/cJnnlc8ZvqrD/yklpUXhx2rIFWUFOlLD75Be9Yt129/+ZB++9GDujzJGVzAQnA4DYCc1DVwWf/l7w/ohVcv6mduXa0/euetWlZB0QpTbaJUf/fA6/TJ77yiTz3TrkOdg/r0+/doc21V2NGAnMbMFoCc4u76xsHzuu8Tz+nlrkH96bt36f/+8m6KVo6Ix0z/9W1b9fAHb1PfyITe8anv6ZEXzmp8ilku4Fosl280unfvXm9raws7BoAAnOkf1Vf3d+kf93fpdP+YWptq9IlfatWGVZVhR8M1XBi8rP/8yH69ePqSlpUX65271+o9e5u0Y0112NGAUJjZPnff+28eD7Jsmdm9kj4pKS7ps+7+v673fMoWkL8mpmd0qm9Ubacv6h8PnNe+M5dkJr1+40q9a89a/fzutSqOM/me65JJ1w9O9evvXzynbx3p1uRMUjsbl+n+1rXa1bhM21ZXczslFIzQy5aZxSW9IukeSZ2SXpT0Pnc/eq2/Q9kCosXdNTXjmpxJamxiWpfGpnRpbFIDY5O6NDal7sFxnewd1onuYZ3uH9NMcvb7z5b6Kr1zd6Pub13DQaURdml0Uv94oEt//+I5He8evvL4uhUV2r46oS31CdUmSrW8okQrK0u0oqpEyytKVFYcV1lxTCXxmMzYaYroulbZCvLXjdsktbv7qVSgL0m6X9I1y1a2/ae/a9OpvtGwhgeybqG/Sl39S5fP+cBTX5/9p5R0VzLpmnHXTDJdsJKamE5qciap6/3+ZiatX1GhLfUJ3XfLam1pSGjH6oQ211bxQzYPLK8s0Qdv36hf/ekNujA4rmMXhlJ/hnXswpC+fbTnuv99SFJpUUylRTEVxWOKmSkek+JmisVMMTOZSSbJzGTS7Ccp6Q8X898S/9UVhtLimB7/jTtCGz/IsrVW0tWHs3RKet3cJ5nZg5IelKR169ZlNdD6lZWKc14P8pwt9MfJNX5ozf5gu+oHnGn2h2DqB2A8JhXFYld+SJYWx1USj6m8JK7lFSVaXlGsmooSLa8svjKLgfxmZlpTU641NeW6a3v9lcenZ5IauDyli6OT6h+Z1KWx2T/jU0mNT81oYjqpiekZTUwlNZ1MaiapK8U+mXQlryr96V8C0ub+krAQvpgnI9LCXpKQcxfS3f0hSQ9Js5cRsznWf3/79my+PADgKkXxmFZVlWpVValUf+PnA/kiyKrXJanpqs8bU48BAADkrSDL1ouSWsxso5mVSHqvpK8HOD4AAEDgAruM6O7TZvYRSU9q9uiHz7n7y0GNDwAAEIZA12y5+zclfTPIMQEAAMLEiYEAAABZRNkCAADIIsoWAABAFlG2AAAAsijQG1Evlpn1STqT5WFWSXoty2MUEt7PzOM9zTze08zi/cw83tPMCur9XO/utXMfzOmyFQQza5vvppG4Obyfmcd7mnm8p5nF+5l5vKeZFfb7yWVEAACALKJsAQAAZBFlK3XTa2QM72fm8Z5mHu9pZvF+Zh7vaWaF+n4W/JotAACAbGJmCwAAIIsoWwAAAFlE2ZJkZu82s5fNLGlmbLW9SWZ2r5mdMLN2M/t42Hmizsw+Z2a9ZnYk7Cz5wMyazOwZMzua+v/9o2FnijozKzOzF8zsYOo9/f2wM+UDM4ub2X4zezzsLPnAzE6b2WEzO2BmbWFkoGzNOiLpXZKeCztIVJlZXNJfSLpP0g5J7zOzHeGmiryHJd0bdog8Mi3pt9x9h6TXS/ow/40u2YSkO919l6RWSfea2evDjZQXPirpWNgh8sxb3b01rLO2KFuS3P2Yu58IO0fE3Sap3d1PufukpC9Juj/kTJHm7s9Juhh2jnzh7hfc/aXUx8Oa/WG2NtxU0eazRlKfFqf+sOtqCcysUdLPSPps2FmQOZQtZMpaSeeu+rxT/CBDjjKzDZJ2S/phyFEiL3XJ64CkXklPuTvv6dJ8QtLHJCVDzpFPXNK3zWyfmT0YRoCiMAYNg5l9R1LDPF/6XXf/WtB5AITDzKokfUXSb7r7UNh5os7dZyS1mlmNpK+a2S3uzjrDm2BmPyup1933mdlbQo6TT97o7l1mVifpKTM7nrpyEJiCKVvufnfYGfJcl6Smqz5vTD0G5AwzK9Zs0fqiuz8Wdp584u4DZvaMZtcZUrZuzu2S3mFmb5dUJqnazL7g7u8POVekuXtX6p+9ZvZVzS57CbRscRkRmfKipBYz22hmJZLeK+nrIWcCrjAzk/Q3ko65+5+FnScfmFltakZLZlYu6R5Jx0MNFWHu/jvu3ujuGzT7PfRpitbSmFmlmSXSH0t6m0L4ZYCyJcnM3mlmnZLeIOkJM3sy7ExR4+7Tkj4i6UnNLjz+B3d/OdxU0WZmj0j6gaStZtZpZg+EnSnibpf0AUl3praAH0jNIODmrZb0jJkd0uwvXE+5O8cVIJfUS/qemR2U9IKkJ9z9W0GH4HY9AAAAWcTMFgAAQBZRtgAAALKIsgUAAJBFlC0AAIAsomwBAABkEWULAAAgiyhbAAAAWUTZAlAwzOxZM9uW+nilmXFbGQBZR9kCUEiaJb2S+ninpMMhZgFQIChbAAqCma2X1OXuydRDOyUdCjESgAJB2QJQKHbpx8vVT4qyBSAAlC0AhaJVUpkkmVmLpPvFZUQAAaBsASgUuyTFzOygpP8h6aikfx9uJACFwNw97AwAkHVmdlLSHncfDjsLgMLCzBaAvGdmCUlO0QIQBma2AAAAsoiZLQAAgCyibAEAAGQRZQsAACCLKFsAAABZRNkCAADIIsoWAABAFlG2AAAAsuj/Aw2GO2E/q1XhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mus = np.linspace(-1, 5, 100)\n",
    "l_mu = [likelihood_norm(mu=m, data=sample) for m in mus]\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.plot(mus, l_mu);\n",
    "ax.set_xlabel(\"$\\mu$\");\n",
    "ax.set_ylabel(\"$L(\\mu)$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e839b56e-d062-4032-81ed-4c3b6fc9662e",
   "metadata": {},
   "source": [
    "## 事前分布\n",
    "\n",
    "$$\n",
    "P(\\mu) = \\mathcal{N}\\left(\\mu | \\lambda, \\tau^2 \\right) = \\frac{1}{\\sqrt{2 \\pi}\\tau} \\exp \\left[ -\\frac{(\\mu-\\lambda)^2}{2 \\tau^2} \\right] \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8d0224e8-6781-41c9-aa5e-4d11887d030a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAEICAYAAAAwSKKsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwaElEQVR4nO3deXxV1b3//9cnOZnJQEISIIQZmWUKMyKKA9pWxKoFax1qaxVxuLW3tfd+e2/rr7fW3qrVK7TifK11nlCrOIATIhAmmSECkoQhCQkJmaf1+yMH7hFBGZKzc5L38/E4j+yz99onn32K5c3ae61lzjlEREREJLSEeV2AiIiIiJw4hTgRERGREKQQJyIiIhKCFOJEREREQpBCnIiIiEgI8nldgBc6derkevbs6XUZIiIiIt9q5cqVRc651CP3t8sQ17NnT7Kzs70uQ0RERORbmdmXR9uv26kiIiIiISioIc7MppnZFjPLMbM7jnI8ysye8x9fZmY9jzje3czKzewXx/uZIiIiIm1R0EKcmYUDc4ELgEHALDMbdESz64AS51xf4D7g7iOO3wu8dYKfKSIiItLmBLMnbgyQ45zb7pyrBZ4Fph/RZjrwpH/7RWCqmRmAmV0M7AA2nOBnioiIiLQ5wQxxGUBuwPs8/76jtnHO1QOlQIqZdQB+BfzuJD4TADO73syyzSy7sLDwpC9CREREpDUIlYENvwXuc86Vn+wHOOfmO+eynHNZqalfG6UrIiIiElKCOcVIPpAZ8L6bf9/R2uSZmQ9IBPYDY4FLzexPQBLQaGbVwMrj+EwRERGRNieYIW4F0M/MetEUtGYCVxzRZgFwNbAUuBRY5JxzwBmHGpjZb4Fy59yD/qD3bZ8pQeaco6qugYqaBipr66mqa6C6rpHqugb/q5Hahkbq6hupa2jarq1vpKHR0eAcjY2O+samn8cSHhZGeBiEhRnhZoSHGVG+MCLC/S9fGJHhYURHhBEdEe5/hRETEU5spI8OUT6iI8LwP3IpIiIScoIW4pxz9WY2B1gIhAOPOec2mNmdQLZzbgHwKPCUmeUAxTSFshP+zBa9kHaisdFRWlXH/opaiitqOVBZy4HKOg5U1VJSWUdpVR1lVXUcrK6nrLrp58HqOipqGqiorccdO3+dkKNlrOb67DCDuEgfsVHhxEdHkBDtIz46gvhoHwkxESTFRJAUG0FSbCRJMRF0jIukY2wknTpEkhAdQViYAqCIiHjHXHP9jRhCsrKyXHtdsaG6roF9ZdXsKa1mX1k1hQdrDr8KDtZQVF5zOLg1HKMnzBdmJMVGkBAQeOKjfcRHRRAX5SMuKrzpZ2Q4MZE+YiObesGifeFE+XvEonxhRIaHE+Gzw71nvrCmHrUwM3xhdsyQ5Jyj0UFDo6PRORoaHfUNjrpGf8+ev4evpr6R6rpGauoaDvcGVtU1UFVbT7m/l7CipoGKmnoO1tRRVtUURMuq6ymragqq9d/wHXSMiyQlLpLU+KjDr7T4aFLjo0iPj6JLYgzpiVFE+cKb7X8/ERFpf8xspXMu68j97XLZrbbKOUfhwRryDlSRX1LF7gNV5Pu395RWs7esmuKK2q+dF+kLI7VDFGkJUWQmxzKiexLJcZGkxEWR0iGSZH8PVKK/NyouMtzT25BmRrhBeAv3hDnnKK+pb+qBrKyjpLKWkspaisprKa6oYX9503ZReQ3bCysoPFhDbUPj1z4nJS6SzonRdEmMJiMphq5JMWR0jCEjKYZuHWPp1CFSt3VFROSEKcSFmLqGRnKLK9m5v4KdRZXsKq4kt9j/s6SS6rqvhoiEaB8ZHWPpmhjNiO5JdEmMpnNiDJ0ToklPaOo5SojxKUQchZn5b69GkJn87e2da7oFXXCwhn1l1ewtbXrt8W/nlVSxbEcxB6vrv3JeTEQ4mckxdE+OJTM5lu7JsfTsFEevlDi6dYzBFx4qg8hFRCSYFOJaIeccReW1fFFY3vQqqGB7UTk7iirIK6n6ym3OuMhwuqfE0atTHGeelkpmcizdOv5fT098dISHV9K+mFnT83OxkZyWHn/MdmXVdeSXNPWQ5pVUkltSdTiMf/rFfiprGw639YUZmcmx9EyJpU9qB/qkdWj6mRpHcpx68ERE2jOFOA8duv25dV85W/cdZFvBQbbuK2fbvoOUBfTWxESE0zs1jiEZiXzv9K5NvTSdYumZor/IQ1FCdAQJXSIY2CXha8ecc+yvqGVnUQU7iioO97huL6pg6fb9X+lpTYyJoF9aB/qlx3NaegdOS4+nX3oHUjtE6c+EiEg7oIENQVJd18C2feVs2lvG5j0H2bSnjM17yyiprDvcpmNsBP3S4+mX1oG+/h6Xvmkd6JwQrZGQQmOjY3dpFV8UVvBFQTk5heXk7Ctna8FBDgT8OUqJi2RAl3gGdE5gQOd4BnZJoF96Bw2wEBEJURrYEESlVXVs3F3Ght2l/p9l5BSWH74NGhMRTv/O8Uwb0pn+6fH+HpR4PeAu3ygszOjWMZZuHWM587T/W3XEOUdheQ3b9pWzZe9Btuw9yOa9ZTy97MvDPXe+MKNvWgcGd01kcNcEBndNYFDXBN1uFxEJYQpxLWDG3CVsL6oAIC0+iiEZiZw3OJ2BXRIY2CWBHsmx6lmTZmNmpMVHkxYfzcS+nQ7vb2h0fLm/gk17DrJxTykbdpfx4dZCXlqV5z8PeneK4/RuSQzNSOT0bokM7ppITKR67EREQoFup7aAt9fvJToijMFdE0mNj2qx3yNyMgrKqtmwu4x1+aV8nlfKuvwD7CurAZqmbemfHs/w7kmMyExiRPckenfqoH90iIh46Fi3UxXiRIR9ZdV8nlfK2twDrMk9wNrcAxysaRpcEx/tY2T3jozq0ZGsHh0ZlplEXJQ68UVEgkUhLoBCnMg3a2x0bC8qZ/WuA6zadYBVX5awteAgzjUtVzawSwKjeyYzplcyo3smq8dZRKQFKcQFUIgTOXGlVXWs3lXCqi9LWLGzhNW5JYcHTvTqFMeYnsmM65PMuN4pdEmM8bhaEZG2QyEugEKcyKmrrW9k/e5SVuwoZsXOYlbsLKG0qmmqkx4psYzvncK43ilM6JNCWkK0x9WKiIQuhbgACnEiza+x0bFpbxmfbS9m6Rf7Wb5j/+FJq/umdWBS305M6JPCuD4pJGhqExGR46YQF0AhTqTlNTQ6Nu0pY0lOEUv8oa66rpEwg+GZSZzRL5XJp3ViWLckrQ8rIvINFOICKMSJBF9NfQOrdx1gSU4RH28rYm3eAZxrGv06sU8nzuyfypT+qXqeTkTkCApxARTiRLx3oLKWJTn7+XhbIR9tLWR3aTUA/dPjmdI/lSn908jq2ZEI9dKJSDunEBdAIU6kdXHOkVNQzgdbClm8pYAVO4upa3DER/mYfFoqUwemMaV/GslxkV6XKiISdApxARTiRFq38pp6Ps0pYvGWAt7fVEDBwRrMYGT3jpw9II3zB6fTJ7WD1hoWkXZBIS6AQpxI6GhsdGzYXcZ7m/bx/uZ9rM8vA5rmpjt3UDrnDkpnZPeOhGtpMBFpoxTiAijEiYSuPaVVvLdxH+9s3Mdn2/dT1+BIiYvkvMHpnD+4MxP6dCLSp+foRKTtUIgLoBAn0jaUVdfx4ZZCFm7Yy+LNBVTUNhAf7eOcgU2Bbkr/VKIjwr0uU0TklCjEBVCIE2l7qusaWJJTxNvr9/Lupn0cqKwjLjKcswem852hnZnSP02BTkRC0rFCnM+LYkREmlt0RDhTB6YzdWA69Q2NfLa9mDfX7WHhhr28vnY3sZHhnD0gje8N68qZp6mHTkRCn3riRKRNOzLQFVfUEh/l4/whnfnesK5M6JOiuehEpFXT7dQACnEi7VN9QyOffrGfBWt3s3D9Xg7W1JMcF8l3hnbh4hFdGdm9o6YtEZFWRyEugEKciNTUN/DhlkIWrN3Ne5v2UV3XSGZyDNOHZXDxiK70TYv3ukQREUAh7isU4kQkUHlNPQvX7+XVNfksySmi0cHQjEQuGZnBRcO6ktIhyusSRaQdU4gLoBAnIsdScLCa19fu4ZXVeazPL8MXZkzpn8olI7sxdWAaUT4NiBCR4FKIC6AQJyLHY8veg7y8Oo9XV+ezr6yGxJgIpg/vymWjMhmSkaDn50QkKBTiAijEiciJaGh0fJJTxIsr81i4YS+19Y0M6BzPpaO6cfGIDDrpdquItCCFuAAKcSJyskor61jw+W5ezM5lbV4pEeHGuYPSuTwrkzP6pWoNVxFpdgpxARTiRKQ5bN13kOdX5PLSqjxKKuvomhjNZVmZXD46k4ykGK/LE5E2QiEugEKciDSnmvoG3ttYwLMrdvFJThEAU05L5YqxPTirfyo+TSYsIqegVYQ4M5sG3A+EA4845/54xPEo4H+BUcB+4AfOuZ1mNgaYf6gZ8Fvn3Cv+c24Ffurf/7Bz7i/fVodCnIi0lNziSp7PzuW5FbkUHKyhc0I0l4/OZOboTLqqd05EToLnIc7MwoGtwLlAHrACmOWc2xjQZjZwunPuBjObCcxwzv3AzGKBWudcvZl1AdYCXYEBwLPAGKAWeBu4wTmX8021KMSJSEura2jk/U0FPLN8Fx9tK8SAqQPTuXJcD87o24kwPTsnIsfpWCHOF8QaxgA5zrnt/oKeBaYDGwPaTAd+699+EXjQzMw5VxnQJho4lDwHAssOHTezD4FLgD+11EWIiByPiPAwpg3pzLQhncktruSZ5bt4bkUu727cR4+UWH44tjuXjcqkY1yk16WKSIgK5oMaGUBuwPs8/76jtnHO1QOlQAqAmY01sw3AOpp62+qB9cAZZpbi7627EMg82i83s+vNLNvMsgsLC5vxskREvllmciy/nDaAT399NvfPHE5afBR/+Odmxt71Pv/6wlrW55d6XaKIhKBg9sSdEufcMmCwmQ0EnjSzt5xzm8zsbuAdoAJYAzQc4/z5+J+ry8rKan+jOUTEc1G+cKYPz2D68Aw27y3jqaVf8vKqfF5YmceoHh25ekJPpg3uTKRPAyFE5NsF8/8p8vlqL1k3/76jtjEzH5BI0wCHw5xzm4ByYIj//aPOuVHOuclACU3P3YmItGoDOifwXzOG8tm/TeU33x1EUXkNtzyzmkl3L+KB97dRVF7jdYki0soFM8StAPqZWS8ziwRmAguOaLMAuNq/fSmwyDnn/Of4AMysB00DGnb636f5f3an6Xm4f7T0hYiINJfEmAium9SLxbdP4fFrRjOgSwL3vruVCX9cxC9eWMuG3brVKiJHF7Tbqf6RpXOAhTRNMfKYc26Dmd0JZDvnFgCPAk+ZWQ5QTFPQA5gE3GFmdUAjMNs5V+Q/9pKZpQB1wE3OuQPBuiYRkeYSFmacNSCNswakkVNQzhOf7uCllfm8uDKPsb2SuW5SL6YOTNeKECJymCb7FRFppUor63guexdPfvol+Qeq6JkSy3WTevH9Ud2IjQyZR5pF5BR5Pk9ca6IQJyKhpL6hkbc37OXhj3ewNvcAiTER/HBsd66e0JP0hGivyxORFqYQF0AhTkRCkXOOVbtKePijHSzcuBdfmHHx8Ayun9ybfunxXpcnIi2kNUz2KyIip8DMGNUjmVE/SubL/RU8+skOns/O5YWVeUwdkMb1k3szplcyZnpuTqQ9UE+ciEgIK66o5X+X7uR/l35JcUUtwzOTuHFKH84dmK6lvUTaCN1ODaAQJyJtTVVtAy+uzGX+x9vJLa6iT2ocN5zZh+nDMzR5sEiIU4gLoBAnIm1VfUMjb67bw98+3M6mPWV0SYzmJ2f0ZtaYTI1oFQlRCnEBFOJEpK1zzvHh1kLmffAFy3cUkxwXyY8n9uRH43uSGBPhdXkicgIU4gIoxIlIe5K9s5i5i3NYvKWQ+CgfV03owY8n9iKlQ5TXpYnIcVCIC6AQJyLt0fr8UuZ9kMNb6/cS7QvnynHd+enk3qTFa645kdZMIS6AQpyItGc5BeXMW5zDq2vyiQgPY9aY7txwZh86JyrMibRGCnEBFOJERGBnUQXzPsjh5VX5hJlx+ehuzJ7Sl65JMV6XJiIBFOICKMSJiPyf3OJK5n3wBS+uzAXgB6MzFeZEWhGFuAAKcSIiX5dX0hTmXsjOxTB+MDqTG6f0UZgT8ZhCXACFOBGRYzsyzM0ak8lNZ/UlLUHPzIl4QSEugEKciMi3yyup5MFFObywMg9fmHHV+B7ccGYfTU0iEmQKcQEU4kREjt+X+yu4//1tvLo6n+iIcK6Z0JOfTe5DYqwmDRYJBoW4AApxIiInLqegnPvf38Ybn++mQ5SPn03uzbUTexEXpeW8RFqSQlwAhTgRkZO3aU8Z97yzlfc27SMlLpLZZ/Xlh2O7Ex0R7nVpIm2SQlwAhTgRkVO3alcJf164hU+/2E+XxGhuO6cf3x/ZDV94mNelibQpxwpx+i9NREROysjuHfnHT8fx9E/GkpYQza9eWsf5f/mIt9fvoT12EIgEm0KciIickol9O/Hq7An87cpRANzw91XMmPcpS7/Y73FlIm2bQpyIiJwyM2PakM4svG0yf/r+6ewrq2bWw59x9WPL2bSnzOvyRNokPRMnIiLNrrqugf9dupMHF+VwsKaeS0Z04/bzTtPqDyInQQMbAijEiYgEx4HKWuZ98AVPfLoTgGsn9mT2lL4kxmiOOZHjpRAXQCFORCS48g9Ucc87W3hldT6JMRHccnY/rhzXg0ifnuoR+TYanSoiIp7JSIrh3suH8+bNZzCkayJ3vrGRc+/7kLfWaSSryMlSiBMRkaAZ1DWBp64bwxPXjibKF8aNT6/i0r8tZdWuEq9LEwk5CnEiIhJUZsaU/mn885YzuOuSoXy5v5JL5n3KLc+sJv9AldfliYQMhTgREfGELzyMWWO68+G/TuHms/uycMNezv7zB/x54RYqauq9Lk+k1VOIExERT8VF+bj9vP4s+sUULhjSmQcX5zDlzx/w/IpcGhv1vJzIsSjEiYhIq5CRFMNfZo7gldkT6J4cyy9f+pyL5n7Cip3FXpcm0iopxImISKsyontHXrxhPA/MGsH+8lou+9tS5vxjlZ6XEzmCQpyIiLQ6ZsZFw7qy6PYp3Dq1H+9t2sfZf/6Ae9/dSlVtg9flibQKQQ1xZjbNzLaYWY6Z3XGU41Fm9pz/+DIz6+nfP8bM1vhfa81sRsA5/2JmG8xsvZk9Y2bRQbwkERFpQTGR4fzLuafx/u1TOG9wZx54fxtT7/mANz/X/HIiQQtxZhYOzAUuAAYBs8xs0BHNrgNKnHN9gfuAu/371wNZzrnhwDTgITPzmVkGcIv/2BAgHJjZ4hcjIiJBlZEUw//MGsHzPxtPYmwkN/1jFbMe/ozNe8u8Lk3EM8HsiRsD5DjntjvnaoFngelHtJkOPOnffhGYambmnKt0zh0abx4NBP7zywfEmJkPiAV2t9gViIiIp8b0SuaNmyfx+4uHsHnvQS68/2P+87X1lFbWeV2aSNAFM8RlALkB7/P8+47axh/aSoEUADMba2YbgHXADc65eudcPvBnYBewByh1zr1ztF9uZtebWbaZZRcWFjbjZYmISDCFhxlXjuvBB7+YwpXjevDUZ19y9j2akkTan5AZ2OCcW+acGwyMBn5tZtFm1pGm3rteQFcgzsyuPMb5851zWc65rNTU1OAVLiIiLSIpNpI7pw/hjZvPoHdqHL986XMu+eunfJ53wOvSRIIimCEuH8gMeN/Nv++obfy3RxOB/YENnHObgHJgCHAOsMM5V+icqwNeBia0SPUiItIqDeqawPM/G8+9lw8jr6SK6XOX8OuX11FSUet1aSItKpghbgXQz8x6mVkkTQMQFhzRZgFwtX/7UmCRc875z/EBmFkPYACwk6bbqOPMLNbMDJgKbGr5SxERkdbEzLhkZDcW/eJMfjyxF89n53L2PR/w3IpdusUqbVbQQpz/Gbc5wEKagtbzzrkNZnanmV3kb/YokGJmOcDPgUPTkEwC1prZGuAVYLZzrsg5t4ymARCraHpWLgyYH6xrEhGR1iUhOoLffHcQb94yib5pHfjVS+u49G+fsmF3qdeliTQ7a4/z7GRlZbns7GyvyxARkRbknOOlVfnc9c9NlFTWctX4ntx+3mnER0d4XZrICTGzlc65rCP3h8zABhERkRNhZlw6qhuLbp/CFWO78+TSnUy950NNFCxthkKciIi0aYmxEfz+4qG8OnsiqfFR3PSPVVzz+Aq+3F/hdWkip0QhTkRE2oVhmUm8dtNE/vN7g1j5ZQnn3fcRDy7aRk291mKV0KQQJyIi7YYvPIxrJ/bivZ+fyTkD0/nzO1v5zgOfsHxHsdeliZwwhTgREWl3OidGM/eHI3n8mtFU1TZw+UNLueOlzzlQqbnlJHSccIgzszj/YvYiIiIh7awBabz788n8bHJvXliZxzn3fshra/I18EFCwreGODMLM7MrzOxNMysANgN7zGyjmf23mfVt+TJFRERaRmykj19fOJDX50wio2Mstz67hqseW05ucaXXpYl8o+PpiVsM9AF+DXR2zmU659JomoD3M+DuY61XKiIiEioGdU3g5Rsn8LuLBrPKP/Dh4Y+2U9/Q6HVpIkf1rZP9mlmEf13SU2rTmmiyXxER+Sa7D1Txm1fX8/7mAoZmJHLXJUMZkpHodVnSTp30ZL/HE85CKcCJiIh8m65JMTxydRYPXjGCPaVVTJ+7hLve2kR1naYjkdbjuAc2mNk6M3vazH5lZheYWTcz+/eWLE5ERMQrZsZ3T+/Kez8/k++PzOChD7dzwf0f89n2/V6XJgKc2OjUM4GHgSpgJrAeuLAlihIREWktkmIj+dOlw/j7dWOpb2xk5vzP+PdX1nGwWjehxFvHHeKcc8XOuQ+ccw84564GRgM5LVeaiIhI6zGpXycW3jaZn0zqxTPLd3HefR+xaPM+r8uSduxEbqeeFvjeObcNOL3ZKxIREWmlYiN9/L/vDuLl2RNJiI7gx09k8y/PraGkQpMES/CdyO3Uh8xsl5ktNbOHzOxJYL2ZxbVUcSIiIq3R8MwkXr95ErdO7cfra3dz7n0f8fb6PV6XJe3M8Uz2awDOubOcc92BHwBv0HQrNQZYZWZbWrRKERGRVibSF8a/nHsaC+ZMIj0hihv+voqbnl5FUXmN16VJO3Fck/2a2c1m1h3AObfLOfc6cDfwV2AZMLcFaxQREWm1BnVN4NWbJvKv5/fn3Y37OPfeD3l97W4t3SUt7nhC3DSgAXjGzA4tt7Ud2EbTKNV7nXMPtGSRIiIirVlEeBg3ndWXN2+ZRPeUOG5+ZjWz1SsnLexbV2z4SmOzCKATUOWcO9BSRbU0rdggIiItpb6hkYc/3sF9724lLiqcO6cP4bund8H/dJLICTvpFRvM7GozKzKzYuARoDyUA5yIiEhL8oWHceOUPk29csmx6pWTFnM8t1N/A5wLDAB2AX9o0YpERETagH7p8bx04wR+Oa0/728q4HyNYJVmdjwhrsw5t9o5V+Cc+w0wpqWLEhERaQt84WHMntKX12+eRJekaG74+ypufXY1Byo1r5ycuuMJcV3M7Hozm2xmqUBESxclIiLSlvTvHM8rsydy2zn9ePPzPZx330cs3lzgdVkS4o4nxP0nMBT4/4AtwBAz+6eZ3WVms1q0OhERkTYiIjyM2845jVdvmkhSbATXPrGCO176nPKaeq9LkxB1QqNTAcysG02h7nRgiHPuRy1RWEvS6FQREfFSTX0D9727jYc++oKMpBjuuWwYY3uneF2WtFLHGp16wiGuLVCIExGR1mDFzmJuf34tuSWV/PSM3vz83NOIjgj3uixpZU56ihERERFpGaN7JvPWrWcwa0x35n+0nYse/IT1+aVelyUhQiFORETEQ3FRPv4wYyiPXzuaA5V1zJi3hHkf5NDQ2P7ulMmJUYgTERFpBc7qn8bC2yZz7qB0/vT2FmbOX0pucaXXZUkrphAnIiLSSnSMi2TuFSO59/JhbN5zkGl/+Yjns3Npj8+vy7dTiBMREWlFzIxLRnbjrdvOYEhGIr988XNu+PtKSio0QbB8lUKciIhIK9StYyzP/HQc/3bhABZtLuD8v3zEx9sKvS5LWhGFOBERkVYqLMy4fnIfXr1pIgkxEfzo0eXc+fpGqusavC5NWoGghjgzm2ZmW8wsx8zuOMrxKDN7zn98mZn19O8fY2Zr/K+1ZjbDv79/wP41ZlZmZrcF85pERERa2uCuibxx8ySumdCTx5bsYPqDS9i8t8zrssRjQQtxZhYOzAUuAAYBs8xs0BHNrgNKnHN9gfuAu/371wNZzrnhwDTgITPzOee2OOeG+/ePAiqBV1r8YkRERIIsOiKc3140mMevHc3+ilouenAJjy/ZoUEP7Vgwe+LGADnOue3OuVrgWWD6EW2mA0/6t18EppqZOecqnXOHFpeLBo72J3Yq8IVz7ssWqF1ERKRVOKt/Gm/fdgaT+nbid69v5MdPrKCovMbrssQDwQxxGUBuwPs8/76jtvGHtlIgBcDMxprZBmAdcENAqDtkJvDMsX65mV1vZtlmll1YqAdDRUQkdHXqEMWjV2dx5/TBLPliP9P+8hGLtxR4XZYEWcgMbHDOLXPODQZGA782s+hDx8wsErgIeOEbzp/vnMtyzmWlpqa2fMEiIiItyMy4anxPXp8ziZS4KK59fAW/e30DNfUa9NBeBDPE5QOZAe+7+fcdtY2Z+YBEYH9gA+fcJqAcGBKw+wJglXNuXzPXLCIi0qr17xzPa3Mmcs2Enjy+ZCcXz/2UnIJyr8uSIAhmiFsB9DOzXv6es5nAgiPaLACu9m9fCixyzjn/OT4AM+sBDAB2Bpw3i2+4lSoiItKWHRr08OjVWewrq+Z7//MJz63YpUEPbVzQQpz/GbY5wEJgE/C8c26Dmd1pZhf5mz0KpJhZDvBz4NA0JJOAtWa2hqbRp7Odc0UAZhYHnAu8HKxrERERaY2mDkznrVvPYGSPJH710jrm/GM1pVV1XpclLcTaY0rPyspy2dnZXpchIiLSIhobHQ99tJ173tlCekI0D8wawageHb0uS06Sma10zmUduT9kBjaIiIjI8QkLM26c0ocXbhiPGVz+0FLmfZBDY2P767hpyxTiRERE2qgR3Tvy5i1nMG1wZ/709haufnw5hQc1p1xboRAnIiLShiXGRPDgFSO465KhLN9RzAX3f8zH2zRfalugECciItLGmRmzxnRnwZxJJMdFcNVjy/nvhZupb2j0ujQ5BQpxIiIi7UT/zvG8dtMkfpCVydzFXzDr4c/YU1rldVlykhTiRERE2pGYyHD++P3TuX/mcDbuLuPC+z9m0WbNlR+KFOJERETaoenDM3j95kl0Tozhx09k819vbqROt1dDikKciIhIO9U7tQOvzJ7AleO68/DHO7j8oaXkH9Dt1VChECciItKORUeE8/uLh/I/s0awbV8533lAt1dDhUKciIiI8L1hXXn95kl09d9eveutTbq92sopxImIiAgAvTrF8fLsCVwxtjsPfbidWfM1erU1U4gTERGRw6IjwvnDjKHcP3M4m/aU8Z0HPtHkwK2UQpyIiIh8zfThGbw2ZxKdOkRy1WPL+ct7W2nQ2qutikKciIiIHFXftA68etNEZozI4C/vbeOax5ezv1xrr7YWCnEiIiJyTLGRPu65bBh/vGQoy3YU850HPmHll8VelyUoxImIiMi3MDNmjunOK7MnEBURxg8e+ozHPtmBc7q96iWFOBERETkug7smsmDOJKb0T+PONzYy55nVlNfUe11Wu6UQJyIiIsctMSaC+T8axa+mDeCtdXu46MFP2LbvoNdltUsKcSIiInJCwsKMG6f04emfjKOsqp6LHlzCa2vyvS6r3VGIExERkZMyvk8Kb94yiSEZCdz67Bp+9/oGrfIQRApxIiIictLSE6L5x0/H8eOJvXh8yU5mzf+MfWXVXpfVLijEiYiIyCmJCA/jP743iAdmjWDD7qZVHpZt3+91WW2eQpyIiIg0i4uGdeW1ORNJiPZxxSPLeOTj7ZqGpAUpxImIiEizOS09ntfmTOScgWn8/s1N3PLsGiprNQ1JS1CIExERkWYVHx3B364cxb+e3583Pt/NjLmfsrOowuuy2hyFOBEREWl2ZsZNZ/XlyWvHsO9gNd978BMWbd7ndVltikKciIiItJjJp6Xy+pxJZHaM5bons7n/vW00Nuo5ueagECciIiItKjM5lpdunMCM4Rnc995Wrn9qJQer67wuK+QpxImIiEiLi4kM557Lh/Gf3xvE4i0FTJ+7hJyCcq/LCmkKcSIiIhIUZsa1E3vx9E/GUlpZx8Vzl/DOhr1elxWyFOJEREQkqMb1TmHBzZPo1SmO659ayb3vbtVzcidBIU5ERESCLiMphhduGM/3R3bjgfe36Tm5k6AQJyIiIp6Ijgjnz5edfvg5uRnzPmV7oZ6TO15BDXFmNs3MtphZjpndcZTjUWb2nP/4MjPr6d8/xszW+F9rzWxGwDlJZvaimW02s01mNj6IlyQiIiKn4NBzck9dN4b95TVMn7uExVsKvC4rJAQtxJlZODAXuAAYBMwys0FHNLsOKHHO9QXuA+72718PZDnnhgPTgIfMzOc/dj/wtnNuADAM2NSiFyIiIiLNbkKfTiyYM4luHWP58RMrmPdBjtZd/RbB7IkbA+Q457Y752qBZ4HpR7SZDjzp334RmGpm5pyrdM4dWngtGnAAZpYITAYeBXDO1TrnDrTsZYiIiEhLaJpPbjzfGdqFP729hVueXUNVbYPXZbVawQxxGUBuwPs8/76jtvGHtlIgBcDMxprZBmAdcIP/eC+gEHjczFab2SNmFne0X25m15tZtpllFxYWNud1iYiISDOJjfTxP7NGHF539bKHPmX3gSqvy2qVQmZgg3NumXNuMDAa+LWZRQM+YCTwV+fcCKAC+Nqzdv7z5zvnspxzWampqUGrW0RERE7MoXVXH7kqi51FlVz04Cdk7yz2uqxWJ5ghLh/IDHjfzb/vqG38z7wlAvsDGzjnNgHlwBCaevPynHPL/IdfpCnUiYiISIibOjCdV2+aQIcoH7Me/oxnl+/yuqRWJZghbgXQz8x6mVkkMBNYcESbBcDV/u1LgUXOOec/xwdgZj2AAcBO59xeINfM+vvPmQpsbOkLERERkeDomxbPazdNYlzvFO54eR2/XbCB+oZGr8tqFXzf3qR5OOfqzWwOsBAIBx5zzm0wszuBbOfcApoGKDxlZjlAMU1BD2AScIeZ1QGNwGznXJH/2M3A0/5guB24NljXJCIiIi0vMTaCx68ZzV1vbebRT3aQU1DO3CtGkhgb4XVpnrL2OHw3KyvLZWdne12GiIiInKDnV+Ty76+uIyMphkeuHk3ftA5el9TizGylcy7ryP0hM7BBRERE5PLRmTzz03GU19Qzo51PDKwQJyIiIiElq2cyr82ZRLfkWK57YgWPfLy9XU4MrBAnIiIiIScjKYaXbhzPeYM68/s3N3HHS+uorW9fAx4U4kRERCQkxUb6mPfDkdx8dl+ey87lykeXUVxR63VZQaMQJyIiIiErLMy4/bz+3D9zOGtyDzB97ids3XfQ67KCQiFOREREQt704Rk8d/04qusauWTepyze3PYHPCjEiYiISJswontHFsyZSI+UWK57cgWPfbKjTQ94UIgTERGRNqNLYgwv3DCecwamc+cbG/l/r66nro2u8KAQJyIiIm1KbKSPv105ihun9OHpZbu45vHllFbWeV1Ws1OIExERkTYnLMz41bQB/Pelp7N8RzEz5i1hZ1GF12U1K4U4ERERabMuy8rk6Z+Mo6SylovnLWHZ9v1el9RsFOJERESkTRvTK5lXb5pISlwkVz66jBdX5nldUrNQiBMREZE2r0dKHC/fOJHRPZP5xQtr+fPCLTQ2hvbIVYU4ERERaRcSYyN48sdjmDk6kwcX53DzM6uprmvwuqyT5vO6ABEREZFgiQgP465LhtI7NY673tpM3oEqHrkqi9T4KK9LO2HqiRMREZF2xcy4fnIf/nblKLbsLWPGvCVsC8GluhTiREREpF06f3Bnnrt+PDX1jVzy10/5ZFuR1yWdEIU4ERERabeGZSbxyuwJdE2M4ZrHl/Pcil1el3TcFOJERESkXevWMZYXbhzP+D4p/Oqlddz99uaQGLmqECciIiLtXkJ0BI9dM5pZY7rz1w++4JZnW//IVY1OFREREaFp5OofZgyhR0osf3xrM3tLq5l/VRbJcZFel3ZU6okTERER8TMzbjizDw9eMYLP80u5ZN4SdrTSNVcV4kRERESO8N3Tu/LMT8dSWlXHJfOWkL2z2OuSvkYhTkREROQoRvVI5pXZE0mKjeSKR5bx5ud7vC7pKxTiRERERI6hZ6c4Xr5xAkMzErnpH6t4+KPtONc6Rq4qxImIiIh8g45xkTz9k7FcOLQz//XPTfx2wQYaWsEUJApxIiIiIt8iOiKcB2eN5PrJvXly6Zf87KmVVNbWe1qTQpyIiIjIcQgLM/7twoH87qLBLNq8j1nzP6O0ss6zejRPnIiIiMgJuHpCT7okRvPmuj10iPYuSinEiYiIiJyg8wZ35rzBnT2tQbdTRUREREKQQpyIiIhICFKIExEREQlBQQ1xZjbNzLaYWY6Z3XGU41Fm9pz/+DIz6+nfP8bM1vhfa81sRsA5O81snf9YdhAvR0RERMQzQRvYYGbhwFzgXCAPWGFmC5xzGwOaXQeUOOf6mtlM4G7gB8B6IMs5V29mXYC1Zva6c+7QBC1nOeeKgnUtIiIiIl4LZk/cGCDHObfdOVcLPAtMP6LNdOBJ//aLwFQzM+dcZUBgiwa8nyZZRERExEPBDHEZQG7A+zz/vqO28Ye2UiAFwMzGmtkGYB1wQ0Coc8A7ZrbSzK5vwfpFREREWo2QmSfOObcMGGxmA4Enzewt51w1MMk5l29macC7ZrbZOffRkef7A971AN27dw9q7SIiIiLNLZghLh/IDHjfzb/vaG3yzMwHJAL7Axs45zaZWTkwBMh2zuX79xeY2Ss03bb9Wohzzs0H5gOYWaGZfdksV3VsnQA9p9d89H02P32nzU/fafPS99n89J02r2B9nz2OtjOYIW4F0M/MetEU1mYCVxzRZgFwNbAUuBRY5Jxz/nNy/QMbegADgJ1mFgeEOecO+rfPA+78tkKcc6nNdlXHYGbZzrmslv497YW+z+an77T56TttXvo+m5++0+bl9fcZtBDnD2BzgIVAOPCYc26Dmd1JU4/aAuBR4CkzywGKaQp6AJOAO8ysDmgEZjvnisysN/CKmR26ln84594O1jWJiIiIeCWoz8Q55/4J/POIff8RsF0NXHaU854CnjrK/u3AsOavVERERKR104oNLWe+1wW0Mfo+m5++0+an77R56ftsfvpOm5en36c5pynXREREREKNeuJEREREQpBCnIiIiEgIUohrQWZ2mZltMLNGM9OQ7pNkZtPMbIuZ5ZjZHV7XE+rM7DEzKzCz9V7X0haYWaaZLTazjf7/3m/1uqZQZ2bRZrbczNb6v9PfeV1TW2Bm4Wa22sze8LqWtsDMdprZOjNbY2bZXtSgENey1gOXcJTJh+X4mFk4MBe4ABgEzDKzQd5WFfKeAKZ5XUQbUg/c7pwbBIwDbtKf0VNWA5ztnBsGDAemmdk4b0tqE24FNnldRBtzlnNuuFdzxSnEtSDn3Cbn3Bav6whxY4Ac59x251wt8Cww3eOaQpp/Wbpir+toK5xze5xzq/zbB2n6S/LIdaHlBLgm5f63Ef6XRuGdAjPrBnwHeMTrWqT5KMRJa5cB5Aa8z0N/QUorZWY9gRHAMo9LCXn+W39rgALgXf/62XLy/gL8kqYJ86V5OOAdM1vpX5896II62W9bZGbvAZ2PcujfnXOvBbseEfGGmXUAXgJuc86VeV1PqHPONQDDzSyJppV5hjjn9BznSTCz7wIFzrmVZjbF43LakknOuXwzSwPeNbPN/jsdQaMQd4qcc+d4XUMblw9kBrzv5t8n0mqYWQRNAe5p59zLXtfTljjnDpjZYpqe41SIOzkTgYvM7EIgGkgws7875670uK6Q5pzL9/8sMLNXaHr8J6ghTrdTpbVbAfQzs15mFknTeroLPK5J5DBrWrz5UWCTc+5er+tpC8ws1d8Dh5nFAOcCmz0tKoQ5537tnOvmnOtJ0/+HLlKAOzVmFmdm8Ye2gfPw4B8ZCnEtyMxmmFkeMB5408wWel1TqHHO1QNzgIU0PTD+vHNug7dVhTYzewZYCvQ3szwzu87rmkLcROBHwNn+qQbW+Hs85OR1ARab2ec0/UPuXeecpsWQ1iQd+MTM1gLLgTedc28HuwgtuyUiIiISgtQTJyIiIhKCFOJEREREQpBCnIiIiEgIUogTERERCUEKcSIiIiIhSCFOREREJAQpxImIiIiEIIU4EZFTZGYfmNkA/3aKmWl5KBFpcQpxIiKnri+w1b99OrDOw1pEpJ1QiBMROQVm1gPId841+nedDnzuYUki0k4oxImInJphfDW0jUIhTkSCQCFOROTUDAeiAcysHzAd3U4VkSBQiBMROTXDgDAzWwv8B7ARuNrbkkSkPTDnnNc1iIiELDPbBox0zh30uhYRaV/UEycicpLMLB5wCnAi4gX1xImIiIiEIPXEiYiIiIQghTgRERGREKQQJyIiIhKCFOJEREREQpBCnIiIiEgIUogTERERCUEKcSIiIiIh6P8HJ2YRYFhBj6AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lam = 0.0\n",
    "tau2 = 100.0\n",
    "p_mu = stats.norm(lam, np.sqrt(tau2)).pdf(mus)\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.plot(mus, p_mu);\n",
    "ax.set_xlabel(\"$\\mu$\")\n",
    "ax.set_ylabel(\"$P(\\mu)$\");\n",
    "#ax.set_ylim((0.0, 0.1));   # y軸を固定"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8765366-5bb7-4b4b-a538-dd47875f6ec1",
   "metadata": {},
   "source": [
    "## 事後分布\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\mu|X=x_1,…,x_n) & \\propto P(\\mu)f(x_1,…,x_n|\\mu) \\\\\n",
    "& \\propto \\exp \\left[ -\\frac{n \\tau^2 + \\sigma^2}{2 \\sigma^2 \\tau^2} \\left(\\mu - \\frac{n \\tau^2 \\bar{x} + \\sigma^2 \\lambda}{n \\tau^2 + \\sigma^2} \\right)^2 \\right] \\\\\n",
    "&= \\mathcal{N} \\left( \\frac{n \\tau^2 \\bar{x} + \\sigma^2 \\lambda}{n \\tau^2 + \\sigma^2} , \\frac{\\sigma^2 \\tau^2}{n \\tau^2 + \\sigma^2} \\right)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "993542f1-2702-47bd-90f8-447d054cd39e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.3666353635392703\n",
      "2.3654526372206597 0.04997501249375312\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAEICAYAAAAAxLONAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAApn0lEQVR4nO3deXzc9X3n8ddndNqSLAnrsC35tvARbMugNU0I4SghpmlxkiYp5CjpIyk9oM12+2ibbLfJLtnHNi2P7aZpTkpYmrbAJiEkNCEhkEAgIQSEbWzjUzY+JGxpLNk6bV3z2T804wjhQ5Y1v99vZt7Px0MPzfyOmY/nAdJb39PcHRERERGJpljYBYiIiIjI2SmsiYiIiESYwpqIiIhIhCmsiYiIiESYwpqIiIhIhOWHXUA6VVVV+aJFi8IuQ0REROS8XnrppWPuXj3xeFaHtUWLFtHc3Bx2GSIiIiLnZWYHz3Rc3aAiIiIiEaawJiIiIhJhCmsiIiIiEaawJiIiIhJhCmsiIiIiEaawJiIiIhJhCmsiIiIiERbYOmtmdh/wm0CHu192hvN/AXxwXF0rgWp37zKzA0AvMAqMuHtTMFWLiITr1PAoR7tP0dk/SGffEJ39Q3T1D3Hd8hpWzZsVdnkiEoAgF8W9H/gC8PUznXT3u4G7Aczst4A/c/eucZdc5+7H0l2kiEhUnBoe5eq/f4p47+Abzj25s51H/viqEKoSkaAFFtbc/RkzWzTJy28FHkxjOSIikffkznbivYP8+dsvZXV9ObNLiphdWsgjm9u4+/HdtHT0saymNOwyRSTNIjdmzcxmAhuAh8cdduBHZvaSmd1+nvtvN7NmM2uOx+PpLFVEJK2+s/k1amcV8cfXLePa5TWsri9nXsUM3tdUT17M+NZLrWGXKCIBiFxYA34L+PmELtC3uvvlwE3AHWb2trPd7O73uHuTuzdVV79hL1QRkYxwvH+In+7p4Oa188iL2evO1ZQVc93yah7Z3MpowkOqUESCEsWwdgsTukDdvS35vQN4BFgfQl0iIoH5/rYjDI86Gxvrznj+vVfU094zyDN71YMgku0iFdbMrBy4BvjuuGMlZlaWegzcCGwPp0IRkWB8d0sbDTWlvOksMz6vX1FL5cwCdYWK5IAgl+54ELgWqDKzVuDTQAGAu38ledm7gR+5e/+4W2uBR8wsVe8D7v7DoOoWEQla6/EBXjxwnL94x3KSP/veoDA/xsbGOh745SFODAxRMbMw4CpFJChBzga9dRLX3M/YEh/jj+0H1qanKhGR6PnultcAuHntvHNe994r6rn/uQP8x8uv8eE3LwqgMhEJQ6S6QUVEcp27853NbTQtrGT+JTPPee1ldeWsnDuLb6orVCSrKayJiETIjiM97O3o413rzjyxYKL3XlHP1tZudh/tTXNlIhIWhTURkQj57pbXyI8Z71w9d1LXv6txHvkx4+FNal0TyVYKayIiETGacB7d8hrXLq+msmRyEwZmlxZx/Yoavr2pjeHRRJorFJEwKKyJiETEL1/t5GjPqbOurXY2772inmN9gzyzR2uuiWQjhTURkYj4zuY2SovyuWFl7QXdd92KGmaXFPLtTW1pqkxEwqSwJiISAYMjo/xg+1He8aY5zCjMu6B7C/JiXHNpNS8c6MJd20+JZBuFNRGRCNje1kPvqRFufNOFtaqlrKkvJ947yNGeU9NcmYiETWFNRCQCtrWeAKBxfsWU7l+TvO/lw93TU5CIRIbCmohIBGxt66a6rIjaWcVTun/V3Fnkx4ytydAnItlDYU1EJAK2t3Wzpq58yvcXF+SxfE4ZW1vVsiaSbRTWRERC1j84QktHH6vrpx7WANbUV7C19YQmGYhkGYU1EZGQ7TjSQ8LHJglcjLX15fScGuFA58A0VSYiUaCwJiISslTX5WUX0Q0KYy1rY6934iIrEpEoUVgTEQnZttYTzJlVTE3Z1CYXpFxaW0pRfkwzQkWyjMKaiEjItrV1X/R4NYD8vBhvmjeLbW0nLr4oEYkMhTURkRD1nhpm/7H+i5oJOt6a+gq2t/Uwok3dRbKGwpqISIheea0Hd7hsGlrWANbOL+fk8Cgt8b5peT0RCZ/CmohIiLYlJxesnsaWNYCtGrcmkjUCC2tmdp+ZdZjZ9rOcv9bMus1sS/LrU+PObTCz3WbWYmafCKpmEZF029bWTV3FDKpKi6bl9RbPLqGsKJ+XNSNUJGsE2bJ2P7DhPNc86+6Nya+7AMwsD/gicBOwCrjVzFaltVIRkYBsa+vmsrpZ0/Z6sZixur5cOxmIZJHAwpq7PwN0TeHW9UCLu+939yHgIWDjtBYnIhKC7pPDvHqs/3TX5XRZU1/BrqM9DI6MTuvrikg4ojZm7c1m9rKZ/cDM3pQ8VgccHndNa/LYGZnZ7WbWbGbN8Xg8nbWKiFyUV9qmd7xaypr6coZHnZ1Heqf1dUUkHFEKa5uAhe6+Fvgn4DtTeRF3v8fdm9y9qbq6ejrrExGZVlvTGNZAOxmIZIvIhDV373H3vuTjx4ACM6sC2oD54y6tTx4TEclo29q6qa+cQWVJ4bS+bl3FDGaXFGrcmkiWiExYM7M5ZmbJx+sZq60TeBFoMLPFZlYI3AI8Gl6lIiLTY1tr90Vv3n4mZsaa+nK1rIlkifyg3sjMHgSuBarMrBX4NFAA4O5fAd4L/JGZjQAngVvc3YERM7sTeBzIA+5z91eCqltEJB1ODAxxqGuAW9cvSMvrr6mv4Kd74vQPjlBSFNiPehFJg8D+D3b3W89z/gvAF85y7jHgsXTUJSIShm3J8WrpaFmDsZ0MEg7b27q5csnstLyHiAQjMt2gIiK5JBXWLpuXnrB2eicDjVsTyXgKayIiIdjW2s3C2TMpn1mQltevKi2irmKGdjIQyQIKayIiIdja2j3tS3ZMtLqu/HQLnohkLoU1EZGAdfYN0nbiZNrGq6WsmFvGoa4BTg1rJwORTKawJiISsB1HeoD0jVdLaagpwx32xfvS+j4ikl4KayIiAdvbPhaeGmrL0vo+DbWlALR0KKyJZDKFNRGRgO3t6KNiZgFVpdO7c8FEi2aXkBez0+FQRDKTwpqISMD2dfTRUFNKctOWtCnMj7Fo9kz2dmhDd5FMprAmIhIgd2dPRy/LatLbBZrSUFPGXnWDimQ0hTURkQB19g9xYmCYZTWlgbxfQ20pBzsHGBzRjFCRTKWwJiISoNRg/4aAwtqymlJGE86BYwOBvJ+ITD+FNRGRAKW6JFMzNdOtIdndqnFrIplLYU1EJEAt7b2UFuUzZ1ZxIO+3pLqEmKEZoSIZTGFNRCRAezv6WBrATNCU4oI8FlwyU2utiWQwhTURkQC1JJftCNKymjJ1g4pkMIU1EZGAdA8M09E7GHhYa6gt5dVj/QyPJgJ9XxGZHgprIiIBaYmPtW4FtWxHSkNNKcOjzsFOzQgVyUQKayIiAfnVsh3BLIibknq/FnWFimQkhTURkYDsbe+juCBGXeWMQN93aU3J6fcXkcwTWFgzs/vMrMPMtp/l/AfNbKuZbTOz58xs7bhzB5LHt5hZc1A1i4hMp70dfSypKiUvFsxM0JSZhfnUV87QtlMiGSrIlrX7gQ3nOP8qcI27rwY+A9wz4fx17t7o7k1pqk9EJK1aOvoCWwx3ooaaUoU1kQwVWFhz92eArnOcf87djyefPg/UB1KYiEgA+gdHaDtxMvCZoCkNtWXsi/cxmvBQ3l9Epi6qY9Y+Cvxg3HMHfmRmL5nZ7ee60cxuN7NmM2uOx+NpLVJEZLL2xcdatZYFPLkgZVlNKUMjCQ53aUaoSKaJXFgzs+sYC2t/Ne7wW939cuAm4A4ze9vZ7nf3e9y9yd2bqqur01ytiMjkpAb3B71sR0qqRU9doSKZJ1JhzczWAPcCG929M3Xc3duS3zuAR4D14VQoIjI1LfE+CvKMhbNnhvL+y06HNS3fIZJpIhPWzGwB8G3gw+6+Z9zxEjMrSz0GbgTOOKNURCSq9rb3sbiqhIK8cH7slhUXMLe8mBYt3yGScfKDeiMzexC4Fqgys1bg00ABgLt/BfgUMBv4UnKD45HkzM9a4JHksXzgAXf/YVB1i4hMh5aOXlbNmxVqDcs0I1QkIwUW1tz91vOc/xjwsTMc3w+sfeMdIiKZ4dTwKIe6Bri5sS7UOhpqynjwhUMkEk4s4LXeRGTqItMNKiKSrV491k/CCW3ZjpSG2lJODo/SduJkqHWIyIVRWBMRSbNU12NYC+KmpMJii7pCRTKKwpqISJq1tPcSM1hcVRJqHZoRKpKZFNZERNKsJd7HwtklFOXnhVpHxcxCqsuKtKG7SIZRWBMRSbO97X2hLYY7kfYIFck8CmsiImk0PJrg1WP9kQprLR19uGuPUJFMobAmIpJGBzsHGEl46DNBU5bVlNI3OEJ7z2DYpYjIJCmsiYikUUtyMH9UWtaWakaoSMZRWBMRSaNUKFpaHY2wtixZx764wppIplBYExFJo5aOPuaVF1NSFNiGMedUXVZEWVG+WtZEMojCmohIGu2L95/ueowCM2NpTala1kQyiMKaiEiaJBLOvnh0lu1IWVpdqpY1kQyisCYikiZHe04xMDQamfFqKctqSunoHaTn1HDYpYjIJCisiYikSar1Knota2PbXu1T65pIRlBYExFJk6jNBE1Jhcd98f6QKxGRyVBYExFJk33xPspnFFBVWhh2Ka+z4JKZFOSZxq2JZAiFNRGRNGnp6GNpdQlmFnYpr5OfF2PR7BLNCBXJEAprIiJpsi8enT1BJ1paXaoxayIZQmFNRCQNugeGOdY3GNmwtqymlINdAwyNJMIuRUTOI9CwZmb3mVmHmW0/y3kzs8+bWYuZbTWzy8edu83M9ia/bguuahGRC9cSH9sTNGqTC1KW1pQwmnAOdmqSgUjUBd2ydj+w4RznbwIakl+3A18GMLNLgE8DVwLrgU+bWWVaKxURuQj7OsZCUGRb1qrLAO0RKpIJAg1r7v4M0HWOSzYCX/cxzwMVZjYXeAfwhLt3uftx4AnOHfpERELVEu+jMD9GfeXMsEs5oyXJtdY0I1Qk+qI2Zq0OODzueWvy2NmOv4GZ3W5mzWbWHI/H01aoiMi57OvoY0lVCXmxaM0ETSkpymdeebHWWhPJAFELaxfN3e9x9yZ3b6qurg67HBHJUS3xvkht4H4mS2u0R6hIJohaWGsD5o97Xp88drbjIiKRc2p4lMNdA5GdXJCytLqUffE+3D3sUkTkHKIW1h4Ffjc5K/TXgG53PwI8DtxoZpXJiQU3Jo+JiETOgc5+Eh7dyQUpS2tKGRga5Uj3qbBLEZFzyA/yzczsQeBaoMrMWhmb4VkA4O5fAR4DfgNoAQaA30ue6zKzzwAvJl/qLnc/10QFEZHQ/GpP0JKQKzm3ZdWpPUL7mFcxI+RqRORsLjismVkJcMrdRy/0Xne/9TznHbjjLOfuA+670PcUEQlaS0cfZtFdYy1lac2vZoRe3aAxviJRdd5uUDOLmdkHzOz7ZtYB7AKOmNkOM7vbzJalv0wRkcyxL95PfeUMigvywi7lnKpLi5hVnK+11kQibjJj1p4ClgKfBOa4+3x3rwHeCjwP/J2ZfSiNNYqIZJSxDdyj3aoGYGaaESqSASbTDXqDuw9PPJgcM/Yw8LCZFUx7ZSIiGSiRcPbH+7hq6eywS5mUZdWlPL1Ha1KKRNl5W9bOFNSmco2ISC5oO3GSwZFE5NdYS1laU0q8d5Duk/oxLhJVk166w8y2mdm/m9lfmdlNZlZvZn+dzuJERDJNqksx6st2pIyfESoi0XQh66xdA/wzcBK4BdjO2DIbIiKSlAo9mTBmDTjdAqhxayLRNemlO5Jj1J5OfmFmDcB/S0tVIiIZqqWjj0tKCrmkpDDsUiZlfuUMCvNialkTibAL6Qa9dPxzd98LrJn2ikREMti+eN/prsVMkJ8XY1HVTPapZU0ksi5kUdyvmtlSxvbk3AoUA9vNrMTd+9NSnYhIhmnp6GPDZXPCLuOCLKspZeeR3rDLEJGzmMyiuAbg7te5+wLgd4DvMbYl1Axgk5ntTmuVIiIZoLNvkOMDwxkzXi1laXUpBzv7GRy54I1pRCQAk2lZe8rMHga+6+6H3P0QcMjMHgeuBm4DmtNZpIhIJtgXH+tkyJRlO1KW1ZSScHj1WD8r5swKuxwRmWAyY9Y2AKPAg2aW2mZqP7CXsVmh/+Dun09nkSIimeD0sh0Z1rK2fE4ZALuPqitUJIrO27Lm7qeALwFfSu5UUAWcdPcTaa5NRCSj7GnvZUZBHnUVM8Iu5YIsqSolP2YKayIRNZkxa7eZ2TEz6wLuBfoU1ERE3mj30V4unVNGLGZhl3JBCvNjLKkuUVgTiajJdIP+DfB2YAVwCPhfaa1IRCQDuTu723tZUVsWdilTsnzOLHYprIlE0mTCWo+7b3b3Dnf/G2B9uosSEck08b5BuvqHTo//yjQr5pTRduIkvae0R6hI1EwmrM01s9vN7G1mVg0UpLsoEZFMsyu5TtmKuZkZ1pYnWwT3tKt1TSRqJhPWPg2sBj4D7AYuM7PHzOxvzezWtFYnIpIhUuO9MnXpi1SLoLpCRaJnMrNB7xn/3MzqGQtvaxjbyP3B9JQmIpI5dh3tpbqsKGP2BJ2ovnIGJYV57FFYE4mcC9luCgB3bwVagR9c6L1mtgH4RyAPuNfdPzvh/P8Brks+nQnUuHtF8twosC157pC733yh7y8iki6723tYkaHj1QDMjEvnlKllTSSCLjisTZWZ5QFfZGxmaSvwopk96u47Ute4+5+Nu/5PgHXjXuKkuzcGVK6IyKSNjCbY097HbW9eGHYpF2XFnDJ+sP0o7k5yp0ERiYDJjFmbLuuBFnff7+5DwEPAxnNcfyvqYhWRDHCgc4ChkQTLM3S8Wsry2jJODAzT0TsYdikiMk6QYa0OODzueWvy2BuY2UJgMfCTcYeLzazZzJ43s3ed7U2SM1ebzaw5Ho9PQ9kiIuf2q8kFmdsNCpwOm+oKFYmWIMPahbgF+Ja7j447ttDdm4APAJ8zs6VnutHd73H3Jndvqq6uDqJWEclxu4/2ELOxDdEz2YrTe4T2hFyJiIwXZFhrA+aPe16fPHYmtzChC9Td25Lf9wNP8/rxbCIiodl5tJfFVSUUF+SFXcpFqSwppKasiN1H+8IuRUTGCTKsvQg0mNliMytkLJA9OvEiM1sBVAK/GHes0syKko+rgKuAHRPvFREJw+6jvRm7vtpEy+eUsbtdLWsiURJYWHP3EeBO4HFgJ/ANd3/FzO4ys/HLcNwCPOTuPu7YSqDZzF4GngI+O34WqYhIWPoHRzjUNZCx20xNtLy2jL3tfYwm/PwXi0ggAlu6A8DdHwMem3DsUxOe//cz3PccYwvxiohESmp7pkyfXJCyfE4ZgyMJDnT2s7Q6s8fgiWSLqE4wEBHJCJm+zdREqX/Hbs0IFYkMhTURkYuw62gvMwvzqK+cEXYp06KhtpSYafkOkShRWBMRuQi7jvawfE4ZsVh2rPhfXJDHotkl2iNUJEIU1kREpsjdkzNBs2O8WsqltWXsbldYE4kKhTURkSmK9w5yfGCY5bXZFdaWzynjQGc/J4dGz3+xiKSdwpqIyBTtTHYVZvqeoBOtmFOGO+ztUOuaSBQorImITFFqW6Zs6wZNrRmnSQYi0aCwJiIyRbuO9lI7q4jKksKwS5lWC2eXUFwQ0/IdIhGhsCYiMkW7j/ZmXRcoQF7MaKgpO73gr4iES2FNRGQKRkYT7O3oY2WWdYGmXFpbpm5QkYhQWBMRmYIDnf0MjSSyZk/QiVbMKSPeO0hX/1DYpYjkPIU1EZEp2HV6JmiWhrW5Y/+unUd6Qq5ERBTWRESmYPfRXvJixrKa7NzsfHVdOQBbDp8ItxARUVgTEZmKbW3dLKsupSg/L+xS0qJiZiFLqkvYfOhE2KWI5DyFNRGRC5RIOJsPneDyhRVhl5JWjfMr2HL4OO4edikiOU1hTUTkAu0/1k/3yWHWza8Mu5S0WregkmN9Q7QePxl2KSI5TWFNROQCbT50HCDrW9bWza8AYLPGrYmESmFNROQCbTp0glnF+Sypys7JBSnL55RRXBA7HU5FJBwKayIiF2jzoeM0LqgkFrOwS0mrgrwYa+oqNCNUJGSBhjUz22Bmu82sxcw+cYbzHzGzuJltSX59bNy528xsb/LrtiDrFhFJ6RscYU977+kuwmzXuKCCV9p6GBwZDbsUkZwVWFgzszzgi8BNwCrgVjNbdYZL/5+7Nya/7k3eewnwaeBKYD3waTPL7pG9IhJJLx8+QcLh8oW58SNo3fwKhkYT7DyiradEwhJky9p6oMXd97v7EPAQsHGS974DeMLdu9z9OPAEsCFNdYqInFVq/FZjjrSsrVswFko1bk0kPEGGtTrg8LjnrcljE/22mW01s2+Z2fwLvBczu93Mms2sOR6PT0fdIiKnbTp0gmU1pZTPKAi7lEDMKS9mzqxiLY4rEqKoTTD4D2CRu69hrPXsXy70Bdz9Hndvcvem6urqaS9QRHKXu7P50HEuX1ARdimBWrdAkwxEwhRkWGsD5o97Xp88dpq7d7r7YPLpvcAVk71XRCTdDnQOcHxg+HTXYK5Yt6CCQ10DHOsbPP/FIjLtggxrLwINZrbYzAqBW4BHx19gZnPHPb0Z2Jl8/Dhwo5lVJicW3Jg8JiISmE0Hk4vh5lhYa0zu1LBFXaEioQgsrLn7CHAnYyFrJ/ANd3/FzO4ys5uTl/2pmb1iZi8Dfwp8JHlvF/AZxgLfi8BdyWMiIoHZfPg4pUX5LKvJ7sVwJ1pdV05ezNQVKhKS/CDfzN0fAx6bcOxT4x5/EvjkWe69D7gvrQWKiJzDpoMnaJxfQV6WL4Y70YzCPFbOLWPzYc0IFQlD1CYYiIhE0sDQCLuO9rAuxyYXpDTOr+Dlw92MJjzsUkRyjsKaiMgkvHy4e2wx3Bwbr5aybn4lfYMj7Iv3hV2KSM5RWBMRmYRNObYY7kSpFkUtjisSPIU1EZFJ2HzoBEuqSqgsKQy7lFAsriqhfEaBJhmIhEBhTUTkPFKL4eba+mrjmRmN8yu0k4FICBTWRETO43DXSTr7h3J2ckFK4/wKdrf30jc4EnYpIjlFYU1E5DxS49VydXJByroFFbjDVnWFigRKYU1E5Dw2HTrOzMI8Lq3NrcVwJ1q3oJK8mPHzfcfCLkUkpyisiYicx/P7O1m3oIL8vNz+kVk+o4CmhZX8eGdH2KWI5JTc/skjInIeBzv72dPex/UrasMuJRJuWFnLrqO9tB4fCLsUkZyhsCYicg5PJluRblhZE3Il0XB98nN4apda10SCorAmInIOT+5o59LaUhbOLgm7lEhYWl3K4qqS0yFWRNJPYU1E5Cy6B4Z54UAXN6xUF+h416+o4Rf7OunXEh4igVBYExE5i6f3dDCacG5YpbA23q+vrGFoNMHPWjQrVCQICmsiImfxxI52qkqLaKyvCLuUSPlPiy6hrCifn6grVCQQCmsiImcwNJLgp7vj3LCyhljMwi4nUgryYrxteTU/3tVBIuFhlyOS9RTWRETO4JevdtI7OKLxamdxw8oajvUNsq2tO+xSRLKewpqIyBk8uaOd4oIYVy2rCruUSLr20hpiBj/WEh4iaaewJiIygbvz5M4O3rqsmhmFeWGXE0mVJYVcsbCSH+9sD7sUkawXaFgzsw1mttvMWszsE2c4/1/MbIeZbTWzH5vZwnHnRs1sS/Lr0SDrFpHcsvNIL20nTvL2VVoI91yuX1HLK6/1cLT7VNiliGS1wMKameUBXwRuAlYBt5rZqgmXbQaa3H0N8C3g78edO+nujcmvmwMpWkRy0pM72zFDW0ydR2pXhx/vUuuaSDoF2bK2Hmhx9/3uPgQ8BGwcf4G7P+XuqQ3nngfqA6xPRAQYC2vr5ldQXVYUdimRtqymlPmXzNASHiJpFmRYqwMOj3vemjx2Nh8FfjDuebGZNZvZ82b2rrPdZGa3J69rjsfjF1WwiOSeo92n2NrarYVwJ8HM+PUVtfys5Rgnh0bDLkcka0VygoGZfQhoAu4ed3ihuzcBHwA+Z2ZLz3Svu9/j7k3u3lRdXR1AtSKSTZ5MDph/u5bsmJRfX1nD4EiC5/ZpNwORdAkyrLUB88c9r08eex0zuwH4a+Bmdx9MHXf3tuT3/cDTwLp0FisiuemJHe0snD2TZTWlYZeSEa5cPJvyGQV866XWsEsRyVpBhrUXgQYzW2xmhcAtwOtmdZrZOuCrjAW1jnHHK82sKPm4CrgK2BFY5SKSEw529vPs3jjvXD0XM+1aMBmF+TE+eOUCHn/lKAc7+8MuRyQrBRbW3H0EuBN4HNgJfMPdXzGzu8wsNbvzbqAU+OaEJTpWAs1m9jLwFPBZd1dYE5Fpde+zr5Ifi3HbWxaFXUpGue0ti8iLGff97NWwSxHJSvlBvpm7PwY8NuHYp8Y9vuEs9z0HrE5vdSKSy471DfKN5sO8e10dtbOKwy4no9TOKubmtXV8o7mVP3v7pVTMLAy7JJGsEskJBiIiQfuX5w4wNJrg9muWhF1KRvr9ty3m5PAo//7LQ2GXIpJ1FNZEJOf1D47w9V8c5MZVtSyt1sSCqVgxZxZXN1Rx/3MHGBzRMh4i00lhTURy3kMvHqb75DB/cM0ZVwSSSfr9q5cQ7x3k0S2vhV2KSFZRWBORnDY8muBrz+5n/eJLuHxBZdjlZLSrG6pYMaeMe599FXcPuxyRrKGwJiI57T9efo3Xuk/xhxqrdtHMjI9dvYTd7b08s1eL5IpMF4U1EclZ7s5Xf7qf5bVlXLe8JuxyssLNa+dRU1bEvc/uD7sUkayhsCYiOevp3XF2t/fyB9cs0SK406QwP8ZHrlrEs3uPseO1nrDLEckKCmsikpPcnS8/vY955cX81tp5YZeTVT64fiEzC/O4+/FdGrsmMg0U1kQkJ9338wO8cKCLP7puGQV5+lE4ncpnFvAX71jOU7vjfE27GohcNP2EEpGc03ygi799bCc3rqrlQ1cuCLucrPSRtyzixlW1/N0Pd7Hl8ImwyxHJaAprIpJTjvUNcscDm6irnMHd71ursWppYmbc/d611JQVc+cDm+g+ORx2SSIZS2FNRHLGaML5+EObOTEwzJc/eAXlMwrCLimrlc8s4J8+sI6j3af4xMNbNX5NZIoU1kQkZ3zuyT38vKWTz7zrMlbNmxV2OTnh8gWV/OWG5fxg+1H+7fmDYZcjkpEU1kQkJzy1q4N/+kkL72+q5/1N88MuJ6d87K1LuG55NZ/53k62t3WHXY5IxlFYE5Gs99M9cT7+0GZWzZ3FXRsvC7ucnBOLGf/7/Y1cUlLIh7/2S57a3RF2SSIZRWFNRLLW4Mgo//N7O7jtvheYU17MVz98BcUFeWGXlZMuKSnkgd+/ktpZxfze/32Rv/vhLkZGE2GXJZIRFNZEJCvti/fxni89x70/e5XfffNCHr3zrcy/ZGbYZeW0JdWlfOeOq7h1/QK+/PQ+brnneV47cTLsskQiT2FNRLLKaMJ56IVD/Obnf8ZrJ07yz7/bxF0bL1OLWkQUF+Txt+9ZzT/e0sjOIz288/PP8v2tR9TKJnIO+WEXICIyHfa29/Lwpja+s7mNoz2neMvS2fzD+xuZU14cdmlyBhsb61hdV84dD2zmjgc2UTuriN++vJ73Nc1ncVVJ2OWJRIoFue6NmW0A/hHIA+51989OOF8EfB24AugEfsfdDyTPfRL4KDAK/Km7P36+92tqavLm5uZp/TeISDScHBqlpaOPFw908cjmNra1dZMXM669tJr3XF7PTZfNIRbTgrdRNzya4Mc7O/hm82Ge2t1BwmH9oku4uXEeq+vKubS2jBmFahWV3GBmL7l70xuOBxXWzCwP2AO8HWgFXgRudfcd4675Y2CNu/+hmd0CvNvdf8fMVgEPAuuBecCTwKXuPnqu91RYE8ks7s7wqDM4Mkr/4CjHB4Y4PjDEiYFhjg8M0d59it3tvexp7+NAZz+pH1+X1c3iPevqublxHlWlReH+I2TK2ntO8e1NbXyz+TD7j/UDEDNYVFXCyrmzaKgppaq0iEtKCk9/VcwooLgwj6L8GIV5Me1IIRntbGEtyG7Q9UCLu+9PFvQQsBHYMe6ajcB/Tz7+FvAFG/s/byPwkLsPAq+aWUvy9X4RUO1n9Af/2sz+eH+YJYik1WT/lBv/R5+Pe+DJc2PfwXESibFjo+6MJlIBLcHgyNjXufzqF3cZGxvnsby2jFXzZrFwtrrNskHtrGL+6Nql/OE1SzjYOcCuoz3sPNLLrqM9bGvt5vtbj5z3NYryYxTlx8jPixEzIy8GeWbEYkbMDDMwxrbDMhh7kpR6eCGBT9EwNxQVxPjen1wd2vsHGdbqgMPjnrcCV57tGncfMbNuYHby+PMT7q0705uY2e3A7QALFqR3g+aFs0vIUzeLZDmb7K+js/zSG/vFOO4XpDH2SzT5CzQvBvmxGEUFMYry807/sp1RmEflzEIqZhZQObNw7KukgKJ8dYllOzNjUVUJi6pK2HDZ3NPHh0YSHB8Yoqt/iOP9Q3QNDHF8YJjB4dGxsJ/6PpJgJJFgNAGJxNgfBomEk3jdHw5n/yNjsvxCLpaMVpAX7nzMrJtg4O73APfAWDdoOt/rv/7GynS+vIiIjFOYH6N2VjG1szRpRHJLkFGxDRi/x0t98tgZrzGzfKCcsYkGk7lXREREJOsEGdZeBBrMbLGZFQK3AI9OuOZR4Lbk4/cCP/GxdupHgVvMrMjMFgMNwAsB1S0iIiISmsC6QZNj0O4EHmds6Y773P0VM7sLaHb3R4GvAf+anEDQxVigI3ndNxibjDAC3HG+maAiIiIi2SDQddaCpqU7REREJFOcbekObTclIiIiEmEKayIiIiIRprAmIiIiEmEKayIiIiIRltUTDMwsDhxM89tUAcfS/B65RJ/n9NNnOv30mU4vfZ7TT5/p9Arq81zo7tUTD2Z1WAuCmTWfaeaGTI0+z+mnz3T66TOdXvo8p58+0+kV9uepblARERGRCFNYExEREYkwhbWLd0/YBWQZfZ7TT5/p9NNnOr30eU4/fabTK9TPU2PWRERERCJMLWsiIiIiEaawJiIiIhJhCmvTwMzeZ2avmFnCzDRVeorMbIOZ7TazFjP7RNj1ZDozu8/MOsxse9i1ZAMzm29mT5nZjuT/7x8Pu6ZMZ2bFZvaCmb2c/Ez/R9g1ZQMzyzOzzWb2vbBryQZmdsDMtpnZFjNrDqMGhbXpsR14D/BM2IVkKjPLA74I3ASsAm41s1XhVpXx7gc2hF1EFhkB/tzdVwG/Btyh/0Yv2iBwvbuvBRqBDWb2a+GWlBU+DuwMu4gsc527N4a11prC2jRw953uvjvsOjLceqDF3fe7+xDwELAx5Joymrs/A3SFXUe2cPcj7r4p+biXsV+GdeFWldl8TF/yaUHyS7PeLoKZ1QPvBO4NuxaZPgprEhV1wOFxz1vRL0KJKDNbBKwDfhlyKRkv2WW3BegAnnB3faYX53PAXwKJkOvIJg78yMxeMrPbwyggP4w3zURm9iQw5wyn/trdvxt0PSISDjMrBR4G/rO794RdT6Zz91Gg0cwqgEfM7DJ31zjLKTCz3wQ63P0lM7s25HKyyVvdvc3MaoAnzGxXsuciMAprk+TuN4RdQ5ZrA+aPe16fPCYSGWZWwFhQ+3d3/3bY9WQTdz9hZk8xNs5SYW1qrgJuNrPfAIqBWWb2b+7+oZDrymju3pb83mFmjzA2bCfQsKZuUImKF4EGM1tsZoXALcCjIdckcpqZGfA1YKe7/0PY9WQDM6tOtqhhZjOAtwO7Qi0qg7n7J9293t0XMfYz9CcKahfHzErMrCz1GLiREP6YUFibBmb2bjNrBd4MfN/MHg+7pkzj7iPAncDjjA3c/oa7vxJuVZnNzB4EfgEsN7NWM/to2DVluKuADwPXJ6fwb0m2YMjUzQWeMrOtjP3B9oS7a7kJiZJa4Gdm9jLwAvB9d/9h0EVouykRERGRCFPLmoiIiEiEKayJiIiIRJjCmoiIiEiEKayJiIiIRJjCmoiIiEiEKayJiIiIRJjCmoiIiEiEKayJiEySmT1tZiuSj2ebmbZFEpG0U1gTEZm8ZcCe5OM1wLYQaxGRHKGwJiIyCWa2EGhz90Ty0Bpga4gliUiOUFgTEZmctbw+nF2BwpqIBEBhTURkchqBYgAzawA2om5QEQmAwpqIyOSsBWJm9jLwKWAHcFu4JYlILjB3D7sGEZHIM7O9wOXu3ht2LSKSW9SyJiJyHmZWBriCmoiEQS1rIiIiIhGmljURERGRCFNYExEREYkwhTURERGRCFNYExEREYkwhTURERGRCFNYExEREYkwhTURERGRCPv/2JIHOSlEUMsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_bar = sample.mean()\n",
    "post_mu = (n * tau2 * x_bar + var_set * lam) / (n * tau2 + var_set)\n",
    "post_var = (var_set * tau2) / (n * tau2 + var_set)\n",
    "print(x_bar)\n",
    "print(post_mu, post_var)\n",
    "\n",
    "p_mu_post = stats.norm(post_mu, np.sqrt(post_var)).pdf(mus)\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.plot(mus, p_mu_post);\n",
    "ax.set_xlabel(\"$\\mu$\")\n",
    "ax.set_ylabel(\"$P(\\mu)$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5d10e13-c6eb-4a67-9ff0-30e06d40dd72",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1d2a605a-81cd-4499-b0c9-bc356e67f97f",
   "metadata": {},
   "source": [
    "# ポアソン分布のパラメータの事後分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dd040c01-5098-4122-b7a3-b56433651693",
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "167ab2d9-d19d-464f-a588-a1cc87bc0d6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42f97bb6-08bc-446a-be3c-57f13a744b3a",
   "metadata": {},
   "source": [
    "## トイデータ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3ebe1a99-e3ec-4b3a-b453-079027289bf5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAEGCAYAAABIAbBLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPKElEQVR4nO3dbaxlZ1nH4f9Np02hVBrpCamUYUo0jYQEiscqhTRYAra2oh/4QAMYiWaiAVKihhRjYkhM5BOCBjGTAqK8NFjahLSAEKlBVF5mSnlpCwabSlvADiKBkmil3n6YXTMM0/ZAn+esfc5cV3Iye++zu587K037m7XWXqu6OwAAjPGopQcAANhNxBUAwEDiCgBgIHEFADCQuAIAGGjP0gMc7cwzz+x9+/YtPQYAwMM6dOjQN7p749jX1yqu9u3bl4MHDy49BgDAw6qqfzve6w4LAgAMJK4AAAYSVwAAA4krAICBxBUAwEDiCgBgoGlxVVXnVtXNR/18u6pePWs9AIB1MO06V939pSTPSJKqOinJ3Umum7UeAMA62K7Dgs9L8q/dfdyLbQEA7BbbdYX2Fyd5z/F+UVX7k+xPkr17927TODDOvitvWHqEbXXH6y9degSAtTZ9z1VVnZLkhUn+5ni/7+4D3b3Z3ZsbGz9wex4AgB1lOw4LXpLkpu7+921YCwBgUdsRV5fnQQ4JAgDsNlPjqqpOS/L8JNfOXAcAYF1MPaG9u7+b5PEz1wAAWCeu0A4AMJC4AgAYSFwBAAwkrgAABhJXAAADiSsAgIHEFQDAQOIKAGAgcQUAMJC4AgAYSFwBAAwkrgAABhJXAAADiSsAgIHEFQDAQOIKAGAgcQUAMJC4AgAYSFwBAAwkrgAABpoaV1V1RlVdU1VfrKrbqupZM9cDAFjansmf/6YkH+ruF1XVKUkeM3k9AIBFTYurqnpckguT/HqSdPd9Se6btR4AwDqYeVjwnCSHk7y9qj5TVVdV1WnHvqmq9lfVwao6ePjw4YnjAADMNzOu9iR5ZpK3dPd5Sb6b5Mpj39TdB7p7s7s3NzY2Jo4DADDfzLi6K8ld3f3J1fNrciS2AAB2rWlx1d1fT3JnVZ27eul5SW6dtR4AwDqY/W3BVyV51+qbgrcnefnk9QAAFjU1rrr75iSbM9cAAFgnrtAOADCQuAIAGEhcAQAMJK4AAAYSVwAAA4krAICBxBUAwEDiCgBgIHEFADCQuAIAGEhcAQAMJK4AAAYSVwAAA4krAICBxBUAwEDiCgBgIHEFADCQuAIAGEhcAQAMJK4AAAYSVwAAA+2Z+eFVdUeS7yS5P8n3untz5noAAEubGlcrv9Dd39iGdQAAFuewIADAQLPjqpN8uKoOVdX+472hqvZX1cGqOnj48OHJ4wAAzDU7rp7T3c9MckmSV1TVhce+obsPdPdmd29ubGxMHgcAYK6pcdXdd6/+vCfJdUnOn7keAMDSpsVVVZ1WVac/8DjJC5J8YdZ6AADrYOa3BZ+Q5LqqemCdd3f3hyauBwCwuGlx1d23J3n6rM8HAFhHLsUAADCQuAIAGEhcAQAMJK4AAAYSVwAAA4krAICBxBUAwEDiCgBgIHEFADCQuAIAGEhcAQAMJK4AAAYSVwAAA4krAICBxBUAwEDiCgBgIHEFADCQuAIAGEhcAQAMJK4AAAaaHldVdVJVfaaqrp+9FgDA0rZjz9UVSW7bhnUAABY3Na6q6uwklya5auY6AADrYs/kz39jktckOf3B3lBV+5PsT5K9e/dOHgdY0r4rb1h6hG1zx+svXXoEYCHT9lxV1WVJ7unuQw/1vu4+0N2b3b25sbExaxwAgG0x87Dgs5O8sKruSHJ1kouq6p0T1wMAWNzDHhasqh9/qN939zcf5PXXJnnt6jOem+T3uvulP/yIAAA7x1bOubopyZOS/GeSSnJGkq+sftdJnjJlMgCAHWgrhwU/kuSXu/vM7n58ksuSfLi7z+nuLYVVd/99d1/2SAYFANgJthJXP9/dH3jgSXd/MMkF80YCANi5tnJY8KtV9QdJHjgZ/SVJvjpvJACAnWsre64uT7KR5Lok164eXz5zKACAneph91ytvg14RVWd1t3f3YaZAAB2rIfdc1VVF1TVrVndH7Cqnl5Vfz59MgCAHWgrhwX/JMkvJvmPJOnuzya5cOZQAAA71Zau0N7ddx7z0v0TZgEA2PG28m3BO6vqgiRdVScnuSKrQ4QAAHy/rey5+q0kr0jyxCR3J3nG6jkAAMd4yD1XVXVSkjd190u2aR4AgB3tIfdcdff9SZ5cVads0zwAADvaVs65uj3JP1bV+5P8/3WuuvsN06YCANihHnTPVVX99erhC5Ncv3rv6Uf9AABwjIfac/UzVfUTSb6S5M+2aR4AgB3toeLqL5L8XZJzkhw86vVK0kmeMnEuAIAd6UEPC3b3n3b3Tyd5e3c/5aifc7pbWAEAHMfDXuequ397OwYBANgNtnT7GwAAtkZcAQAMJK4AAAYSVwAAA02Lq6o6tao+VVWfrapbqup1s9YCAFgXW7n9zY/qv5Nc1N33VtXJST5eVR/s7k9MXBMAYFHT4qq7O8m9q6cnr3561noAAOtg5p6rVNVJSQ4l+ckkb+7uTx7nPfuT7E+SvXv3zhwnSbLvyhumr7FO7nj9pUuPAAAnlKkntHf3/d39jCRnJzm/qp52nPcc6O7N7t7c2NiYOQ4AwHTb8m3B7v5WkhuTXLwd6wEALGXmtwU3quqM1eNHJ3l+ki/OWg8AYB3MPOfqrCTvWJ139agk7+3u6yeuBwCwuJnfFvxckvNmfT4AwDpyhXYAgIHEFQDAQOIKAGAgcQUAMJC4AgAYSFwBAAwkrgAABhJXAAADiSsAgIHEFQDAQOIKAGAgcQUAMJC4AgAYSFwBAAwkrgAABhJXAAADiSsAgIHEFQDAQOIKAGAgcQUAMNC0uKqqJ1XVjVV1a1XdUlVXzFoLAGBd7Jn42d9L8rvdfVNVnZ7kUFV9pLtvnbgmAMCipu256u6vdfdNq8ffSXJbkifOWg8AYB1syzlXVbUvyXlJPrkd6wEALGXmYcEkSVU9Nsn7kry6u799nN/vT7I/Sfbu3Tt7HIC1t+/KG5YeYVvd8fpLf+R/1rbauhNpWz2S7TTC1D1XVXVyjoTVu7r72uO9p7sPdPdmd29ubGzMHAcAYLqZ3xasJG9Nclt3v2HWOgAA62TmnqtnJ3lZkouq6ubVzy9NXA8AYHHTzrnq7o8nqVmfDwCwjlyhHQBgIHEFADCQuAIAGEhcAQAMJK4AAAYSVwAAA4krAICBxBUAwEDiCgBgIHEFADCQuAIAGEhcAQAMJK4AAAYSVwAAA4krAICBxBUAwEDiCgBgIHEFADCQuAIAGEhcAQAMJK4AAAaaFldV9baquqeqvjBrDQCAdTNzz9VfJrl44ucDAKydaXHV3R9L8s1Znw8AsI4WP+eqqvZX1cGqOnj48OGlxwEAeEQWj6vuPtDdm929ubGxsfQ4AACPyOJxBQCwm4grAICBZl6K4T1J/jnJuVV1V1X9xqy1AADWxZ5ZH9zdl8/6bACAdeWwIADAQOIKAGAgcQUAMJC4AgAYSFwBAAwkrgAABhJXAAADiSsAgIHEFQDAQOIKAGAgcQUAMJC4AgAYSFwBAAwkrgAABhJXAAADiSsAgIHEFQDAQOIKAGAgcQUAMJC4AgAYSFwBAAw0Na6q6uKq+lJVfbmqrpy5FgDAOpgWV1V1UpI3J7kkyVOTXF5VT521HgDAOpi55+r8JF/u7tu7+74kVyf5lYnrAQAsrrp7zgdXvSjJxd39m6vnL0vyc939ymPetz/J/tXTc5N8acpAyzszyTeWHmIHsJ22zrbaOttqa2ynrbOttma3b6cnd/fGsS/uWWKSo3X3gSQHlp5jtqo62N2bS8+x7mynrbOtts622hrbaetsq605UbfTzMOCdyd50lHPz169BgCwa82Mq08n+amqOqeqTkny4iTvn7geAMDiph0W7O7vVdUrk/xtkpOSvK27b5m13g6w6w99DmI7bZ1ttXW21dbYTltnW23NCbmdpp3QDgBwInKFdgCAgcQVAMBA4moytwDamqp6W1XdU1VfWHqWdVdVT6qqG6vq1qq6paquWHqmdVRVp1bVp6rqs6vt9LqlZ1pnVXVSVX2mqq5fepZ1VlV3VNXnq+rmqjq49DzrrKrOqKprquqLVXVbVT1r6Zm2i3OuJlrdAuhfkjw/yV058g3Ky7v71kUHW0NVdWGSe5P8VXc/bel51llVnZXkrO6+qapOT3Ioya/69+r7VVUlOa27762qk5N8PMkV3f2JhUdbS1X1O0k2k/xYd1+29DzrqqruSLLZ3bv5wphDVNU7kvxDd1+1umrAY7r7WwuPtS3suZrLLYC2qLs/luSbS8+xE3T317r7ptXj7yS5LckTl51q/fQR966enrz68bfJ46iqs5NcmuSqpWdhd6iqxyW5MMlbk6S77ztRwioRV7M9McmdRz2/K/4nyEBVtS/JeUk+ufAoa2l1qOvmJPck+Uh3207H98Ykr0nyvwvPsRN0kg9X1aHV7ds4vnOSHE7y9tXh5quq6rSlh9ou4gp2qKp6bJL3JXl1d3976XnWUXff393PyJE7RJxfVQ45H6OqLktyT3cfWnqWHeI53f3MJJckecXqlAZ+0J4kz0zylu4+L8l3k5ww5x2Lq7ncAogpVucQvS/Ju7r72qXnWXerwxE3Jrl44VHW0bOTvHB1LtHVSS6qqncuO9L66u67V3/ek+S6HDn9gx90V5K7jtpbfE2OxNYJQVzN5RZADLc6UfutSW7r7jcsPc+6qqqNqjpj9fjROfLFki8uOtQa6u7XdvfZ3b0vR/4b9dHufunCY62lqjpt9SWSrA5xvSCJbzgfR3d/PcmdVXXu6qXnJTlhvnQz7fY3uAXQD6Oq3pPkuUnOrKq7kvxhd7912anW1rOTvCzJ51fnEyXJ73f3B5YbaS2dleQdq2/tPirJe7vbZQZ4JJ6Q5Lojf7/JniTv7u4PLTvSWntVknetdi7cnuTlC8+zbVyKAQBgIIcFAQAGElcAAAOJKwCAgcQVAMBA4goAYCBxBQAwkLgCABhIXAG7UlX9bFV9rqpOXV1Z+xb3FgS2g4uIArtWVf1RklOTPDpH7nP2xwuPBJwAxBWwa61uu/HpJP+V5ILuvn/hkYATgMOCwG72+CSPTXJ6juzBApjOnitg16qq9ye5Osk5Sc7q7lcuPBJwAtiz9AAAM1TVryX5n+5+d1WdlOSfquqi7v7o0rMBu5s9VwAAAznnCgBgIHEFADCQuAIAGEhcAQAMJK4AAAYSVwAAA4krAICB/g9NjWYsv5dyOgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu_set = 3.5\n",
    "n = 20\n",
    "sample = stats.poisson(mu_set).rvs(n)\n",
    "max_sample = sample.max()\n",
    "count = collections.Counter(sample)\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.bar(count.keys(), count.values());\n",
    "ax.set_xlabel(\"x\");\n",
    "ax.set_ylabel(\"freq\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b0e306c-9a5f-4577-962a-f5bb2211353a",
   "metadata": {},
   "source": [
    "## 尤度関数\n",
    "\n",
    "$$\n",
    "L(\\lambda) = \\prod^n_{i=1} \\mathrm{Poi}(x | \\lambda) = \\prod_{i=1}^{n} \\frac{\\lambda^{x_i} e^{-\\lambda}}{x_i!}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "20c009ef-80de-4598-9427-9e1bbae1b98b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood_poi(lam, data):\n",
    "    ps = [stats.poisson(lam).pmf(x) for x in data]\n",
    "    return np.prod(ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "85b524f3-d553-49ad-83e5-1192f6242ca9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAloAAAETCAYAAADu5KjTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAs4UlEQVR4nO3deXiV5Z0+8Pt7tuwL2ckCYY8EJECMIC6AomBF0FErVtufXbTt2LH7MtOZTjvTuabTmVanrdNibdXaWvcqilpWFWULeyAkQAiQfSN7Qs7y/f2RxCIGCJDzPme5P9eVKycnh/PcXkcO93me931eUVUQERER0cizmQ5AREREFKpYtIiIiIj8hEWLiIiIyE9YtIiIiIj8hEWLiIiIyE9YtIiIiIj8JKCLloj8TkQaRKRkhJ7vLRFpFZHXz7j/ehHZKSK7RWSTiEwcifGIiIgovAV00QLwJIDFI/h8PwVw3xD3/x+AT6lqAYA/Afj+CI5JREREYSqgi5aqvgug5fT7RGTCwMzUDhF5T0TyLuD51gHoGOpXAOIHbicAqLnYzERERESDHKYDXISVAL6oqodE5EoAjwFYeInP+XkAq0WkB0A7gDmX+HxEREREwVW0RCQWwFUAXhCRwbsjBn53O4AfDfHHqlX1pvM89dcA3KyqW0XkWwB+hv7yRURERHTRgqpooX+ps3XgWKqPUNWXAbx8oU8oIqkAZqjq1oG7ngPw1qWEJCIiIgIC/BitM6lqO4CjInInAEi/GZf4tCcBJIjI5IGfFwEovcTnJCIiIoKoqukMZyUizwKYDyAFQD2AHwBYj/6zBEcDcAL4s6oOtWQ41PO9ByAPQCyAZgCfU9W3ReQ29C87+tBfvD6rqhUj+19DRERE4SagixYRERFRMAuqpUMiIiKiYBKwB8OnpKRobm6u6RhERERE57Vjx44mVU098/6ALVq5ubkoLi42HYOIiIjovETk2FD3c+mQiIiIyE9YtIiIiIj8hEWLiIiIyE9YtIiIiIj8hEWLiIiIyE9YtIiIiIj8hEWLiIiIyE8Cdh8tIqJg5vb68EJxFU529yHSaUeU045olx2RTjsKc0chJTbCdEQisgCLFhHRCCutbcc3X9iD/TXtQ/5+VLQT/3XHDCyamm5xMiKyGosWEdEIcXt9+L+NR/CL9YeQEOXEr++djYV5aejp86LH7UV3nwdNnX344ar9+MLTxbhvzlj80ycuQ6TTbjo6EfkJixYR0Qg4fRZr6YxM/PDWfCTFuAAALocNCXACAManAi9/+Sr89K0y/HbTUWw72oL/XTETUzLiTMYnIj/hwfBERJdoz4lWLPvV+6hv78Wv752FX6yY+WHJGkqEw47v3zIVT95/BZq7TmHpLzfhpR1VFiYmIquwaBERXYKOXje+8uwupMS48ObD12LxtNHD/rPzp6ThzYevxewxo/Cdl/Zi5/GTfkxKRCawaBERXSRVxT+9UoKqk914dMVMpMZd+JmEqXER+PV9s5GREImv/GkX2rrdfkhKRKawaBERXaQXdlThtT01+NoNk3FFbtJFP09ClBO/vGcW6tt78a0X90BVRzAlEZlkSdESkSkisvu0r3YR+aoVYxMR+cPhhg784NX9mDs+GV9eMPGSn68gJxHfXZKHvx6ox5MfVF56QCIKCJacdaiqZQAKAEBE7ACqAbxixdhERCOt1+3FQ3/ahSiXHY/cXQC7TUbkeT939ThsPtKM/1hditljR+Hy7MQReV4iMsfE0uH1AI6o6jEDYxMRXbL/WF2Kg3Ud+O87L0d6fOSIPa+I4L/vnIHU2Ag89KddaO/l8VpEwc5E0bobwLND/UJEHhCRYhEpbmxstDgWEdH5bT7SjKc3H8Pnrx6HhXkjv7P7qBgX/nfFTFS39uD7r5SM+PMTkbUsLVoi4gJwK4AXhvq9qq5U1UJVLUxNTbUyGhHRsPx8bTnS4iLwzZum+G2MwtwkPLRgIl7bU4Mdx1r8Ng4R+Z/VM1pLAOxU1XqLxyUiumSbjzRj29EWfGn+BL9fNufB68YjJTYCP3mzjGchEgUxq4vWCpxl2ZCIKNA9uq5/NmtF0Ri/jxXtcuDh6ydiW2ULNpbxUAqiYGVZ0RKRGACLALxs1ZhERCNlS0UztlS04IvX+X82a9DdRWMwNjkaP3nrIHw+zmoRBSPLipaqdqlqsqq2WTUmEdFIeXTtIaTGReCeK/0/mzXIabfh64sm42BdB17bU2PZuEQ0crgzPBHReWw72oLNFc2WzmYNWnp5JvIz4/E/a8rQ5/FZOjYRXToWLSKi83h0XTlSYiPwKQtnswbZbIJvL87DiZYePLvtuOXjE9GlYdEiIjqH7ZUteP9wM7543XjLZ7MGXTspBXPGJ+EX6w+h65THSAYiujgsWkRE5/Do2kNIiXXhU1eONZZBRPCdxXlo6uzDE5uOGstBRBeORYuI6Cx2Hj+JTYeb8OC1ExDlMjObNWjmmFG4KT8dK9+t4KV5iIIIixYR0Vk8s/kY4iIc+NQc64/NGspDCyah85QHz28/YToKEQ0TixYR0RDaetxYXVKLZTMzEe1ymI4DAJienYDCsaPw1OZKeLmvFlFQYNEiIhrCa3tq0Ov24ZOFgTGbNej+eeNwoqUH60p5JTOiYMCiRUQ0hOe2H8fU0fGYlhVvOspH3JSfjsyESPz+/UrTUYhoGFi0iIjOUFLdhpLqdtxdlAMRMR3nIxx2Gz59VS42VzSjtLbddBwiOg8WLSKiMzxffAIuhw3LZmSZjjKku6/IQaTThic5q0UU8Fi0iIhO0+v24pVd1bh5WgYSop2m4wwpMdqF22dl4y+7q9HS1Wc6DhGdA4sWEdFp3iqpQ0evB5+8IrAOgj/T/Vfl4pTHx8vyEAU4Fi0iotP8eftxjE2OxpzxSaajnNOk9DhcMykFf9h8DG4vLzZNFKhYtIiIBlQ2dWFLRQvuKgy8g+CHcv+8XNS19+LNkjrTUYjoLFi0iIgGPF98AjYB7pidbTrKsMyfnIZxKTH4/fu8/iFRoGLRIiIC4PH68MKOKizMS0N6fKTpOMNiswk+M3csdh1vRUl1m+k4RDQEy4qWiCSKyIsiclBESkVkrlVjExGdz4ayRjR2nAr4g+DPdNusbEQ4bHi+mNc/JApEVs5oPQrgLVXNAzADQKmFYxMRndOru6uRHOPCgimppqNckIQoJxZPy8Cru2vQ6/aajkNEZ7CkaIlIAoBrATwBAKrap6qtVoxNRHQ+PX1erCttwE3TMuCwB98RFXcV5qCtx401B3j9Q6JAY9U7yjgAjQB+LyK7ROS3IhJj0dhEROe0sawBPW4vPjF9tOkoF2Xu+GRkJUZx+ZAoAFlVtBwAZgH4P1WdCaALwHfPfJCIPCAixSJS3NjYaFE0Igp3b+yrRXKMC1eOC+y9s87GZhPcMTsbmw43obq1x3QcIjqNVUWrCkCVqm4d+PlF9Bevj1DVlapaqKqFqanBdZwEEQWnXrcX6w8G77LhoDtmZ0MVeGlHlekoRHQaS95VVLUOwAkRmTJw1/UADlgxNhHRuWwsa0B3X/AuGw7KSYrGvInJeGHHCfh8ajoOEQ2w8uPbVwD8UUT2AigA8B8Wjk1ENKQ39tUhKYiXDU93V2EOTrT0YMvRZtNRiGiAw6qBVHU3gEKrxiMiOp9etxfrSuuxrCArqJcNB92Un4G4SAdeKK7CVRNSTMchInBneCIKYxvLGkNi2XBQpNOOZQWZWL2vFu29btNxiAgsWkQUxlbvq0VSjAtzxgf/suGguwpzcMrjw6o9NaajEBFYtIgoTA0uG96Unx4Sy4aDpmclIC8jDi8U8+xDokAQOu8uREQXYGNZI7r6vLg5RJYNB4n076m1+0Qryus7TMchCnssWkQUllbvq8WoaCfmjk82HWXELZ+ZBbtN8Jdd1aajEIU9Fi0iCjt/WzYM7k1KzyYlNgJXT0zBq7truKcWkWGh9w5DRHQe75SH5rLh6ZbPzER1aw+Kj500HYUorLFoEVHYeaukDonRTsydEHrLhoNunJqBKKcdf9nN5UMik1i0iCiseLw+rD/YgIV5aXCG4LLhoJgIB27MT8fqfbXo8/hMxyEKW6H7LkNENIQdx06irceNRZelm47id8sLstDa7cY75Y2moxCFLRYtIgora0vr4bLbcM3kVNNR/O7qSSlIinFx+ZDIIBYtIgor60obcOX4JMRGWHapV2OcdhtuuXw01h6oRwcvyUNkBIsWEYWNI42dqGjqwqKpob9sOGhZQRZOeXx4e3+96ShEYYlFi4jCxrrS/rKxMC/NcBLrzBqTiDFJ0XiVy4dERrBoEVHYWFvagMtGxyN7VLTpKJYRESwryMT7h5vQ0NFrOg5R2GHRIqKwcLKrD8WVLbjhsvCZzRq0rCALPgVW7ak1HYUo7LBoEVFY2FjeAJ8CN4TBtg5nmpgWi2lZ8Vw+JDKARYuIwsLaAw1IjYvA9KwE01GMWF6Qhb1Vbaho7DQdhSisWFa0RKRSRPaJyG4RKbZqXCKiPo8P75Q34vq8NNhsYjqOEUtnZEIEeHV3jekoRGHF6hmtBapaoKqFFo9LRGFs29EWdJ7yhOWy4aD0+EjMGZeMVXtroKqm4xCFDS4dElHIW1tajwiHDfMmppiOYtTSGZmoaOzC/pp201GIwoaVRUsB/FVEdojIA0M9QEQeEJFiESlubOS1uYjo0qkq1pbW45pJKYhy2U3HMWrJtAw4bIJVe7l8SGQVK4vW1ao6C8ASAH8vItee+QBVXamqhapamJoa+tchIyL/K6vvQNXJHlwfxsuGg0bFuHDNpBS8vqeWy4dEFrGsaKlq9cD3BgCvACiyamwiCl/rShsAANeH0W7w57J0RiaqW3uw8/hJ01GIwoIlRUtEYkQkbvA2gBsBlFgxNhGFt7Wl9ZiRnYC0+EjTUQLCoqnpiHDYuHkpkUWsmtFKB7BJRPYA2AbgDVV9y6KxiShMtXT1YfeJVizgbNaH4iKdWJiXhtf31sLj9ZmOQxTyHFYMoqoVAGZYMRYR0aB3yxuhCiyYwqJ1ultnZOLNkjpsPdoS9mdiEvkbt3cgopC1oawByTGusN0N/mwW5KUhxmXHqj08+5DI31i0iCgkeX2Kd8obcd2U1LDdDf5sIp123JifgTdL6tDn4fIhkT+xaBFRSNp9ohWt3W4s5PFZQ7p1Ribaetx47xD3LCTyJxYtIgpJG8saYLcJrpnIPfmGMm9iChKjnXiNy4dEfsWiRUQhaf3BBsweMwoJ0U7TUQKSy2HDkmkZWHOgHj19XtNxiEIWixYRhZz69l7sr2nH/DzOZp3L0hmZ6O7zYt3BetNRiEIWixYRhZx3yvqPO+K2Dud25bhkpMVF8OxDIj9i0SKikLOhrAEZ8ZHIy4gzHSWg2W2CT1w+GhvKGtHe6zYdhygksWgRUUhxe31471ATFuSlQoTbOpzPrTMy0efx4e2SOtNRiEISixYRhZTiypPoPOXBfC4bDktBTiJykqKwai+vfUjkDyxaRBRSNpY1wGkXXM1LywyLiGDp5Zl4/3ATmjtPmY5DFHJYtIgopGwoa8CV45IRE2HJpVxDwq0FmfD6FKv3cVaLaKSxaBFRyKg62Y3y+k7Mn8JtHS5EXkY8JqfHcvNSIj9g0SKikLFhcFsHXnbngi29PBPbK0+iprXHdBSikMKiRUQhY+PBBoxJisb4lBjTUYLO0hmZAIDX93JWi2gksWgRUUjodXvx/pEmzJ/CbR0uRm5KDGZkJ3D5kGiEsWgRUUjYerQFvW4flw0vwdIZmSipbkdFY6fpKEQhg0WLiELChoMNiHTaMHd8sukoQeuWyzMhAqzaw7MPiUaKpUVLROwisktEXrdyXCIKbaqKDWUNuGpCCiKddtNxglZGQiSKcpPw2p5qqKrpOEQhweoZrYcBlFo8JhGFuKNNXTjW3I0F3Nbhkt1akIkjjV04UNtuOgpRSLCsaIlINoBPAPitVWMSUXgY3NaBl925dEumjYbDJjwonmiEWDmj9QiAbwPwne0BIvKAiBSLSHFjY6NlwYgouG0sa8DEtFjkJEWbjhL0kmJcuHZyKlbtroHPx+VDokt1wUVLRGJE5IIOghCRWwA0qOqOcz1OVVeqaqGqFqamcgmAiM6v65QHWytasJBnG46Y5TOzUNPWi22VLaajEAW98xYtEbGJyD0i8oaINAA4CKBWRA6IyE9FZOIwxpkH4FYRqQTwZwALReSZS0pORATg/cNN6PP6eNmdEbTosnTEuOz4y65q01GIgt5wZrQ2AJgA4HsAMlQ1R1XTAFwNYAuAn4jIved6AlX9nqpmq2ougLsBrFfVc/4ZIqLh2FDWiNgIBwrHJpmOEjKiXHbcNC0Db+yrRa/bazoOUVAbTtG6QVX/TVX3quqHx1epaouqvqSqfwfgOf9FJCIamqpiY1kDrp6YApeD2wKOpNtmZqGj14MNBxtMRyEKaud9Z1JV90g85rTHblTVW4b7eCKisymr70BtWy8W5HHZcKRdNSEFqXEReIXLh0SXxHEhDxaRHAD5AKYBmA4gX1UL/RGMiOh8Nhzktg7+YrcJls3IxFObK9Ha3YfEaJfpSERBaTgHwz8oIh+ISCuAcgCfBxAL4DUA9/g3HhHR2W042ID8zHikx0eajhKSls/MgturWL2vznQUoqA1nIMavgfgawBmA3gdQCSA3w0cn1Xuz3BERGfT1u3GjuMnsYCzWX6TnxmPiWmxPPuQ6BIMp2jdoqpbVfWIqt4J4FcAVonI10SER58SkRHvHW6E16c8PsuPRAS3zczCtsoWnGjpNh2HKCgN52D4kjN+fhNAEYAkAO/7KRcR0TltONiIxGgnCnJGmY4S0m6dkQkAvCQP0UUazjFacuZ9qnpKVf8ZwGfO9hgiIn/x+RTvlDfgusmpsNv49uNPOUnRKMpNwiu7qqHKS/IQXahhbVgqIl8RkTGn3ykiLgDZIvIUBgoXEZEVdle1oqmzj5fdscjymVk43NCJ/TXtpqMQBZ3hFK3FALwAnhWRmoFL71QAOARgBYBHVPVJP2YkIvqIdaX1sNsE8yezaFnh5ukZcNqFe2oRXYTz7qOlqr0AHgPwmIg4AaQA6FHVVj9nIyIa0rrSBlyROwoJ0U7TUcJCYrQLC/PS8Oruanx3SR6cdp4HRTRcF/S3RVXdqlo7WLJEhAfDE5GlTrR042BdB264LN10lLBy5+wcNHX2YWNZo+koREHlUj+WZI5ICiKiYVpXWg8AuJ5Fy1Lzp6QiJTYCLxSfMB2FKKgM56zDX4jIAyIyV0Tizvg1T0EhIkutO9iACakxGJcSYzpKWHHYbbh9VhbWH2xAU+cp03GIgsZwZrT2of+6hv8JoFJEjorIayLyYwBnFi8iIr/p6HVjS0Uzlw0NuXN2Njw+5U7xRBdgOAfDrzz9ZxHJRn/xuhzA237KRUT0Me+WN8HtVdwwlUXLhEnpcSjIScTzxSfwuavHgVsoEp3fBR+jpapVqvqmqv5EVe/1RygioqGsK63HqGgnZo3hbvCm3FmYjfL6TuytajMdhSgo8BxdIgoKHq8P68sasGBKGneDN2jpjExEOGx4YQcPiicaDhYtIgoKO4+3orXbzWVDw+IjnVgyLQOv7q5Br9trOg5RwGPRIqKgsK60Hk674JpJKaajhL07C3PQ0evB2/vrTEchCniWFC0RiRSRbSKyR0T2i8gPrRiXiELHmtJ6zBmfjLhI7gZv2tzxychKjMKLO6pMRyEKeFbNaJ0CsFBVZwAoALBYROZYNDYRBbmjTV2oaOzitg4BwmYT3DE7G5sON6G6tcd0HKKAZknR0n6dAz86B7642SkRDcvfdoPnRaQDxR2zs6EKvMRZLaJzsuwYLRGxi8huAA0A1qjq1iEe84CIFItIcWMjr6dFRP3WHKhHXkYcskdFm45CA3KSonHVhGQ8X3wCPh8/NxOdjWVFS1W9qloAIBtAkYhMG+IxK1W1UFULU1NTrYpGRAGstbsPxcdOctkwAK0oGoOqkz145xA/GBOdjeVnHapqK4ANABZbPTYRBZ81B+rh9SluzGfRCjQ35WcgJTYCz2w+ZjoKUcCy6qzDVBFJHLgdBWARgINWjE1Ewe3NkjpkJUZhelaC6Sh0BpfDhhVFOVhf1oATLd2m4xAFJKtmtEYD2CAiewFsR/8xWq9bNDYRBan2XjfeO9SIm6dn8Lp6AWpF0RgIgGe3HTcdhSggnfei0iNBVfcCmGnFWEQUOtaXNsDtVSyeNtp0FDqLzMQoXH9ZOp7bfgIP3zAJEQ676UhEAYU7wxNRwFq9rxYZ8ZGYmZNoOgqdw71zxqK5qw9vlXCneKIzsWgRUUDqPOXBxvJGLJ6WARsvIh3QrpmYgrHJ0fgDD4on+hgWLSIKSBsONqDP48OSaRmmo9B52GyCe68ci+JjJ1Fa2246DlFAYdEiooD0VkkdUmIjUJibZDoKDcMds7PhctjwzBbOahGdjkWLiAJOT58X6w82YPG0dNi5bBgURsW4sPTyTLyyqxodvW7TcYgCBosWEQWcd8ob0OP24maebRhU7ps7Ft19XvxlV7XpKEQBg0WLiALO6n11GBXtRNE4LhsGkxnZCZiWFY+nNx+DKq9/SASwaBFRgOl19y8b3pSfAYedb1HBRETw6bm5ONTQiXcPNZmOQxQQ+C5GRAFl06EmdJ7yYMl0LhsGo2UFmUiLi8DKd4+YjkIUEFi0iCigrC6pRXykA3PHJ5uOQhchwmHH/fPG4f3DzSipbjMdh8g4Fi0iChh9Hh/WHqjHoqkZcDn49hSs7rlyDGJcdqx8t8J0FCLj+E5GRAFj0+FGtPd6cPN0blIazBKinFhRNAZv7KtF1clu03GIjGLRIqKA8dLOaiTFuHDNpFTTUegSffbqcRAAT2w6ajoKkVEsWkQUENp63FhzoB63zsjksmEIyEyMwtIZmXhu+wm0dXMDUwpffDcjooCwel8t+jw+3D4ry3QUGiFfuGY8uvu8eGYrL8tD4YtFi4gCwss7qzAxLRbTsxJMR6ERMjUzHtdOTsXv369Er9trOg6RESxaRGTcseYubK88idtnZUGE1zYMJQ9eOx5Nnad4WR4KWyxaRGTcK7uqIQIsL+CyYai5akIy8jPjsfK9Cvh8vCwPhR9LipaI5IjIBhE5ICL7ReRhK8YlosCnqnh5ZzWumpCMzMQo03FohIkIHrxuAioau7C6pNZ0HCLLWTWj5QHwDVWdCmAOgL8XkakWjU1EAWzHsZM43tKN22dmm45CfvKJ6aMxOT0WP1tTDo/XZzoOkaUsKVqqWquqOwdudwAoBcA1AiLCSzurEeW0Y/E0blIaquw2wdcXTUZFYxde3V1jOg6RpSw/RktEcgHMBLB1iN89ICLFIlLc2NhodTQisliv24vX99ZgybQMxEQ4TMchP7opPwPTsuLxyLpy9Hk4q0Xhw9KiJSKxAF4C8FVVbT/z96q6UlULVbUwNZU7QxOFunWlDejo9eD2WVw2DHUigm/cOAUnWnrwwo4TpuMQWcayoiUiTvSXrD+q6stWjUtEgevlnVXIiI/E3AnJpqOQBeZPTsXssaPwi3WHua8WhQ2rzjoUAE8AKFXVn1kxJhEFtqbOU9hY3ojlM7Ngt3HvrHDQP6s1GXXtvfjT1uOm4xBZwqoZrXkA7gOwUER2D3zdbNHYRBSAXtpRBa9P8Xe85E5YuWpCCq6akIzHNh5Gd5/HdBwiv7PqrMNNqiqqermqFgx8rbZibCIKPF6f4unNx3DluCRMSo8zHYcs9o0bp6Cpsw9PflBpOgqR33FneCKy3NrSelS39uD/XZVrOgoZMHvsKCzMS8Nv3qlAe6/bdBwiv2LRIiLLPfVBJTITIrFoarrpKGTI1xdNRnuvG/+79pDpKER+xaJFRJYqq+vAB0eace/csXDY+RYUrqZlJeDuK8bg9x9Uory+w3QcIr/huxwRWeqpzZWIcNhw9xVjTEchw7510xTERjjwr6/thyovOE2hiUWLiCzT1u3GKzursawgE0kxLtNxyLCkGBe+eeNkfHCkGav31ZmOQ+QXLFpEZJnni0+gx+3FZ3gQPA2458qxmDo6Hv/+xgFu90AhiUWLiCzh9Sme3lKJotwk5GcmmI5DAcJuE/xoWT5q23rx2IYjpuMQjTgWLSKyxPqDDTjR0sPZLPqYwtwk3D4zCyvfrUBlU5fpOEQjikWLiCzx1AeVGJ0QiRvzuaUDfdx3l+TB5bDhR68fMB2FaESxaBGR3x2q78Cmw024d85YOLmlAw0hLT4SX71hEtYfbMBbJTwwnkIH3/GIyO9+827FwJYOOaajUAD7zFW5yM+Mxz+9sg9NnadMxyEaESxaRORXRxo78fLOKtw3ZyySYyNMx6EA5rTb8PNPFqDjlAfffWkf99aikMCiRUR+9cjaQ4h02vHF+RNMR6EgMDk9Dt++aQrWltbj+eITpuMQXTIWLSLym9LadqzaU4P75+UihbNZNEyfnTcOc8cn40erDuB4c7fpOESXhEWLiPzm52vKERfpwAPXcDaLhs9mE/z3XTNgE8HXn98Nr49LiBS8WLSIyC/2VrXirwfq8YVrxiMh2mk6DgWZrMQo/HBZPoqPncTKdytMxyG6aCxaROQX//PXcoyKduL+ebmmo1CQum1mFm6enoGfrSlDSXWb6ThEF4VFi4hG3PbKFrxT3ogvXjcBcZGczaKLIyL48fLpSIpx4YvP7EAzt3ygIGRJ0RKR34lIg4iUWDEeEZmjqvjp22VIjYvAp+fmmo5DQW5UjAu/ua8QDR2n8KU/7kSfx2c6EtEFsWpG60kAiy0ai4gMev9wM7YdbcFDCyYiymU3HYdCQEFOIv7r7y7HtqMt+NdV+7m/FgUVS4qWqr4LoMWKsYjInD6PDz96fT+yEqNwdxF3gaeRs3xmFr543QT8aetxPLPlmOk4RMMWUMdoicgDIlIsIsWNjY2m4xDRBXps42GU13fi35bnI8LB2SwaWd+6aQquz0vDv646gA8ON5mOQzQsAVW0VHWlqhaqamFqaqrpOER0AcrrO/CrDYexrCATC/PSTcehEGS3CR65uwDjUmLw5T/txLHmLtORiM4roIoWEQUnr0/x7Rf3IjbCgX+5ZarpOBTC4iKd+O2nCwEA9z2xDTWtPYYTEZ0bixYRXbKnPqjE7hOt+MHSfF44mvwuNyUGT95fhJNdfVjx+BbUtrFsUeCyanuHZwFsBjBFRKpE5HNWjEtE/neipRs/fbsM86ekYllBpuk4FCYKchLx9OeK0NzZh3se34q6tl7TkYiGZNVZhytUdbSqOlU1W1WfsGJcIvIvVcU/vrIPNgF+fNt0iIjpSBRGZo4Zhac+W4SG9l7c8/gWNLSzbFHg4dIhEV20F3dU4b1DTfjOkjxkJUaZjkNhaPbY/rJV196Lux/fgoYOli0KLCxaRHRRSqrb8C+v7kdRbhLuvXKs6TgUxgpzk/Dk/UWoa+vFnb/ejMMNHaYjEX2IRYuILlhDey++8HQxRkU78ctPzYTNxiVDMqtoXBKe+fyV6DrlxW2/+gAbDjaYjkQEgEWLiC5Qr9uLLzxdjLYeNx7/TCHS4iJNRyICAMwaMwqvPjQPOUnR+OxT2/H4uxW8XA8Zx6JFRMOmqvjGC3uwt7oNj3yyAPmZCaYjEX1EVmIUXvzSXCzOz8CPV5fiWy/uxSmP13QsCmMsWkQ0bI+sPYQ39tbiO4vzcGN+huk4REOKdjnwq3tm4R+un4QXd1Thrt9swZHGTtOxKEyxaBHRsLy2pwaPrjuEO2Zn48Frx5uOQ3RONpvg64sm47FPzUJlUxdufvQ9PP5uBbw+LiWStVi0iOi8ntt+HF97bjeKcpPw49umcb8sCho3Tx+NNV+7FtdMSsGPV5firt9sRgVnt8hCLFpEdFaqip+vKcd3XtqHeRNT8Lv7r0CEw246FtEFSYuPxOOfLsTP7pqBQ/UdWDIwu9Xn8ZmORmGARYuIhuT2+vCdl/Z+uFz4xGcKERvhMB2L6KKICG6flY01X78OV0/sn926/mcb8Zdd1fBxOZH8iEWLiD6m65QHn3+qGM8XV+Efrp+En95xOZx2vl1Q8EuPj8RvP1OIJ++/ArERTnz1ud34xC82YUNZA7eCIL+QQP0fq7CwUIuLi03HIAo7+2va8M0X9qK8vgP/vnwaVhSNMR2JyC98PsWqvTX4n7+W43hLN4pyk/D5a8bh+svSYecmvHSBRGSHqhZ+7H4WLSICgO4+Dx5ZewhPbDqKUdFO/PTOGVgwJc10LCK/6/P48Nz243hs4xHUtvUiKzEKn5ozBndfMQZJMS7T8ShIsGgR0VltLGvA9/9SgqqTPVhRlIPvLr4MCdFO07GILOXx+rC2tB5PfXAMmyua4XLYcMvlo7G8IAtzJyRz+ZzO6WxFi0e2EoWx0tp2/HLDYbyxtxYTUmPw/INzUTQuyXQsIiMcdhsWTxuNxdNGo7y+A3/YfAyv7KrGyzurkRDlxI1T03Hz9NGYNzEFLgdLFw0PZ7SIwozXp1hzoB6/f/8oth5tQaTThgevnYAvL5jArRuIztDr9uK9Q014c18t1hyoR8cpD+IiHJgzIRlzxydj7oRkTEmP44XViTNaROHuWHMX3iqpw9Obj6G6tQdZiVH43pI8fPKKHCRG8zgUoqFEOu1YNDUdi6am45THiw8ON+Pt/XX44Egz1hyoBwCMinZizvhkFOQkIj8zAfmZ8RjFY7toAIsWUYjqdXuxuaIZ75Q14p3yRhxt6gIAFI1Lwj/fchluuCwdDh5zQjRsEQ47FuSlYUFe/0kiVSe7saWiBZuPNGPr0Wa8WVL34WOzEqMwNTMeU9LjkJsSg3Ep0chNjkFSjItXVggzli0dishiAI8CsAP4rar+57kez6VDouFRVbR09aGsrgOldR04WNuOsvoOHKzrQJ/HhwiHDXMnJGP+5FTMn5KG3JQY05GJQlJLVx8O1LSjpKYN+2vasb+mDceauz9yfcX4SAdykqIxOiESGQmRGJ0QhYz4SKTFR2BUtAvJsS4kxbi4jB+EjJ51KCJ2AOUAFgGoArAdwApVPXC2P8OiReHG7fWh1+3FKU//9163D919HnT0etDe4+7/3utGS1cf6tp7UdfWi7r2XtS39aKrz/vh86TERiAvIw5TM+Mxb2IKrhyXhEgn37SJTHB7fTjR0o3K5i4cberG0aZOVJ/sQe3A39/WbveQfy7GZUditAvxUU7ERToQH+lAfKQTsZEORLnsiHY6EO2y99922RHhsCPCYUOE04YIhx0uhw1Ou8Blt8Fpt8ExcNtuEzhsNtjtAodNYLcJ7CI8xmwEmD5GqwjAYVWtGAjzZwDLAJy1aPnbg38oRkVjl6nhKQAM9yPGUB9G9GM3+m+q6sB3QKH93xXwqQ589T/G41N4vf3fPT4fPL7+xw6HwyZIj+//NHxZRjzmT05D1qgoTEmPw5SMOKTGRQzzv4yI/M1pt2F8aizGp8YO+fuePi/q2nvR0N6Lk919aOly42R3H5o7+9Da3Yf2Xg86et2oae3Fwd4OdJ7yoLvP65frNNoEsNsENhn8AmwiEAFsA/cLAJH+Sxp9eBsy8B0fLosOro4O/v70+wBg8Obpy6jysRsff/xwnLk0G+Gw4Y1/uOYCnmFkWVW0sgCcOO3nKgBXnvkgEXkAwAMAMGaMf3ejHpscw51/6cM3gGE88Kx3nflGcfobzuD3wTcwGXjTcgx8qnTYZeATZv+nzUinHRFOGyId/d+jXQOfZAc/1UY5Eety8NMnUYiIctkxLiUG4y5wSd/j9aHH7UVPnxc9AzPhp9w+nPIM3PZ40edRuL0+eHw+uD2KPq8PXl//BzzfwHeP1wfvwIdAn0/7b/tO/2DY/0FRtf93/R8iBz5MfvjzGfcPfgLVv30WPf0D69/uwxD3neOD7XAM8WCn3ez7ZUAdDK+qKwGsBPqXDv051j/efJk/n56IiMhvHHYb4uw2xEVyY+FAZ9UpR9UAck77OXvgPiIiIqKQZVXR2g5gkoiMExEXgLsBvGbR2ERERERGWLJ0qKoeEXkIwNvo397hd6q634qxiYiIiEyx7BgtVV0NYLVV4xERERGZxm2hiYiIiPyERYuIiIjIT1i0iIiIiPyERYuIiIjITyy7qPSFEpFGAMf8PEwKgCY/j0EXhq9JYOLrEnj4mgQmvi6Bx6rXZKyqpp55Z8AWLSuISPFQF4Akc/iaBCa+LoGHr0lg4usSeEy/Jlw6JCIiIvITFi0iIiIiPwn3orXSdAD6GL4mgYmvS+DhaxKY+LoEHqOvSVgfo0VERETkT+E+o0VERETkNyxaRERERH4SlkVLRBaLSJmIHBaR75rOQ4CI/E5EGkSkxHQW6iciOSKyQUQOiMh+EXnYdCYCRCRSRLaJyJ6B1+WHpjNRPxGxi8guEXnddBbqJyKVIrJPRHaLSLGRDOF2jJaI2AGUA1gEoArAdgArVPWA0WBhTkSuBdAJ4GlVnWY6DwEiMhrAaFXdKSJxAHYAWM6/K2aJiACIUdVOEXEC2ATgYVXdYjha2BORrwMoBBCvqreYzkP9RQtAoaoa20Q2HGe0igAcVtUKVe0D8GcAywxnCnuq+i6AFtM56G9UtVZVdw7c7gBQCiDLbCrSfp0DPzoHvsLrE3MAEpFsAJ8A8FvTWSiwhGPRygJw4rSfq8B/PIjOSURyAcwEsNVwFMKHS1S7ATQAWKOqfF3MewTAtwH4DOegj1IAfxWRHSLygIkA4Vi0iOgCiEgsgJcAfFVV203nIUBVvapaACAbQJGIcLndIBG5BUCDqu4wnYU+5mpVnQVgCYC/HzhMxVLhWLSqAeSc9nP2wH1EdIaBY4BeAvBHVX3ZdB76KFVtBbABwGLDUcLdPAC3DhwP9GcAC0XkGbORCABUtXrgewOAV9B/+JClwrFobQcwSUTGiYgLwN0AXjOciSjgDBx0/QSAUlX9mek81E9EUkUkceB2FPpP7DloNFSYU9XvqWq2quai/9+U9ap6r+FYYU9EYgZO5IGIxAC4EYDlZ7aHXdFSVQ+AhwC8jf6De59X1f1mU5GIPAtgM4ApIlIlIp8znYkwD8B96P90vnvg62bToQijAWwQkb3o/+C4RlW5nQDRx6UD2CQiewBsA/CGqr5ldYiw296BiIiIyCphN6NFREREZBUWLSIiIiI/YdEiIiIi8hMWLSIiIiI/YdEiIiIi8hMWLSIiIiI/YdEiIiIi8hMWLSIKCyIyXUSOiciXTGchovDBokVEYUFV96H/8iifNp2FiMIHixYRhZMGAPmmQxBR+GDRIqJw8p8AIkRkrOkgRBQeWLSIKCyIyBIAMQDeAGe1iMgiLFpEFPJEJBLATwB8GcA+ANPMJiKicMGiRUTh4PsAnlbVSrBoEZGFWLSIKKSJyBQAiwA8MnAXixYRWUZU1XQGIiIiopDEGS0iIiIiP2HRIiIiIvITFi0iIiIiP2HRIiIiIvITFi0iIiIiP2HRIiIiIvITFi0iIiIiP/n/Pigc4tYBYYoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lams = np.linspace(0, 5, 100)\n",
    "l_lam = [likelihood_poi(lam=l, data=sample) for l in lams]\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.plot(lams, l_lam);\n",
    "ax.set_xlabel(\"$\\lambda$\");\n",
    "ax.set_ylabel(\"$L(\\lambda)$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "219327af-2033-439e-a49c-377c7a8493ef",
   "metadata": {},
   "source": [
    "## 事前分布\n",
    "\n",
    "$$\n",
    "P(\\lambda) = \\mathrm{Gam}(\\lambda | a, 1) = \\frac{1}{\\Gamma(a)} \\lambda^{a-1} e^{-\\lambda}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "01e6f6a1-d18a-485f-9c3b-ddb824d35dd3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAEMCAYAAACr2eC6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjLklEQVR4nO3de5Bmd33f+fdn+jYzPfcbEpIAeaWYCEjh0AinHLOOHeHBS5CqVjbC2IiEjZIQZZOi7Fgk8XpXa6dMpXZxvMtmrZirA8gElngSQmSyiGSdWLJ6sJCQZMEgLppB0tzvM91z+e4fz3m6n+npmemWZp7uM/N+VT11zvmd3znP7+gp1Xz69zu/c1JVSJIkqV2WLHQDJEmSNH+GOEmSpBYyxEmSJLWQIU6SJKmFDHGSJEktZIiTJElqob6FuCSbkzydZFuSe2bZ//4kTyZ5LMn/m+SVPfvuTPKt5nNnT/kbkjzenPO3k6Rf1yNJkrSQ0o/nxCUZAL4J3AJsBx4B3llVT/bU+SvAw1V1NMnfAX6iqt6RZB0wDowBBWwF3lBV+5L8CfA/Ag8D/x747ar60iW/IEmSpAXWr564m4FtVfVMVU0C9wO39laoqger6miz+RBwbbP+08CXq2pvVe0DvgxsTnI1sKqqHqpOEv0kcFsfrkWSJGnBDfbpe64Bnu3Z3g686Tz13wt0e9RmO/aa5rN9lvKzJLkLuAtgdHT0Da9+9avn03ZJkqQFsXXr1t1VtXG2ff0KcXOW5BfoDJ3+txfrnFV1H3AfwNjYWI2Pj1+sU0uSJF0ySb53rn39Gk7dAVzXs31tU3aGJH8V+MfA26tq4gLH7mB6yPWc55QkSboc9SvEPQLcmOT6JMPAHcCW3gpJfgT4HToBbmfPrgeAtyRZm2Qt8Bbggap6DjiY5EebWanvBv6gHxcjSZK00PoynFpVJ5PcTSeQDQAfraonktwLjFfVFuCfASuAf908KeT7VfX2qtqb5H+lEwQB7q2qvc36+4CPA8vo3EPnzFRJknRF6MsjRhYT74mTJEltkWRrVY3Nts83NkiSJLWQIU6SJKmFDHGSJEktZIiTJElqIUOcJElSCxniJEmSWsgQJ0mS1EKGOEmSpBYyxEmSJLWQIU6SJKmFDHGSJEktZIiTJElqIUOcJElSCxniJEmSWsgQJ0mS1EKGOEmSpBYyxEmSJLWQIU6SJKmF+hbikmxO8nSSbUnumWX/m5N8LcnJJLf3lP+VJI/2fI4nua3Z9/Ek3+nZ9/p+XY8kSdJCGuzHlyQZAD4M3AJsBx5JsqWqnuyp9n3gPcAv9R5bVQ8Cr2/Osw7YBvxhT5VfrqrPXbLGS5IkLUJ9CXHAzcC2qnoGIMn9wK3AVIirqu82+06f5zy3A1+qqqOXrqmSJEmLX7+GU68Bnu3Z3t6UzdcdwGdmlP1GkseSfCjJyIttoCRJUpu0ZmJDkquB1wEP9BR/AHg18EZgHfAr5zj2riTjScZ37dp1ydsqSZJ0qfUrxO0AruvZvrYpm4+fA75QVSe6BVX1XHVMAB+jM2x7lqq6r6rGqmps48aN8/xaSZKkxadfIe4R4MYk1ycZpjMsumWe53gnM4ZSm945kgS4DfjGS2+qJEnS4teXEFdVJ4G76QyFPgV8tqqeSHJvkrcDJHljku3AzwK/k+SJ7vFJXkWnJ+8/zTj1p5I8DjwObAB+/ZJfjCRJ0iKQqlroNvTV2NhYjY+PL3QzJEmSLijJ1qoam21fayY2SJIkaZohTpIkqYUMcZIkSS1kiJMkSWohQ5wkSVILGeIkSZJayBAnSZLUQoY4SZKkFjLESZIktZAhTpIkqYUMcZIkSS1kiJMkSWohQ5wkSVILGeIkSZJayBAnSZLUQoY4SZKkFjLESZIktZAhTpIkqYUMcZIkSS3UtxCXZHOSp5NsS3LPLPvfnORrSU4muX3GvlNJHm0+W3rKr0/ycHPO308y3I9rkSRJWmh9CXFJBoAPA28FbgLemeSmGdW+D7wH+PQspzhWVa9vPm/vKf8g8KGqugHYB7z3ojdekiRpEepXT9zNwLaqeqaqJoH7gVt7K1TVd6vqMeD0XE6YJMBPAp9rij4B3HbRWixJkrSI9SvEXQM827O9vSmbq6VJxpM8lOS2pmw9sL+qTl7onEnuao4f37Vr1zybLkmStPgMLnQD5uiVVbUjyQ8BX0nyOHBgrgdX1X3AfQBjY2N1idooSZLUN/3qidsBXNezfW1TNidVtaNZPgN8FfgRYA+wJkk3iM7rnJIkSW3WrxD3CHBjM5t0GLgD2HKBYwBIsjbJSLO+Afgx4MmqKuBBoDuT9U7gDy56yyVJkhahvoS45r61u4EHgKeAz1bVE0nuTfJ2gCRvTLId+Fngd5I80Rz+54HxJF+nE9p+s6qebPb9CvD+JNvo3CP3kX5cjyRJ0kJLp0PryjE2Nlbj4+ML3QxJkqQLSrK1qsZm2+cbGyRJklrIECdJktRChjhJkqQWMsRJkiS1kCFOkiSphQxxkiRJLWSIkyRJaiFDnCRJUgsZ4iRJklrIECdJktRChjhJkqQWMsRJkiS1kCFOkiSphQxxkiRJLWSIkyRJaiFDnCRJUgsZ4iRJklrIECdJktRCfQtxSTYneTrJtiT3zLL/zUm+luRkktt7yl+f5I+TPJHksSTv6Nn38STfSfJo83l9ny5HkiRpQQ3240uSDAAfBm4BtgOPJNlSVU/2VPs+8B7gl2YcfhR4d1V9K8nLga1JHqiq/c3+X66qz13SC5AkSVpk+hLigJuBbVX1DECS+4FbgakQV1Xfbfad7j2wqr7Zs/6DJDuBjcD+S95qSZKkRapfw6nXAM/2bG9vyuYlyc3AMPDtnuLfaIZZP5Rk5BzH3ZVkPMn4rl275vu1kiRJi05rJjYkuRr4PeCvV1W3t+4DwKuBNwLrgF+Z7diquq+qxqpqbOPGjX1pryRJ0qXUrxC3A7iuZ/vapmxOkqwCvgj846p6qFteVc9VxwTwMTrDtpIkSZe9foW4R4Abk1yfZBi4A9gylwOb+l8APjlzAkPTO0eSALcB37iYjZYkSVqs+hLiquokcDfwAPAU8NmqeiLJvUneDpDkjUm2Az8L/E6SJ5rDfw54M/CeWR4l8qkkjwOPAxuAX+/H9UiSJC20VNVCt6GvxsbGanx8fKGbIUmSdEFJtlbV2Gz7WjOxQZIkSdMMcZIkSS1kiJMkSWohQ5wkSVILGeIkSZJayBAnSZLUQoY4SZKkFjLESZIktZAhTpIkqYUMcZIkSS1kiJMkSWohQ5wkSVILGeIkSZJaaN4hLslokoFL0RhJkiTNzQVDXJIlSX4+yReT7AT+DHguyZNJ/lmSGy59MyVJktRrLj1xDwL/DfAB4Kqquq6qNgF/GXgI+GCSX7iEbZQkSdIMg3Oo81er6sTMwqraC3we+HySoYveMkmSJJ3TBXviZgtwL6aOJEmSLp55TWxIcl2SzUl+KcknkozP49jNSZ5Osi3JPbPsf3OSryU5meT2GfvuTPKt5nNnT/kbkjzenPO3k2Q+1yNJktRWc5nY8LeS/Nck+4FvAv8DsALYAvz8XL6kmc36YeCtwE3AO5PcNKPa94H3AJ+ecew64NeANwE3A7+WZG2z+18AfxO4sflsnkt7JEmS2m4u98R9AHgHsBv4TWAZ8NGq+v48vudmYFtVPQOQ5H7gVuDJboWq+m6z7/SMY38a+HJzDx5JvgxsTvJVYFVVPdSUfxK4DfjSPNp10f3r8Wf5o227uWrVUq5avXR6uXopG1eMMDjgo/kkSdJLN5cQ97aq+kaz/rNJ3gr82yQfB/55Vc0MXbO5Bni2Z3s7nZ61uZjt2Guaz/ZZys+S5C7gLoBXvOIVc/zaF2fX4Qm+9v19vHBggslTZ/6nWRLYsGLkjHD3slUz1lcvZcXIXH4WSZJ0JZtLWniid6OqvpTkQeCfAP8F+EtJUlV1KRp4MVTVfcB9AGNjY5e0ne/7iRt430/cQFWx98gkzx88zvMHjvP8weO80CyfO3Cc7+45wkPP7OHg8ZNnnWPFyCAvWzUye8hr1jesGGFgibcASpJ0pZpLiHswyeeBP+gOoVbV8ST3At9P8gk6z5L7+HnOsQO4rmf72qZsLnYAPzHj2K825de+yHNecklYv2KE9StGeM3LV5+z3rHJU1NB74Um4L1wsPN5/uBxHvr2HnYemuDk6TOz58CSsHHFCC9bNcLLVk334m1aOR3+XrZqKauWDuJ8D0mSLj9zCXGbgb8BfCbJDwH76NwXtwT4Q+C3qupPL3COR4Abk1xPJ2jdwRwnRQAPAP+0ZzLDW4APVNXeJAeT/CjwMPBu4P+Y4zkXjWXDA1y/YZTrN4yes87p08XuIxO8cGBiKtxNB70JvrfnKA9/Zy8Hjp39pJdlQwPnD3orl7Jp1QhLh3yTmiRJbZL5jII2D/XdAByrqv3z+qLkZ4DfAgboTIz4jaY3b7yqtiR5I/AFYC1wHHi+ql7THPs3gH/UnOo3qupjTfkYnR7AZXQmNPy9Cw3rjo2N1fj4nJ+M0irHT5xqwt0Ezx88zs6eodydTdkLB48zcfLs2xhXLxviqlWdQNcJfJ3lppVLp4Z2N6wYYciJGZIk9U2SrVU1Nuu+C4W45rls/xudnrd/C9xdVYcueiv75HIOcXNRVRw4dmIq2L1w8Dg7D01Mh71DE+xslqdmDOEmsH60dwh3hI1NyHvZyqVTZeu9X0+SpIvifCFuLsOpvwrcQmcY9O8B/7RZqoWSsGb5MGuWD/Pqq85d79TpzsSMF6aGbruBb/r+vce2H2DPkQlm/h3QnYXb6ckbYVOz7Ia8bu+eYU+SpBdvLiHuYM89b7+a5OFL2SAtDgNLwsaVI2xcOcJrrzn3xIwTp06z5/Dk1PDtC4cm2NUzpPuDA8f5+vb97D48edax3bC3qenJ29QEvE09QW/TyqVsWDHs8/UkSZphLiHu6uY5a38GPAX4sntNGRpYMvUw4/OZPHma3Yenh293Ti0neOFQZ1bu17fvZ8+RybN69jrDuMNsXNnt0esNeyNT5ZtWjTAy6AQNSdKVYS4h7teA1wHvapYrkvx74OvAY1X1mUvYPl0mhgeX8PI1y3j5mmXnrdft2ZsKe4c69+5NLyf4s+cPsvvw5Fn37EFngkY30G1qwl23R7E3+K0Y8dErkqR2u2CIax6UOyXJtXTC3F8AfgYwxOmimWvPXu89e7sOnRnydh7qBMA/+c5edh06+80Z0Hn0SifYdQLfxhWde/emA19nuX7U+/YkSYvTvN/vVFXb6bziakHfUaorW+89e+fTnY3bHbrddbgT9jrBrxP4nn7+EP/fod0cmuXtGQNL0gzlTge7TSvPDnsbV46wfNjXpUmS+sd/dXRZ652N++detvK8dY+fODXVq9cNebua8Lfz0HF2HZ7gyefOPZQ7OjzQ6c1bMR3sNq4cOWt7/agTNSRJL50hTmosHRrgunXLuW7d8vPWO3262Ht08oygN/U53Jm08dTzB/nP35qYtXcvgXXLh88Z8jauGGFDs1yzfMh79yRJszLESfO0ZEnYsGKEDStG+PNXn7/udO/eBLsPzwx7nbJndh1h1+EJJmd5k8Zg810bV46wYcVws5wOfL3rK52sIUlXFEOcdAnNtXevqjh4/CS7Zgl73bKdh84/nDs8uKSnF+/MwNcNnd0g6OxcSWo/Q5y0CCRh9bIhVi8b4oZNK85b9/TpYt/RSXYfnpwKeFPBr1nu2H+cR589wN4jE8yS9xgZXNIJdTMCX2/Y29CUrVpq4JOkxcgQJ7XMkiVh/YrOa8t++KrzT9boPoqlG/R2H55g96FJdh2eYHcT+nbsP87Xtx9gz+HZA9/wwBI2rBhmfTfcNeFvQ+92s752+TBLfCSLJPWFIU66jM31USzQCXydHr5O0Jvq4evZ7g7p7jk8yclZEt+SwLrR6WHb9aPToW/9aNO7NzrChpXDrBsd9g0bkvQSGOIkAZ3A1+1V46rz1z19uvP8vT1HJtjVE/j2HO4Nf5N8Z/cRdh+e4PiJsydtAKxaOjj1neubXr31U7173d6/TpkTNyTpTIY4SfO2ZElYOzrM2tFhbth0/rpVxdHJUz1DupNTYW9Ps73r8ATf2nmYh57Zw76jJ2Y9z/DAEtavGO58Rkd6gl5nuzcE2ssn6UpgiJN0SSVhdGSQ0ZFBXrl+9IL1T5w6zb4jnWDXvZ+vE/qmg9+eI5Ns23n4nI9mAVjZ9PKtGx1m/ej0PX3d9W7oWzfauZfP16tJahtDnKRFZWhgCZtWLWXTqvO/Pxc6vXyHJ06y5/Ake46c2cvXDYB7j0zyvT1H2fq9few7Ojnr5I0E1i7vBrzpnr31oyOsWzHMhib4rRsdZsOKYVYtHXICh6QFZ4iT1FpJWLl0iJVLh3jVhgv38p06Xew/Osmenh6+vUcmO8O6RybZ2wTAp547yJ4jkxw4NvvQ7sCSTPXwrWs+3V69btBbNzrdC7h6maFP0sXXtxCXZDPwz4EB4Her6jdn7B8BPgm8AdgDvKOqvpvkXcAv91T9C8BfrKpHk3wVuBo41ux7S1XtvLRXIqmtBnoez3Khd+nC9NDuniOTU719U8Gv6fnbe2SSJ35wkN2HZ3/NWvd7uz1966Z6+5qg16yvXT48dT+fw7uS5qIvIS7JAPBh4BZgO/BIki1V9WRPtfcC+6rqhiR3AB+kE+Q+BXyqOc/rgH9TVY/2HPeuqhrvx3VIurLMZ2gXYPLk6amAt/fIZDOkO8neZrsbAJ/8wfl7+rrDu92evfXNJJLenr/pfSOsHR1yIod0BepXT9zNwLaqegYgyf3ArUBviLsV+J+b9c8B/2eSVFXvHSzvBO6/9M2VpPkbHlzCVauXctXquYW+3p6+vd1lcx/f3qOTU8Fv287D7D0yec57+gBWjgyyrnngcm/oW9sNfMuHWbdieukjW6T261eIuwZ4tmd7O/Cmc9WpqpNJDgDrgd09dd5BJ+z1+liSU8DngV+fEfoASHIXcBfAK17xipdwGZJ08cy3p+9U83y+Ts9eZ7nnyOQZQXDvkUmeP3i881DmI5PnnL07NJAzevvWdoNez3Z3mLezbW+ftNi0ZmJDkjcBR6vqGz3F76qqHUlW0glxv0jnvrozVNV9wH0AY2Nj5/g7VpIWt+6EinWjw3OqX1UcmTzFvm7AO9qZvNFd74a/fUcmeeoHB9l7dJL953hOH8CKkUHWjg6xbvmZPXzd9anAt3yItaPDrFk2xODAkot1+ZJm6FeI2wFc17N9bVM2W53tSQaB1XQmOHTdAXym94Cq2tEsDyX5NJ1h27NCnCRdiZKwYmSQFSODXLdu+ZyOOXnqNPuPnZgKePuPTk71+nWX+46eYO+RSb71wmH2HZ3k6OSpc55v9bKhqWB3RtBrAuCabnmzvWrZkJM6pDnqV4h7BLgxyfV0wtodwM/PqLMFuBP4Y+B24CvdodEkS4CfA368W7kJemuqaneSIeBtwH+81BciSZezwYElU69Cu3GOxxw/cYp9zT18+46caHr8OmFv39HJzvLIJD/Yf5xv7Oj0+J1rmDeBNcs6PXlrl3c/naC3Zvkw60aHmmXT47e88wgXe/x0JepLiGvucbsbeIDOI0Y+WlVPJLkXGK+qLcBHgN9Lsg3YSyfodb0ZeLY7MaIxAjzQBLgBOgHuX/bhciRJPZYODXD16mVcvXrZnOpXFcdOnJoKfXuOTLC/G/imhno72zv2H+MbOw6cN/hBp8evO4y7ttvD1wz1rmnC3pqe3sA1y73HT+2XWeYBXNbGxsZqfNwnkkhSm3SDX7dXrztbd9+RTk/f/qOT7O0uj0xOhcLzDfWODg+wZnln0kYn2HV697rLbtjr9giuGR1yVq/6LsnWqhqbbV9rJjZIkq5cSVg+PMjy4UGuWTO3Hj/oDPVO9fL19PB17/Xb3y0/eoJn9x5l39ET53x+H8DgkrCmJ+itWT48Nfw7Hfi6+4ebuvb66dIwxEmSLltLhwa4avXAnJ/dB9OPcumGvengN32PXzf8Pbv3KI81Zecb7l0+PMCaZU3o697Lt7wzBLxm2ZlDvt06zu7VhRjiJEnqMd9HuXQdmzw11eO3/+iJqdC3fyoAnuDAsU7ge+r5gxw4eoL9x05w6lxPcKbzEOc1o9NBb6rnb/kQq6d6AYdYvWw6+DnR48phiJMk6SJYNjzAsuFlvHwew71VxaGJk+zv9vYdOzEV+rohcLpX8ATb9x2bKjvfLe0rlw42oa4T7lYvG5plezr4rWlm+Q4PGv7axBAnSdICScKqpUOsWjrEK9bP7Vl+AKdPF4eOn5wl+HW3T0yFvwPHTrBj37Gpeufp+Jsa9l29fJjVywanQ19PAOz29q3uCYOjwwNO+FgAhjhJklpmyZJ0QtTyoXkdd/p0p+fvYBP09jfDuweasLe/GeLd3wz9PrP78FSv4OSpc9/z153wsWrZmT173R7A3uXqZdP77P17aQxxkiRdIZYsyVR4um7d3I+rKo6fOD3Vs9ft6TtwbHJqff+xpuzoCXYeOs43XzjEgWMnOHT85HnPvXx44IxQd2bw67zFY/Wynh7A5uPbPQxxkiTpApK8qHv+oPMqt0PHT04N5x7ohr0m8O2fsf29PUf5+vZOveMnzt37B52JH6tmhLtuCJyt/HILgIY4SZJ0yQwOLOm8SWN0GBid17ETJ09Nhbtu0Nvfs37g2AkO9qx/e9fhqfWJ8zzyBWDFyOBUoFu9bPCcYa83DK5auriGgA1xkiRpURoZHGDTygE2rZz7c/66jp841bn3r6eX74xewCYAHjzeWf/O7iNT5RfqAVw21BkC/t07x3jtNatf7OW9ZIY4SZJ02Vk6NMDSoQE2rZp/AJw4eYqDx07OGvZ6e/7WzvNZghebIU6SJKnHyOAAG1cOsHHlyEI35bwWx6CuJEmS5sUQJ0mS1EKGOEmSpBYyxEmSJLWQIU6SJKmFDHGSJEkt1LcQl2RzkqeTbEtyzyz7R5L8frP/4SSvaspfleRYkkebz//dc8wbkjzeHPPbSdr/Dg1JkqQ56EuISzIAfBh4K3AT8M4kN82o9l5gX1XdAHwI+GDPvm9X1eubz9/uKf8XwN8Ebmw+my/VNUiSJC0m/eqJuxnYVlXPVNUkcD9w64w6twKfaNY/B/zU+XrWklwNrKqqh6qqgE8Ct130lkuSJC1C/Qpx1wDP9mxvb8pmrVNVJ4EDwPpm3/VJ/jTJf0ry4z31t1/gnAAkuSvJeJLxXbt2vbQrkSRJWgTaMLHhOeAVVfUjwPuBTydZNZ8TVNV9VTVWVWMbN268JI2UJEnqp36FuB3AdT3b1zZls9ZJMgisBvZU1URV7QGoqq3At4E/19S/9gLnlCRJuiz1K8Q9AtyY5Pokw8AdwJYZdbYAdzbrtwNfqapKsrGZGEGSH6IzgeGZqnoOOJjkR5t7594N/EE/LkaSJGmhDfbjS6rqZJK7gQeAAeCjVfVEknuB8araAnwE+L0k24C9dIIewJuBe5OcAE4Df7uq9jb73gd8HFgGfKn5SJIkXfbSmdh55RgbG6vx8fGFboYkSdIFJdlaVWOz7WvDxAZJkiTNYIiTJElqIUOcJElSCxniJEmSWsgQJ0mS1EKGOEmSpBYyxEmSJLWQIU6SJKmFDHGSJEktZIiTJElqIUOcJElSCxniJEmSWsgQJ0mS1EKGOEmSpBYyxEmSJLWQIU6SJKmFDHGSJEktZIiTJElqob6FuCSbkzydZFuSe2bZP5Lk95v9Dyd5VVN+S5KtSR5vlj/Zc8xXm3M+2nw29et6JEmSFtJgP74kyQDwYeAWYDvwSJItVfVkT7X3Avuq6oYkdwAfBN4B7Ab+WlX9IMlrgQeAa3qOe1dVjffjOiRJkhaLfvXE3Qxsq6pnqmoSuB+4dUadW4FPNOufA34qSarqT6vqB035E8CyJCN9abUkSdIi1a8Qdw3wbM/2ds7sTTujTlWdBA4A62fU+e+Br1XVRE/Zx5qh1F9NkovbbEmSpMWpNRMbkryGzhDr3+opfldVvQ748ebzi+c49q4k40nGd+3adekbK0mSdIn1K8TtAK7r2b62KZu1TpJBYDWwp9m+FvgC8O6q+nb3gKra0SwPAZ+mM2x7lqq6r6rGqmps48aNF+WCJEmSFlK/QtwjwI1Jrk8yDNwBbJlRZwtwZ7N+O/CVqqoka4AvAvdU1X/pVk4ymGRDsz4EvA34xqW9DEmSpMWhLyGuucftbjozS58CPltVTyS5N8nbm2ofAdYn2Qa8H+g+huRu4Abgf5rxKJER4IEkjwGP0unJ+5f9uB5JkqSFlqpa6Db01djYWI2P+0QSSZK0+CXZWlVjs+1rzcQGSZIkTTPESZIktZAhTpIkqYUMcZIkSS1kiJMkSWohQ5wkSVILGeIkSZJayBAnSZLUQoY4SZKkFjLESZIktZAhTpIkqYUMcZIkSS1kiJMkSWohQ5wkSVILGeIkSZJayBAnSZLUQoY4SZKkFjLESZIktZAhTpIkqYX6FuKSbE7ydJJtSe6ZZf9Ikt9v9j+c5FU9+z7QlD+d5Kfnek5JkqTLVV9CXJIB4MPAW4GbgHcmuWlGtfcC+6rqBuBDwAebY28C7gBeA2wG/q8kA3M8pyRJ0mWpXz1xNwPbquqZqpoE7gdunVHnVuATzfrngJ9Kkqb8/qqaqKrvANua883lnJIkSZelwT59zzXAsz3b24E3natOVZ1McgBY35Q/NOPYa5r1C50TgCR3AXc1m4eTPP0irmE+NgC7L/F3aP78XRYff5PFyd9l8fE3WXz69Zu88lw7+hXiFlRV3Qfc16/vSzJeVWP9+j7Njb/L4uNvsjj5uyw+/iaLz2L4Tfo1nLoDuK5n+9qmbNY6SQaB1cCe8xw7l3NKkiRdlvoV4h4BbkxyfZJhOhMVtsyoswW4s1m/HfhKVVVTfkcze/V64EbgT+Z4TkmSpMtSX4ZTm3vc7gYeAAaAj1bVE0nuBcaragvwEeD3kmwD9tIJZTT1Pgs8CZwE/m5VnQKY7Zz9uJ456NvQrebF32Xx8TdZnPxdFh9/k8VnwX+TdDq7JEmS1Ca+sUGSJKmFDHGSJEktZIi7yHwV2OKT5KNJdib5xkK3RR1JrkvyYJInkzyR5O8vdJuudEmWJvmTJF9vfpP/ZaHbpI7mLUV/muTfLXRb1JHku0keT/JokvEFa4f3xF08zavAvgncQufhw48A76yqJxe0YVe4JG8GDgOfrKrXLnR7BEmuBq6uqq8lWQlsBW7z/5WF07whZ7SqDicZAv4I+PtV9dAFDtUlluT9wBiwqqrettDtUSfEAWNVtaAPYLYn7uLyVWCLUFX9ZzoznrVIVNVzVfW1Zv0Q8BTTb2LRAqiOw83mUPPxr/wFluRa4L8Dfneh26LFxxB3cc32ejH/YZLOI8mrgB8BHl7gplzxmmG7R4GdwJeryt9k4f0W8A+B0wvcDp2pgD9MsrV5teeCMMRJWjBJVgCfB/5BVR1c6PZc6arqVFW9ns4bcG5O4u0HCyjJ24CdVbV1oduis/zlqvqLwFuBv9vcttN3hriLy1eBSXPU3Hf1eeBTVfX/LHR7NK2q9gMPApsXuClXuh8D3t7cf3U/8JNJ/tXCNkkAVbWjWe4EvkDndqq+M8RdXL4KTJqD5ib6jwBPVdX/vtDtESTZmGRNs76MzgStP1vQRl3hquoDVXVtVb2Kzr8nX6mqX1jgZl3xkow2E7JIMgq8BViQpx8Y4i6iqjoJdF8F9hTw2UX0KrArVpLPAH8M/HCS7Uneu9BtEj8G/CKdnoVHm8/PLHSjrnBXAw8meYzOH6RfriofaSGd7WXAHyX5Op13uX+xqv7DQjTER4xIkiS1kD1xkiRJLWSIkyRJaiFDnCRJUgsZ4iRJklrIECdJktRChjhJkqQWMsRJkiS1kCFOkl6CJK9L8r0kf2eh2yLpymKIk6SXoKoep/NKpHcvdFskXVkMcZL00u0EXrPQjZB0ZTHESdJL95vASJJXLnRDJF05DHGS9BIkeSswCnwRe+Mk9ZEhTpJepCRLgQ8C7wMeB167sC2SdCUxxEnSi/dPgE9W1XcxxEnqM0OcJL0ISX4YuAX4rabIECepr1JVC90GSZIkzZM9cZIkSS1kiJMkSWohQ5wkSVILGeIkSZJayBAnSZLUQoY4SZKkFjLESZIktdD/D/+obhQsn4dkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = 1\n",
    "b = 0.1\n",
    "p_lam = stats.gamma(a=a, scale=1./b).pdf(lams)\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.plot(lams, p_lam);\n",
    "ax.set_xlabel(\"$\\lambda$\")\n",
    "ax.set_ylabel(\"$P(\\lambda)$\");\n",
    "ax.set_ylim((0.0, 0.2));   # y軸を固定"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea8e5f73-aa0b-46ee-aba6-5f44f7d22715",
   "metadata": {},
   "source": [
    "## 事後分布\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\lambda|X=x_1,…,x_n) & \\propto P(\\lambda)f(x_1,…,x_n|\\lambda) \\\\\n",
    "& \\propto \\lambda^{a + n \\bar{x} - 1} e^{-(n+1)\\lambda} \\\\\n",
    "&= \\mathrm{Gam}(a + n \\bar{x}, n+1)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a3d21b83-efa9-44cd-ac71-1835b8971c74",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.8\n",
      "2.7142857142857144\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAEICAYAAADxz+gAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAArzUlEQVR4nO3deXyVZ53+8c/3ZCUbAZJASMIetgJtgdLW7puytEXbqqXWVm2tOtbRccZRx9FxHJ1R6zg6rlMdtdrN/qwdodDS2o3WliXQskMIW0iALIRAFrKdc//+yKGmNEAoybnPcr1fr7xylofzXHro4TrPfT/3Y845RERERMSfgO8AIiIiIolOhUxERETEMxUyEREREc9UyEREREQ8UyETERER8SzZd4CzkZeX58aMGeM7hoiIiMhprV27tt45l9/bczFdyMaMGUNZWZnvGCIiIiKnZWZ7T/achixFREREPFMhExEREfFMhUxERETEMxUyEREREc9UyEREREQ8UyETERER8UyFTERERMSzmF6HTEQk2oRCjvLaJsr2HCYvK5V3Tx1BIGC+Y4lIlFMhExE5SztqmnhmSw1lexoo23uYprauN5+bWpjDF+ZO4sqJ+ZipmIlI71TIRETOwis76rnrgTW0d4UoLcji+hmFzB49lFmjh7C+qpH/fKacj/56DXPGDuWLcycxa/RQ35FFJAqpkImIvEOv7qzn7t+uYWxeJr/56BxGDE5/y/Nj8jKZN62Q36+p5L+fr+Dmn73G+84v4nvvP5ckDWOKSA+a1C8i8g6s3HWIu35TxqihGTx094VvK2PHpSYH+PDFY3jpC1fy6avG88Tr1dy3fHuE04pItNMRMhGRM7RmTwMf+80aioYM4qG7L2JYVtpp/0xGajJfeM9kDrd28vOXdjJ1ZA43njsyAmlFJBboCJmIyBlYu7eBj/xqNSMGp/Pwxy8kP/v0Zaynr99wDheMGcI//mE9m6qPDFBKEYk1KmQiIn1U39zOR3+9hoKcdB75+EUUZPc+THkqqckBfvqhWQzJSOUTv1vLoeb2AUgqIrFGhUxEpI++/2w5rR1BfnHHbIbnnHkZOy4/O43/+fAs6pvb+ZuH1tEZDPVjShGJRSpkIiJ9sO3gUR5dXcntF41mQkHWWb/ejOJcvn3zdFbtbuCbT27ph4QiEstUyERETsM5xzef3Ep2egqfu7a03173fecX89FLxvDAa3tZV3m4315XRGKPCpmIyGk8v62WVyrq+dy1peRmpPbra//DuyeRl5XGt5dtwznXr68tIrFDhUxE5BQ6gyG+tWwr4/Izuf2i0f3++plpyXz22lJW72ng+W21/f76IhIbVMhERE7hwZV72VXXwlfmTyElaWA+Mm+9oISxeZl85+ltBEM6SiaSiFTIREROorG1gx/8eQeXTsjj6skFA7aflKQAX3jPJMprmvnjuqoB24+IRC8VMhGRk/jBn3fQ1NbJP18/BbOBvfbkvGkjOLckl+8/W05bZ3BA9yUi0UeFTESkF9WNx3hw5V5unTOKySNyBnx/ZsaX5k7mwJE2Hnh1z4DvT0SiiwqZiEgvHl61l5BzfPqqCRHb58Xjh3HVpHx+8kIFR1o7I7ZfEfFPhUxE5AQdXSF+v2YfV08eTlHuoIju+x/nTqapvYufvlgR0f2KiF8RKWRm9iszqzWzTSd53szsv82swsw2mNnMSOQSEenN05sPUt/cwe0XjYr4vqcU5nDT+cX8+tU91B5ti/j+RcSPSB0h+w0w9xTPzwNKwz/3AD+LQCYRkV49uHIvo4ZmcHlpvpf9f+bqCXQGQzy0qtLL/kUk8iJSyJxzK4CGU2yyEPit67YSyDWzwkhkExHpqbymidW7G/jQhaMIBAb2zMqTGZOXyZUT83loVSUdXbrwuEgiiJY5ZEXAvh73q8KPvY2Z3WNmZWZWVldXF5FwIpI4Hly5l9TkAO+fXeI1x0cuGUt9czvLNh7wmkNEIiNaClmfOefud87Nds7Nzs/3M5wgIvGppb2LP66r5vrphQzN7N9rVp6pyybkMS4vk99oCQyRhBAthawa6Pl1tDj8mIhIxPzfG9U0t3fxoQG4ZuWZCgSMOy4ezRv7Glm/r9F3HBEZYNFSyBYDd4TPtrwIOOKc03F6EYkY5xy/e20vUwpzmDkq13ccAG6eVUxmapIWihVJAJFa9uIR4DVgkplVmdldZvZJM/tkeJNlwC6gAvgF8DeRyCUicty6ysNsO9jEhy8aPeCXSeqr7PQUbplVzJIN+6lravcdR0QGUHIkduKcW3Sa5x3w6UhkERHpzYMrK8lKS2bheSN9R3mLO941hgde28sjqyv522tKfccRkQESLUOWIiLeNLR0sHTDAW6eWURmWkS+p/bZ+PwsLp+Yz0Or9tIZ1BIYIvFKhUxEEt7SjQfoCIa4dU7kV+bvi4+8azQ1R9t5etNB31FEZICokIlIwlu6YT8TCrKYPCLbd5ReXTmxgNHDMrQEhkgcUyETkYRW29TGqt0NLJheGDWT+U8UCBgfvmg0a/ceZlP1Ed9xRGQAqJCJSEJ7etNBnIMFM6L7am3vn1VCalKAx9dV+Y4iIgNAhUxEEtqTGw4wcXgWE4dH53DlcYMzUrh2agGL39ivyf0icUiFTEQSVs3RNtbsaWDB9Oha6uJkbp5ZzKGWDl7cruv4isQbFTIRSVhPbTwQHq4c4TtKn1w+MZ9hman8UcOWInFHhUxEEtbSjQeYPCKbCQXRPVx5XEpSgIXnFfHc1loaWzt8xxGRfqRCJiIJ6eCRNtbsOcyC6dE9mf9EN80soiMYYskGXe5XJJ6okIlIQlq2sbvQzI/ysytPdM7IHCaPyObxtRq2FIknKmQikpCWbjzAlMIcxudn+Y5yRsyMm2YW8ca+RnbWNfuOIyL9RIVMRBLO/sZjrN17mOtj7OjYce89r4iAocn9InFEhUxEEs6bw5UxNn/suIKcdC4rzeeJddWEQs53HBHpBypkIpJwlm48wDkjcxibl+k7yjt286xi9h9pY+WuQ76jiEg/UCETkYRSdbiV1ysbo/5SSafz7qnDyU5L5vF11b6jiEg/UCETkYTyzOYaAOZPi+1Clp6SxIIZhTy16QAt7V2+44jIWVIhE5GE8ty2GkoLshgTw8OVx900s5jWjiDLNx/0HUVEzpIKmYgkjKa2TlbtauDqKQW+o/SL2aOHUJQ7iCe1SKxIzFMhE5GE8fKOerpCjmsmD/cdpV8EAsaCGYW8vKNOl1ISiXEqZCKSMJ7bWsvgQSnMHJXrO0q/uX5GIZ1Bp2FLkRinQiYiCSEYcry4vZYrJ+WTnBQ/H33TiwYzamiGhi1FYlz8fCqJiJzC+qpGDrV0cPXk+Jg/dpyZcf2MQl7deYhDze2+44jIO6RCJiIJ4fmttSQFjCsm5vuO0u9uOHckwZDjqU0athSJVSpkIpIQnttWy6zRQ8jNSPUdpd9NHpHN+PxMntyw33cUEXmHVMhEJO7tbzzG1gNHuSbOhiuP6x62HMmq3Q3UHG3zHUdE3oGIFTIzm2tm282swsy+1Mvzo8zsBTN73cw2mNn8SGUTkfj2/LZaAK6Jk/XHenPDuYU499cLp4tIbIlIITOzJOAnwDxgKrDIzKaesNk/A485584HbgV+GolsIhL/nt9Wy+hhGYzPz/IdZcBMKMhm8ohsnW0pEqMidYRsDlDhnNvlnOsAHgUWnrCNA3LCtwcDmgwhImftWEeQv1TUc/XkAszMd5wBdcO5I1m79zD7G4/5jiIiZyhShawI2NfjflX4sZ6+DtxuZlXAMuAzvb2Qmd1jZmVmVlZXVzcQWUUkjry6s572rlDcrM5/KtfP6L5g+lIdJROJOdE0qX8R8BvnXDEwH/idmb0tn3PufufcbOfc7Pz8+Dt9XUT613PbaslMTWLO2KG+owy40cMymV40mCU621Ik5kSqkFUDJT3uF4cf6+ku4DEA59xrQDqQF5F0IhKXnHM8v7WWyyfmk5ocTd8/B871MwrZUHWEvYdafEcRkTMQqU+oNUCpmY01s1S6J+0vPmGbSuAaADObQnch05ikiLxjWw4c5eDRtrhbnf9UFoSHLTW5XyS2RKSQOee6gHuB5cBWus+m3Gxm3zCzG8Ob/T3wcTNbDzwCfMQ55yKRT0Ti0/Nbu5e7uHJS4hSy4iEZnFeSy1ObVMhEYklypHbknFtG92T9no99rcftLcAlkcojIvHvpfI6phcNJj87zXeUiFowvZBvLdtK5aFWRg3L8B1HRPogMSZViEjCOdrWyev7GuPy2pWnM3faCACW6SiZSMxQIRORuPRqRT3BkOPyBCxkJUMzmFE8mKe0ar9IzFAhE5G49FJ5HVlpyZw/Ktd3FC/mTy9kfdURqg63+o4iIn2gQiYiccc5x4ryei6ZMIyUpMT8mJsXHrZ8auNBz0lEpC8S85NKROLazroWqhuPJeRw5XGjh2VyzsgczSMTiREqZCISd14q717C8PLSxC1k0D1s+Xplo65tKRIDVMhEJO6sKK9jXH4mJUMTe8mHN4ctN2nYUiTaqZCJSFxp6wyyavehhD86BjAuP4vJI7J1tqVIDFAhE5G4snp3A22dIa6YpEIG3cOWZXsPc/BIm+8oInIKKmQiEldWlNeRmhzgorHDfEeJCvOnd1/b8mlN7heJaipkIhJXVuyoY86YoQxKTfIdJSpMKMhi4vAslmkemUhUUyETkbixv/EY5TXNCXm5pFOZN62QNXsaqG3SsKVItFIhE5G48fKO8HIXKmRvsWBGIc7Bch0lE4laKmQiEjdWlNczIiedicOzfEeJKqUFWYzPz2SpzrYUiVoqZCISF7qCIV7eUcflE/MwM99xooqZsWB6Iat3N1Df3O47joj0QoVMROLC+qojHG3r0nDlScybXkjIwdMathSJSipkIhIXVpTXETC4dEKe7yhRafKIbMblZfKUlr8QiUoqZCISF1bsqGNGcS65Gam+o0QlM2P+9EJe23mIQxq2FIk6KmQiEvOOtHayfl+jhitPY970EYQcLN9c4zuKiJxAhUxEYt5fdtYTcnB5qYYrT2VqYQ5jhmVo2FIkCqmQiUjMW1FeR3ZaMueV5PqOEtWOD1u+uvMQDS0dvuOISA8qZCIS05xzvLyjnndNGEZykj7STmf+9EKCIcczm3W2pUg00aeXiMS0nXUtVDce0/yxPjpnZA6jhmZokViRKKNCJiIx7c3LJZWqkPVFz2HLwxq2FIkaKmQiEtNWlNcxNi+TkqEZvqPEjAXhYctnt+hsS5FooUImIjGrvSvIyl0NXKazK8/ItKIciocM0rClSBRRIRORmLV2z2GOdQY1XHmGjl/b8i8V9Rxp7fQdR0SIYCEzs7lmtt3MKszsSyfZ5gNmtsXMNpvZw5HKJiKxacWOepIDxkXjh/mOEnPmTy+kK+R4ZovOthSJBhEpZGaWBPwEmAdMBRaZ2dQTtikFvgxc4pw7B/hcJLKJSOxaUV7HrNFDyEpL9h0l5swoHkxR7iCWadhSJCqccSEzs8xwwToTc4AK59wu51wH8Ciw8IRtPg78xDl3GMA5V3um2UQkcdQ1tbPlwFEtd/EOmRnXzyjk5R31NLbqbEsR305byMwsYGa3mdlSM6sFtgEHwkOL95nZhD7spwjY1+N+VfixniYCE83sL2a20szmniTPPWZWZmZldXV1fdi1iMSjVyq03MXZuuHckXSFHE9v0rCliG99OUL2AjCe7uHEEc65EudcAXApsBL4jpnd3g9ZkoFS4EpgEfALM8s9cSPn3P3OudnOudn5+fogFklUL5fXMzQzlXNG5viOErPOGZnD2LxMlmzY7zuKSMLry8SLa51zbzsNxznXADwOPG5mKad5jWqgpMf94vBjPVUBq8L72m1m5XQXtDV9yCgiCSQUcqzYUc+lE/IIBMx3nJhlZtwwo5Afv1BBbVMbBdnpviOJJKzTHiHrrYy9g23WAKVmNtbMUoFbgcUnbPN/dB8dw8zy6B7C3HW6fYtI4tl2sIn65natP9YPbjh3JCEHT23UsKWIT2c0qd/MSsLLV/yDmT1gZmV9+XPOuS7gXmA5sBV4zDm32cy+YWY3hjdbDhwysy10D5N+wTl36EzyiUhiWHH8ckma0H/WSodnM2l4NkvWa9hSxKfTDlma2SeAO+leriINWApsovsI17f6uiPn3DJg2QmPfa3HbQd8PvwjInJSL++oY9LwbIbnaIitP9xwbiHfe6ac6sZjFOUO8h1HJCH15QjZl4G/A2YBTwLpwK+cc48758oHMpyIyIla2rtYs/swV0zS0bH+cv2MkQAs1eR+EW/6Usiud86tcs7tdM69n+4FXpeY2d+ZmS69JCIR9erOQ3QEQ1yp4cp+MyYvkxnFg3lygxaJFfGlL4Vqc887zrmngAuBYcBfAMxMpzmJSES8uL2WzNQkZo8Z6jtKXLlhxkg2VB1hT32L7ygiCalP65CZ2WfMbNTxB5xzbcA3gF+b2QN0zzETERlQzjle3F7HuybkkZqsA/T9acGMQgCe1LCliBd9+USbCwSBR8zs+Ar9u4EdwAXAD5xzvxnAjCIiAOysa6G68RhXaLiy343MHcQFY4awZL2GLUV86Ms6ZG3OuZ865y4BRgHXAOc750Y75z7unHt9wFOKiNA9XAlwpSb0D4gbzh3J9pomth9s8h1FJOH05VqWd5pZvZk1AL8Emp1zjQOeTETkBC+V1zGhIIviIRm+o8SledMKCZiGLUV86MuQ5VeB64DJQCXw7wOaSESkF60dXaza1aCzKwdQfnYaF48fxuL1++leGlJEIqUvheyoc+5151ytc+6rwJyBDiUicqLXji93ManAd5S4tvC8IvYeamVdZaPvKCIJpS+FrNDM7jGzy80sHzjdhcRFRPrdS+V1DEpJ4oKxQ3xHiWvzpo0gPSXAH9dV+Y4iklD6Usj+BZgO/BuwHZhmZsvM7D/MbNGAphMRocdyF+OHkZac5DtOXMtOT+E954zgyQ0HaO8K+o4jkjD6cpbl/c65zzjnrnDODQXGAT8CGoH5A5xPRITd9S1UNrTq7MoIuWlmMUeOdfLCtlrfUUQSxmkvLn4i51wVUAU81f9xRETe7sXtdQCaPxYhl4wfRn52Go+vq2butELfcUQSgpa6FpGo92J5HePyMykZquUuIiE5KcB7zxvJi9traWjp8B1HJCGokIlIVGvrDLJq1yGtzh9hN80spjPotCaZSISokIlIVHtt1yHau7TcRaRNKcxh8ohsHl9X7TuKSEJQIRORqPbS9jrSUwJcOHao7ygJ5+aZxazf18jOumbfUUTingqZiEQt5xwvbK/l4nHDSE/RcheRtvC8kQQMntBRMpEBp0ImIlGroraZvYdauWbKcN9RElJBTjqXlebzxOvVhEK6lJLIQFIhE5Go9ezWGgCuVSHz5qaZRVQ3HmPV7gbfUUTimgqZiEStZ7fUMKN4MCMGp/uOkrDePXUEmalJPPG6LqUkMpBUyEQkKtU2tfHGvkYdHfNsUGoS86cXsmzjQVo7unzHEYlbKmQiEpWe31qLc3DdVBUy3z5wQQnN7V08uf6A7ygicUuFTESi0rNbaijKHcTkEdm+oyS82aOHUFqQxUOrK31HEYlbKmQiEnVaO7p4paKe66YOx8x8x0l4ZsaiOaNYv6+RzfuP+I4jEpdUyEQk6ry8o572rpCGK6PITTOLSEsO8IiOkokMiIgVMjOba2bbzazCzL50iu1uNjNnZrMjlU1EosuzW2rISU9mjlbnjxq5GaksmF7I/72+X5P7RQZARAqZmSUBPwHmAVOBRWY2tZftsoHPAqsikUtEok8w5Hh+Wy1XTS4gJUkH8aPJbReOorm9iyXrdcFxkf4WqU+7OUCFc26Xc64DeBRY2Mt2/wZ8B2iLUC4RiTLrKg/T0NKh5S6i0Kzw5P6HV+/zHUUk7kSqkBUBPf8Lrgo/9iYzmwmUOOeWnuqFzOweMyszs7K6urr+TyoiXv15Sw0pScYVk/J9R5ETmBm3XajJ/SIDISrGA8wsAHwf+PvTbeucu985N9s5Nzs/Xx/YIvHm2S01XDRuGDnpKb6jSC9uOr9Yk/tFBkCkClk1UNLjfnH4seOygWnAi2a2B7gIWKyJ/SKJpaK2mV31LTq7MooNzkhhwQxN7hfpb5EqZGuAUjMba2apwK3A4uNPOueOOOfynHNjnHNjgJXAjc65sgjlE5Eo8OfwxcSv0fyxqHbbHE3uF+lvESlkzrku4F5gObAVeMw5t9nMvmFmN0Yig4hEv2e31HDOyByKcgf5jiKnMGv0ECYO1+R+kf4UsTlkzrllzrmJzrnxzrlvhR/7mnNucS/bXqmjYyKJpeZoG+sqD2u4Mgb0XLl/Q1Wj7zgicSEqJvWLiCzdcADn4PoZhb6jSB/cMquY7LRkfvnybt9RROKCCpmIRIUlG/YzpTCHCQW6mHgsyE5P4dY5JSzdeIDqxmO+44jEPBUyEfFuX0Mrr1c2cuO5I31HkTPwkUvGAvDAq3v8BhGJAypkIuLdkg3dZ+tpuDK2FOUOYv70Qh5ZVUlTW6fvOCIxTYVMRLxbsv4A54/KpWRohu8ocobuvnQsTe1dPFZW5TuKSExTIRMRrypqm9h64Cg3zNBwZSw6tySXOWOG8qtXdtMVDPmOIxKzVMhExKsl6w9gpuHKWHbXZWOpbjzG8s01vqOIxCwVMhHxxjnHkvX7uWjsMApy0n3HkXfo2inDGTMsg1+8vAvnnO84IjFJhUxEvNm8/yi76lu4QWdXxrSkgPGxS8fyxr5G1lUe9h1HJCapkImIN0s27Cc5YMybNsJ3FDlLt8wqZvCgFH6xQgvFirwTKmQi4oVzjifXH+Cy0jyGZKb6jiNnKSM1mQ9dOIrlWw6yp77FdxyRmKNCJiJerKtspLrxmIYr48hHLhlDalKAHz1f4TuKSMxRIRMRL5as309qckAXE48jBdnp3HHxaJ54vYqddc2+44jEFBUyEYm4YMjx5IYDXD2pgOz0FN9xpB994orxpKck8cM/7/AdRSSmqJCJSMS9VF5LfXM7C8/TcGW8yctK4853jWHJhv1sP9jkO45IzFAhE5GIe3hVJXlZaVyr4cq4dM9l48hMTeaHz5X7jiISM1TIRCSi9jce4/lttXxgdjEpSfoIikdDMlP52KVjWbbxIJv3H/EdRyQm6NNQRCLq92v24YBFc0b5jiID6K5Lx5KTnsx/Pau5ZCJ9oUImIhHTFQzx+zX7uKw0n5KhGb7jyAAaPCiFj182jj9vrWH9vkbfcUSingqZiETMC9vrOHi0jdt0dCwhfPTSsQzJSOG//qy5ZCKno0ImIhHz8Kq9FGSncc2UAt9RJAKy0pL5xBXjeXF7HWV7GnzHEYlqKmQiEhFVh1t5sbyOD15Qosn8CeSOi0czIiedry/ZTDDkfMcRiVr6VBSRiPj9mn0AfPCCEs9JJJIyUpP5yoIpbKo+yqNrKn3HEYlaKmQiMuA6w5P5r5yYT/EQTeZPNNfPKOSicUO5b/l2Drd0+I4jEpVUyERkwD23tZbapnZuu3C07yjigZnxrzdOo6mti+89s913HJGopEImIgPu4dWVjMhJ56pJ+b6jiCeTRmRz58VjeHh1JZuqtVisyIkiVsjMbK6ZbTezCjP7Ui/Pf97MtpjZBjN7zsz0VVokDuw91MLLO+q4dU4JyZrMn9A+d10pwzJT+dqfNhHSBH+Rt4jIp6OZJQE/AeYBU4FFZjb1hM1eB2Y752YAfwC+G4lsIjKwfv7STlKSAlqZX8hJT+FL86awrrKRP75e7TuOSFSJ1NfVOUCFc26Xc64DeBRY2HMD59wLzrnW8N2VQHGEsonIAKluPMYf1lbxwdklDM9J9x1HosBN5xcxc1Qu335qK0eOdfqOIxI1IlXIioB9Pe5XhR87mbuApwY0kYgMuJ+/uBOAT1453nMSiRaBgPGNhdM41NLBt5/a6juOSNSIugkdZnY7MBu47yTP32NmZWZWVldXF9lwItJnB4+08fs1+7hlVjFFuYN8x5EoMq1oMPdcPo5HVu/jz1tqfMcRiQqRKmTVQM/VIIvDj72FmV0LfAW40TnX3tsLOefud87Nds7Nzs/XGVsi0ep/Vuwk6ByfumKC7ygShT5/3USmFubwxcc3UNfU68e9SEKJVCFbA5Sa2VgzSwVuBRb33MDMzgf+h+4yVhuhXCIyAOqa2nl4VSXvO7+IUcO0EKy8XVpyEj+89Tya27v44uMbcE5nXUpii0ghc851AfcCy4GtwGPOuc1m9g0zuzG82X1AFvD/zOwNM1t8kpcTkSj3y5d30RkM8emrdHRMTq50eDZfnjeZ57fV8tAqXVZJEltypHbknFsGLDvhsa/1uH1tpLKIyMBpaOngdyv3csO5Ixmbl+k7jkS5Oy4ew/Pb6/jm0i1cPH4Y4/OzfEcS8SLqJvWLSGz731d2cawzyL06OiZ9EAgY990yg0EpSXzu0TfoDIZ8RxLxQoVMRPpNY2sHD7y6l/nTCikdnu07jsSI4Tnp/MdN09lYfYT7lutal5KYVMhEpN/88LkdtHR0ce/VOjomZ2butEI+fNFo7l+xi8fK9p3+D4jEGRUyEekXm6qP8MCre/jQhaOYUpjjO47EoK/dMJXLSvP4yhMbWbnrkO84IhGlQiYiZy0YcvzTExsZmpnGF94z2XcciVEpSQF+fNtMRg/L5JMPrmV3fYvvSCIRo0ImImftwZV72VB1hK9eP4XBg1J8x5EYNnhQCr+68wICZnzsN2tobO3wHUkkIlTIROSs1Bxt477l27msNI8bzx3pO47EgVHDMrj/w7OoPnyMTz64lo4unXkp8U+FTETOyjee3EJHMMS/LZyGmfmOI3Fi9pihfPeWGazc1cAXH99AMKSV/CW+RWxhWBGJPy9ur2XphgN8/rqJjNEisNLP3nt+EVWHW/neM+UEQ47//MC5pCTpOILEJxUyEXlH2jqDfPVPmxiXn8knrhjnO47EqXuvLiUpEOA7T2+jvSvIjxbNJDVZpUzij/5Wi8g78q2lW9nXcIxvvncaaclJvuNIHPvUleP5lxumsnxzDZ98cC1tnUHfkUT6nQqZiJyxR1ZX8ruVe7n70rG8a3ye7ziSAD56yVj+/X3TeWF7LXc/UEZrR5fvSCL9SoVMRM7Imj0NfO1Pm7isNI8vzdOaYxI5t104iu/dci6v7qzn9l+uovZom+9IIv1GhUxE+qy68Rif/N1aiodk8ONFM0nWBGuJsJtnFfOT22ay9UATC370Cqt3N/iOJNIv9GkqIn1yrCPIPb8to70rxC/umMXgDC0AK37Mm17In+69hOy0ZBb9YiW/fHkXzmlZDIltKmQiclrOOb7wh/VsOXCU/150HhMKsn1HkgQ3cXg2f7r3Eq6bMpxvLt3Kpx9eR3O75pVJ7FIhE5FTcs7x3eXbeXLDAb7wnklcPXm470giAGSnp/Cz22fyT/Mns3xzDTf86BVe26mLkktsUiETkZPqDIb44uMb+NmLO1k0p4RPXTHedySRtzAz7rl8PA/dfSFdoRCLfrGSzz/2BvXN7b6jiZwRFTIR6VVrRxf3/LaMx8qq+NurJ/Dv75uuSyNJ1Lpo3DCe+dwV3HvVBJas3881//kSj6yuJKRLLkmMUCETkbepb25n0f0ream8jm+9bxqff/cklTGJeoNSk/iH90ziqc9exuQR2Xz5jxu5+eev8sqOek36l6hnsfyXdPbs2a6srMx3DJG4sqe+hTt/vZqao238aNFMrpuqOWMSe5xzPL6umvuWb6PmaDvnleTymasncPXkAn25EG/MbK1zbnavz6mQiQh0zxf79V9281/P7iA9JcAv77yAWaOH+I4lclbau4L8YW0VP3txJ1WHjzGlMId7r5rAu88ZrguVS8SpkInIKa2rPMw//XEj2w42cc3kAv514TkUD8nwHUuk33QGQ/zpjf389IUKdtW3kJeVysLzirhlVjFTCnN8x5MEoUImIr06cqyT7z69jYdXVzI8O52v3ziV95wzQkM6EreCIccL22r5w9oqnttWQ2fQcc7IHG6eWcx7po2gKHeQ74gSx1TIROQtKmqbeXDlXh5fV0VLexd3vmsMf//uSWSlJfuOJhIxDS0dLH6jmj+sq2JT9VEASguyuHJSPldOKuCCMUNJTdawpvQfFTIRoTMY4tktNTy4ci+v7jxESpIxf3ohH79sHNOKBvuOJ+JVRW0TL2yr48XyWlbvbqAz6MhMTeL8UUM4ryS3+2dULnlZab6jSgxTIRNJUPsaWnmlop6/VNTz6s5DNLR0UJQ7iNsuHMUHLyjRPy4ivWhp7+LVnYdYUV7H2r2H2V7TRDC8nlnJ0EFMLcxhQkFW909+NuMLMslI1dFlOb1TFbKI/Q0ys7nAD4Ek4JfOuW+f8Hwa8FtgFnAI+KBzbk+k8onEslDIUd14jIq6ZnbWNrP9YBMrdx9iX8MxAAqy07hiYj4Lphdy1eQCkgKaIyZyMplpyVw3dfibS760dnSxqfoob+w7zBv7Gtl+sInnttbS1WPR2eE5aYzMHUTR8Z8hgxiek05eVirDMtMYlpVKVlqy5mfKSUWkkJlZEvAT4DqgClhjZoudc1t6bHYXcNg5N8HMbgW+A3wwEvlEokUw5OgMhmjvCtHeFaSlPUhLe1f3T0cXTW1d1DW1U9/cQX1zO3VN7dQcbWPPoRbaOkNvvs7QzFRmjR7C3ZeO45IJwxifn6V/CETeoYzUZOaMHcqcsUPffKyjK0RlQwsVtc1U1Daz91Ar+48cY1P1EZ7ZXENHMPS210lNDjA0I5XBg1LIGZRMTnoKOYNSyE5PJiM1mYzUpPBP9+205ABpKQHSksO3k5NISTZSkgKkJgVITuq+nRwwkgJGciAQ/m0E9KUr5kTqCNkcoMI5twvAzB4FFgI9C9lC4Ovh238Afmxm5jyOqf6lop6vL97sa/cSYX39i9bbX0l3wg0X3q77N4Scw7nux0IOukKOkHMEQ45QyNEVcnQEQ28Oi5xOSpKRl5VGXlb3t/JLJuS9OYQyPj+LoZmpffxfIyLvRGpygAkF2UwoyH7bc6GQo76lnZoj7RxqaedQc0eP3x00tXVy9FgXB4+2UV7bRFNbF60dQTq63l7izkZSwAgYBMzCtw0L3z/+uIUfMwj/7nn/r6XOrPsH/rpN9216bGNve+yvT/bpoYh9cextL3dfNpYPXjAqIvvvTaQKWRGwr8f9KuDCk23jnOsysyPAMKC+50Zmdg9wD8CoUQP7f1xmWjKlw7MGdB8SXaz3j5LeNjzpQz0/lI5/+B3/cOu+D0mBAEkBSLLub7LJASM1OUBqUvc34NSkAGnJATLTkslMSyarx++8rO5v2DriJRKdAgGjIDudguz0M/pzXcEQrZ1BWtuDtHZ0hY+Uh2jvDL55uzN4/Me9ebsr2P3lLhj+ktcZDBEKdX/5C7ruL33B8P1Q+AtlyB3/Ughw/AsjuOO36b4ffvZtXzaP3z7O9Xj+RKf8EnvaB/ufO8mOhmT4/SIbc7MQnXP3A/dD96T+gdzXeSW5/PRDswZyFyIiIgAkJwXISQqQk57iO4p4EKkFVqqBkh73i8OP9bqNmSUDg+me3C8iIiIS1yJVyNYApWY21sxSgVuBxSdssxi4M3z7FuB5n/PHRERERCIlIkOW4Tlh9wLL6V724lfOuc1m9g2gzDm3GPhf4HdmVgE00F3aREREROJexOaQOeeWActOeOxrPW63Ae+PVB4RERGRaKGLdImIiIh4pkImIiIi4pkKmYiIiIhnKmQiIiIinlksryxhZnXA3gHeTR4nXC1AooLel+ij9yQ66X2JPnpPolMk3pfRzrn83p6I6UIWCWZW5pyb7TuHvJXel+ij9yQ66X2JPnpPopPv90VDliIiIiKeqZCJiIiIeKZCdnr3+w4gvdL7En30nkQnvS/RR+9JdPL6vmgOmYiIiIhnOkImIiIi4pkKmYiIiIhnKmSnYGZzzWy7mVWY2Zd85xEws1+ZWa2ZbfKdRbqZWYmZvWBmW8xss5l91nemRGdm6Wa22szWh9+Tf/WdSf7KzJLM7HUze9J3FgEz22NmG83sDTMr85ZDc8h6Z2ZJQDlwHVAFrAEWOee2eA2W4MzscqAZ+K1zbprvPAJmVggUOufWmVk2sBZ4r/5b8cfMDMh0zjWbWQrwCvBZ59xKz9EEMLPPA7OBHOfc9b7zJDoz2wPMds55XaxXR8hObg5Q4Zzb5ZzrAB4FFnrOlPCccyuABt855K+ccwecc+vCt5uArUCR31SJzXVrDt9NCf/o23cUMLNiYAHwS99ZJLqokJ1cEbCvx/0q9I+MyCmZ2RjgfGCV5ygJLzws9gZQCzzrnNN7Eh1+APwjEPKcQ/7KAc+Y2Vozu8dXCBUyEekXZpYFPA58zjl31HeeROecCzrnzgOKgTlmpiF+z8zseqDWObfWdxZ5i0udczOBecCnw1NjIk6F7OSqgZIe94vDj4nICcLzlB4HHnLO/dF3Hvkr51wj8AIw13MUgUuAG8Nzlh4FrjazB/1GEudcdfh3LfAE3VOWIk6F7OTWAKVmNtbMUoFbgcWeM4lEnfAE8v8Ftjrnvu87j4CZ5ZtZbvj2ILpPTtrmNZTgnPuyc67YOTeG7n9TnnfO3e45VkIzs8zwyUiYWSbwbsDLWfwqZCfhnOsC7gWW0z1J+THn3Ga/qcTMHgFeAyaZWZWZ3eU7k3AJ8GG6v+2/Ef6Z7ztUgisEXjCzDXR/uXzWOaclFkTebjjwipmtB1YDS51zT/sIomUvRERERDzTETIRERERz1TIRERERDxTIRMRERHxTIVMRERExDMVMhERERHPVMhEREREPFMhExEREfFMhUxEJMzMppvZXjP7lO8sIpJYVMhERMKccxvpvqTNHb6ziEhiUSETEXmrWuAc3yFEJLGokImIvNW3gTQzG+07iIgkDhUyEZEwM5sHZAJL0VEyEYkgFTIREcDM0oHvAH8DbASm+U0kIolEhUxEpNs/A791zu1BhUxEIkyFTEQSnplNAq4DfhB+SIVMRCLKnHO+M4iIiIgkNB0hExEREfFMhUxERETEMxUyEREREc9UyEREREQ8UyETERER8UyFTERERMQzFTIRERERz/4/yNTbBJ3tDrQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_bar = sample.mean()\n",
    "a_post = a + n * x_bar\n",
    "b_post = n + 1\n",
    "p_lam_post = stats.gamma(a=a_post, scale=1./b_post).pdf(lams)\n",
    "\n",
    "print(x_bar)\n",
    "print(a_post / b_post)\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.subplots(1, 1)\n",
    "ax.plot(lams, p_lam_post);\n",
    "ax.set_xlabel(\"$\\lambda$\")\n",
    "ax.set_ylabel(\"$P(\\lambda)$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae742279-4847-4a55-a8e0-612eea07f02f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}