{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Intitutive Explanation of Expectation Maximization Algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Probabilistic models, such as hidden Markov models or Bayesian networks, are commonly used to model biological data. Much\n", "of their popularity can be attributed to the existence of efficient and robust procedures for learning parameters from observations. \n", "\n", "Often, however, the only data available for training a probabilistic model are incomplete. Missing values can occur, for example, in medical diagnosis, where patient histories generally include results from a limited battery of tests.\n", "\n", "Alternatively, in gene expression clustering, incomplete data arise from the intentional omission of gene-to-cluster assignments in the probabilistic model. \n", "\n", "**The expectation maximization algorithm enables parameter estimation in probabilistic models with incomplete data.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What is EM?\n", "\n", "- The Expectation-Maximization (EM) algorithm is a way to find **maximum-likelihood estimates** for model parameters when your data is incomplete, has missing data points, or has unobserved (hidden) latent variables.\n", "\n", "- It is an **iterative way to approximate the maximum likelihood function.**\n", "\n", "- While maximum likelihood estimation can find the “best fit” model for a set of data, it doesn’t work particularly well for incomplete data sets.\n", "\n", "- The more complex EM algorithm can find model parameters even if you have missing data. It works by choosing random values for the missing data points, and using those guesses to estimate a second set of data. The new values are used to create a better guess for the first set, and the process continues until the algorithm converges on a fixed point." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How EM is different from MLE?\n", "\n", "Although Maximum Likelihood Estimation (MLE) and EM can both find “best-fit” parameters, how they find the models are very different.\n", "\n", "- MLE accumulates all of the data first and then uses that data to construct the most likely model.\n", "- EM takes a guess at the parameters first — accounting for the missing data — then tweaks the model to fit the guesses and the observed data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How EM works?\n", "\n", "The basic steps for the algorithm are:\n", "\n", "1. An initial guess is made for the model’s parameters and a probability distribution is created. This is sometimes called the “E-Step” for the “Expected” distribution.\n", "2. Newly observed data is fed into the model.\n", "3. The probability distribution from the E-step is tweaked to include the new data. This is sometimes called the “M-step.”\n", "4. Steps 2 through 4 are repeated until stability (i.e. a distribution that doesn’t change from the E-step to the M-step) is reached." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This will try to use very less math and just give intution, how EM converges." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Expectation Maximization using Two Coin example: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- We have two biased coin and we dont know probability of getting head on each coins. Assume $\\theta_A$ and $\\theta{_B}$\n", "- We randomly choose one coin and toss it 10 times and repeat this experiment five times.\n", "- A simple way to estimate $\\theta_A$ and $\\theta_B$ is to return the observed proportions of heads for each coin:\n", "\n", " $$\\hat{\\theta_A} = \\frac{Number ~of~ heads ~using ~coin~ A}{total ~number~ flips~ using~ coin ~A }$$ \n", " \n", " $$\\hat{\\theta_B} = \\frac{Number~ of~ heads ~using ~coin ~B}{total number~ flips ~using~ coin~ B }$$\n", " \n", " \n", "This intuitive guess is, in fact, known in the statistical literature as Maximum Likelihood Estimation . **As shown in below figure (a)**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**More Complicated Problem**:\n", "- If in experiment we just know count of head, but which coin is used for experiment is not known.\n", "- How to estimated probablities of head for each coin? This problem can be solved using Expectation Maximization.\n", "\n", "Steps are shown in above **figure (b)**:\n", " 1. Assume some probablity of head for both the coins.\n", " 2. Create the probability distribution (using binomial distribution). (Expectation)\n", " Example: for second experiment we have 9 heads and 1 tail.\n", " \n", " $$p(A) = \\theta_A^9(1-\\theta_A)^{10-9} \\approx 0.004$$\n", " \n", " $$p(B) = \\theta_B^9(1-\\theta_B)^{10-9} \\approx 0.001$$\n", " \n", " This give us:\n", " $$\\frac{0.004}{0.004+0.001}=0.8$$\n", " \n", " $$\\frac{0.001}{0.004+0.001}=0.2$$\n", " \n", " Similiary we can calculate the distribution for Coin A and Coin B for each experiment.\n", " \n", " 3a. Based on distribution calculate the new estimated probability. (Maximization)\n", " 3b. Repeat above two steps.\n", " 4. Give the estimated probablity after convergence.\n", " \n", "Theoritically it is proven the EM always converges [ref](https://www.cmi.ac.in/~madhavan/courses/dmml2018/literature/EM_algorithm_2coin_example.pdf)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Expecation maximization for Nomral distribution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this experiment we will try to find mean and std of two normal distribution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Lets generate some data" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "import warnings\n", "\n", "warnings.filterwarnings('ignore')\n", "\n", "np.random.seed(110) # for reproducible random results\n", "\n", "# set parameters\n", "red_mean = 3\n", "red_std = 0.8\n", "\n", "blue_mean = 7\n", "blue_std = 2\n", "\n", "# draw 20 samples from normal distributions with red/blue parameters\n", "red = np.random.normal(red_mean, red_std, size=20)\n", "blue = np.random.normal(blue_mean, blue_std, size=20)\n", "\n", "both_colours = np.sort(np.concatenate((red, blue))) " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import matplotlib.mlab as mlab\n", "\n", "plt.rcParams[\"figure.figsize\"] = (20,3)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5,1,'Data(With Label)')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABIsAAADSCAYAAADKQUBCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGMpJREFUeJzt3X2UbXdZH/Dvk/sihgQDuVeKSciFikpEKGFAkAKpuly8qPGlKhQVqF2xCbbapVaQtlg0Yit1WZYUTDUgJgYRaRsxVlhRiIiANwTCSwxGDOSSSK5gYiJUCXn6xzmzOUzmzpyZOXPPmcnns9ZZc85v7/3bz/6dvc+c+71776nuDgAAAAAkyQnzLgAAAACAxSEsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAO41quqlVfUjm1z22VX15jWmn1NVRzZf3dR1vKaqfmary1bVI6vqHbOtDgDYDYRFAMBxV1U3VtVnquqOqrqtqt5RVf+6qqb6blJVh6qqq2rvBtZ5MMn3J/nl8evrq+q7J6Y/cdznyrY7q2pvd1/a3d80Ma2r6sunXf8q9by1qv7VZpffqu6+NsltVfUt86oBAFhMwiIAYF6+pbtPTnJmkp9L8hNJfnUb1/fcJFd092fGr69K8pSJ6U9O8mertL2ju+/axrrm6dIkPzjvIgCAxSIsAgDmqrtv7+7Lk3xPkudU1SOSpKqeUVXXVNXfVtVNVfVTE4tdNf552/jMnydU1T+uqj+oqk9W1V9X1aVVdcrEMk9L8rYVfTx54vWTkvyXVdquGtfz3Kp6+/j58vrfN17/9ywvUFU/WlW3VtUtVfW8zYxJVf1WVf1VVd1eVVdV1VevmOVAVb1lfGbW26rqzIllv2o87VMrz55axVuTfENVfdFm6gQAdidhEQCwELr73UmOZBTQJMnfZXTZ2ClJnpHk/Kr6tvG05UDnlO4+qbv/JEkleWmSL0vy8CRnJPmpiVV8TZLrJ16/LclXV9UDxpe/LSX5zSSnTLR9XT4fTE3Wurz+R43X/5vj1/8oyZckOS3JDyR5RVXdf8ODkfxekocl+dIk78noDKBJz07y00kOJHnv8vSqum+StyT5jfGyz0ryP1YJm5a34+NJPpvkKzdRIwCwSwmLAIBFcnOSByRJd7+1u9/f3XeP769zWb7wErEv0N03dPdbuvvvu/tokl9YMf8pSe6YmP9jST6WUTj1qCR/Pr5E7Y8n2u6T5F0bqP+zSV7S3Z/t7iuS3JlNBDHdfXF339Hdf59R4PWoqvqSiVl+t7uvGk9/UZInVNUZSb45yY3d/eruvqu735Pkt5P88zVWd0dGYwMAkCSZ+qaQAADHwWlJPpUkVfW1Gd3L6BFJ9if5oiS/dawFq+pLk7w8o6Dn5Iz+U+xvJmb5m3H7pOVL0T6W5I/GbW+faHvXOJCZ1idX3N/o00lO2sDyqao9SS5M8l1JDia5ezzpQJLbx89vWp6/u++sqk9ldEbVmUm+tqpum+hyb5JfX2OVJye5bY3pAMC9jDOLAICFUFWPzSgsevu46TeSXJ7kjO7+kiSvyuhSsyTpVbp46bj9kd19vyTfOzF/klyb5CtWLLMcFj0pnw+L/mii7R6XoB0H/yLJuUm+MaNL2g6N2ye35YzlJ1V1UkZnY92cUYj0tu4+ZeJxUnefv9qKqurLMgrirl9tOgBw7yQsAgDmqqruV1XfnOR1SS7p7vePJ52c5FPd/f+q6nEZhSjLjmZ0xs1DJ9pOzuiyr9uq6rQkP75iVVfknpexXZXk0eP2Px63vT/JQ5L8s6wdFn1ixfo3Y29V3WfisS+j7fj7JJ9McmKSn11luadX1T+tqv0Z3bvoXd19U5I3JfmKqvq+qto3fjy2qh5+jPWfk+QPNnj2FACwywmLAIB5+Z2quiOjs2FelNE9hib/etgFSV4ynuc/JXn98oTu/nRGl2r9cVXdVlWPT/Kfk5yd0aVav5vkjSvW99qMQpYvnujnw0luTXJLd982brs7ybuT3C/JO9ao/6eS/Np4/Wv9xbG1vDLJZyYerx7X+dEkH0/yoSTvXGW530jy4owu2XtMRje8TnffkeSbkjwzozON/iqjv/B2rL929uyMztgCABhU92pncQMA7D5V9bNJbu3uX5x3LfNWVV+T5KLufsK8awEAFouwCAAAAICBy9AAAAAAGAiLAAAAABgIiwAAAAAYCIsAAAAAGOyddwErHThwoA8dOjTvMgAAAAB2jauvvvqvu/vgNPMuXFh06NChHD58eN5lAAAAAOwaVfXRaed1GRoAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAg3XDoqq6uKpuraoPHGN6VdXLq+qGqrq2qs5eMf1+VfXxqvqlWRUNAAAAwPaY5syi1yR56hrTn5bkYePHeUleuWL6Tyd522aKAwAAAOD4Wjcs6u6rknxqjVnOTfLaHnlnklOq6kFJUlWPSfLAJG+eRbEAAAAAbK9Z3LPotCQ3Tbw+kuS0qjohyX9L8uMzWAcAAAAAx8EswqJapa2TXJDkiu6+aZXpX9hB1XlVdbiqDh89enQGJQEAAACwGXtn0MeRJGdMvD49yc1JnpDkSVV1QZKTkuyvqju7+wUrO+jui5JclCRLS0s9g5oAAAAA2IRZhEWXJ/mhqnpdkq9Ncnt335Lk2cszVNVzkyytFhQBAAAAsDjWDYuq6rIk5yQ5UFVHkrw4yb4k6e5XJbkiydOT3JDk00met13FAgAAALC91g2LuvtZ60zvJM9fZ57XJHnNRgoDAAAA4PibxQ2uAQAAANglhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAzWDYuq6uKqurWqPnCM6VVVL6+qG6rq2qo6e9z+T6rqT6rqg+P275l18QAAAADM1jRnFr0myVPXmP60JA8bP85L8spx+6eTfH93f/V4+V+sqlM2XyoAAAAA223vejN091VVdWiNWc5N8tru7iTvrKpTqupB3f3hiT5urqpbkxxMctsWawYAAABgm8zinkWnJblp4vWRcdugqh6XZH+Sv5jB+gAAAADYJrMIi2qVth4mVj0oya8neV53371qB1XnVdXhqjp89OjRGZQEAAAAwGbMIiw6kuSMidenJ7k5Sarqfkl+N8l/6O53HquD7r6ou5e6e+ngwYMzKAkAAACAzZhFWHR5ku8f/1W0xye5vbtvqar9Sf5XRvcz+q0ZrAcAAACAbbbuDa6r6rIk5yQ5UFVHkrw4yb4k6e5XJbkiydOT3JDRX0B73njR707y5CSnVtVzx23P7e73zrB+AAAAAGZomr+G9qx1pneS56/SfkmSSzZfGgAAAADH2ywuQwMAAABglxAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwWDcsqqqLq+rWqvrAMaZXVb28qm6oqmur6uyJac+pqj8fP54zy8IX2qWXJocOJSecMPp56aXzrmj3Wh7rqmTv3tHP3TbmK/enCy5YfZsvuCA5cGD0umr0/HiNw2b2+dWW2eqxc+mla4/BLI7NafpYkM+AmZaxVmfrjftmul2QMVzLDihxR5nF4b/a8hvpd73dfFHe70WqZbdalDGe1f67kfkXZduPt+34WgKM/omy/E+WPXuSk07aWcfUvfpzoLvXfCR5cpKzk3zgGNOfnuT3klSSxyd517j9AUk+Mv55//Hz+6+3vsc85jG9o11ySfeJJ3Ynn3+ceOKondlabax325ivtY3TPPbv3/5x2Mw+v9oy+/aN6t3s+3jJJaM+jjUGszg2p+ljQT4DZlrGWp1dcsk937fl93OdlR2z2/P/aCHGcC0L8jbvGlsdz2Mtf/750/e73m6+KO/3ItWyWy3KGG+kjo3WPItjZjfZjq8lwOgzZa1/qiz6MbUovw9mKcnhXieTWX7UaP61VdWhJG/q7kesMu2Xk7y1uy8bv74+yTnLj+7+wdXmO5alpaU+fPjwujUtrEOHko9+9J7tZ56Z3Hjj8a5mdzvWWC/bDWO+3jZOY7vHYTP7/Ea2a9r61+rzzDNHP7d6bE6zrQvyGTDTMtbqLFl73NdY2TG73XMkN37ujA33dzwtyNu8a2x1PI+1/J49yec+N12/m9nN5/F+2/e236KM8Ubq2GjNszhmdpPt+FoCjM4oWu0zZdIiH1OL8vtglqrq6u5emmreGYRFb0ryc9399vHrK5P8REZh0X26+2fG7f8xyWe6+2Wr9HFekvOS5MEPfvBjPrrVfxzP0wknjELHlaqSu+8+/vXsZsca62W7YczX28ZpbPc4bGaf38h2TVv/Wn1WjX5u9dicZlsX5DNgpmWs1Vmy9rivsbJjdpu7c3f2bLi/42lB3uZdY6vjudGPytX63cxuPo/32763/RZljDdSx0ZrnsUxs5tsx9cS4PO/Q9ebZ1GPqUX5fTBLGwmLZnGD69V2gV6j/Z6N3Rd191J3Lx08eHAGJc3Rgx+8sXY2b70x3Q1jPott2O5x2Mw+v5Gapp13vfXN4ticpo8F+QyYaRlrdbaF9/mY3e65eVP9HU8L8jbvGlsdz2PNt2eVzPFY829mN5/H+71ItexWizLGG6ljozXP4pjZTbbjawlw7M+USYt8TC3K74N5mUVYdCTJ5PUCpye5eY323e3CC5MTT/zCthNPHLUzW6uN9bLdMuZrbeM09u/f/nHYzD6/2jL79o3q3Ug/K/vct++e7ctjMItjc5o+FuQzYKZlrNXZhRfe831LRu/FOis7Zrfn3bgQY7iWBXmbd42tjuexlj/vvOn7XW83X5T3e5Fq2a0WZYw3UsdGa57FMbObbMfXEmD0mbKWRT+mFuX3wdxMc2OjJIdy7BtcPyNfeIPrd4/bH5DkLzO6ufX9x88fsN66dvwNrrtHd7w688zuqtHPnXwHrEW3PNZJ9549o5+7bcxX7k/nn7/6Np9/fvepp37+7munnnr8xmEz+/xqy2z12LnkkrXHYBbH5jR9LMhnwEzLWKuz9cZ9M90uyBiuZQeUuKPM4vBfbfmN9Lvebr4o7/ci1bJbLcoYz2r/3cj8i7Ltx9t2fC0BRv9EWf4nywkndN/3vjvrmNptnwOZ5Q2uq+qyjO4/dCDJJ5K8OMm+cdD0qqqqJL+U5KlJPp3ked19eLzsv0zyk+OuLuzuV68XXu34G1wDAAAALJiN3LNo73ozdPez1pneSZ5/jGkXJ7l4mkIAAAAAmL9Z3LMIAAAAgF1CWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAYKqwqKqeWlXXV9UNVfWCVaafWVVXVtW1VfXWqjp9Ytp/raoPVtV1VfXyqqpZbgAAAAAAs7NuWFRVe5K8IsnTkpyV5FlVddaK2V6W5LXd/cgkL0ny0vGyX5fkiUkemeQRSR6b5Ckzqx4AAACAmZrmzKLHJbmhuz/S3f+Q5HVJzl0xz1lJrhw//8OJ6Z3kPkn2J/miJPuSfGKrRQMAAACwPaYJi05LctPE6yPjtknvS/Kd4+ffnuTkqjq1u/8ko/DolvHj97v7uq2VDAAAAMB2mSYsWu0eQ73i9Y8leUpVXZPRZWYfT3JXVX15kocnOT2jgOnrq+rJ91hB1XlVdbiqDh89enRDGwAAAADA7EwTFh1JcsbE69OT3Dw5Q3ff3N3f0d2PTvKicdvtGZ1l9M7uvrO770zye0kev3IF3X1Rdy9199LBgwc3uSkAAAAAbNU0YdGfJnlYVT2kqvYneWaSyydnqKoDVbXc1wuTXDx+/rGMzjjaW1X7MjrryGVoAAAAAAtq3bCou+9K8kNJfj+joOf13f3BqnpJVX3reLZzklxfVR9O8sAkF47b35DkL5K8P6P7Gr2vu39ntpsAAAAAwKxU98rbD83X0tJSHz58eN5lAAAAAOwaVXV1dy9NM+80l6EBAAAAcC8hLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYDBVWFRVT62q66vqhqp6wSrTz6yqK6vq2qp6a1WdPjHtwVX15qq6rqo+VFWHZlc+AAAAALO0blhUVXuSvCLJ05KcleRZVXXWitleluS13f3IJC9J8tKJaa9N8vPd/fAkj0ty6ywKBwAAAGD2pjmz6HFJbujuj3T3PyR5XZJzV8xzVpIrx8//cHn6OFTa291vSZLuvrO7Pz2TygEAAACYuWnCotOS3DTx+si4bdL7knzn+Pm3Jzm5qk5N8hVJbquqN1bVNVX18+MzlQAAAABYQNOERbVKW694/WNJnlJV1yR5SpKPJ7kryd4kTxpPf2yShyZ57j1WUHVeVR2uqsNHjx6dvnoAAAAAZmqasOhIkjMmXp+e5ObJGbr75u7+ju5+dJIXjdtuHy97zfgStruS/O8kZ69cQXdf1N1L3b108ODBTW4KAAAAAFs1TVj0p0keVlUPqar9SZ6Z5PLJGarqQFUt9/XCJBdPLHv/qlpOgL4+yYe2XjYAAAAA22HdsGh8RtAPJfn9JNcleX13f7CqXlJV3zqe7Zwk11fVh5M8MMmF42U/l9ElaFdW1fszuqTtf858KwAAAACYiepeefuh+VpaWurDhw/PuwwAAACAXaOqru7upWnmneYyNAAAAADuJYRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyqu+ddwxeoqqNJPjrvOmAVB5L89byLgB3MMQRb4xiCrXMcwdY4hna2M7v74DQzLlxYBIuqqg5399K864CdyjEEW+MYgq1zHMHWOIbuPVyGBgAAAMBAWAQAAADAQFgE07to3gXADucYgq1xDMHWOY5gaxxD9xLuWQQAAADAwJlFAAAAAAyERbCGqjqjqv6wqq6rqg9W1Q/PuybYiapqT1VdU1VvmnctsBNV1SlV9Yaq+rPx76QnzLsm2Emq6t+Nv8t9oKouq6r7zLsmWHRVdXFV3VpVH5hoe0BVvaWq/nz88/7zrJHtIyyCtd2V5Ee7++FJHp/k+VV11pxrgp3oh5NcN+8iYAf770n+b3d/VZJHxfEEU6uq05L82yRL3f2IJHuSPHO+VcGO8JokT13R9oIkV3b3w5JcOX7NLiQsgjV09y3d/Z7x8zsy+nJ+2nyrgp2lqk5P8owkvzLvWmAnqqr7JXlykl9Nku7+h+6+bb5VwY6zN8kXV9XeJCcmuXnO9cDC6+6rknxqRfO5SX5t/PzXknzbcS2K40ZYBFOqqkNJHp3kXfOtBHacX0zy75PcPe9CYId6aJKjSV49vpzzV6rqvvMuCnaK7v54kpcl+ViSW5Lc3t1vnm9VsGM9sLtvSUb/sZ7kS+dcD9tEWARTqKqTkvx2kh/p7r+ddz2wU1TVNye5tbuvnnctsIPtTXJ2kld296OT/F2c9g9TG99T5dwkD0nyZUnuW1XfO9+qABabsAjWUVX7MgqKLu3uN867HthhnpjkW6vqxiSvS/L1VXXJfEuCHedIkiPdvXxm6xsyCo+A6Xxjkr/s7qPd/dkkb0zydXOuCXaqT1TVg5Jk/PPWOdfDNhEWwRqqqjK6R8R13f0L864HdprufmF3n97dhzK6megfdLf/zYUN6O6/SnJTVX3luOkbknxojiXBTvOxJI+vqhPH3+2+IW4SD5t1eZLnjJ8/J8n/mWMtbKO98y4AFtwTk3xfkvdX1XvHbT/Z3VfMsSYA7n3+TZJLq2p/ko8ked6c64Edo7vfVVVvSPKejP7S7TVJLppvVbD4quqyJOckOVBVR5K8OMnPJXl9Vf1ARkHsd82vQrZTdfe8awAAAABgQbgMDQAAAICBsAgAAACAgbAIAAAAgIGwCAAAAICBsAgAAACAgbAIAAAAgIGwCAAAAICBsAgAAACAwf8HujtWMn5VZnUAAAAASUVORK5CYII=\n", "text/plain": [ "