Skip to content

Instantly share code, notes, and snippets.

@andiac
Last active June 7, 2021 08:46
Show Gist options
  • Select an option

  • Save andiac/c9973bbc4e21b1027f712c0970cf3aba to your computer and use it in GitHub Desktop.

Select an option

Save andiac/c9973bbc4e21b1027f712c0970cf3aba to your computer and use it in GitHub Desktop.
KL divergence between a gaussian and a mixture of gaussian
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.stats import norm, bernoulli\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"p_mu = 2\n",
"p_sigma = 0.5\n",
"\n",
"q_mu1 = 0\n",
"q_sigma1 = 1\n",
"q_mu2 = 5\n",
"q_sigma2 = 1\n",
"w1 = 0.3\n",
"w2 = 0.7"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"N_MONTE = 30000"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvFUlEQVR4nO3dd5yU5bn/8c812/vCFgSWJoKKMTZERQ0aS7D8RBNNLIkpJoYkaso5OdFUE0+OyTE9khAs0USj8RhjSMQae9QIVkAEV6Qs7LIF2Mb2uX5/3DMwLLvs7O7MPM/MXO/XixdTnpnnUna/c8/93EVUFWOMMckv4HUBxhhjYsMC3RhjUoQFujHGpAgLdGOMSREW6MYYkyIyvTpxeXm5Tp061avTG2NMUnrllVcaVbVioOc8C/SpU6eyYsUKr05vjDFJSUQ2DvacdbkYY0yKsEA3xpgUYYFujDEpwgLdGGNShAW6McakiKgCXUTmi8haEakWkWsHeL5ERP4uIm+IyGoR+XTsSzXGGLM/Qwa6iGQAi4CzgFnAJSIyq99hXwLeUtUjgFOAn4pIdoxrNcYYsx/RtNDnANWqul5Vu4F7gQX9jlGgSEQEKAS2A70xrdSYCN29Qf740kaeeGub16UY4xvRTCyaCGyOuF8DHNfvmJuBpcBWoAj4mKoG+7+RiFwJXAkwefLkkdRrDADffnAl962oAWDRpUdzzvvHe1yRMd6LpoUuAzzWf1eMDwGvAxOAI4GbRaR4nxepLlHV2ao6u6JiwJmrxgxpbV0r962o4VNzp3L4xBL+Z9ka+oK2UYsx0QR6DTAp4n4VriUe6dPAA+pUA+8Bh8SmRGP2ds/Lm8jODPCV02fwhVOms2VnB8+ua/C6LGM8F02gLwdmiMi00IXOi3HdK5E2AacBiMg44GBgfSwLNQYgGFQeWVXHvJkVlOZnc/qh4yjMyeSxt+q8Ls0Yzw0Z6KraC1wFPAqsAe5T1dUislBEFoYOuwGYKyIrgX8C31DVxngVbdJXdUMbdS2dnHHoOACyMwOcPKOcJ9+ux/bHNekuqtUWVXUZsKzfY4sjbm8Fzoxtacbsa/mG7QDMmTZ292MnzSjn4VV1bN7eweSyfK9KM8ZzNlPUJJUVG3ZQXpjDlIjgPnJSKQCvbd7hUVXG+IMFukkqyzds59ipY3BTHpyDxxWRl5XBa5t2eleYMT5ggW6Sxo72bmp2dOxukYdlZgQ4vKqE1zbv9KQuY/zCAt0kjbfrWgE4ZPw+Uxw4oqqENbUtNh7dpDULdJM03q5rAeDQA4r2eW7muCK6e4NsbGpPdFnG+IYFukkaa+taGZOfRUVRzj7PzRznQn7dtrZEl2WMb1igm6Sxpq6VQw4o3uuCaNhBlYUArNvWmuiyjPENC3STFFSVd+vbmDmucMDnC3IymTQ2zwLdpDULdJMUGtu6aevqZWp5waDHzKwsorreulxM+rJAN0khfLFzatnggT6lrICNTbtsCQCTtizQTVLY0LQLYL8t9Cll+XT09NHQ2pWosozxFQt0kxQ2NLaTERCqxuQNekx4HZeN23clqixjfMUC3SSFDU3tVI3JIytj8B/ZKWNDgd5kgW7SkwW6SQobm3YxZT/95wBVY/IJCGyyyUUmTVmgG99TVTY0tTNtiKVxszMDjC/Jsy4Xk7aiCnQRmS8ia0WkWkSuHeD5r4vI66E/q0SkT0TGDvRexgzXjl09tHb2DtlCB3dh1LpcTLoaMtBFJANYBJwFzAIuEZFZkceo6k2qeqSqHglcBzyjqtvjUK9JQzU7XEBPGjv05hVTyvLZZC10k6aiaaHPAapVdb2qdgP3Agv2c/wlwD2xKM4YgK07OwEYX5I75LGTxxawvd1NQjIm3UQT6BOBzRH3a0KP7UNE8oH5wF8Gef5KEVkhIisaGmyXdhOdrTs7AJhYOviQxbAJpS70a0OvMSadRBPo+66EBINNxft/wL8G625R1SWqOltVZ1dUVERbo0lzW3d2kJeVQWl+1pDHTgiF/hYLdJOGogn0GmBSxP0qYOsgx16MdbeYGNva3MGE0twBV1nsL9wtU9vcGe+yjPGdaAJ9OTBDRKaJSDYutJf2P0hESoB5wN9iW6JJd1t2du5ueQ9lXHEuItblYtLTkIGuqr3AVcCjwBrgPlVdLSILRWRhxKEXAI+pqs3qMDG1dWcHE0qiC/SsjADjinLZai10k4YyozlIVZcBy/o9trjf/TuAO2JVmDEAXb1usa1oW+gA40tzqW22FrpJPzZT1Pjatma3cmJ49Eo0JpTkUbvTWugm/VigG1/bMowhi2HjS3LZ2txh66KbtGOBbnwtPAZ9/LC6XPLo7AmyY1dPvMoyxpcs0I2v7Q70KGaJhk0IHbvVRrqYNGOBbnxta3MHZQXZ5GZlRP2acGvexqKbdGOBbnxtW0sXBwyjdQ57Wug20sWkGwt042sNrV1UFOUM6zXlhTlkZcjuRb2MSRcW6MbXGlq7qCgcXqAHAkJlUS7bWizQTXqxQDe+FQwqjW3Db6EDVBbnUN9qgW7SiwW68a2dHT30BnVEgT6uKJf6lq44VGWMf1mgG99qaHWBPNIWunW5mHRjgW58KxzolUXDG+XiXpNDS2cvnT19sS7LGN+yQDe+Fe4DH1kL3X0IhD8UjEkHFujGt0bV5RJ6jXW7mHRigW58q6G1i7ysDAqyo58lGhbupqm3FrpJIxboxrcaQkMWo9l6rr9xxa6FXm8tdJNGogp0EZkvImtFpFpErh3kmFNE5HURWS0iz8S2TJOORjJLNGxMfjaZAWGbtdBNGhlyxyIRyQAWAWfgNoxeLiJLVfWtiGNKgd8A81V1k4hUxqlek0YaWruYXlE4otcGAkJFUY6NRTdpJZoW+hygWlXXq2o3cC+woN8xlwIPqOomAFWtj22ZJh01jHCWaFhlca7NFjVpJZpAnwhsjrhfE3os0kxgjIg8LSKviMjlA72RiFwpIitEZEVDQ8PIKjZpoau3j527enaPVhmJSmuhmzQTTaAPdEWq/95emcAxwDnAh4DviMjMfV6kukRVZ6vq7IqKimEXa9JHY1s3MLIhi2GVRbaei0kvQ/ah41rkkyLuVwFbBzimUVXbgXYReRY4AlgXkypN2hnNGPSwccW57NjVQ1dvHzmZwx/6aEyyiaaFvhyYISLTRCQbuBhY2u+YvwEni0imiOQDxwFrYluqSSexCPRwd43NFjXpYsgWuqr2ishVwKNABnC7qq4WkYWh5xer6hoReQR4EwgCt6rqqngWblJbTAI9PBa9tYuqMfkxqcsYP4umywVVXQYs6/fY4n73bwJuil1pJp2FA72sYDQt9NBsUbswatKEzRQ1vtTQ1smY/CyyM0f+I7qnhW4XRk16sEA3vtTQ2jWiZXMjlRXkEBBotD50kyYs0I0vjWbaf1hGQCgrzLEFukzasEA3vlQfg0AHqCjMsVEuJm1YoBvfUdWYtNDB9aM3tFmgm/RggW58p7Wrl67eIBWFsWmh2ygXky4s0I3vxGIMelhFUQ6NbV0Eg/1XqzAm9VigG9+JdaD3BpWdHT2jfi9j/M4C3fhOLAN9z1Z0NhbdpD4LdOM74UAfzdK5YRW2notJIxboxnca2rrIyhBK8rJG/V4W6CadWKAb36lv6aKicGSbQ/cXbuXb5CKTDizQje+Mduu5SAU5meRnZ1gL3aQFC3TjO7GaVBRWUWSzRU16sEA3vhPrQLet6Ey6sEA3vtIXVLa3d8VklmiYtdBNuogq0EVkvoisFZFqEbl2gOdPEZFmEXk99Oe7sS/VpIOm9i6CChXFo1s6N5It0GXSxZA7FolIBrAIOAO3GfRyEVmqqm/1O/Q5VT03DjWaNLJ7UlEMW+iVxbm0dPbS2dNHbpZtFm1SVzQt9DlAtaquV9Vu4F5gQXzLMumqPoazRMPCHw7WSjepLppAnwhsjrhfE3qsvxNE5A0ReVhEDhvojUTkShFZISIrGhoaRlCuSXWxnCUatntykS2ja1JcNIE+0OyO/kvXvQpMUdUjgF8DDw70Rqq6RFVnq+rsioqKYRVq0kM40MtjfFE08r2NSVXRBHoNMCnifhWwNfIAVW1R1bbQ7WVAloiUx6xKkzYaWrsoyskkLzt2fd02W9Ski2gCfTkwQ0SmiUg2cDGwNPIAETlAQvO0RWRO6H2bYl2sSX2xnCUaNrYgGxFroZvUN+QoF1XtFZGrgEeBDOB2VV0tIgtDzy8GLgS+ICK9QAdwsarajgJm2BpauyiPcaBnZgQoK8i2QDcpb8hAh93dKMv6PbY44vbNwM2xLc2ko8bWLmZNKI75+1YU5dJgs0VNirOZosZXYj3tP8xmi5p0YIFufKOju4/Wrt74BLrNFjVpwALd+EY8ZomGVRbn0NDWhV3aManMAt34RkOb6+OOVwu9p0/Zucs2izapywLd+EYsN4fuz2aLmnRggW58I56BvntyUYsFukldFujGNxpauwgIlBXEs4VuQxdN6rJAN77R0NZFWWEOGYHRbw7dn63nYtKBBbrxjYbW2O5UFKkwJ5O8rAzrcjEpzQLd+EZ9nCYVAYiIm1xkF0VNCrNAN74Rr1miYTZb1KQ6C3TjC8Gg0hiHlRYjVRTm2BK6JqVZoBtfaO7ooadP49aHDqHZohboJoVZoBtfCPdtx7uF3tzRQ1dvX9zOYYyXLNCNL8RjL9H+wh8WjW3dcTuHMV6yQDe+EM9ZomGVxeHZoja5yKSmqAJdROaLyFoRqRaRa/dz3LEi0iciF8auRJMOEhHoFYW5e53LmFQzZKCLSAawCDgLmAVcIiKzBjnux7it6owZlvrWTnKzAhTmRLWJ1ojYAl0m1UXTQp8DVKvqelXtBu4FFgxw3NXAX4D6GNZn0kR4DHpor/G4KCt0m0XbbFGTqqIJ9InA5oj7NaHHdhORicAFwGL2Q0SuFJEVIrKioaFhuLWaFNbQFr9p/2FZGQHG5mdbC92krGgCfaAmU/9tX34BfENV9zseTFWXqOpsVZ1dUVERZYkmHdS3dFFZlBv389hsUZPKoumwrAEmRdyvArb2O2Y2cG/o63I5cLaI9Krqg7Eo0qS++tYuTpheFvfzVBTZbFGTuqIJ9OXADBGZBmwBLgYujTxAVaeFb4vIHcA/LMxNtDp7+mju6InrGPSwiqIc1je0x/08xnhhyEBX1V4RuQo3eiUDuF1VV4vIwtDz++03N2YoeyYVJa7LRVXjegHWGC9ENUZMVZcBy/o9NmCQq+qnRl+WSSfhLpCK4vi30CuLcunuC9Lc0UNpfnbcz2dMItlMUeO5hlY3czNRXS7unNaPblKPBbrxXH0iu1wKLdBN6rJAN56rb+kiIyCUFcS/CyTcQreRLiYVWaAbz9W3dlJemE0gDptD9xdeoMta6CYVWaAbz9W3JmZSEUBRTiY5mQGbLWpSkgW68ZybJRr/C6KwZ7NoW0LXpCILdOO5+tau3V0hiVBZlGMtdJOSLNCNp3r7gjS1d1GRoC4XsPVcTOqyQDeeamrvRjUxY9DDLNBNqrJAN54Kr02eyECvLMplx64eunuDCTunMYlggW48VR+eJVqc2C4XgEbrRzcpxgLdeGrPLNEEdrnYbFGToizQjafCXS7lcd6tKFJ4RI3NFjWpxgLdeKq+tZOxBdlkZybuR9EW6DKpygLdeMrNEk1c6xygrMAC3aSmqAJdROaLyFoRqRaRawd4foGIvCkir4c2gT4p9qWaVFTf2rW7xZwo2ZkBxhZk774ga0yqGDLQRSQDWAScBcwCLhGRWf0O+ydwhKoeCXwGuDXGdZoUVd/SmbB1XCJVFNpYdJN6ommhzwGqVXW9qnYD9wILIg9Q1TZV1dDdAkAxZgi9fUHqW7sYX+JBoNv0f5OCogn0icDmiPs1ocf2IiIXiMjbwEO4Vrox+9XY1k1fUDnAo0APj7AxJlVEE+gDLVK9TwtcVf+qqocA5wM3DPhGIleG+thXNDQ0DKtQk3pqmzsAPGmhhxfo2vPF0pjkF02g1wCTIu5XAVsHO1hVnwWmi0j5AM8tUdXZqjq7oqJi2MWa1FLX7C5KetVC7+4N0tLZm/BzGxMv0QT6cmCGiEwTkWzgYmBp5AEicpCISOj20UA20BTrYk1qqQ0F+viSvISfe89YdBvpYlJH5lAHqGqviFwFPApkALer6moRWRh6fjHwEeByEekBOoCPqX2XNUOoa+kkJzPAmPyshJ87cm/RgyqLEn5+Y+JhyEAHUNVlwLJ+jy2OuP1j4MexLc2kutrmTsaX5BL6cpdQlck8W7R+Dbz6B3j3SWiqhmAfjJkCU0+Goz4Ok4/3ukLjkagC3Zh4qGvu8KT/HPas7hjux08KHTvh0W/C63+CjCyYNg9mfggkwwX76r/Ca3+Eg06Hc38BpZOGekeTYizQjWdqmzs5dupYT85dnJtFUU7m7n5839vyKtz3SWjZAnOvghO/CgVlex/TvQuW3wrP/Bh+eyJcsBgOOdubeo0nbC0X44lgUNnW0ulZCx1gfGkuW3d2eHb+qL37FNxxrrt9xWNw5n/vG+YA2flw4jWw8DkYOw3+fBm8cmdiazWeskA3nmhq76anTz0Zgx42oTSPrc0+D/QNz8OfPgpjpsJnH4eq2UO/ZuyB8OllMP2D8Pdr4LW74l6m8QcLdOOJ3WPQE7hTUX8TSvPYutPHXS71a+DeS2HMNPjUP6DogOhfm10AF98DB54KS6+BtY/Er07jGxboxhN7Zokmfgx62MTSPLa3d9PR3edZDYPq2AF3fxQy8+Dj90P+CK41ZGbDx+6CA94HD3wOmt6NfZ3GVyzQjSfqWrybJRo2odSd23fdLqrw4JegtRYu/hOUTh75e+UUulAPZMKfPw7d7bGr0/iOBbrxRG1zJ1kZQllBtmc1TAh9O/DdhdF/L4a1D8EZ34eqY0b/fqWT4cLbXBfOo98c/fsZ37JAN56oa+5kXHEugUDiJxWFTSh1gV7rp370+rfh8e/CzLPg+C/G7n2nfxDmXg2v3AHVT8TufY2vWKAbT2zZ2bG7heyVA0pyEXG1+EKwD5ZeBdmFcN6vIdYzaE/9FlQcAn+72k1SMinHAt14YsuODqrGeBvoWRkBKoty/NPl8vISqFkOZ/0YCuOwGmlWLpz/W2jbBo99O/bvbzxngW4SrqcvSG1zB1Vj870uxT9j0Xduhn/+AGacCYdfFL/zTDwaTviiWyJg8/L4ncd4wgLdJFztzk6CiuctdPDRWPTHv+NGt5zz09h3tfQ37xtQNB4e+prr5jEpwwLdJFzNjl2APwJ9YmkeW3Z2eLtz0XvPuYW1TvrK6IYoRiunCD70Q6h7E1bcHv/zmYSxQDcJV7PDdXFMGuODLpeSXLp7gzS1d3tTQF8vPHItlEyCudck7ryHfdit1vjkDbBre+LOa+LKAt0kXM2OXQTE20lFYeGhi55dGH31Dti2Cs68wS2ulSgiMP9H0NUKz96UuPOauIoq0EVkvoisFZFqEbl2gOcvE5E3Q39eEJEjYl+qSRU1OzoYX5JHVob37Ymq0LeEzds9CPTOFnjyhzDlJJh1fuLPP24WHHkZvHwLbF+f+PObmBvyN0pEMoBFwFnALOASEZnV77D3gHmq+n7gBmBJrAs1qWPzjl1M9EH/OcDkMhfoG7d7MCX+xZuhY7trnXuwaxPgxqZnZLkRNibpRdNEmgNUq+p6Ve0G7gUWRB6gqi+o6o7Q3ZeAqtiWaVJJjQ/GoIcV5mRSXpjNpqZdiT1xeyO8uAhmLXBDCb1SPN7NIF39V6hZ4V0dJiaiCfSJwOaI+zWhxwZzBfDwQE+IyJUiskJEVjQ0NERfpUkZ3b1B6lo6d3d1+MHksflsTHSgP/dT6NkFp/pggs/ca6Cg0k02sr3dk1o0gT7Qd8EB/9VF5FRcoH9joOdVdYmqzlbV2RUVcZgJZ3yvtrkDVZjkkxY6uEDftD2Bgb5zk9sq7shLoWJm4s47mJxCOPU62PQirLN105NZNIFeA0TuNlsFbO1/kIi8H7gVWKCqTbEpz6SacEt4kg9miYZNLitga3MHXb0JmmTz9I8BgXn7jC/wzlGfgLHT4Ynv22SjJBZNoC8HZojINBHJBi4GlkYeICKTgQeAT6jqutiXaVLFhiZ38XFaeYHHlewxZWw+qnvGx8dVw1p4409w7GehdNLQxydKRhac9h1oWANv/tnraswIDRnoqtoLXAU8CqwB7lPV1SKyUEQWhg77LlAG/EZEXhcRu7piBrS+oZ387Awqi3K8LmW3KaGRLgnpdnnyBsjKh5O/Fv9zDdes82HC0W4oZY8PlkMww5YZzUGqugxY1u+xxRG3Pwt8NralmVS0oamdqWUFiFfD9AYQHroY95EuW16BNX+HU66DgvL4nmskROD06+EP57k+/rlXeV2RGSbvZ3aYtPJeYzvTKvzT3QJQUZhDXlZG/Ee6/PMHkF8GJ3wpvucZjQPnuc0wnvsJdDZ7XY0ZJgt0kzA9fUFqdnQwrcxfgS4ioZEucZxctP5p9+fk/3CLY/nZ6de7Tar/9SuvKzHDZIFuEmbz9l30BdVXF0TDppUXsL4xToGu6lrnxVUw+4r4nCOWxh8B7/sIvPQbaK3zuhozDBboJmHeCwXmVB8G+vTKAjY27aK7Nxj7N3/7H67//JRr3a5ByeCD34a+bnjmx15XYobBAt0kTDjQD/RhoM+oLKIvqLuHVcZMsA/+eQOUz4QjLonte8fT2APhmE/DK3dC07teV2OiZIFuEua9xnZK8rIYU5DtdSn7OKiyEIDq+rbYvvEb90LjWtfizYhqUJl/zPsvyMx1Qy1NUrBANwnzXmO7L7tbAA4MjbyJaaD3dMLTN8KEo+DQ82L3volSWOlG5Kz+K2x51etqTBQs0E3CrNvWxsHjCr0uY0D52ZlMLM2LbaCvuB2aN7tRIz4adz8sc692Qy2fuN7rSkwULNBNQjS1ddHY1sXMcf4dsndQZWHsAr2r1Y3lnjYPDjwlNu/phdxi+MDX4b1n4N0nva7GDMEC3STEum0uKP0e6Osb2wgGY7CE7IuLYFcTnP690b+X12Z/xm1e/cT1EIzDKCATMxboJiHWbWsF4OAD/B3onT1Btox2f9H2Rnjh167ffOIxsSnOS5k5bmej2jfgrb96XY3ZDwt0kxDrtrVSkpflq0W5+psRGukS/vAZsWdvcptXfNAHm1fEyuEXQeVhbghmb7fX1ZhBWKCbhFi3rZWDxxX5alGu/sLfHt7a2jLyN2l8xy1sdfTlUHFwjCrzgUCG6z7a8R68eqfX1ZhBJNnAWBNzPR1Q+ybs2ADt9W4iTGYOFB0AZTOg8lD3yzwKqsrbda0sOHJCbGqOk6LcLKaW5bN6NIH+2LchM891UaSaGWfC5LnwzP+6SVI5/hyxlM4s0NNR6zZYdb9byrVmBQR7Bj82u9CN1DjsfDj47BH9Em/avovWzl4Om1Ay8poT5LAJJby5ZefIXvzuU24Lt9Ovd2O4U40InPF9uO0MeOm3MO/rXldk+rFATyd1K+G5n8FbfwPtgwPeDyd8ESYdD2UHQdE4CGS5VntrLdS/BZtegrUPw9qHIKcEjvkkHP8FKI6+tb1yi1uG9fCJ/g/0WROKeWhlLc0dPZTkZUX/wr5eePSbUDoFjvtC/Ar02qQ5cMi58K9fwuxP+3Nd9zQWVR+6iMwXkbUiUi0i+2yEKCKHiMiLItIlIv8Z+zLNqOzYCPddDotPgnced7P/vvQyLHwOzvgBHHK226w4twSy86GgDA54H7z/o3Duz+Crq+HTD8NBp8GLN8OvjnJ7T3ZG1zWxsqaZ7IyAr4cshh02oRgYQT/6a39wH4Bn/CB5FuAaqdO+Cz3t8NxPva7E9DNkoItIBrAIOAuYBVwiIrP6HbYduAb4ScwrNCPX2wVP3QiL5rggP+U6+OoqOPOG4V2wCwRgyly46PdwzWtuON7zP4NfHw0r73fLw+7Hyi3NHDq+iOxM/1+DD3cLrd46jM0d2hrcB9yUE2HWgjhV5iMVB8ORl8HLt7g9Uo1vRPMbNgeoVtX1qtoN3Avs9VOrqvWquhzYT2esSai6VXDLB+GZH7m+76uWu+Vb80pH975jpsJHboHPPQklk+AvV8B9n4C2+gEPV1VWbmnmfUnQ3QJQUZTDhJJcXt+8M/oXPf4d6G6Hc3+evFP8h+u070F2Afzja0N+oJvEiSbQJwKbI+7XhB4bNhG5UkRWiMiKhoaGkbyFGUowCM//HJac4kL2kj+7lnVJVWzPM/EYuOJxOP37sO4xWHSc+xbQz4Ymd0E0GfrPw46eMoZXN+6I7uD3noU37oETr0mtYYpDKaxwF0g3Pu/++40vRBPoAzU5RvSRrKpLVHW2qs6uqKgYyVuY/elsgT9f5qZoH3wWfPElOHh+/M6XkQknfcX1xRdPgLsvdOfu6919yPIN2wE4ZsqY+NURY8dMGcPW5k5qm4eYMdrb5VqoY6a69U7SzVGXw6Tj3FDNXdu9rsYQXaDXAJMi7lcBW+NTjhmxxnfg1tNg3aNw1k3w0T+4i5uJUHEwfPYJOOZT7tvBnedCSy0AL7+3nTH5WbvXG08GR092Hz6vbty5/wOf+V9oegfO/ilk5cW/ML8JBFw3U8dOF+rGc9EE+nJghohME5Fs4GJgaXzLMsOy9hHXX76rCS7/Gxx3ZeL7crPy4P/9Ej58i5uotOQUqFnB8g3bOXbqWF/PEO1v1oRicrMCvLK/bpfNy92F4SMvgxmnJ644vxl3mPuW9vrdbnir8dSQga6qvcBVwKPAGuA+VV0tIgtFZCGAiBwgIjXA14Bvi0iNiBTHs3CD6y9/5ia452L3tf/KZ2Dayd7W9P6Pwmcfh8wc9PdnMXvHI8yZNtbbmoYpKyPAEVWl/Pu9poEP6N4FDy6E4okw/8bEFudH866FcYfD0qvdwmTGM1GNI1PVZao6U1Wnq+oPQ48tVtXFodt1qlqlqsWqWhq6PYr502ZIXa3wf5fDU//tFk76zKNQOmno1yXCuMPgc0/ROOYofpq9mAV1N+/Vr54MTp5RzuqtLTS1de375OPfhaZqWLDIjd1Pd5nZcMFi6GyGf3zVRr14yP8Dg82+mt6FW8+Atx+CM38IH17iJgT5SUEZN1XeyN2cRcXq29wF044oR474wMkz3EX756v7tThXPQDLb4HjvwQHzvOgMp864H1u/Zo1S+GV33tdTdqyQE821U/ALadCWx18/AGYe5Uvxz4Hg8qT63bw0sHfgPNuhg3Pw5JTYdtbXpcWlfdNLKE0P4vn3okI9IZ1rluhao5br8Xsbe41MP00ePgbsPU1r6tJSxboyULVjSC5+yI3oefKp2H6qV5XNajVW1tobOvig4dUwNGfgE895NYIv/U0WPUXr8sbUkZAOPGgcp57p8HtYNTV5iZQZebARXe4bgazt0DAXRQvqHRLTSTRN7JUYYGeDLpa3S/IE9e7qeVXPOYugvrYk2/XIwIfCHVdMPk4+PyzbkGw+z8Dj37L9/3qpx9aybaWLl7b2Aj/9yk3NPQjt0HJiObVpYeCMveB11Lr/p/12eTxRLJA97vGd+CW0+Dtf8CZ/w0X/t5NufYxVeWhlVs5ZvIYygojdigqOgA++Xc49nNuka+7LvD1qIjTDh1HdoYgD/0HVD8O5/zU19+KfGPSsW4I6/qn4e9fsYukCWSB7mdvPxQaX94In3gQ5l7ty/7y/tbUtrJuWxsLjhqgJZuZDef8BM7/LWx+GX43Dzb9O/FFRqE4J5Ofly/l6Ma/oSd+zS0Xa6Jz1GUw7xvw+l3wzI+9riZtWKD7UU8nLPs63HsplE1348uTaETF317fQmZAOOfw8YMfdOSlbqhlIAN+fxY8/SN/dcGowpM3cE7zPfyp91RemvZFrytKPqdc53Y2evpGt2m2iTsLdL+pf9u1yl9e4obG+Wl8eRS6evt44LUtnHJwBWMLhrhwOOFIWPi8G0f/9I1wx9mw/b2E1LlfwT545Fp47qf0HvlJ/jdzIXe9vHno15m9ibgRTrPOd0sDWKjHnQW6X/T1wr9+BUvmub09L/sLzP8fN6oiifzjjVoaWru4/ISp0b0gtxg+/Dt3sbF+DfzmBLcbjlcX07ra4N7L4N+L4fgvknneL7hw9mQeXVVHfUunNzUls4xM928bDvUnrncznE1cWKD7Qd1KN5zv8e+4cbwL/5WU64OoKrc+/x4zxxVy8oxhbk12+IXwxRdh+gfdTMwlp7jt7xKp9k337eidR+Hsn7hp/YEAlx0/hT5VbvuXD749JKNwqB/zaTf09i9XuG0OTcxZoHupvQmW/Ze7MNiyBS66Ey6+2+3tmYQeXlXHmtoWPnfygSNbjKukCi75E3zsbrcc6+0fgnsudd1Q8dTbDc//wn2odjbDJ/4Kcz63++lp5QUsOGICd76wgfpWa6WPSEamW5nxjB/A6gfcyK2GdV5XlXIs0L3Q0+H6E391lJtGfswn3R6fh52fFKNYBtLV28eND6/hkAOK+PDRo9xM49Bz4eoV8MFvw4bn4LcnuDHNm1+OSa27qbrlhn87F574Hsw4E77wAhx4yj6HfuX0mfT0KT97zEJoxETgxC+77sS2Ote9+PIt7pqFiQlRj8aIzp49W1esWOHJuT3T2QIrbocXF7l+8oPOcPt7Vh7qdWWjduOyNfzu2fX88Yo5u9dBiYn2JnjhV7Di99DV7HZKOuIS1ydbOMLzdLfDmr/DCzfDtpUw9kCY/yOY+aH9vux/lq1hybPrueuK4zhpuF1KZm8tW+HBL8L6p2DibNd6H/9+r6tKCiLyiqrOHvA5C/QEqH8bXr3TrRnd2QwHngof+E+YepLXlcXEU2vr+fTvl3PZcZP54QWHx+ckXW1uq7Plt0HDGpCA2y1n6sluA+vKQ6Fw3MDfcPp63IJmm150W8ate8QtQ1A+07UYD/9oVFP5O3v6OPuXz9He3cuDXzqR8SVpuKlFLKnCyv+DR65zcy0Ou8ANdUynrfxGwALdCy1b3cSglffD5pcgkAWHnOP2npx4jNfVxczyDdu5/LaXmVpewANfmEtedkb8T7pttVsP5t2noPZ10NCoiZxiKKyEnCLIzHWh3bETmmtAQ1/rCyrdv8PhF8LkuW79kWFYU9vCRYtfpGpMHnd99jjKC5NrFJIvdexw35b+vdh9e5pxBhz7WTjodDdPwezFAj0RejpcH+/Gf0H1P2FL6L+tfCYc9XE44tKRdxH4kKpy/ys1fOvBVUwszeO+z59ARZEH4dbZDFtedUskNK5zuzZ1tbj9PrMLXLiXTnETtKrmuL9HeZ3iuXca+NwfVlBZlMuiS4/m8CpbEz0m2pvg5d/BK3dA2zb3jeuQc+DQ89y3sCQbwhsvow50EZkP/BLIAG5V1R/1e15Cz58N7AI+paqv7u89kzbQVd36I9vfdcMNt62CulVQ9yb0dbuugPFHwiFnux/EFPv6qKr8q7qJ3zxdzQvvNnHctLH85rKj916zJQ28umkHC//4Ck3t3Xx09iSu/MCBTCv39xo7SaOvx327Xf0AvPO4+6aVkQNVs2HyCTD+CNfFNmaaGz2TZkYV6CKSAawDzsBtGL0cuERV34o45mzgalygHwf8UlWP29/7eh7ofT2uVd3TAb0de273dLhWX8d219rbtd3dbq2DHRth5yZ3fFhuidt+a+JRrj938vEpsYtNX1Bp7exh564eaps7qW5oY2XNTp57p5Ha5k7KCrL58ukzuOy4KWQEknNkzmg17+rhJ4+t5c/LN9PdF+SwCcUcf2AZh00oZmp5AZVFOZQX5pCbZd0GI9bT4Rb52vA8bHwBat/Y032Wke1CvWQiFE9wWwIWVkJuqfsdDP/Jynet+8wc98GQmZO0o8lg9IF+AnC9qn4odP86AFW9MeKY3wFPq+o9oftrgVNUtXaw9x1poK98+n5Knv0uAYIEVBGCBAj/Hb6tBDTo/g497m674zJCf6LRQybNFNEkY6iVSmqlkjqpZIuM493AVOopI/L/4F63NfJxHeTxgY/v/+zgr9FBHh/6GPbznrt6+vappzg3k5NmlHP6oeM4+/DxFlQh9S2d/PW1LTz+1jZWbmmmq3fvn62MgJCTGSA3K4OczMDuD0ARECT0N4gIAtD/vtktR7uYEtzElL5NTA1uYnywjopgIxXaxBjdQYDoupC7yaSHLPrIQEUikkOITJU99wWV8GODGfhfSwd4vHb6RRx/2fei+4/uf5b9BHo031cmApELWdTgWuFDHTMR2CvQReRK4EqAyZMnR3HqfWUVjKE+f6b7HysZBGVPVCsBVFx8I6F4F4Fw1If/QSRAj2TTIzn0BHLd7UDodiCHjkAhuzJKaM8ooSuQBxEt0MhfsekC03c/HvnfyYDH7/X4Xv/GAx8zovfd6/iBf8CieZ+CnExK87Iozc+ioiiH6RWFjC/JHdmEoRRXWZzL5+dN5/PzptPbF2R9Yzs1O3bR0NpFQ2sXHT19dPUE6eoN0tnTR1/QRYLiPjzd33vfR/duBJiwQpQyNnAUG/o9E9BeCvtayAu2kdfXRn6wjfxgK9nBLjK1m0ztIUt7Qrfd/Qztg3CChBqIe+4HIyJ9z/2BDB7zAz+eGafJg9EE+kC/wf2rjOYYVHUJsARcCz2Kc+/jkGNPg2NPG8lLjYm7zIwAM8cVMXNckdelmDQUzZitGiByub8qYOsIjjHGGBNH0QT6cmCGiEwTkWzgYmBpv2OWApeLczzQvL/+c2OMMbE3ZJeLqvaKyFXAo7hhi7er6moRWRh6fjGwDDfCpRo3bNG2djHGmASLahCnqi7DhXbkY4sjbivwpdiWZowxZjhstUVjjEkRFujGGJMiLNCNMSZFWKAbY0yK8Gy1RRFpADaO8OXlQGMMy4m3ZKo3mWqF5Ko3mWqF5Ko3mWqF0dU7RVUHXLrVs0AfDRFZMdhaBn6UTPUmU62QXPUmU62QXPUmU60Qv3qty8UYY1KEBboxxqSIZA30JV4XMEzJVG8y1QrJVW8y1QrJVW8y1Qpxqjcp+9CNMcbsK1lb6MYYY/qxQDfGmBSR9IEuIv8pIioi5V7Xsj8icpOIvC0ib4rIX0Wk1Oua+hOR+SKyVkSqReRar+sZjIhMEpGnRGSNiKwWkS97XVM0RCRDRF4TkX94Xcv+iEipiNwf+nldE9qG0rdE5Kuhn4NVInKPiOR6XVMkEbldROpFZFXEY2NF5HEReSf095hYnCupA11EJuE2r97kdS1ReBx4n6q+H7fp9nUe17OX0Gbgi4CzgFnAJSIyy9uqBtUL/IeqHgocD3zJx7VG+jKwxusiovBL4BFVPQQ4Ah/XLCITgWuA2ar6PtwS3xd7W9U+7gDm93vsWuCfqjoD+Gfo/qgldaADPwf+i8E27vMRVX1MVXtDd1/C7erkJ3OAalVdr6rdwL3AAo9rGpCq1qrqq6HbrbjAmehtVfsnIlXAOcCtXteyPyJSDHwAuA1AVbtVdaenRQ0tE8gTkUwgH5/tlqaqzwLb+z28ALgzdPtO4PxYnCtpA11EzgO2qOobXtcyAp8BHva6iH4G2+jb10RkKnAU8G+PSxnKL3CNj6DHdQzlQKAB+H2oe+hWESnwuqjBqOoW4Ce4b+m1uN3SHvO2qqiMC+/qFvq7MhZv6utAF5EnQv1i/f8sAL4FfNfrGiMNUW/4mG/hugzu9q7SAUW10befiEgh8BfgK6ra4nU9gxGRc4F6VX3F61qikAkcDfxWVY8C2olRd0A8hPqeFwDTgAlAgYh83NuqvBPVjkVeUdXTB3pcRA7H/QO+ISLgui9eFZE5qlqXwBL3Mli9YSLySeBc4DT13wSApNroW0SycGF+t6o+4HU9QzgROE9EzgZygWIRuUtV/Rg8NUCNqoa/8dyPjwMdOB14T1UbAETkAWAucJenVQ1tm4iMV9VaERkP1MfiTX3dQh+Mqq5U1UpVnaqqU3E/hEd7GeZDEZH5wDeA81R1l9f1DCCazcB9Qdyn+G3AGlX9mdf1DEVVr1PVqtDP6sXAkz4Nc0K/Q5tF5ODQQ6cBb3lY0lA2AceLSH7o5+I0fHwRN8JS4JOh258E/haLN/V1Cz3F3AzkAI+HvlW8pKoLvS1pj8E2A/e4rMGcCHwCWCkir4ce+2Zo71szelcDd4c+2Nfj403fVfXfInI/8CquK/M1fLYMgIjcA5wClItIDfA94EfAfSJyBe5D6aKYnMt/3/yNMcaMRFJ2uRhjjNmXBboxxqQIC3RjjEkRFujGGJMiLNCNMSZFWKAbY0yKsEA3xpgU8f8BAbgztM0Orv0AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-4, 10, 1000)\n",
"plt.plot(x, norm.pdf(x, p_mu, p_sigma))\n",
"plt.plot(x, w1*norm.pdf(x, q_mu1, q_sigma1) + w2*norm.pdf(x, q_mu2, q_sigma2))\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"KL(q, p): 15.2689387314446\n"
]
}
],
"source": [
"# calc KL\n",
"z = bernoulli.rvs(p=w1, size=N_MONTE)\n",
"x_q = z * norm.rvs(loc=q_mu1, scale=q_sigma1, size=N_MONTE) + (1-z)* norm.rvs(loc=q_mu2, scale=q_sigma2, size=N_MONTE)\n",
"KL_div = np.sum((np.log(w1*norm.pdf(x_q, q_mu1, q_sigma1) + w2* norm.pdf(x_q, q_mu2, q_sigma2)) - np.log(norm.pdf(x_q, p_mu, p_sigma)))) / N_MONTE\n",
"print(f\"KL(q, p): {KL_div}\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATa0lEQVR4nO3dfayc5Znf8e9vHUrYZFGIMMixnZquvO0CbcxiuW6RqnRJixtWMfkjkiM1oG4kp4i0SZWqhazUzaqyRNW8dFELkhNYTJcGWXkRVgLbODSrVSRecmAJxjg0VqBwsIvPbpWGtBJbO1f/mNvV1IzPmXOOPePj+/uRRvPMNfczc404/Pyc+9zPM6kqJEl9+KVpNyBJmhxDX5I6YuhLUkcMfUnqiKEvSR1527QbWMill15aGzZsmHYbkrSiPP30039WVatPrZ/zob9hwwZmZmam3YYkrShJ/tuoutM7ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUkXP+jFxJZ9+G27891riX77zxLHeis23BI/0kb0/yVJIfJjmY5Pda/XNJXkvybLt9cGifO5IcTvJikhuG6tcmOdCeuytJzs7HkiSNMs6R/pvAb1bVz5NcAHw/yaPtuS9V1eeHBye5EtgBXAW8B/hukl+rqhPAPcBO4AngEWAb8CiSpIlY8Ei/Bn7eHl7QbvN9se524KGqerOqXgIOA1uSrAEurqrHa/DFvA8ANy2re0nSooz1h9wkq5I8CxwD9lfVk+2pTyZ5Lsl9SS5ptbXAq0O7z7ba2rZ9an3U++1MMpNkZm5ubvxPI0ma11ihX1UnqmoTsI7BUfvVDKZqfhXYBBwFvtCGj5qnr3nqo95vd1VtrqrNq1e/5XLQkqQlWtSSzar6KfDHwLaqer39Y/AL4MvAljZsFlg/tNs64EirrxtRlyRNyDird1YneVfbvgj4APCjNkd/0oeB59v2PmBHkguTXAFsBJ6qqqPAG0m2tlU7NwMPn7mPIklayDird9YAe5KsYvCPxN6q+laS/5hkE4MpmpeBTwBU1cEke4EXgOPAbW3lDsCtwP3ARQxW7bhyR5ImaMHQr6rngGtG1D82zz67gF0j6jPA1YvsUdISjXvSlfrhZRgkqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1JFxrrIpScD4F3B7+c4bz3InWiqP9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6siCoZ/k7UmeSvLDJAeT/F6rvzvJ/iQ/bveXDO1zR5LDSV5McsNQ/dokB9pzdyXJ2flYkqRRxjnSfxP4zap6H7AJ2JZkK3A78FhVbQQea49JciWwA7gK2AbcnWRVe617gJ3AxnbbduY+iiRpIQuGfg38vD28oN0K2A7safU9wE1tezvwUFW9WVUvAYeBLUnWABdX1eNVVcADQ/tIkiZgrDn9JKuSPAscA/ZX1ZPA5VV1FKDdX9aGrwVeHdp9ttXWtu1T66Peb2eSmSQzc3Nzi/g4kqT5jBX6VXWiqjYB6xgctV89z/BR8/Q1T33U++2uqs1VtXn16tXjtChJGsOiVu9U1U+BP2YwF/96m7Kh3R9rw2aB9UO7rQOOtPq6EXVJ0oSMs3pndZJ3te2LgA8APwL2Abe0YbcAD7ftfcCOJBcmuYLBH2yfalNAbyTZ2lbt3Dy0jyRpAsa5tPIaYE9bgfNLwN6q+laSx4G9ST4OvAJ8BKCqDibZC7wAHAduq6oT7bVuBe4HLgIebTdJ0oQsGPpV9RxwzYj6nwPXn2afXcCuEfUZYL6/B0iSziLPyJWkjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSPjnJwl6Ryz4fZvT7sFrVAe6UtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpIwuGfpL1Sb6X5FCSg0k+1eqfS/Jakmfb7YND+9yR5HCSF5PcMFS/NsmB9txdSXJ2PpYkaZRxrrJ5HPhMVT2T5FeAp5Psb899qao+Pzw4yZXADuAq4D3Ad5P8WlWdAO4BdgJPAI8A24BHz8xHkSQtZMEj/ao6WlXPtO03gEPA2nl22Q48VFVvVtVLwGFgS5I1wMVV9XhVFfAAcNNyP4AkaXyLmtNPsgG4BniylT6Z5Lkk9yW5pNXWAq8O7Tbbamvb9qn1Ue+zM8lMkpm5ubnFtChJmsfYoZ/kncDXgU9X1c8YTNX8KrAJOAp84eTQEbvXPPW3Fqt2V9Xmqtq8evXqcVuUJC1grNBPcgGDwH+wqr4BUFWvV9WJqvoF8GVgSxs+C6wf2n0dcKTV142oS5ImZJzVOwHuBQ5V1ReH6muGhn0YeL5t7wN2JLkwyRXARuCpqjoKvJFka3vNm4GHz9DnkCSNYZzVO9cBHwMOJHm21T4LfDTJJgZTNC8DnwCoqoNJ9gIvMFj5c1tbuQNwK3A/cBGDVTuu3JGkCVow9Kvq+4yej39knn12AbtG1GeAqxfToCTpzPGMXEnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6sg4l2GQpEXZcPu3xxr38p03nuVOdCqP9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1ZMHQT7I+yfeSHEpyMMmnWv3dSfYn+XG7v2RonzuSHE7yYpIbhurXJjnQnrsryagvXJcknSXjHOkfBz5TVb8ObAVuS3IlcDvwWFVtBB5rj2nP7QCuArYBdydZ1V7rHmAnsLHdtp3BzyJJWsCCoV9VR6vqmbb9BnAIWAtsB/a0YXuAm9r2duChqnqzql4CDgNbkqwBLq6qx6uqgAeG9pEkTcCi5vSTbACuAZ4ELq+qozD4hwG4rA1bC7w6tNtsq61t26fWR73PziQzSWbm5uYW06IkaR5jh36SdwJfBz5dVT+bb+iIWs1Tf2uxandVba6qzatXrx63RUnSAsYK/SQXMAj8B6vqG638epuyod0fa/VZYP3Q7uuAI62+bkRdkjQh46zeCXAvcKiqvjj01D7glrZ9C/DwUH1HkguTXMHgD7ZPtSmgN5Jsba9589A+kqQJGOebs64DPgYcSPJsq30WuBPYm+TjwCvARwCq6mCSvcALDFb+3FZVJ9p+twL3AxcBj7abJGlCFgz9qvo+o+fjAa4/zT67gF0j6jPA1YtpUJJ05nhGriR1xNCXpI4Y+pLUEUNfkjpi6EtSR8ZZsilpAjbc/u1pt6AOeKQvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjqyYOgnuS/JsSTPD9U+l+S1JM+22weHnrsjyeEkLya5Yah+bZID7bm7kpzuy9YlSWfJOEf69wPbRtS/VFWb2u0RgCRXAjuAq9o+dydZ1cbfA+wENrbbqNeUJJ1FC36JSlX9SZINY77eduChqnoTeCnJYWBLkpeBi6vqcYAkDwA3AY8upWlNzrhf7PHynTee5U4knQnLmdP/ZJLn2vTPJa22Fnh1aMxsq61t26fWR0qyM8lMkpm5ublltChJGrbUr0u8B/jXQLX7LwC/DYyap6956iNV1W5gN8DmzZtPO07nDn8jkFaGJR3pV9XrVXWiqn4BfBnY0p6aBdYPDV0HHGn1dSPqkqQJWlLoJ1kz9PDDwMmVPfuAHUkuTHIFgz/YPlVVR4E3kmxtq3ZuBh5eRt+SpCVYcHonyVeB9wOXJpkFfhd4f5JNDKZoXgY+AVBVB5PsBV4AjgO3VdWJ9lK3MlgJdBGDP+D6R1xJmrBxVu98dET53nnG7wJ2jajPAFcvqjtJ0hnlGbmS1BFDX5I6YuhLUkeWuk5fK9y46+olnV880pekjnikL2lqPJN78jzSl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOuK1dzRRi7m6p9dbkc48j/QlqSOGviR1ZMHQT3JfkmNJnh+qvTvJ/iQ/bveXDD13R5LDSV5McsNQ/dokB9pzdyXJmf84kqT5jHOkfz+w7ZTa7cBjVbUReKw9JsmVwA7gqrbP3UlWtX3uAXYCG9vt1NeUJJ1lC4Z+Vf0J8D9OKW8H9rTtPcBNQ/WHqurNqnoJOAxsSbIGuLiqHq+qAh4Y2keSNCFLXb1zeVUdBaiqo0kua/W1wBND42Zb7f+07VPrIyXZyeC3At773vcusUXp3OD3EetccqaXbI6ap6956iNV1W5gN8DmzZtPO05vZcBIms9SV++83qZsaPfHWn0WWD80bh1wpNXXjahLkiZoqaG/D7ilbd8CPDxU35HkwiRXMPiD7VNtKuiNJFvbqp2bh/aRJE3IgtM7Sb4KvB+4NMks8LvAncDeJB8HXgE+AlBVB5PsBV4AjgO3VdWJ9lK3MlgJdBHwaLtJkiZowdCvqo+e5qnrTzN+F7BrRH0GuHpR3UmSzijPyJWkjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kd8ZuzdM4a95ISfsOWND6P9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOekSvpnDfu2dngGdoLMfRXiMX80EvS6Ti9I0kdMfQlqSPLCv0kLyc5kOTZJDOt9u4k+5P8uN1fMjT+jiSHk7yY5IblNi9JWpwzcaT/d6tqU1Vtbo9vBx6rqo3AY+0xSa4EdgBXAduAu5OsOgPvL0ka09mY3tkO7Gnbe4CbhuoPVdWbVfUScBjYchbeX5J0GssN/QK+k+TpJDtb7fKqOgrQ7i9r9bXAq0P7zrbaWyTZmWQmyczc3NwyW5QknbTcJZvXVdWRJJcB+5P8aJ6xGVGrUQOrajewG2Dz5s0jx0jT5jJarUTLOtKvqiPt/hjwTQbTNa8nWQPQ7o+14bPA+qHd1wFHlvP+kqTFWXLoJ3lHkl85uQ38feB5YB9wSxt2C/Bw294H7EhyYZIrgI3AU0t9f0nS4i1neudy4JtJTr7Of6qqP0ryA2Bvko8DrwAfAaiqg0n2Ai8Ax4HbqurEsrqX8AvUpcVYcuhX1U+A942o/zlw/Wn22QXsWup7SpKWxzNyJakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkf85qwp81R+SZNk6Es6r3iy3vwMfekU/val85mhr24Y5pJ/yJWkrhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdccnmWeDSQEnnKkNfUpd6PXPX6R1J6oihL0kdmfj0TpJtwO8Dq4CvVNWdk+5hqZyrl7TSTTT0k6wC/gPw94BZ4AdJ9lXVC5PsQ5LGdb7N/U/6SH8LcLiqfgKQ5CFgOzDV0PcIXtJyrZR/HCYd+muBV4cezwJ/89RBSXYCO9vDnyd5cczXvxT4s2V1OD0rtXf7nryV2vtK7RvOYO/5N2fiVcbyl0cVJx36GVGrtxSqdgO7F/3iyUxVbV5KY9O2Unu378lbqb2v1L5hZfd+qkmv3pkF1g89XgccmXAPktStSYf+D4CNSa5I8peAHcC+CfcgSd2a6PROVR1P8kngPzNYsnlfVR08g2+x6Cmhc8hK7d2+J2+l9r5S+4aV3fv/J1VvmVKXJJ2nPCNXkjpi6EtSR87b0E/yz5NUkkun3cs4kvzbJD9K8lySbyZ517R7mk+SbUleTHI4ye3T7mdcSdYn+V6SQ0kOJvnUtHtajCSrkvxpkm9Nu5fFSPKuJF9rP+OHkvytafc0jiT/rP2cPJ/kq0nePu2eluu8DP0k6xlc6uGVafeyCPuBq6vqbwD/Fbhjyv2c1tDlNP4BcCXw0SRXTrersR0HPlNVvw5sBW5bQb0DfAo4NO0mluD3gT+qqr8GvI8V8BmSrAX+KbC5qq5msPhkx3S7Wr7zMvSBLwH/ghEnfp2rquo7VXW8PXyCwTkM56r/dzmNqvoL4OTlNM55VXW0qp5p228wCJ+10+1qPEnWATcCX5l2L4uR5GLg7wD3AlTVX1TVT6fa1PjeBlyU5G3AL3MenFd03oV+kg8Br1XVD6fdyzL8NvDotJuYx6jLaayI4ByWZANwDfDklFsZ179jcDDziyn3sVh/BZgD/qBNTX0lyTum3dRCquo14PMMZgyOAv+zqr4z3a6Wb0WGfpLvtjm2U2/bgd8B/tW0exxlgb5PjvkdBlMQD06v0wWNdTmNc1mSdwJfBz5dVT+bdj8LSfJbwLGqenravSzB24DfAO6pqmuA/wWc838HSnIJg99grwDeA7wjyT+cblfLtyK/LrGqPjCqnuSvM/gP9MMkMJgieSbJlqr67xNscaTT9X1SkluA3wKur3P7BIoVfTmNJBcwCPwHq+ob0+5nTNcBH0ryQeDtwMVJ/rCqVkIIzQKzVXXyN6qvsQJCH/gA8FJVzQEk+Qbwt4E/nGpXy7Qij/RPp6oOVNVlVbWhqjYw+GH7jXMh8BfSvlzmXwIfqqr/Pe1+FrBiL6eRwdHAvcChqvritPsZV1XdUVXr2s/1DuC/rJDAp/3/92qSv9pK1zPly6mP6RVga5Jfbj8317MC/gC9kBV5pH+e+vfAhcD+9lvKE1X1j6fb0mgTuJzG2XQd8DHgQJJnW+2zVfXI9Frqwj8BHmwHCT8B/tGU+1lQVT2Z5GvAMwymXP+U8+ByDF6GQZI6cl5N70iS5mfoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI78Xxboi03qOrrwAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(x_q, bins=30)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# setting 2\n",
"p_mu = 2\n",
"p_sigma = 2"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8ZklEQVR4nO3dd3yUVfb48c9JDwlJKGmQUAJJKKGHIiiKoIKN1bVhr6xrXcvuuruu+123+dNd197dta6ua0VFUREbRTqETmhJCClAAqmkzP398SQQY0Imycw8M5Pzfr3ymswzTzmB5Myde+9zrhhjUEop5b8C7A5AKaWUe2miV0opP6eJXiml/JwmeqWU8nOa6JVSys8F2R1AS3r37m0GDBhgdxhKKeUzVq1atd8YE9vSa16Z6AcMGMDKlSvtDkMppXyGiOxp7TXtulFKKT+niV4ppfycJnqllPJzmuiVUsrPaaJXSik/p4leKaX8nCZ6pZTyc145j14p1Qn1tbDlY9i/HeKHQdpMCAi0OyplI030SvmTgzvhjTlQvOXYtqQJcPGr0D3BvriUrbTrRil/UVYIL50N5YVw8WvwuwI471ko3AivzIaqUrsjVDbRRK+UPzAG3psLVSVw5TwYeg4Eh8OoS+DSN61unPl32x2lsokmeqX8wYZ3YOdXcNr9kDjyh68NnAon/wqy/gfZC20JT9lLE71Svq7uCHz2e0gcDZnXtrzPiXdATH9rP4fDo+Ep+2miV8rXrXsDyvJhxv+1PrsmKBRO/T0UbYRtn3o0PGU/TfRK+TJHPSx+FPqMgZRTjr/v8PMguh8secwjoSnvoYleKV+W/YU1pXLybSBy/H0Dg+CEmyBnKexd5Zn4lFfQRK+UL1v9CnTrDUPOdm7/0ZdCUDised29cSmvooleKV9VXmT1t4+6BIJCnDsmLNqaernhbaitdm98ymtoolfKV214Fxx1MOaK9h03+lKoPgRbP3ZPXMrraKJXyldt+gDihkHckPYdN/BkiEyAje+7JSzlfTTRK+WLygqtQdVhs9t/bEAADDnLGsitqXR9bMrraKJXyhdt+RAwMPTcjh0/9ByorYQdX7o0LOWdNNEr5Ys2zYNeqRA3tGPHDzgRwmJg84cuDUt5J030Svma6sOwZ7HV/dLW3PnWBAZD+pmw7ROor3NtfMrraKJXytfs/taabTN4RufOk3a6NftGb57ye5rolfI12QshOAKSJ3buPCmngATADq1o6e+cSvQiMlNEtopItojc08Lrl4nI+oavJSIyqslru0UkS0TWishKVwavVJe0YyEMPMn5m6RaE94D+o7T0sVdQJuJXkQCgSeBWcAwYI6IDGu22y7gZGPMSOBPwHPNXp9mjBltjMl0QcxKdV0HdkDJbhg03TXnG3Qq5K+GyoOuOZ/ySs606CcA2caYncaYGuBN4AeTd40xS4wxJQ1PlwFJrg1TKQUcmw452FWJfjoYB+z62jXnU17JmUTfF8ht8jyvYVtrrgM+afLcAJ+JyCoRmdvaQSIyV0RWisjK4uJiJ8JSqgvaschaQKRnimvO13cchEZZq1MpvxXkxD4tzd8yLe4oMg0r0Z/YZPMUY0y+iMQBn4vIFmPMNz86oTHP0dDlk5mZ2eL5lerSHA5rWuXQczo+rbK5wCDoNwn2LHXN+ZRXcqZFnwckN3meBOQ330lERgIvALONMQcatxtj8hsei4D3sLqClFLtVbwZqkuh/xTXnrf/ZNi/Fcr1k7S/cibRrwBSRWSgiIQAlwDzmu4gIv2Ad4ErjDHbmmyPEJHujd8DpwMbXBW8Ul3KniXWY//Jrj1v4xtHzhLXnld5jTa7bowxdSJyC7AACAT+ZYzZKCI3Nrz+DHAf0At4SqyPlHUNM2zigfcatgUB/zHG6IKVSnXEnsUQlQQx/Vx73sTRENwNdi/uWJE05fWc6aPHGDMfmN9s2zNNvr8euL6F43YCo5pvV0q1kzFWi37gya7rn28UFAJJ4499YlB+R++MVcoXHNwJ5YWu77ZpNOBEKNwAVSVt76t8jiZ6pXzBnsXWo6sHYhv1nwwYyFnmnvMrW2miV8oX5HwP3XpB71T3nL/PWAgIgrwV7jm/spUmeqV8wd6V0DfT9f3zjUK6QXyGJno/pYleKW9XfQiKt0KSm0tFJWXC3jXgqHfvdZTHaaJXytvtXQ0Yq1yBOyWNh5oy601F+RVN9Ep5u70N1b09kehBu2/8kCZ6pbxd3ironQbhMe69Ts8Uq0a9Jnq/o4leKW9mzLGBWHcTsVr1ebo+kL/RRK+UNyvNgYpiSHJzt02jpPFQvMUaAFZ+QxO9Ut7saP+8hxZn6zsOMA0DwMpfaKJXypvlrYKgMIgf7pnrNQ747tXuG3+iiV4pb5a3wqouGRjsmeuFx1iDsvlrPXM95RGa6JXyVvW1ULDe/TdKNZc4Gvat8+w1lVtpolfKW+3fDnXVkOjhSt99RsOhXKg40OauyjdoolfKWxVkWY8JIzx73cTR1uO+NZ69rnIbTfRKeavCLAgMhV5uqljZmsZPENpP7zc00SvlrQqyIH4YBDq1EJzrhMdAj4Gwb61nr6vcRhO9Ut7ImIZEn2HP9fuMhnwdkPUXmuiV8kZl+6DyACSMtOf6iaPhUA5UHrTn+sqlNNEr5Y0KNliPnh6IbdRntPWYrwOy/kATvVLeqGC99eipO2KbaxyQ1X56v6CJXilvVJBlDYiGRdlz/fAe0GOAzrzxE5rolfJGBVmQYNNAbKPE0dqi9xOa6JXyNkfK4eBO+wZiGyWOssokV5XaG4fqNE30Snmbok2AsW8gtlHj9Qs32huH6jSnEr2IzBSRrSKSLSL3tPD6ZSKyvuFriYiMcvZYpVQzjQOx3pLoG0sxKJ/VZqIXkUDgSWAWMAyYIyLDmu22CzjZGDMS+BPwXDuOVUo1VZBlDYZG9bU3jsh4iIi1SjEon+ZMi34CkG2M2WmMqQHeBGY33cEYs8QYU9LwdBmQ5OyxSrmaw2EwxtgdRsc13hErYm8cIlYc2qL3ec4U0egL5DZ5ngdMPM7+1wGftPdYEZkLzAXo16+fE2EpBcYYVueU8vXWIr7fdZDdByooKjuCMRAWHMCAXhEMSejOyemxTEuPI6ZbiN0hH5+jHgo3Qea1dkdiSRgB3z9r1cb31OInyuWcSfQtNStabC6JyDSsRH9ie481xjxHQ5dPZmamDzfHlCfU1Dn478pcXlq8ix3FFQQIZPSNZmpqLAnRYYgI5dV17D5Qwbfb9/P+2nxCggI4d1Qf5k5NIS2+u90/QssO7IC6Kvv75xsljID6I1Zt/HjtdfVVziT6PCC5yfMkIL/5TiIyEngBmGWMOdCeY5Vqj0837OPPH28mr6SK0ckxPPjTkcwckUBUWMstTofDsH7vIf63Mpd3V+/l3dV5XDy+H3efnkavyFAPR98GbxmIbdR0QFYTvc9yJtGvAFJFZCCwF7gEuLTpDiLSD3gXuMIYs609xyrlrAPlR7jvg418nLWPoYlRvHLtCE5K7Y200ZcdECCMTo5hdHIMd5+ezmNfbufVpXv4fFMBfz1vBKcPT/DQT+CEgiwICIbeaXZHYumVatXEL8wCLrY7GtVBbSZ6Y0ydiNwCLAACgX8ZYzaKyI0Nrz8D3Af0Ap5q+KOrM8Zktnasm34W5cfW55Xys1dXcaC8hl+ekc7cqSkEB7b/NpAeESH84ZzhXDK+H3f8dy1zX13F9ScO5J5ZQwjqwPlcriAL4oZAkJeMJQQGQdxQHZD1cU6taGCMmQ/Mb7btmSbfXw9c7+yxSrXHB2v38su31xMbGcq7N00mo290p8+ZntCd92+ewl/nb+aF73axraicx+eMITrc5gHHgixIPc3eGJpLGAFb51s18u2eCaQ6xAuaMEq17tWlu7n9zbWMSY5h3i1TXJLkG4UEBfB/5w7nb+ePYOmO/Vz6/DIOVtS47PztVlYIFUXe0z/fKGGkVRu/rMDuSFQHaaJXXuuZr3fw+w82MmNoPC9fO8FtA6dzJvTjuSszyS4q55LnllJUVu2W67Sp0KbFwNvSWFxNu298liZ65ZVeXbaHBz7Zwjmj+vD05WMJCw506/Wmpcfx76vHk1dSxRUvLOdQVa1br9eixkRqVw361jTG0zgjSPkcTfTK63y0Pp/7PtjA9CFxPHzRqA4NunbE5MG9ee6KTHbuL+eGl1dSXVvvkeseVZAF0f2s8gfeJCzaqk1fuMHuSFQHaaJXXmXZzgPc8d+1ZPbvwROXjvVYkm90Ympv/nHRaJbvPsjtb67B4fDgvXsFG7yv26aRlkLwaZroldfIK6nkptdX069nN164cjzhIe7trmnNuaP68Puzh7FgYyGPfLGt7QNcoaYSDmz33kSfMNK6a7emwu5IVAdooldeoaqmnrmvrKK23sHzV2YS3c3eaY7XThnARZlJPPZlNp9u8MBsk6LNYBz2ryrVmoQRgLHq8Cifo4le2c4Ywz3vrmdzwWEeu2QMKbGRdoeEiHD/7AxGJcdw11tryS4qc+8Fva30QXNHZ97ogKwv0kSvbPe/VXl8sDafO2ekMW1InN3hHBUWHMizl48jPCSQW/6zxr2DswVZEBoFMf3dd43OiE62BmW1n94naaJXttpRXM4fPtjICSm9uGnaYLvD+ZGE6DD+fuEothSU8cAnW9x3ocKGgVhvvfNUBOJH6MwbH6WJXtnmSF09t/5nDWHBAfzz4tEEBnhnkjslPY5rpwzkpSW7Wbi50PUXcDi8e8ZNo4QMq4/e4bA7EtVOmuiVbf7x2TY27TvMQxeMIiE6zO5wjuvXs9IZmhjFL99e7/o7Z0t2QW2FNYXRmyWMsOIs2WV3JKqdNNErW6zJKeGFb3cyZ0I/ZgyLtzucNoUGBfL4nNGUH6njDx+4uACrtw/ENorXAVlfpYleeVx1bT2/fHs9CVFh/PbMIXaH47TBcd35xYxUPtlQwPysfa47cUEWBARBrJf/W8QOAQm0upmUT9FErzzusYXbyS4q528/HUn3VlaF8lZzT0phRN9o7vtgAyWuqnRZsAF6p0Owd3dfERwGsek6IOuDNNErj8rKO8Sz3+zkwnFJnJwWa3c47RYUGMCDF4yktLKW+z9y0c1DBVne323TSEsh+CRN9Mpj6h2G37y3nl4RIdx7lu+uPzo0MYqbpw3mvTV7+WprUedOVrEfyvK9947Y5hIy4PBeqDxodySqHTTRK4/5z/d72LD3MPeePcz2EgeddfO0waTERvB/8zZ27kaqAi+tQd+axgFZ7b7xKZrolUfsLz/CQwu2MnlQL84ZmWh3OJ0WEhTAH88dzu4DlTz/zc6On6gxYcb7SKJvfEPSAVmfooleecTf5m+hqrae+2dnIN5692c7nZQay1kjEnliUTa5Bys7dpKCLIjqCxG9XBucu0TGQWS89tP7GE30yu1W7D7IO6vzuOGkFAbH2V+wzJXuPXsogQHCHz/s4MBsQZb33yjVXHzGsWUPlU/QRK/cqt5h+P37G+gbE84tp3pfLZvOSowO5/bpqXyxubD95RFqq6F4q+/0zzdKGGHFXWfjQuqqXTTRK7d6a2UuWwrK+N1ZQ+kWEmR3OG5xzZSBDIqN4M8fb6amrh11YIo3g6n3zURfXwP7PbQoi+o0TfTKbcqqa/nHZ1sZP6AHszIS7A7HbUKCArj3rGHs2l/B69/vcf7AxgFNX0v0OvPG52iiV27z9Fc72F9ew71nDfObAdjWnJIey0mpvXnki+2UVjrZpVGQBSGR0GOge4NztV6DIShMB2R9iCZ65RZ5JZW88N0uzh/Tl1HJMXaH43Yiwu/OGkpZdS2Pf5nt3EEFWRA/HAJ87M8wMAjihmqi9yE+9humfMX/+3QrAQJ3n5FudygeMyQhiovHJ/PK0t3s2t/GItoOh2+VPmguPsPqujHG7kiUE5xK9CIyU0S2iki2iNzTwutDRGSpiBwRkbubvbZbRLJEZK2IrHRV4Mp7rdpTwofr8pl7Ugp9YsLtDsej7jgtjZDAAB74ZPPxdyzdAzVlvpvoE0ZA5QEo88DC6arT2kz0IhIIPAnMAoYBc0SkeaGSg8BtwN9bOc00Y8xoY0xmZ4JV3s8Yw18+3kRs91B+dvIgu8PxuLjuYdw0bTALNhby/c4Dre9Y6KMDsY2O3iGr3Te+wJkW/QQg2xiz0xhTA7wJzG66gzGmyBizAqh1Q4zKhyzaWsTqnFLumJFGRKh/Tqdsy3UnDiQ+KpQHF2zFtNa1UZAFEgBxPlrcLX649ag3TvkEZxJ9XyC3yfO8hm3OMsBnIrJKROa2tpOIzBWRlSKysri4uB2nV97C4TD8fcE2+vXsxoWZSXaHY5uw4EBun57Gqj0lfLG5leqWBVnQKxWCfbRrKywaYvppzRsf4Uyib2leXHtGYKYYY8Zidf3cLCJTW9rJGPOcMSbTGJMZG+t7dcoVfLqxgE37DvOLGakEB3btcf6LMpNI6R3BQwu2UO9o4c/FlwdiGyWM1K4bH+HMX2MekNzkeRKQ7+wFjDH5DY9FwHtYXUHKz9Q7DA9/vo3BcZHMHt2eD3z+KSgwgLtOT2dbYTnvr9n7wxcrD8KhXN+pQd+a+Aw4uANqOljQTXmMM4l+BZAqIgNFJAS4BJjnzMlFJEJEujd+D5wO6Gc9PzRv3V6yi8q587Q0AgP8++YoZ83KSGBE32ge/nwbR+qa1Kw/WoN+pD2BuUpCBhgHFLUxw0jZrs1Eb4ypA24BFgCbgbeMMRtF5EYRuRFARBJEJA+4E7hXRPJEJAqIB74TkXXAcuBjY8yn7vphlD1q6x388/PtDEuMYuZw/y110F4BAcKvZqazt7SK15flHHuhMdEnjrInMFdp7HrSAVmv59S0CGPMfGB+s23PNPm+AKtLp7nDgI//Nqu2vL0qj5yDlbx4VSYB2pr/gZNSY5kyuBdPLMrmovHJRIYGQcF66N4HInrbHV7nxPSH0Cjtp/cBXXvETHVadW09jy3czujkGE4dEmd3OF7pV2cM4WBFzbGVqPat9/2BWAARa5qlzrzxeproVae8uTyHfYeq+eUZ6X5fuKyjRiXHMHN4Ai9+t4uS0kNWed9EH++fb5QwAgo3WiUdlNfSRK86rLKmjicW7WBSSk8mD/KRpfBscufpaVTU1PHB5wt9swZ9a+IzrFIOpbvtjkQdhyZ61WGvLN3D/vIj3H26tubbkhbfnXNG9mFn1hJrg6/PuGnUOEVUu2+8miZ61SFl1bU88/UOTkmPJXNAT7vD8Qm3z0glzeyiOjACegywOxzXiBtmlXLQRUi8miZ61SEvfreL0spa7jqt65Qh7qxBsZGcGLmP9XX9KCw7Ync4rhEcbpVy0Jk3Xk0TvWq3kooaXvx2F2cMj2dEUrTd4fgORz39aney0dGfpxY5uTiJL0jI0K4bL6eJXrXbc9/upLymjju1Nd8+B3YQUFdFZP+xvLE8l72lVXZH5BrxGXAoB6pK7Y5EtUITvWqX4rIjvLR4N+eO6kN6Qne7w/EtBesBmHrydAyGJ5xdctDbNQ4saz+919JEr9rlqa+yqal3cPv0VLtD8T371kFgCPEpo7hkfD/+tzKXnAN+UBBMZ954PU30ymn5DTVbfjq2LymxkXaH43sKsqxFtQODuXnaYAIChMe+3G53VJ0XGQ/demvNGy+miV457YlF2RgMt56qrfl2M8bqumm4USohOozLJ/bn3dV57Cwutzm4ThKxfi6deeO1NNErp+QcqOStFbnMmdCP5J7d7A7H9xzOtxbTTjhW4+/npwwiNCiQRxf6Qas+caRVrriuxu5IVAs00SunPLJwG4EBws3TBtsdim9qGIhtWuMmtnsoV07uz7x1+WwvLLMpMBfpMwbqa6Bok92RqBZoolfHV7KHoi+fInX9P3i2/yLiC7+DOj+52ceT9q0D5Nii2g1+NnUQ3YIDeeQLH2/VJ462HvPX2BqGaplT9ehVF3RgB3x2L2ydTxxwfWAgwXvr4fVnISIOTv4VZF4HAdpWcEr+GohNh9AfTkntGRHCtScO5PEvs7k5/zDD+kTZFGAn9RgAYTENif4am4NRzelfqfqxDe/AMyfBrm8pHnsHJx95mMcmL4bf5sOlb0HvNJh/N7x2vrX+qTo+Y2Dvaqt7owXXn5RCVFgQD3++zcOBuZCI9fNpi94raaJXP7TiBXj7WmsWxS3L+U3JWZSEJnH91MEQEgFpZ8DVH8E5j8GeJfDvWXB4n91Re7fD+VBR1Gqijw4PZu7UFL7YXMja3FLPxuZKfUZbA7K11XZHoprRRK+OWfMafHwXpM2Eq+axpjScLzYXMXdqCtHhwcf2E4FxV8Hlb8OhPHj9Qqg+bF/c3q6xldtnbKu7XD1lID26BfOPz7Z6KCg36DMGHLVQtNHuSFQzmuiVJed7+PAXkHIKXPQKBIXy8Ofb6BkRwtVTBrZ8zMCpcNHL1kyL/10NjnoPBuxD8teABB67g7QFkaFB/PyUQXy7fT/Ld/lod1jjJxbtvvE6mugVVOyH/14OMclw4UsQFMqynQf4dvt+bjplkLWgdWsGz4Cz/g47FsJ3//RYyD4lf7VVtz04/Li7XTFpALHdQ/n7Z1sxxngoOBeKToZuvSB/rd2RqGY00Xd1xsBHv4DqUrj4NQjvgTGGf3y2lbjuoVw+qX/b5xh3DWT8FBb9FXKXuzti32KM1cLt23L/fFPhIYHcMm0wy3cdZHH2AQ8E52Ii1jRLTfReRxN9V5f1Nmz+EKb99ugc72+272fF7hJuPXUwYcGBbZ9DBM7+J0T1hQ9u0bsjmyrdA1UlrQ7ENnfJhGT6RIf5bqu+zxirK6/WT0ow+wlN9F1Z5UH45FeQNB4m3wZwtDXfNyaci8YnO3+usGg46x+wfyssedRNAfugowOxziX60KBAbpueytrcUr7cUuTGwNykzxhr8fNCHZD1Jprou7Kv/5/VZXP2PyHAarkv2FjI+rxD/GJGKqFBTrTmm0o7HYb9BL5+yLrhSlnz5wNDIG542/s2+Om4JPr17MY/PtuGw+Fjrfo+o61HHZD1Kprou6qizbD8eRh39dGKivUOqzU/KDaC88b07dh5Zz4AgcGw8I+ui9WX5a+xVmAKCnH6kODAAH4xI5VN+w7z6cYCNwbnBlF9ISJWE72XcSrRi8hMEdkqItkick8Lrw8RkaUickRE7m7PscoGxsCnv4HQSJh279HNH6zdy/aicu48LZ2gwA62AaISYfKtsOkDyF3hooB9lMNh1bhxstumqdmj+zIoNoKHP99GvS+16vUOWa/U5l+ziAQCTwKzgGHAHBEZ1my3g8BtwN87cKzytF3fwM5FcPKvIaIXADV1Dh75YjvD+0QxKyOhc+c/4RarHs7n91lvKl3VgWw4chj6tn6jVGsCA4Q7T0snu6icD9fluyE4N+ozFoq3wBEfr8jpR5xptk0Aso0xO40xNcCbwOymOxhjiowxK4Da9h6rPMwYaxpk9z5WUbIGb63MJedgJXefkU5AgHTuGqGRMO03kLMEtn7SyYB9WF7DVNOk8R06fFZGAkMTo3jki23U1jtcGJibJY0H49BWvRdxJtH3BXKbPM9r2OYMp48VkbkislJEVhYXFzt5etVuO76E3GUw9S4IDgOguraexxZuZ/yAHpySFuua64y5EnqmwNcPdN1Wfe5yazZSr46tyBUQINx1Whq7D1Ty7uo8FwfnRknjrEe9p8JrOJPoW2reOfuX6/SxxpjnjDGZxpjM2FgXJRv1Q42t+ehkGHPF0c2vLN1NUdkRfnnGEEQ62ZpvFBgEJ95p9VFv/9w15/Q1eSus1m0nSjlPHxrHqOQYHluYzZE6HykxEd7DqnCat9LuSFQDZ34D84CmE6qTAGc7DTtzrHK17C9g70qYejcEhQJQVl3L01/tYGpaLBMG9nTt9UZdAtH94JsHu16rvvqQNbMpaUKnTiNiter3llbxn+9zXBScBySNt97outr/u5dyJtGvAFJFZKCIhACXAPOcPH9njlWutuQxa/rb6MuObnrxu12UVNbyy9PTXX+9wGA48RfWH/zOr1x/fm+2dxVgILlj/fNNnZTamxNSevH4l9mUVTcfBvNSSeOhcj+U7LI7EoUTid4YUwfcAiwANgNvGWM2isiNInIjgIgkiEgecCdwr4jkiUhUa8e664dRx1GwwZptM+EGKwEDBytqeOHbXczKSGBEUrR7rjvmcmvg95uH3HN+b5W7AhDom9npU4kI98wawsGKGp7/ZmfnY/OExgFo7b7xCk51Hhpj5htj0owxg4wxf2nY9owx5pmG7wuMMUnGmChjTEzD94dbO1bZ4PtnICgcxl51dNPjX26nsqaOu05Pc991g0Jh8i2wZzHkrXLfdbxN3nKIGwphrlkacFRyDGeNTOT5b3dRdNgHFvaIGwohkdanOWU7vTO2K6jYD+vfgtFzoJvVD7/nQAWvLdvDxeP7MTiuexsn6KQxV0BoFCx9wr3X8RYOx7GBWBf65enp1NY7eHShDywkHhBo3T+gM2+8gib6rmDlv6H+CEy88eimBxdsJTgwgDtmdGzqX7uERVkrUm36AEp9aECxow5stwZjkzs3ENvcgN4RXDqxH2+uyGVHcblLz+0WSeOhcAPUVNodSZenid7f1dVY68AOmg6x1oDrmpwSPl6/jxtOSiEuKswzcUy80bo9ftkznrmenRpbsZ2ccdOS26anEhYUwN8X+MCSg0njwVEH+9baHUmXp4ne3216H8oLYNLPAasM8d8+2ULvyFBumJriuTiik2D4ebD6Fau1689yv4ewGOg12OWnbvx/+2RDAatzSlx+fpdqfKPLWWZvHEoTvV8zBpY9Zd2ZOWg6AF9sLmL5roP8Ykbq8ZcIdIcTboaaMivZ+7M9i6H/5E7dKHU8N5yUQu/IEB6Yv8W7FyeJ6AWxQ2DPErsj6fI00fuz3OVWvZFJN0JAAHX1Dh74ZDMpsRFc3J5FRVylzxjof6LVfVPvI/PB2+vwPji4E/pPcdslIkKDuH16Kst3H+TzTYVuu45L9J9stejr6+yOpEvTRO/Plj1l1VoZNQeAt1bmsaO4gntmDiG4o2WIO+uEm+FwnrV8oT/as9h67D/ZrZe5ZEI/BsdF8pf5m727NEL/KdanuMIsuyPp0jTR+6vSXCuZjr0KQiI4XF3Lw59vZfyAHpw2LN6+uNJmQo+BsOxp+2Jwpz2LIaQ7JIx062WCAwO496yh7DlQyctLdrv1Wp3S+Ian3Te20kTvr1Y8bz1OmAvA4wu3c6CihvvOHu66wmUdERBgDQznLffPuyb3LIF+k6yibm52Snoc09JjeXxhNvvLj7j9eh0S1cd6Y9dEbytN9P6opgJWvQRDz4aYZHYUl/Pvxbu5aFyy+0odtMfoS60bqJY9ZXckrlWx31pww83dNk3de/Ywqmrr+cdn2zx2zXbrP8VK9A4fqqnvZzTR+6N1b1hTGCfdBMCfP9pEeHAgd5/hhsJlHRHaHcZeCRvfh0N77Y7GdRpbrQNO9NglB8VGcuUJA/jvihw25R/22HXbpf9kqDoI+31g7r+f0kTvbxwO+P5Za4ZL8kQWbSli0dZibpueSmz3ULujO2bCXMAc62LyB3sWQ3A3SBzt0cvePj2V6PBg7v9oo3dOtzzaT7/Y3ji6ME30/mbHl7B/G0z8OTX1hj99tImU3hFcNXmA3ZH9UI/+MORsqzyDv9wiv+sbSJ4IQSEevWx0t2DuPC2NZTsP8umGAo9e2yk9BlgVTHdroreLJnp/s+wpiIyH4efx8pLd7Nxfwe/PHkZIkBf+V0+6CapLYf2bdkfSeYf3QdEmGDTNlsvPmdCPIQnd+dNHm6g44mVz1kWs7qzd32o/vU288K9fdVjxVtixEMZfz76Keh75YhvT0mOZNiTO7sha1m+S1c2x7GnfTwCNC6sMOtWWywcFBvDnn2SQf6iax7yxuuWgaVBRDEW6HIUdNNH7k++fgcBQGHcNf5y3iTqH4f7ZGXZH1ToRq1W/f5vV5eTLdnwJEbEQN9y2EDIH9OTizGRe/G4XWwvKbIujRSmnWI++/v/sozTR+4vKg7D2DRh5IQtzHXy6sYDbpqeS3LOb3ZEd3/DzIDLBt6daOhxWiz7lFLfVt3HWPbOG0D0siHvfz/KugdmoPlbdmx2L7I6kS9JE7y9WvwJ1VVSNm8t9H2wkNS6SG07yYHXKjgoKgQnXW11ORVvsjqZjijZCRZFt3TZN9YgI4Z5ZQ1ixu4S3V+XZHc4PpUyDnKVQ6wMrZPkZTfT+oL4Olj8PA07i0axQ9pZW8eefZHjnAGxLxl0DQWFW15MvauyOaOyesNmF45IZ178Hf/tkCyUVNXaHc8ygU6Gu2kr2yqN8JBOo49ryIRzOI2/I1bzw7U4uHJfExJRedkflvIjeMPIiWPem1QXla3YsgtihVveEFwgIEP5yXgaHqmr5y/zNdodzzIApEBAMO7X7xtM00fuDZU9jegzg5uWxRIUH85szh9odUftN/DnUVcGqf9sdSfscKbNuBBo83e5IfmBIQhQ3npzC26vyWLS1yO5wLCER1n0GOiDrcZrofV3eKsj9nu96/pR1+eX8aXYGPSM8e8OOS8QPs7o+lj/vW7Xqd3wJ9TWQfqbdkfzIbdNTSY2L5LfvZlFW7SX/poOnQ0EWHM63O5IuRRO9r1v2JPUh3blty3DOHJHAWSMT7Y6o4ybdDGX7rEXEfcXWT61lA5Mn2h3Jj4QGBfLgBSMpPFzNX+d7yUB3+izrcdun9sbRxWii92WluZiN7zMv8DQkLMq758w7Y/AMa53VpU9ayyB6O0e9lbDSzvBIWeKOGNOvB9eflMIby3NYnL3f7nCsKZY9BlhvkMpjNNH7suXPYgw8VHIK988eTu9ILypa1hEBATDxRshfDXkr7I6mbbnLraqMja1UL3XnaWmk9I7g1++sp9zu8ggikDbLuu+gpsLeWLoQTfS+6kgZ9StfYr5jIqMyMjhrhA932TQ1ao61/OHSJ+2OpG3bPrFmkQzyroHY5sKCA3nowpHkl1bxhw+8oARB+iyoP3KsbIRyO6cSvYjMFJGtIpItIve08LqIyGMNr68XkbFNXtstIlkislZE/HBJIXvUrHiFwJoy/hd8Ln/+SYa9q0a5Umiktfzh5nlQmmN3NK0zxqqnP3AqhEXZHU2bxvXvyS3TBvPO6jw+XGfzQGj/yRAaDVvn2xtHF9JmoheRQOBJYBYwDJgjIsOa7TYLSG34mgs0XxB0mjFmtDEms/MhKxz1VHzzOMsd6Vx/yQX08vUum+YmzAXEmoHjrfJXQ+keyDjf7kicdtv0VMb0i+F372Wxt7TKvkACgyF1Bmxb4PvF7HyEMy36CUC2MWanMaYGeBOY3Wyf2cArxrIMiBERP+lL8D7rPn+VHjX7yEm7hpNSY+0Ox/VikmHYubD6ZThSbnc0LdvwrtVtM+QsuyNxWlBgAI9cPJp6h+GO/66l3mHjgPeQs6xqlnqXrEc4k+j7ArlNnuc1bHN2HwN8JiKrRGRuRwNVloLSKoKXPsLegD6ce9H1dofjPpNuspZDXPeG3ZH8mMNhddsMng7hPeyOpl3694rg/tkZLN91kCe+zLYvkLSZ1mpcG96xL4YuxJlE31Lnb/OmwPH2mWKMGYvVvXOziExt8SIic0VkpYisLC4udiKsrqemzsHzLz3PMHYRfPJdhIQE2x2S+ySNt74WPwp1XlSvBawZQYfzYLjvdNs0df7Yvpw3pi+PLNzGN9ts+lsLibCS/ab3rVpNyq2cSfR5QHKT50lA89GcVvcxxjQ+FgHvYXUF/Ygx5jljTKYxJjM21g+7I1zgr/M3c8bB16gKTyRuypV2h+NeInDyPXAoF9a+Znc0P7T+v1YRNi+fVtkaEasWTlpcd25/c419/fUZP4XKA7Dra3uu34U4k+hXAKkiMlBEQoBLgHnN9pkHXNkw+2YScMgYs09EIkSkO4CIRACnAxtcGH+X8cHavWxc+ikTArYSfsodHl+X1BaDp1ut+m/+4T2t+toqyHobhp7rE7NtWtMtJIinLx9Lbb3hptdXc6Su3vNBDJ4BoVHWeIdyqzYTvTGmDrgFWABsBt4yxmwUkRtF5MaG3eYDO4Fs4Hngpobt8cB3IrIOWA58bIzRW+LaaUvBYe55J4vfdf8YExELY/28Nd9IBE65x+omWfOq3dFYNn8IRw7B2CvsjqTTUmIj+fuFI1mXW8qfPtrk+QCCw6wF4jd/CHVHPH/9LsSp+7aNMfOxknnTbc80+d4AN7dw3E5gVCdj7NIOlB9h7iurmBCyi9E1q2Dq/0FwuN1hec6g6ZA0Ab59GMZcDkE2TyVd/Yp1C3//E+2Nw0VmZiTys6kpPPvNTlLjunPV5AGeDWDET2Hdf6w59cPP8+y1uxC9M9aLVdfW87NXV1F4uJrH4z+Cbr1gvB/PtGlJ01b9SptLGB/cBbu/td5wbF4y0JV+NXMIM4bG88cPN/KVp0sap0yDqL6w2ks+sfkp//lt9TPGGH79znpW7inhpVOqiMr/Dk66C0K72x2a5w061boD9ev/B1Wl9sWx4gUICIJRl9oXgxsEBgiPXjKa9IQobv3PGrYVenBh8YBA641zx5dQmtv2/qpDNNF7qUe+2M4Ha/P55elpnLD7KavVk3md3WHZQwRO/wtUlcC3f7cnhurDsOplq3shuvltJL4vIjSIF6/KJDwkkGtfWkHRYQ+u6zr6Mutx7eueu2YXo4neC728ZDePLtzOBeOSuKnPdmve9sm/tgavuqrEkTD6Uvj+WasLxdPWvAo1ZdaNXH6qT0w4L1yVycGKGq7813JKKz0006lHfxg0Dda8ZpV+Vi6nid7LvL9mL3+Yt5EZQ+N5YPYQ5Is/QM9BVpLr6k691+o6+fw+z163vg6WPQP9p0DfsW3v78NGJsXw/JWZ7Cyu4JqXVlDhqbLG466x7pnY8rFnrtfFaKL3Il9sKuSu/63jhJRePHHpGIJW/wv2b4Mz/mIVgurqovrAiXdalS23LfDcddf/Fw7lwAk/mljml6YM7s1jc8awLreUua+u9Mwc+yFnQUx/3yhP7YM00XuJhZsLuen11QzvE8XzV2USVlMCi/5mDUSmzbQ7PO8x5XZrlaKP7/JMwbO6GmsQOHGUV64L6y4zMxJ48IJRLM4+wI2vrqK61s3JPiAQJv0ccpdZ6yArl9JE7wU+3VDAja+tYmhid169diKRoUHw5Z+hphzO+Js1GKksQSFwzqPWx/xFf3X/9da+ZpUjnnZvl/t/uGBcEn89bwSLthZzwysrqapxc7Ifc7l1p+zSJ9x7nS5IE73NPlqfz83/Wc2IvtG8ev1EorsFWy2a1S/DhBsgbojdIXqffpMg81pY9hTs+tZ916kqtd5MkidC6mnuu44Xu3RiPx66YCTfZe/nmpeWu7fPPrQ7ZF5jFTor3uq+63RBmuht9O/Fu7j1jTWM69eDV66bSFRYsNVVMO8WiIyHab+1O0TvddqfoNcgePcGqDjgnmt89Teo2A+z/l+Xa803dWFmMo9cPJoVu0uY8/wyisvcWK5g8m0QFG51lymX0URvA4fD8JePN/HHDzdx+rB4XrlugtVdA7D4ESjaBGf/01o7VbUsNBIu+JdV/XDeLdbSfq6Uv8Za4SrzGugzxrXn9kGzR/fl2cvHsb2wnPOeWkx2kZvGRyJ6w8SfWYXOCm2ov+OnNNF7WGVNHbe+sYbnv93FVSf056nLxhEWHGi9WLgRvnnIqnPuoyVwPSpxFJx2v1Un5au/ue68NZXwzg0QGQen/t515/VxM4bF8+bcSVTX1vPTp5ewbKebPklNvhVCIq1xKuUSmug9aM+BCs5/agnzN+zjt2cO4f/OHU5gQEOXQE0lvH0thMXArAdtjdOnTLwRRl9ufdRf/1bnz2cMLPgNHNgOP3kauvXs/Dn9yKjkGN67aQq9I0O47IXvefG7XRhXf5rq1hNOuhO2fgzZC1177i5KE72HLNpSxDmPf0fB4WpevmYCc6cOQpr2+y74LRRvgfOfhUhdeMVpIlY314CT4P2fw6bmSyW00/LnYNVLMOUX1t2a6keSe3bj/ZunMGNoHH/6aBO3vbmWyhoXD9KecDP0GAif3gP1ta49dxekid7Nqmvruf/DTVzz0gqSenTjw1tOZGpas0Se9Tas+rc1R3zQqfYE6suCQuCS/0CfsfC/q61/z47IettKLEPOhul/cGmI/qZ7WDDPXD6OX81M5+P1+Zz9+Heszyt13QWCQmHmA9YNg3oTVadponejzfsOM/uJxfxr8S6unjyAd2+aTHLPbj/cae9q+OBm6HeCNVdbdUxYFFzxrjUV8p3rrP5dZ+umGGMNvL5zPfSbDOc/51dliN1FRLjplMG8dt1EKo/Uc/5TS3jiy+3UO1zUlZN2hvWmu+gvOjDbSeLy/jUXyMzMNCtXrrQ7jA6rrq3nqUXZPP31DqLDQ3jowpFMS4/78Y6H98Hz0yAgGOYusmYcqM6pOwIf32kVyOozxurWOd6smfIi+ORXsPE9SD0DLnwJQrq1vr9q0aHKWu79YAMfrstnVHIMf/lJBhl9XTBrrGI/PDnRKn9x/cKusYRmB4nIKmNMZouvaaJ3re+27+fe97PYfaCS88b05d6zhtIrsoVVkSoPwr/PhNIcuO4zSMjwfLD+yhjY8I7VDVNRbHWHZVxg3WjVPQFqq6Foo7WE3ZrXof4ITPud1S+vLflOmbcun/s/3MjBihqumTKQO05LOzZ1uKM2fwT/vQwm/AzO1IkKrdFE7wHZReU8tGALCzYWMqBXN/78kxGcmNpKC736ELx8LhRthsv+ByknezbYrqKqFFa+CMtfgLL8H78eGArDZlsloHsP9nh4/upQZS0PLtjCf5bnENc9lDtmpHHBuCSCAjvxJvrpb2HZkzD7KRhzmeuC9SOa6N2o4FA1jy7czlsrcwkPDuRnU1O4YWrKsbnxzVXsh9cvhIL11gBi2hmeDbgrMsb69y7cCGUF1pq7vQZD8gS9Kc2NVueU8KePNrEmp5RBsRH88owhnDE8/oezzZxVXwevnQ85S+Gyt7Vx1AJN9G6ws7ic577Zybur92IwXDaxP7eeOrjlbppGJbvh1fPhcL7VF5yuVSmVfzPG8NmmQh78dAs7iisYmhjFjSencNaIxPa38Jt2d175vvVGrY7SRO8iDodh6c4DvLZsD59uLCAkMIALM5P42dRBP55N09z2L6y6LBi49C39JVVdSl29g/fX5vPM1zvILiqnb0w410wZwE/HJtEjoh0DrGUF8O9Z1ifji1+FlFPcFrOv0UTfSUVl1by7ei9vLs9h94FKosODuXxSP66ePJDY7sdpwYM1C+SrB+C7f0L8cLjoFasYl1JdkMNh+HJLEc9+s4MVu0sICQzgjIwE5oxPZlJKLwICnOjWObQXXr/AmmN/zqNWeWOlib4jisuO8OmGfXyctY/vdx3EGJgwsCeXTujHzIyE1vvgm8pZBvNug/1brV/GWQ/p1D2lGmwpOMyby3N5d3Ueh6vriI8KZVZGImeOSCSzf4/jJ/2qUnjrStj1NYyaA2c+ZJU57sI00Tuhrt7B2txSvtm+n2+2FbM+rxSHgcFxkZw5IpFzR/VhcFykcycr3GTd5LHlI4juB2c/3GXrmSvVluraehZsLGB+1j6+2lrMkToHvSNDmZram5PSenPi4NiWPznX11lFAL95ELonwul/sgoCdtGS0proW1BWXcva3FJW7ylldU4Jq3NKKKuuI0BgdHIMU9NimZWRSFp8pHOzBBwO2LkIVrwAWz+xWheTb4VJN1kldZVSbSo/UseXW4r4fFMhi7P3c7CiBoAhCd0Z278HY5JjGNu/BwN7RRxr8ecut5aWLFgPSRPgxDus5Te72D0RXTrR19Q5yDlYwdaCcrYWlrGtoIxthWXsOlCBMdabf1qc9Ut0Umpvpgzqba3y5AyHA/JWWItVb55nzQbo1hvGXQUn3KKVD5XqBIfDsGnfYb7ZXszSHQdYm1tKWbVVPC06PJj0hO6kx3cnLaE76bHdGF44j27LH0FKc6DnIBh5MYy4oMuMiXU60YvITOBRIBB4wRjzQLPXpeH1M4FK4GpjzGpnjm1JRxJ9vcPw7uo8ckuqyDtYSV5JFbkllRQcrj66JkWAwIBeEaTFd2doYhRj+8cwKjnGWtmpLcZY0yL3b7UWpchZBrnfWzc/BQRblQ5HXGjdgBPUxgCtUqrdHA7DjuJy1uSUsia3lG0NDbeyJssbxoQKl0Ss5DzHF6RXrwOgPHIAFX2mYAZMISJ5DJF90pDATt6t64U6lehFJBDYBpwG5AErgDnGmE1N9jkTuBUr0U8EHjXGTHTm2JZ0JNEbY8j4wwIqa+tJiAojuUc3knqGk9yjG/16diM9oTuD4yIJCwoA4wBHHdRVw5FyaxHumvJj31eVWNO4ygqgvMAa5d+/HWrKjl2wd7p1S/3AqVb/u954o5THGWMoOFzN1oIysovKySupIq+kktyDVdSW5DCtfilTAjYwIWALkVINQJUJYU9AEgeC4jkUHEdZaAI1Yb0hLIaAbtEEdutBcLdoAkO7ERQSSkhoOCHBIYQGBxEaHEBYUCChwQGEBgUQHBhAYIAQKEJgYMNjgPUVFCAduzmsg46X6J15W5sAZBtjdjac7E1gNtA0Wc8GXjHWu8YyEYkRkURggBPHuoSIsDb2PoLqKhFHPZTVwaE62FEPpt5K7I1fzgqLhsgEiEqE0XOgdxrEpkPccIjo5eofQSnVTiJCYnQ4idHhnNKscKAxhpLKiygqq2bdoQpq8rMIKNpEeMkWosp3MPDIXnpUrSW8sqrN6ziMUEMQNQRRSxAOBBAcCHUINQim8ctY248+F2tf0xAvQGP6NzR5LlARGE3GvUtc9c9zlDOJvi+Q2+R5Hlarva19+jp5LAAiMheYC9CvXz8nwvqx4H7jrX7zgEAICGry2Ox7CbSeB4VaS5aFdm94jLQew6Kt4lfB4R2KQyllPxGhZ0QIPSNCICEK0hOB03+4kzFW92vlAagqxVSXUltewpGKEuprqqivOUJdbTWO2iPWV90RTF0N9Y56HA6DcTgwxmGtsvWDRwfG4QCObTfGupzBHL320f4UAwaoD3HPxA1nEn1Lnz2a9/e0to8zx1objXkOeA6srhsn4vqx2bpAgVKqHUQgPMb6wkpYIQ1f/sSZRJ8HJDd5ngQ0LwXY2j4hThyrlFLKjZyZaLoCSBWRgSISAlwCNF+Ycx5wpVgmAYeMMfucPFYppZQbtdmiN8bUicgtwAKsKZL/MsZsFJEbG15/BpiPNeMmG2t65TXHO9YtP4lSSqkW+f0NU0op1RUcb3pl17pHWCmluiBN9Eop5ec00SullJ/TRK+UUn7OKwdjRaQY2NPBw3sD+10Yjjv5UqzgW/H6UqzgW/H6UqzgW/F2Jtb+xpjYll7wykTfGSKysrWRZ2/jS7GCb8XrS7GCb8XrS7GCb8Xrrli160YppfycJnqllPJz/pjon7M7gHbwpVjBt+L1pVjBt+L1pVjBt+J1S6x+10evlFLqh/yxRa+UUqoJTfRKKeXn/DrRi8jdImJEpLfdsbRGRB4SkS0isl5E3hORGLtjak5EZorIVhHJFpF77I7neEQkWUQWichmEdkoIrfbHVNbRCRQRNaIyEd2x9KWhmVC3274nd0sIifYHVNrROSOht+BDSLyhoiE2R1TUyLyLxEpEpENTbb1FJHPRWR7w2MPV1zLbxO9iCRjLUqeY3csbfgcyDDGjMRaSP03NsfzAw0LvD8JzAKGAXNEZJi9UR1XHXCXMWYoMAm42cvjBbgd2Gx3EE56FPjUGDMEGIWXxi0ifYHbgExjTAZWmfRL7I3qR14CZjbbdg+w0BiTCixseN5pfpvogX8Cv6KVpQu9hTHmM2NM44rly7BW4fImRxeHN8bUAI0LvHslY8w+Y8zqhu/LsBJRX3ujap2IJAFnAS/YHUtbRCQKmAq8CGCMqTHGlNoa1PEFAeEiEgR0w8tWtzPGfAMcbLZ5NvByw/cvAz9xxbX8MtGLyLnAXmPMOrtjaadrgU/sDqKZ1hZ+93oiMgAYA3xvcyjH8whWg8RhcxzOSAGKgX83dDW9ICIRdgfVEmPMXuDvWJ/o92GteveZvVE5Jb5hdT4aHuNccVKfTfQi8kVD31vzr9nA74D77I6xURuxNu7zO6xuh9fti7RFTi/w7k1EJBJ4B/iFMeaw3fG0RETOBoqMMavsjsVJQcBY4GljzBigAhd1LbhaQ9/2bGAg0AeIEJHL7Y3KPs4sDu6VjDEzWtouIiOw/nPXiQhYXSGrRWSCMabAgyEe1VqsjUTkKuBsYLrxvhsbnFkc3quISDBWkn/dGPOu3fEcxxTgXBE5EwgDokTkNWOMtyakPCDPGNP4CeltvDTRAzOAXcaYYgAReReYDLxma1RtKxSRRGPMPhFJBIpccVKfbdG3xhiTZYyJM8YMMMYMwPrlHGtXkm+LiMwEfg2ca4yptDueFvjUAu9ivbu/CGw2xjxsdzzHY4z5jTEmqeH39BLgSy9O8jT8DeWKSHrDpunAJhtDOp4cYJKIdGv4nZiOlw4cNzMPuKrh+6uAD1xxUp9t0fuRJ4BQ4POGTyDLjDE32hvSMT64wPsU4AogS0TWNmz7rTFmvn0h+ZVbgdcb3vR3AtfYHE+LjDHfi8jbwGqsLtE1eFkpBBF5AzgF6C0iecAfgAeAt0TkOqw3qwtdci3v6ylQSinlSn7XdaOUUuqHNNErpZSf00SvlFJ+ThO9Ukr5OU30Sinl5zTRK6WUn9NEr5RSfu7/AwauaxAzau99AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-4, 10, 1000)\n",
"plt.plot(x, norm.pdf(x, p_mu, p_sigma))\n",
"plt.plot(x, w1*norm.pdf(x, q_mu1, q_sigma1) + w2*norm.pdf(x, q_mu2, q_sigma2))\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"KL(q, p): 0.674923650198428\n"
]
}
],
"source": [
"# calc KL\n",
"z = bernoulli.rvs(p=w1, size=N_MONTE)\n",
"x_q = z * norm.rvs(loc=q_mu1, scale=q_sigma1, size=N_MONTE) + (1-z)* norm.rvs(loc=q_mu2, scale=q_sigma2, size=N_MONTE)\n",
"KL_div = np.sum((np.log(w1*norm.pdf(x_q, q_mu1, q_sigma1) + w2* norm.pdf(x_q, q_mu2, q_sigma2)) - np.log(norm.pdf(x_q, p_mu, p_sigma)))) / N_MONTE\n",
"print(f\"KL(q, p): {KL_div}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment