HYChou0515/poisson_mixture.ipynb

## poisson_mixture.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d04a5a8f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [00:08<00:00, 11.75it/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD4CAYAAAC5S3KDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABB4UlEQVR4nO3deXxU5fX48c+Z7AkJkLCTsCiI7MjurmziCgoqirt1qdXa+rVVW6u2Vau/Wq1Wq+KKW4WiVqogIqCo7JssgmUnYQuEbBBCMpnz+2NucIjZgEzuTOa8X6/7yszd5lyMc/I897nnEVXFGGOMCTUetwMwxhhjKmMJyhhjTEiyBGWMMSYkWYIyxhgTkixBGWOMCUnRbgdwPDwejyYkJLgdhjHGhJWioiJV1ZBvoIR1gkpISODAgQNuh2GMMWFFRA66HUNthHwGNcYYU/9EZKSI/CAiG0Tk/kq2x4nIJGf7QhHpUGF7OxHZLyL31vacFVmCMsYYcwQRiQJeAM4HugFXiUi3CrvdDOSqaifgGeDJCtufBqYf5TmPYAnKGGNMRQOBDaq6SVVLgPeBURX2GQVMdF5PAYaKiACIyGhgM7DmKM95hLC+B2WMMXWhtLSUrKwsiouL3Q6lTsXHx5Oenk5MTEzFTdEisiTg/QRVnRDwvi2QGfA+CxhU4RyH91FVr4jkA2kiUgzcBwwH7q1s/2rOeWSQ1W00xphIkJWVRXJyMh06dMBpBIQ9VSUnJ4esrCw6duxYcbNXVfsH6aMfAZ5R1f3H+29pCcoYE/GKi4sbVHICEBHS0tLYs2fPsRy+HcgIeJ/urKtsnywRiQYaAzn4W0VjReT/AU0An9OqWlqLcx7BEpQxxkCDSk7ljuOaFgOdRaQj/iQyDri6wj5TgeuB+cBYYLb6p8c4M+DzHwH2q+rzThKr6ZxHsAQVwQ4Ul/LFumy25RTRqnE8w7u1pElirNthGWNc5txTuhOYAUQBr6vqGhH5E7BEVacCrwFvi8gGYB/+hHPU56zuGAnn+aCSkpLUHtQ9SiUHYP4L7F/8Lgn7t7JPk5lZ1o8XvKPJjW3F/eefzLWD2zfIvyaNqcratWvp2rWr22HUyptvvsmIESNo06ZNrfav7NpEpEhVk4IRX12yYeaRZPcaeOkMmPMYy/KTeX7fqczaWMplzGa651cMyHqXhz5eQ9qwWxERRIS2Ge3cjtoYE+DNN99kx44dbodRL6yLL1LsWAFvjaJIY7nx0B/YkNSbTt2TEBG+9e7i7tzHeb3dJ9wf3ZV/DxhNz/OupnWTeCbddprbkRsTEd555x2ee+45SkpKGDRoEP/85z+5+eabWbJkCSLCTTfdREZGBkuWLGH8+PEkJCQwf/58GnI9UktQkaBwF/xrHCXRSVyUdx9rNmdyzrlJh7vx9kS34tG0J/h9zu94rPQZtiU8zsJ9bUiOt18PE3n++N81fL+joE7P2a1NCg9f3L3K7WvXrmXSpEl8++23xMTEcMcdd/Doo4+yfft2Vq9eDUBeXh5NmjTh+eef56mnnqJ//2CNEg8d1sXX0PnKYPJ1aHEBv+Q+ihLT2Tv1rz+5x1TsSeSp1Ic54Eniac/fSY7ysnHPAYiyJGVMsM2aNYulS5cyYMAA+vTpw6xZs9i3bx+bNm3irrvu4rPPPiMlJcXtMOudffs0dPOfh8yFfNH1z3y2vBmvXteD4b+v/K/D/KimvND0NzyY8zvub/QJv88fTXLfi+o5YGPcVV1LJ1hUleuvv56//OUvR6x/7LHHmDFjBi+99BKTJ0/m9ddfr/fY3GQtqIZs73qY/RiHTjyfu1Z34rzuLRnWrWW1h6yO68vXCUO48tCH9IrPpvGpV5JfVFpPARsTmYYOHcqUKVPIzs4GYN++fWzduhWfz8eYMWN49NFHWbZsGQDJyckUFha6GW69CWqCEpEmIjJFRNaJyFoROVVEUkVkpoisd342dfYVEXnOKcO+UkT6BjO2Bk8Vpv0GYuJ5Jv52vGXw+wuqLRx82Dspt1AisdwfOwlPfBIvzd0Y5GCNiWzdunXj0UcfZcSIEfTq1Yvhw4ezZcsWzjnnHPr06cM111xzuHV1ww03cPvtt9OnTx8OHgyLaZ2OWVCfgxKRicDXqvqqiMQCicDvgH2q+oQzH0hTVb1PRC4A7gIuwF8q41lVrbaQoD0HVY31M+HdseSe+UcGzDqJcQMzeHR0T8D/dPmVL8+r9vAxhe9weeE7DF03huwuV/LtA0NIif9JwUljGoRweg7qaNlzUJUQkcbAWfifNkZVS1Q1jyNLtE8ERjuvRwFvqd8CoImItA5WfA1aWSnM+D2knsjje04nOkq4a0jnozrFtKRLKZRk7k34L4WHvPxr4bYgBWuMMZULZhdfR2AP8IaILBeRV0UkCWipqjudfXYB5TdFKivF3rbiSUXkVhFZIiJLvF5vEMMPYyveg70/sO/0B/lw5R7GD2pPy5T4ozrFQU8Snza6jPPblzCufSGvfbOZQ96yIAVsjDE/FcwEFQ30BV5U1VOAA8ARU/w6hQWPqo9RVSeoan9V7R8dbYMQf6KsFL7+G7TpyyvZXVFVbjitwzGd6oukCzlQovwqaSbZhYf4bPWuuo3VGGOqEcwElQVkqepC5/0U/Alrd3nXnfMz29lem/LupiYrJ0PeVg6dfi/vLcpkRLdWZKQmHtOp9ntSeGNFKS23TuWUpsW8Z918xph6FLQEpaq7gEwR6eKsGgp8z48l2nF+fuy8ngpc54zmGwzkB3QFmtrwlflbT6168kFhD/IPlnLTGT+ZqOyoPLPgEFJWyu9azGfh5n1s3LO/joI1xpjqBfs5qLuAd0VkJdAHeBx4AhguIuuBYc57gGnAJmAD8ApwR5Bja3hWfwj7NqJn/YbX522he5sUBnRoelyn3JSr0GkofXM+Jc7js8ESxph6E9QEpaornPtFvVR1tKrmqmqOqg5V1c6qOkxV9zn7qqr+QlVPVNWeqrokmLE1OD4fzP0rtOjG11GD2ZC9n5tO71g302b0u4Go/Tu4u/1WpizLssESxoSgDh06sHfvXrfDqFNWSaKhWPsx7P0BzrqXN+ZtpVmjOC7qXUej9E8aCUktuNwzm7yiUr764ZimkDbG1JKq4vP53A7DdZagGgKfD+Y+BWmd2dh8GHN+2MM1g9sRFx1VN+ePioFTrqHZjjl0SSzk4xWRMReNMfVpy5YtdOnSheuuu44ePXrw5z//mQEDBtCrVy8efvjhw/uNHj2afv360b17dyZMmOBixMFn47Qbgh+mwe7VcOnLTFyQSWyUh/GD2tftZ/S9FvnmaX7bagl3rG1MQXGpVZYwDdP0+2HXqro9Z6uecP4TNe62fv16Jk6cSEFBAVOmTGHRokWoKpdccglz587lrLPO4vXXXyc1NZWDBw8yYMAAxowZQ1paWt3GGyKsBRXuVGHu/4OmHck/cRT/XpLFxb3b0Dw5rm4/J/UEaH8GpxfN4pC3zJ6JMiYI2rdvz+DBg/n888/5/PPPOeWUU+jbty/r1q1j/fr1ADz33HP07t2bwYMHk5mZeXh9Q2QtqHC3fibs/A4ueZ5Jy3ZwsLSMG0/vEJzP6nU58f+9m+FNdvLximZc0T+j5mOMCTe1aOkES1KSvzyeqvLAAw9w2223HbH9yy+/5IsvvmD+/PkkJiZyzjnnUFxcHJRYRGQk8CwQBbyqqk9U2B4HvAX0A3KAK1V1i4gMBMr7HgV4RFU/co7ZAhQCZYBXVaudddFaUOFMFb56Ehq3w9vjCibO28rAjqn0aNs4OJ/XbRRExXJ706XM25hDdmFw/scwJtKdd955vP766+zf73/ucPv27WRnZ5Ofn0/Tpk1JTExk3bp1LFiwICifLyJRwAvA+UA34CoRqTgdws1Arqp2Ap4BnnTWrwb6q2ofYCTwsogENobOVdU+NSUnsAQV3jbOhu1L4Mx7mPnDPrbnHeSm04/vwdxqJTSFziPonT8LUR8z1uwO3mcZE8FGjBjB1VdfzamnnkrPnj0ZO3YshYWFjBw5Eq/XS9euXbn//vsZPHhwsEIYCGxQ1U2qWgK8j7+gd6DAwt9TgKEiIqpapKrlhVLjOcpydoGsiy9clbeeUtKhz9W8/upS0psmMLyGCQmPW68riF73CWObbmT6quZcO7iOB2MYE6E6dOjA6tWrD7+/++67ufvuu3+y3/Tp0ys9fsuWLUfzcdEiEvis6QRVDRwSWFnx7orTHx3eR1W9IpIPpAF7RWQQ8DrQHrg2IGEp8LmIKPByhc/8aZBHc0UmhGyeC5kL4YKnWLnrIIu35PLghV2J8tTBg7nV6XwexDXm+kYLuXhTZ3L2HyKtUR0PyDDGBFuN93+Oh1ODtbuIdAUmish0VS0GzlDV7SLSApgpIutUdW5V57EuvnBU3npKbg2nXMsrX2+mUVw0Vw6oh0ELMfHQ7RJOzvuKGC3h8++tm8+YBqg2xbsP7+PcY2qMf7DEYaq6FtgP9HDeb3d+ZgMf4e9KrJIlqHC04QvY+i2ccQ/bDyjTVu1k3IAMkuvruaQeY4gqPcAVjdcybZXV8zUNQzBnF3fLcVzTYqCziHR0ZkMfh7+gd6DAwt9jgdmqqs4x0QAi0h44GdgiIkkikuysTwJG4B9QUSVLUOHGVwYzH4amHaHfDbz57WYAbjzOquVHpcOZkJjGVY2WMX9jDnlFJfX32cYEQXx8PDk5OQ0qSakqOTk5xMcf3WSlzrFe4E5gBrAWmKyqa0TkTyJyibPba0CaiGwA7uHH+f7OAL4TkRX4W0l3qOpe/JPTfiMi3wGLgE9V9bPq4rB7UOFm5STIXgNj36DQK7y/KJMLerambZOE+oshKhq6XkyX7yYT5RvPzO93c7k9E2XCWHp6OllZWezZ07DqTMbHx5Oenn5Mx6rqNPyzTASueyjgdTFweSXHvQ28Xcn6TUDvo4nBElQ4KTkAsx+DNn2h+6W889UmCg95ueXMemw9les2mqilbzIm5XumrWprCcqEtZiYGDp2dOH/I1Mt6+ILJ18+AQVZcN7j7C8pY8LcjZx9UnN6pTep/1icbr7xjZbz7YYcCopL6z8GY0yDZgkqXOxaDfNfgL7XQftTmThvC7lFpfx6+EnuxON0851c8C1SVsycddnuxGGMabAsQYWDslL47y/9lRyG/ZHC4lJe+XoTQ05uQZ+MJu7F1W00Ud4iLkn6nhlrrHisMaZuWYIKB3Meh+1L4cKnIDGVf8zeQP7BUu5xq/VUzunmuzZ5KXPW7aG41GbaNcbUHUtQoW7DLPjmGX/XXvdL2ZC9n9e/2cyV/TOCVxS2tpxuvu775+MrPcjc/zWsEVDGGHdZggplu9fAv2+AFl1h5BP4fMrDU1eTEBvFved1cTs6P6eb74KE1TZHlDGmTlmCClU5G+HdyyE2Ccb/G2KTmDh/C99uyOH+80+mWajUv3O6+a5LXs4Xa3dT4vW5HZExpoGwBBWKdiyH10eCt9ifnBqn811mHn+Zvo5hXVtw9cB2bkf4I6ebr2fRfA4VF7FgU07NxxhjTC1YggolPh8snACvjYCoWLhxOrTqyZa9B7h54hJaJMfx5JheiAS5YvnR6jaaaG8R58Wu5DMbzWeMqSOWoOpJUYmXeRv3MvW7Hcxet5utOQd+rPvlK4P1X8CrQ2D6b+CEc+D2r6F5F1Zk5jH2pfl4fT7euGFAaE5t0eFMSGzGdY2X8/ma3ZT5Gk49M2OMe4Ja6qiy+edFJBWYBHQAtgBXqGqu+JsFzwIXAEXADaq6LJjx1YeC4lKe/vx/TFqcycEKw7B7J+7j1sYLOKdoJkmHdqPJbZDLXkF7jGHLvmLemrWGt+ZvpVVKPO/fNIhOLZJduooaREVDt0vos/xfHDhwLcu25TKgQ6rbURljwlx91OI716lkW+5+YJaqPiEi9zvv7wPOBzo7yyDgRX46g2NYydxXxA1vLGLz3gOM6ZvOhb1ak95IiP7hvySsfo+W+xbjyxXmlvVictmVzMsfSPJnSRz4eDb7DpQQ5RHG9k3ndxd2pXFCPU2lcay6X0b0ktcZEfMdn63uagnKGHPc3CgWOwo4x3k9EfgSf4IaBbyl/n6vBSLSRERaq2pYTjiUX1TKNa8tJPdACe/dMpjBbeNg/vOwaAIU5UDTDjDkQTy9r6KTL5Whm/bRYns+eUUlJMRGc3KrZIZ1a1m/VcqPR/vToFFLrvUt4+7V5/DghV1D716ZMSasBDtBVTb/fMuApLML/xwhEDC/vSPLWXdEghKRW4FbAWJjY4MY+rFTVX49eQU78g7y/q2D6XdoMfzjl7B/F3S5AAbeCh3PBo//FmA6kN4vkTH9jq0sflB5omuVaNqkZ7D95Ss5ZclE8or2sXp7AT3TXX6Q2BgT1oKdoH4y/3zgRmf2xaO6o+4kuQkASUlJIXk3/tNVO5m9Lps/XNiVfptfgS8fh5Y94Iq3oF2Y9Vr6vFz58rwad5t022nQ/TKiFk1gePRyPlvTwxKUMea4BHUUXxXzz+8WkdYAzs/yMtiH57d3pDvrwkpxaRmPfrKW7q2TufHQO/7k1Ptq+NkX4ZecjlbGIEhuwzWNllpVCWPMcQtagqpm/vnAeeyvBz52Xk8FrhO/wUB+ON5/+vfSLHYVFPNsl5V4vvkb9L0eRr0AMWFyL+l4eDzQfTSnHFpC9p49bMgudDsiY0wYC2YLqqr5558AhovIemCY8x78UwtvAjYArwB3BDG2oCgt8/HyVxu5rPU+Tlz8RzhxCFz0zOF7TRGh+2VEaSnDPNaKMsYcn6B9c6rqJlXt7SzdVfUxZ32Oqg5V1c6qOkxV9znrVVV/oaonqmpPVV0SrNiCZfa6bHbkHuBhXkbiG8Nlr4Inyu2w6ld6f2icwfhGS6yqhDFhTERGisgPIrLBeSSo4vY4EZnkbF8oIh2c9QNFZIWzfCcil9b2nBVF0J/2wTd5cSa3JH5N49xVMPIJSEpzO6T6J+Lv5itdzrbtO8jcV+R2RMaYoyQiUcAL+J9P7QZcJSLdKux2M5Crqp2AZ4AnnfWrgf6q2gcYCbwsItG1POcRLEHVkeyCYub9kMWdUR9CxmDoMcbtkNzT/VKi1Mt5UUtspl1jwtNAYIPTE1YCvI//WdVAo/A/ywowBRgqIqKqRarqddbH43/cqLbnPIIlqDry6aqdjPPMIrl0Dwx50N+SqANtM9ohIrVa2maESJXzNn2haQfGJSy2+1DGhKZoEVkSsNxaYXtVz6VWuo+TkPKBNAARGSQia4BVwO3O9tqc88ggj+6aTFVmrt7B32JnQMZp0PHMOjvvjqzMWj2HBM6zSKFABHqM5ZSvn2bbtk1kF/SlRUq821EZY37kVdX+wTq5qi4EuotIV2CiiEw/lvNYC6oO5B4oIWnbbFprNgy6ze1wQkPvcXjwMcrzLZ9/v9vtaIwxR6c2z6Ue3kdEooHGwBETwqnqWmA/0KOW5zyCJag6MHtdNuM9MylJbAUnX+h2OKGhWWc0fQBXxX3LZ6vC7nE2YyLdYqCziHQUkVhgHP5nVQMFPtM6FpjtVAfq6CQsRKQ9cDL+mStqc84jWIKqAwtXreOMqFXE9L0aompXdby295bCmfQexwm+reRtXkp2YbHb4Rhjasm5Z3QnMANYC0xW1TUi8icRucTZ7TUgTUQ2APfgn5kC4AzgOxFZgb+C0B2qureqc1YXh92DOk4+n9J4y3Si8UHPsbU+rrb3lkLmvtKx6H4ZOv0BLvV8zdQV5/OzM09wOyJjTC2p6jT8BRQC1z0U8LoYuLyS494G3q7tOatjLajj9P3OAob7vqYguRO07O52OKElMRXpMpIxMfOZumyr29EYY8KMJajj9N33axno+QHpOeaohoRHjN5X0UTzaLb7G37YZbX5jDG1Z118x6lk7WcAJPcezY6s34ffkPBg6zwcX0Izxvm+5MPlY3ng/K5uR2SMCRPWgjoO3jIf6Tnfkh/TAlrYF2+lomLw9L2GoZ5lzFu2ijJfSE7hZYwJQZagjsO6HfsYxCpy25xdZ5UjwpYz825lywlX/JkofAwp+oxGJ5wSOhUvjDEhzbr4jsP2VV/RQw5S0v08t0NxXw0z767Y+3vG6Rz+ddOLLH74/HoMzBgTrqwFdRw8m+bgxUNaz+FuhxLyvki6kNayjz7FC/Ek2lTwxpiaWYI6Di1yl7Et7iQkoYnboYS8ZfGD2OtJ45qomTTqOcztcIwxYcAS1DHKL9zPyd7/kd+sn9uhhAWfRDEr6ULOilpF33598NlgCWNMDSxBHaPNq74lTkqJO/EMt0MJG18kXkgxsdzedDFfb9jrdjjGmBBnCeoYFa3/BoCMPue6HEn4KIxqzFeJIxgd9Q1TvlridjjGmBBnCeoYJWcvZqu0JTm1tduhhJVpjS4jBi9dtv6LtTsL3A7HGBPCLEEdC1XaFa0ms1FPtyMJO7uj2/DR/5Rro77g7Tkr3Q7HGBPC7DmoY1C8ZwuNtZBDzfu4HUpYevTLA4zpIrRa+zo78vrRpkmC2yHVj/ws2PQV7FkLuVuh9CBEx0FyK2jdGzqcAalW8d2YcpagjsHOdfPoCCR0sBF8x2LFLh9FJ17AjRum8c85K7jv0lPdDil4ivbBivf8S7Yz9U1UHDRtD7FJUFoMm7+Gxa/6t6UPgEG3Q/fLwGMdHCayWYI6Bge3LKVUo0jv0t/tUMJW4og/4Ns4ncbLXmT7uX1o29BaUYW7Ye5fYfnb4C2G9IF4h/4J34lDkOYnExMTMLGlKuxdD/+bDsvfgQ9uhq+fhov/DhkDXbsEY9wW9AQlIlHAEmC7ql4kIh2B94E0YClwraqWiEgc8BbQD/+89leq6pZgx3csYvesZAMZdGne1O1QwlfLbhSfNIrrfpjO3z+bx+/GDXU7orpRUoT3q6eQBf8EXwnzG43gbd9I5mW1pnCDF9gGbCMtKZb0pgn0Sm/CoBNSGXLyCSSefjecehes+RC+eAReGwFn/AqG/AE8Ue5elzEuqI8W1N34p/dNcd4/CTyjqu+LyEvAzcCLzs9cVe0kIuOc/a6sh/iOjiot969jYcJpdPVEeIHY45Q48hG866fRdc3TbNwziBObN3I7pGO2/5CXRXM/o/ui39KydDsfl53G096xSOIJnNQymTFNEkhLisXjEUrLfOwuOMTWnAN8uCyLtxdsJSEmigt6tub2s0+gc8+x0HkEfP57+OYZ2LUKxr4O8VYiykSWoCYoEUkHLgQeA+4R/0x9Q4CrnV0mAo/gT1CjnNcAU4DnRURUNaRKDmjeVpK1kKJmPdwOJfylduTQoDu5dMHT/HnKZB68/cawmszR51MWbdnHh4s20XntP7iJqWRLM1474Vk6DDifj9s3pUlibLXn8Jb5WLo1l/+s2MHHK7bzwbIsRvVpw+8v6EqLS/4BbfrCtHvhrdFw7YeQYK12Uz9EZCTwLBAFvKqqT1TYXmmvl4gMB54AYoES4DeqOts55kugNXDQOc0IVc2uKoZgt6D+DvwWSHbepwF5qup13mcBbZ3XbYFMAFX1iki+s/8RJQdE5FbgVoDY2Or/5w+G3E3LSAVi0k+p989uiJKG3Mv+5e8yeuezfLZqBOf3Snc7pBoVl5bx0fLtvPr1JuL3rubvsS/RWTLZc9KVtLrsr9xcTUunbUY7dmRlVrrNk5BCyoDR/Mc7mv8s2shfrxrE5f1vQJJbweTrYOLFcP1/LUmZoHNuzbwADMf/Pb1YRKaq6vcBu1XV67UXuFhVd4hID2AGP37PA4xX1Vo9qR+0BCUiFwHZqrpURM6pq/Oq6gRgAkBSUlK9t64Ktq4kFUjt0Lu+P7phik0i4aIn6PnBjcz9+C/s7/IsjeJCc+zOvgMlvLNgK2/N30Le/iIeafoZV8dPRpKawagpNO9cc1X7HVmZNc66fLC0jIXffM1vP0hgwaYcHr10OInj/gX/GgeTroVrPoTo+v/jzESUgcAGVd0EICLv4+/lCkxQVfV6LQ/YZw2QICJxqnroaIMI5jjW04FLRGQL/kERQ/A3F5uISPk3UDqw3Xm9HcgAcLY3xt9sDCm6+3syfc3p2LaF26E0GFE9LiW3wwX8zPs+L02e6nY4P7F57wEe/M8qTntiFk/P/B8jm+9jRZsnuebge3h6XIbcMR9qkZxqKyEmit2THuRXwzrz0YrtXPXKQva1OQtGPQ9bvoapd/lH/hlz7KJFZEnAcmuF7Yd7tByBvV0/2cfpFSvv9Qo0BlhWITm9ISIrROQPUkOfftASlKo+oKrpqtoBGAfMVtXxwBxgrLPb9cDHzuupznuc7bND7f4TQGL+ejZ5MmjeKM7tUBoOEZpe/g+8sSlcsP4Rpi3f5HZEACzduo/b3l7CkL99yeTFWYzq1ZLF56zi0ew7aVS8G654G8a8Aompdf/h6uNXw07i5Wv6sW5nAWNfnMfODqPgnN/Byvd/fG7KmGPjVdX+AcuEuv4AEemOv9vvtoDV41W1J3Cms1xb3TnceBLwPvwDJjbgz7avOetfA9Kc9fcA97sQW/XKSkkr3kZO4olhdTM/LCQ1I27Mi3TzbKX0P3ezbme+K2GU+ZTpq3Zy6T+/ZcyL81mwaR+/OKcT829tx5P599F8wV/gpPPgjgXQ7ZKgxzOieyve+dkg9hQeYvwrC9nT95f+EX4zfgc7ltd8AmOOzeEeLUdgb9dP9qnY6+UMkPsIuE5VN5YfoKrbnZ+FwHv4uxKrVC8JSlW/VNWLnNebVHWgqnZS1cvLm36qWuy87+RsD40/owPlbCQaL8VNT3I7kgYp+uSR7B98L6NkLtNf+xPZBcX19tlFJV4mztvCuU99yc/fXUbO/hL+eEl35t93NvemzCLt7aGw9we47FV/y6lR83qLbUCHVF6/cQA78g9y7euLKRj5PCS1gH/fAIcK6y0OE1EWA51FpKOIxOLvBavY/15pr5eINAE+Be5X1W/LdxaRaBFp5ryOAS4CVlcXhNVSOQpFO/z/ljGturscScPVaMTvKWg3jLtLX+Olfz7F7iAnqeyCYv46Yx2n/mU2D09dQ1qjWF4c35c5957D9Z1LSHznIpjxAHQ8C+5YCL0uBxdazwM6pPLKdf3ZuGc/d3y0Be9lr0LeNpj5cL3HYho+557SnfhH4K0FJqvqGhH5k4iUdx1U1et1J9AJeMi517RCRFoAccAMEVkJrMDfAnulujhCc7hUiCrYupI4FZq2twQVNB4PKde8TeGrF/PA7qd58IVYbr7p55zUMrnmY2vJ51Pmb8rh3YVb+XzNbspUOa9bK245qyP92qeCtwS+fRq+fBJiE+GyV6CnO4kp0Jmdm/PY6J789oOVPNqiA48MvgPmPw/dRsEJZ7sam2l4VHUaMK3CuocCXhcDl1dy3KPAo1Wc9qgKmFqCOgreXd+zTVvQsXUzt0Np2GITSb7pQ4peuZDHcv7CQy/k0PviX3B5vww8x1G9Y3veQf773Q7eX7SNLTlFNEmM4cbTOzB+UHs6NEvy7/S/GfDZA7BvI3S9BC78GzQKnRGbVwzI4Ifdhbz2zWZ6XXYjl6V9BlPv9N8Ti01yOzxj6pQlqKOQkLuetWRwTmqi26E0fPGNSbx1Oofeu4bHt77EKx9vY9yi27hjWDfOPql5rQapqCqb9x7gq//t4ZOVO1m6NReA/u2b8qthJzGyRyviY5wadzuWw+zHYMNMSOsE46fU6dDxuvTA+SezZkc+D36ygUFj/krbjy6Fr/8GQx+q+WBjwoglqNoq89KkOJM98f2IibJbd/UiLpm466agnz3ALYtf4cw9a/nNmzfzh6bdGXpyS05p14QTmjUiJSEajwj5B0vZs/8Q63cXsnZnIQs35bAj338P6+RWyfzmvC5c1Ks17dOcloYqZC72Vx1fPwPim8CIR2HgbSH9IGx0lIe/X3kKI5+dy61fxTC155VEzfsH9BkPaSe6HZ4xdaZWCUpETg8cjVHVugYtbytRlFGc3MHtSCJLVAxy4VNw4rl0+e/d/JcH+VaG8fii83hzXsXnBn/UIjmOfu2bcse5zTijU7Mfu/AADubC9x/D4tdg10p/6aAhD/oTU3xKlecMJa0ax/P/xvTi1reX8ny7a7k7apq/a3L8ZLdDM6bO1LYF9Q+gby3WNViasxEBf/ePqRcV69alxMHvzozjlwNn8mnMF3y5K5qpu5oxe28aOXHpvPLCc6Q1iqVT80Y0TQpoAXlLIGspZC70t5S2fAM+L7To7r/H1OtKiKu7QRj1ZUT3Vlw9qB3PLtzG2DN+SdvFj/nvoZ10ntuhGVMnqk1QInIqcBrQXETuCdiUgr/CbcQ4sPN/NAISWtkzUPWlsrp1y4FfluUzpGg6Z3WYxdPpmcAuCg6tJmXeOEhs5h8s4PNCcT7kZ0L+dvCV+k+Q1hlOuwtOvhja9nV9ZN7xeuD8k5mzLptbfujPp6knIl88Ap2G2fxRpkGoqQUVCzRy9gv8E7OAH8sVRYSinevwaQItW4d+te2GrjCqMR8nj+Pj5HG09mbSpeR7vJ8/yZ3dUqFoL+RthahYiEmEtv2g+6XQuo9/OvXGVXcLVlRd5fFAbdIz2J657Tiu6Nglx8fw6Oge3DxxCdP63MKF6+6HlZOhz1WuxGNMXao2QanqV8BXIvKmqm6tp5hCkuZsZIu2ol1a+E6q1xDtjM5gZ3QGk6b/kTunfVjj/rVNOuVqqjwOMOm202p9vmAY2rUll/Ruw69XKcPa9CJuzuPQ4zKItnqRJrzV9h5UnIhMADoEHqOqQ4IRVCiKK9jCZu3AyNQEt0Mxx6E2012UczvxHI0HL+rK7HXZ/MMznnvz74Mlr8Pgn7sdljHHpbYJ6t/AS8CrQFnwwglR3kOkHNpFTuwZxEVb374JPS2S4/nl0E48Ps3LTe1PI3XuX+GUayHOWvwmfNX2gR6vqr6oqotUdWn5EtTIQknuFjz4OGhDzOuGJxoRqXExR+eG0zpyQrMkHioYBUU5sPQNt0My5rjUtgX1XxG5A3/59MMTT6nqvqBEFWpy/NXi1R6CrBs+b1jc23GVk8Rro3yQRmy0hz9c3I0b3zjAA60H0vbb52DAzyDGuqVNeKptgiovqf6bgHUKnFC34YSmQ9nricOGmJsqHEUyqbVaJnE4MpGf26UFZ53UnIcyL+A1fQSWvQWDbqv6YGNCWK0SlKp2DHYgoezArvUc1CRatWztdigmFB1jMgmW357XhYv+sYfM5n3I+Obv0O8GG9FnwlKt7kGJyHWVLcEOLlR8N/dTMrU5l4040+6ZmJDXo21jLu7dhkfyLoDCHbDiXbdDMuaY1LaLb0DA63hgKLAMeKvOIwpBbeOL2aCduPiRR4j2VJ3TI/qeiQkp/zf8JIat2kFmclcyvn0W+l5v1SVM2KlVC0pV7wpYbsFfgy8yxq/6fLRLLGY7LapNTsaEkg7Nkhg3sB1P5J8HuVtg3Sduh2TMUTvWb9wDQGTclzqQTbzHx04JnUnrjKmNXw7pzGwZSE5MG5j3vNvhGHPUansP6r8iMtVZPgV+wD/kvOHL9Vd4yo5q6XIgxhydFinxXDmwA/84OAKyFsG2hW6HZMxRqW0L6ingb87yOHCWqt4ftKhCSV55gmrlciDGHL3bzz6Rj/QciqKSYf4/3A7HhBERGSkiP4jIBhH5yfe9iMSJyCRn+0IR6eCsHy4iS0VklfNzSMAx/Zz1G0TkOalhdFlt70F9BazDX9G8KVByFNcZ1kr2bgIgJ9paUCb8tGoczyUDOvNmyVB07SeHHzo3pjoiEgW8AJwPdAOuEpFuFXa7GchV1U7AM8CTzvq9wMWq2hP/M7RvBxzzInAL0NlZRlYXR227+K4AFgGXA1cAC0UkIqbbKM7eTLY2safxTdi6/ZwTecc3gjKJhgUvuh2OCQ8DgQ2quklVS4D3gVEV9hkFTHReTwGGioio6nJV3eGsXwMkOK2t1kCKqi5QVcU/Cnx0dUHUtovv98AAVb1eVa9zgv9DdQeISLyILBKR70RkjYj80Vnf0WkObnCah7HO+kqbi27z5W4hU5sTF20j+Ex4atskgbP79eLjstPwLX/HP+W9iXTRIrIkYLm1wva2QOC8NFnOukr3UVUvkA+kVdhnDLBMVQ85+2fVcM4j1PZb16Oq2QHvc2px7CFgiKr2BvoAI0VkMP5m4DNOszAXfzMRqm4uuiq6IJMsbU6sJSgTxu4450Qmlp2Hx3sQltuDuwavqvYPWCbU9QeISHf83+PHXGurtt+6n4nIDBG5QURuAD4FplV3gPrtd97GOIsCQ/A3B8HfPBztvK60uVjL+IKjzEviwV1k+ppZgjJhLSM1kRN7nc5S7ULZolfA53M7JBPatgMZAe/TnXWV7iMi0UBj/I0XRCQd/0jv61R1Y8D+gVOSV3bOI1T7rSsinUTkdFX9DfAy0MtZ5gM1ZlwRiRKRFUA2MBPYCOQ5zUE4solXm+YiInJrebPU6/VW3Fy3CrLwUMbW4iQ8VsrIhLlbzzqBN0pHEJW3BTZ84XY4JrQtBjo7t2RigXHA1Ar7TOXHQuJjgdmqqiLSBH8j5n5V/bZ8Z1XdCRSIyGCn8XEd8HF1QdTULPg7UOCc/ENVvUdV78GfGf9e0xWqapmq9sGfKQcCJ9d0TC3OOaG8WRodXdtKTccobxsAWwqD/DnG1IOurVM4cOL57KEpZQtfdjscE8KcRsKdwAxgLTBZVdeIyJ9E5BJnt9eANBHZANwDlA9FvxPoBDwkIiucpbzSwR34J77dgL/BMr26OGr65m2pqqsqCX7V0QxiUNU8EZkDnAo0EZFo5x8gsIlX3lzMqthcdE15gsrzYTNBmYbglrO78M4bQ/j1xg/8Q85tjjNTBVWdRoVbOar6UMDrYvwjuyse9yjwaBXnXAL0qG0MNbWgmlSzrdpx1yLS3GnqISIJwHD8mXgO/uYg+JuH5U28SpuLNcQXVJrvH3CSua/YzTCMqTOnnpjGshaj8BKFLnrV7XCMqVZNCWqJiNxScaWI/Ayoacr31sAcEVmJvz9zpqp+AtwH3OM0C9PwNxOh6uaia0r2ZbFHUyjKj4yJg03DJyJcfs4AppUNxLvsbTi0v+aDjHFJTV18vwI+EpHx/JiQ+gOxwKXVHaiqK4FTKlm/Cf/9qIrrK20uuqlkXya7NJWygj1uh2JMnbmgRyvu/PQSLimZD6smQ/+b3A7JmEpV24JS1d2qehrwR2CLs/xRVU9V1V3BD89dUrCdnZqGt3Cv26EYU2eiozwMPmskq30dOPjNS+BuT7oxVaptLb45qvoPZ5kd7KBCRWzRTnZoGt6C7Jp3NiaMXD6gHf/2nEdC3g+QaVXOTWiyp0+rUlxArHc/2ZKG70C+29EYU6eS4qJJ6nsFhZpA0fxX3A7HmEpZgqpKgb/WYXFCK/wFMIxpWK46oxv/8Z1B7LqpUGQDgUzosQRVlQL/EPOy5GprGRoTtjJSE9nS/gqitYSSZVafz4QeS1BVcVpQUU3Sa9jRmPB13tBhLPV1pnj+qzZYwoQcS1BV8OVl4VMhqZklKNNwDejQlC+TLyLlwBZ081y3wzHmCJagqlCcs429NKZl0xS3QzGmTrTNaIeIHLF4PB7+8s5X5GkSU/541eH1bTPauR2uMTU+qBuxvLmZ7NBU2jaxmXRNw7AjK5MrX573k/U+VT7Keo5rOn7OLf+cRkFUEybddpoLERpzJGtBVUEKd7JT02jdJN7tUIwJKo8IMxPOJ4YyTt//mdvhGHOYJajKqBJXtJNdmkrrxmHUgvJE/6QLp7LFmIqKm3Riga8rQw9MR9QmMzShwbr4KlOcT2xZETlRzUmJD6N/Ip+30i6ciqz7xlQUG+1havR5PO77O90OLnM7HGMAa0FVrsA/RVVxYmtrcZiIsbbJ2eRoMmcXfuJ2KMYAlqAq5zwDpSn2kK6JHHHxCXws53KadxGtk6PcDscYS1CVciYqjGlqz0CZyPJVowuIFh+3D23vdijGWIKqjDc3kzJ7SNdEoOLk9nzj68lN3Q6Br8ztcEyEswRViYM5mWTTlNZNk90OxZh6JSJMjz+f9JhCshb9x+1wjItEZKSI/CAiG0TkJzOci0iciExyti8UkQ7O+jQRmSMi+0Xk+QrHfOmcc4WztKguBktQlSjLy2KnptKmsT0DZSLP+iZnsNvXmP3fTnA7FOMSEYkCXgDOB7oBV4lItwq73Qzkqmon4BngSWd9MfAH4N4qTj9eVfs4S7WT7VmCqkRU4Q5/grIqEiYCeaJjeWtHe04qWEj+9vVuh2PcMRDYoKqbVLUEeB8YVWGfUcBE5/UUYKiIiKoeUNVv8Ceq42IJqiJV4g/6q0i0shaUiVAvztqKAhtnvOB2KCY4okVkScBya4XtbYHMgPdZzrpK91FVL5APpNXis99wuvf+IDU8xxNGT6HWk4O5xPgOURDbgvgYG2prItPmLZksTxjBCds+xFvyBNGx9sdaA+NV1f4ufO54Vd0uIsnAB8C1wFtV7WwtqIqch3QPJbZxORBj3OUZ8DNSyWflFzaZYQTaDmQEvE931lW6j4hEA42BnOpOqqrbnZ+FwHv4uxKrZAmqonznv0GKJSgT2Xqfcxk7pCVxy99wOxRT/xYDnUWko4jEAuOAqRX2mQpc77weC8xWrXrWSxGJFpFmzusY4CJgdXVBWIKqQJ2HdGNTbT4cE9mioqLIOuFKupeuYv3qxW6HY+qRc0/pTmAGsBaYrKprRORPInKJs9trQJqIbADuAQ4PRReRLcDTwA0ikuWMAIwDZojISmAF/hbYK9XFEbR7UCKSgb9vsSWgwARVfVZEUoFJQAdgC3CFquY6N8ueBS4AioAbVLXeq1Yeys0iWj2kNLcWlIlgTmX85o0Tybw7ntkv/JKTXlrwk93apGewPXObCwGaYFPVacC0CuseCnhdDFxexbEdqjhtv6OJIZiDJLzA/6nqMueG2FIRmQncAMxS1Sech7/uB+7DP96+s7MMAl50ftar4r3bOEBTWjdtVN8fbUzoCKiMP2fXY1zVYilfPv8FGpN4xG5WGd8EU9C6+FR1Z3kLyLkhthb/sMTAsfMTgdHO61HAW+q3AGgiIq2DFV9VND/LP1FhOM0DZUwQfZV8ESlSRK+8WW6HYiJMvdyDckpgnAIsBFqq6k5n0y78XYBQu3H3iMit5WP3vV5vnccavb/8IV0bVmsMwObEXmwggwsPTcdX9T1wY+pc0BOUiDTCP979V6paELjNGfFxVL/xqjpBVfurav/o6DruoVQl4eBudtGMFsmWoIwBQIQZCRfQy7OJZvlr3I7GRJCgJihnKOEHwLuq+qGzend5153zs7wWU23G3QdXUQ7RWsKBuBZEeWyiQmPKLU4ZTpHGMeTAp1QzktiYOhW0BOWMynsNWKuqTwdsChw7fz3wccD668RvMJAf0BVYP5wh5qWNbASfMYGKoxoxK+YsLuRb9GC+2+GYCBHMFtTp+MtYDAkorX4B8AQwXETWA8Oc9+AfzrgJ2IB/bPwdQYytck4VCUmxeaCMqejrxpeQICUMyJ/hdigmQgRtmLlTzbaqfrKhleyvwC+CFU9t+PK34wHi0uwhXWMq2hbXmTXSmdFlnzO/dCyxVqvSBJlVkghwcO9WSjSKxs3rfXS7MWFhVqML6ezZTtu8JW6HYiKAJagAh/ZlsltTad0kye1QjAlJixudSy7JXHzov/h8NljCBJclqACav50dpNkzUCa8OGWJalrqQqnEMSN+JMNkKdEFmTUfYMxxsPmgAsTs38lObc/ZVkXChJOAskTVqauyRF+mXMKYgx8y4sDH/KtOzmhM5awFVc7nI+nQbvZKM5okxrgdjTla9diKiHR50c2ZG3MalzGH1I7d3Q7HNGDWgip3YA9R6qUooZV9kYWjem5FRLrZKZcydN/X3DJ8sNuhmAbMWlDlCvwP6XrtIV1jarQxrivfSyduabmWtTvy3A7HNFCWoMo5M+lGNflJfVpjTEUizEy+lBM9O5k7/X23ozENlCUohzfP34KKS2vvciTGhIdFSWex81A8Xbe+S1ZukdvhmAbIEpSjaM9WijWG1Gb2kK4xtVEmMby0pJSzPCv5zxdz3A7HNECWoByl+zLZoWm0bmpDzI2prRfn5VEqsbRY9Rq5B0rcDsfUIREZKSI/iMgGZ/bzitvjRGSSs32hM+8fIpImInNEZL+IPF/hmH4isso55jmpYUSaJSiHFNhMusYcrT1FyoGuVzJKvuKDuUvdDsfUERGJAl4Azge6AVeJSLcKu90M5KpqJ+AZ4ElnfTHwB+DeSk79InAL0NlZRlYXhyUoR+yBney0KhLGHB1PNP1v/yfRlFH85TN4YuKqfAatbYYVYQ4jA4ENqrpJVUuA94FRFfYZBUx0Xk8BhoqIqOoBp1h4ceDOzvx/Kaq6wCkO/hYwurog7DkogDIviSV72RfdgsRY+ycxptZ8Xvo/MY95e/7MdXzNtMc+JKVJaqW72jNoISVaRAIr/k5Q1QkB79sCgbWssoBBFc5xeB9V9YpIPpAG7K3iM9s65wk8Z7XDpq0FBVC4Ew8+DiW2cjsSY8LStJQrSJEizir8FJ/NuBsOvKraP2CZUPMh9c8SFByeqFBT7BkoY47F5rgurIjuxQ2eaeQW7Hc7HHP8tgMZAe/TnXWV7iMi0UBjIKeGcwbOBlvZOY9gCQrQPH9LNrppRg17GmOqMi3lClpJLv3yZ1krKvwtBjqLSEcRiQXGAVMr7DMVuN55PRaY7dxbqpSq7gQKRGSwM3rvOuDj6oKwBAUc3LsNgKTm9pCuMcdqZVw/Nno68jPPVPYWHHQ7HHMcVNUL3AnMANYCk1V1jYj8SUQucXZ7DUgTkQ3APcDhoegisgV4GrhBRLICRgDeAbwKbAA2AtOri8NGBOBPUF5NpEWz5m6HYkz4EuGTlCu4O+9JehTMZVfKeXis8HLYUtVpwLQK6x4KeF0MXF7FsR2qWL8E6FHbGKwFBXjz/A/ptrWHdI05LgsSzmKbJ51fyBRrRZnjZgkKiCrcwU5NpU0TS1DGHA+VKD5KGU8XTxY9C+bavShzXCxBAQkHd7JbmpGWFOt2KMaEvfJW1M/lA/bkWyvKHLuITVBtM9ohIiTECEnePLYWRuPxeGwWVmOOU3kr6mRPJj0K5lLms1aUOTYRO0hiR1YmV748j1be7ZB9MwWNO1U5I6s9AW/M0VmQcBaXFbzLLzwfcEvembRNTXI7JBOGgtaCEpHXRSRbRFYHrEsVkZkist752dRZL05l2w0islJE+gYrrorSyvYAkBNlI/iMqSuBrai+hXMoLfO5HZIJQ8Hs4nuTn1aqvR+YpaqdgVn8OG7+fH6sbnsr/oq39SLVmw1AbnSL+vpIYyLC/ISz2BzVgXuiJpOdW+h2OCYMBS1BqepcYF+F1YHVbyfyYyXbUcBb6rcAaOJUvg26JqX+BJUfYy0oY+qSShSTGt9Me082Qw5MIyrF/h8zR6e+B0m0dMpdAOwCWjqvK6ucW2lhPBG5VUSWiMgSr9d73AE19+5itzZBYmyaDWPq2oq4/qyK6cVd0R+Rce7Vbodjwoxro/icmk1HPbxHVSeUV+CNjj7+MR7Ny3aRqS2Ii47YAY3GBI8I/2p8M82kgF/1yGfR5oqdKsZUrb6/lXeXd905P7Od9bWpnBsUrcp2k6nNibUEZUxQbIrtwrz4M/mZ5xOe+/hrG3Zuaq2+v5UDq99ez4+VbKcC1zmj+QYD+QFdgUHj0TKak8MOaWE1w4wJokkpNxKjpVya8wqTl2TWfIAxBHeY+b+A+UAXp5rtzcATwHARWQ8Mc96DvyDhJvwVbl/BX/E26NLK9hCNj2xPy5p3NsYcs93RbXh6fjFjor5mxmdTyT9Y6nZIJgwE7UFdVb2qik1DK9lXgV8EK5aqtCjbBcCeKJtJ15hge3TuIf5vRDvu2f8qf//8LB4e1cvtkEyIi+gbL81K/Qlqb4wlKGOC7UApxIx8jF6ezRQvmsjybbluh2RCXEQnqKalO/Gqh/wYe0jXmHrRcyze9MHcFzOJx6bMswoTploRnaCae3exU9OIibUq5sbUCxGiL3qKFCniin0vM2HuJrcjMiEsohNUi7LdZGlz4m2IuTHB54n2zxDQuhdPzC3iiuivmPfZv4hJbXvE7AFtM9q5HakJERFbzRyglW83X9KL6ChLUMYEnc97eMaAtVrC9uyf80T8RP531z85oU3zw1Pb2OwBplzEfjPHRUFzctllQ8yNqXelEssrTX5Nhuzh52XvsSOv2O2QTAUiMlJEfnBmmbi/ku1xIjLJ2b5QRDoEbHvAWf+DiJwXsH6LiKwSkRUisqSmGCI2QZ3Q1H/puzw2gs8YN6yL68GMxIu4IXoG6flL2H/o+GtrmrohIlHAC/hnmugGXCUi3SrsdjOQq6qdgGeAJ51juwHjgO74Z7T4p3O+cueqah9V7V9THBGboE5q5v/32hmT7nIkxkSu91J+xs6oNjwT8yJ7sndZGaTQMRDYoKqbVLUEeB//rBOBAmenmAIMFX8/7SjgfVU9pKqb8RdgGHgsQURsgjq5bQoAeyxBGeOaQ554nm96P2lSwMNMYOveA26HFCmiy2eFcJZbK2yvzQwTh/dRVS+QD6TVcKwCn4vI0ko+86dB1vZqGpqTWyWyR8EXl+x2KMZEtM2xnZmcfD3jC19jBbNZ7HZAkcFbmy62IDhDVbeLSAtgpoisc+YOrFTEtqA6pwqbtDUJMVE172yMCapPGo1hVWxvRnu+djsU41ebGSYO7yMi0UBjIKe6Y1W1/Gc28BE1dP1FbILq1KiYrbS2IebGhAAVD8+m/p7H0h53OxTjtxjoLCIdRSQW/6CHqRX2CZydYiww26mrOhUY54zy6wh0BhaJSJKIJAOISBIwAlhdXRCR2cVXnE/L2GKyaON2JMYYx35PitshGIeqekXkTmAGEAW8rqprRORPwBJVnQq8BrwtIhuAffiTGM5+k4HvAS/wC1UtE5GWwEfO827RwHuq+ll1cURkgtKcjQiwM6rSWeWNMW5yKk7URpv0DLZnbgtyQJFJVafhnwopcN1DAa+LgcurOPYx4LEK6zYBvY8mhohMUIXb15EC7I7NqHFfY0w9C6g4UROrOtGwReQNmLydG/CpkBvb2u1QjDHGVCEiE9TcltfT59AEouMS3A7FGGNMFSIyQbVIjmPX/1YSayP4jDEmZEXkN/SI7q3Y89Fjtb4Ra4wxpv5FZIIyxhgT+ixBGWOMCUmWoIwxxoQkS1DGmAavbUa7I6aVr26xKedDR0Q+qGuMaSCOouqEPfwbfkIqQYnISOBZ/LWfXlXVJ1wOyRgTympZdcKSTngKmS6+Wk4xbIwxweW0yqwr0H2h1II6PMUwgIiUTzH8vatRGWMiS21bZT8/q9bdi1ExcZSVHqqz/SAyCuWKf/oO94nIWGCkqv7MeX8tMEhV76yw361A+VTBfYGDx/iR0fhLwTd0kXCdkXCNEBnXaddYPxJUNWR60KoSSi2oWlHVCcCE4z2PiCxxacrjehUJ1xkJ1wiRcZ12jSZQKGXQ2kwxbIwxJkKEUoKqzRTDxhhjIkTIdPFVNcVwED/yuLsJw0QkXGckXCNExnXaNZrDQmaQhDHGGBMolLr4jDHGmMMsQRljjAlJEZmgRGSkiPwgIhtE5H6346krIvK6iGSLyOqAdakiMlNE1js/m7oZ4/ESkQwRmSMi34vIGhG521nfYK5TROJFZJGIfOdc4x+d9R1FZKHzezvJGUwU1kQkSkSWi8gnzvuGeI1bRGSViKwQkSXOugbz+xpMEZegGnhJpTeBkRXW3Q/MUtXOwCznfTjzAv+nqt2AwcAvnP9+Dek6DwFDVLU30AcYKSKDgSeBZ1S1E5AL3OxeiHXmbmBtwPuGeI0A56pqn4DnnxrS72vQRFyCIqCkkqqWAOUllcKeqs4F9lVYPQqY6LyeCIyuz5jqmqruVNVlzutC/F9ubWlA16l++523Mc6iwBBgirM+rK8RQETSgQuBV533QgO7xmo0mN/XYIrEBNUWyAx4n+Wsa6haqupO5/UuoKWbwdQlEekAnAIspIFdp9P1tQLIBmYCG4E8VS0vkdMQfm//DvwW8Dnv02h41wj+Py4+F5GlTqk2aGC/r8ESMs9BmeBTVRWRBvFcgYg0Aj4AfqWqBYFFOxvCdapqGdBHRJoAHwEnuxtR3RKRi4BsVV0qIue4HE6wnaGq20WkBTBTRNYFbmwIv6/BEoktqEgrqbRbRFoDOD+zXY7nuIlIDP7k9K6qfuisbnDXCaCqecAc4FSgiYiU/1EZ7r+3pwOXiMgW/N3sQ/DPBdeQrhEAVd3u/MzG/8fGQBro72tdi8QEFWkllaYC1zuvrwc+djGW4+bcp3gNWKuqTwdsajDXKSLNnZYTIpIADMd/r20OMNbZLayvUVUfUNV0Ve2A///B2ao6ngZ0jQAikiQiyeWvgRHAahrQ72swRWQlCRG5AH//d3lJpcfcjahuiMi/gHOAZsBu4GHgP8BkoB2wFbhCVSsOpAgbInIG8DWwih/vXfwO/32oBnGdItIL/43zKPx/RE5W1T+JyAn4WxupwHLgGlWt3eRBIczp4rtXVS9qaNfoXM9Hztto4D1VfUxE0mggv6/BFJEJyhhjTOiLxC4+Y4wxYcASlDHGmJBkCcoYY0xIsgRljDEmJFmCMsYYE5IsQRljjAlJlqCMMcaEpP8P6/XXZlEHoIEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from functools import lru_cache, partial\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "import scipy.special\n",
    "import seaborn as sns\n",
    "from tqdm import tqdm\n",
    "\n",
    "rng = np.random.default_rng(0)\n",
    "\n",
    "real_N = 5\n",
    "Li = np.array([123, 234, 543, 321, 599]) * 4\n",
    "L = sum(Li)\n",
    "real_pi = Li.astype(float) / L\n",
    "real_lam = np.array([3, 6, 12, 20, 30])\n",
    "\n",
    "data = np.concatenate([\n",
    "    rng.poisson(3, Li[0]),\n",
    "    rng.poisson(6, Li[1]),\n",
    "    rng.poisson(12, Li[2]),\n",
    "    rng.poisson(20, Li[3]),\n",
    "    rng.poisson(30, Li[4]),\n",
    "])\n",
    "\n",
    "N = 5\n",
    "p = np.empty((L, N))\n",
    "q = np.empty((L, N))\n",
    "sumpg = np.empty(L)\n",
    "pi = rng.dirichlet(np.ones(N, float))\n",
    "lam = rng.exponential(10, N)\n",
    "\n",
    "\n",
    "def real_pdf(x, k):\n",
    "    return np.power(real_lam[k],\n",
    "                    x) / sp.special.factorial(x) * np.exp(-real_lam[k])\n",
    "\n",
    "\n",
    "@lru_cache\n",
    "def fact(x):\n",
    "    return sp.special.factorial(x)\n",
    "\n",
    "\n",
    "@lru_cache\n",
    "def pdf(x, k):\n",
    "    return np.power(lam[k], x) / fact(x) * np.exp(-lam[k])\n",
    "\n",
    "\n",
    "def em_step():\n",
    "    # E step\n",
    "    for i in range(L):\n",
    "        sumpg[i] = 0\n",
    "        for k in range(N):\n",
    "            sumpg[i] += pi[k] * pdf(data[i], k)\n",
    "\n",
    "    for k in range(N):\n",
    "        for i in range(L):\n",
    "            p[i, k] = pi[k] * pdf(data[i], k) / sumpg[i]\n",
    "    p[:] = p / p.sum(axis=1).reshape(-1, 1)\n",
    "\n",
    "    # M step\n",
    "    sump = p.sum(axis=0)\n",
    "    sumpx = p.T @ data\n",
    "    lam[:] = sumpx / sump\n",
    "    pi[:] = sump / L\n",
    "\n",
    "\n",
    "for t in tqdm(range(100)):\n",
    "    em_step()\n",
    "\n",
    "\n",
    "def mixture(x, pdfs, p):\n",
    "    s = 0\n",
    "    for pp, ppdf in zip(p, pdfs):\n",
    "        s += pp * ppdf(x)\n",
    "    return s\n",
    "\n",
    "\n",
    "est_mixture = np.vectorize(\n",
    "    partial(mixture, pdfs=[partial(pdf, k=i) for i in range(N)], p=pi))\n",
    "real_mixture = np.vectorize(\n",
    "    partial(mixture,\n",
    "            pdfs=[partial(real_pdf, k=i) for i in range(real_N)],\n",
    "            p=real_pi))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "sns.histplot(data, ax=ax)\n",
    "xx = np.linspace(0, 40, 1000)\n",
    "ax2.plot(xx, est_mixture(xx), label='est')\n",
    "ax2.plot(xx, real_mixture(xx), label='real')\n",
    "ax2.legend()\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8e1db19f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>esti_lam</th>\n",
       "      <th>esti_pi</th>\n",
       "      <th>real_lam</th>\n",
       "      <th>real_pi</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.025379</td>\n",
       "      <td>0.061927</td>\n",
       "      <td>3</td>\n",
       "      <td>0.067582</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5.864372</td>\n",
       "      <td>0.134654</td>\n",
       "      <td>6</td>\n",
       "      <td>0.128571</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>12.277207</td>\n",
       "      <td>0.311448</td>\n",
       "      <td>12</td>\n",
       "      <td>0.298352</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>20.969611</td>\n",
       "      <td>0.185048</td>\n",
       "      <td>20</td>\n",
       "      <td>0.176374</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>30.361888</td>\n",
       "      <td>0.306923</td>\n",
       "      <td>30</td>\n",
       "      <td>0.329121</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    esti_lam   esti_pi  real_lam   real_pi\n",
       "0   3.025379  0.061927         3  0.067582\n",
       "1   5.864372  0.134654         6  0.128571\n",
       "2  12.277207  0.311448        12  0.298352\n",
       "3  20.969611  0.185048        20  0.176374\n",
       "4  30.361888  0.306923        30  0.329121"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.concat(\n",
    "    [\n",
    "        pd.DataFrame(list(zip(lam, pi)), columns=[\n",
    "            \"esti_lam\", \"esti_pi\"\n",
    "        ]).sort_values(\"esti_lam\").reset_index(drop=True),\n",
    "        pd.DataFrame(list(zip(real_lam, real_pi)),\n",
    "                     columns=[\n",
    "                         \"real_lam\", \"real_pi\"\n",
    "                     ]).sort_values(\"real_lam\").reset_index(drop=True),\n",
    "    ],\n",
    "    axis=1,\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "jupyter-home",
   "language": "python",
   "name": "jupyter-home"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}