{
 "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 <font color ='red'>Maximum Likelihood Estimation </font>. **As shown in below figure (a)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://i.stack.imgur.com/mj0nb.gif\" width=\"500\" height=\"500\">"
   ]
  },
  {
   "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 <font color='red'>Expectation Maximization.</font>\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": [
       "<Figure size 1440x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = [1 for i in range(len(red))]\n",
    "plt.plot(red, y, 'ro', blue, y, 'bo')\n",
    "plt.title(\"Data(With Label)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**If label is given it is easy problem and we can calucate the mean and std using MLE**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Red Mean is 2.8132116984626867\n",
      "Red Std is 0.8001119672360739\n",
      "\n",
      "Blue Mean is 6.972553898467638\n",
      "Blue Std is 2.0388437926931493\n"
     ]
    }
   ],
   "source": [
    "print(\"Red Mean is {}\".format(np.mean(red)))\n",
    "print(\"Red Std is {}\".format(np.std(red)))\n",
    "print(\"\")\n",
    "print(\"Blue Mean is {}\".format(np.mean(blue)))\n",
    "print(\"Blue Std is {}\".format(np.std(blue)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**If label is not given. USE EM** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Data(Without Label)')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIsAAADSCAYAAADKQUBCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGZhJREFUeJzt3X+0dXVdJ/D3Bx7ATBSFJ5TfMoFCjqbd/JFLYHRWS6mRfkw/HEuhZlH+aKypZvzRGhtnyBqrpY6lUaEyEmaOFRWNuiiiH2o9hKJGEDEKjw/B4w8QpFGBz/xxzrM9Xu5zf57LOffyeq11F2d/93fv/dnfs/e5h/ez977V3QEAAACAJDlg1gUAAAAAMD+ERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAMC2UFWvraofX+eyz6+q9y0z/4yq2r3+6majqrqqvn6jy1bVL1fVj063OgBgXgmLAICpqqpPVNU/V9UdVXVbVf1VVf1oVa3qe0dVnTAOKnasYZs7k7wgya+Np6+tqu+dmP/08ToXt91ZVTu6+6Lu/taJeesOWdaqqi6vqn+/zPw1j8cmeF2SV1XVwTOsAQC4nwiLAIDN8G+6+9Akxyf5+ST/OclvbuL2zk5yaXf/83j6iiSnT8w/LcnfL9H2V9199ybWtS10980Zjd9zZ10LALD5hEUAwKbp7tu7+5Ik35fkhVX1uCSpqm+rqquq6vNVdVNV/ezEYleM/3vb+Mqfp1XVv6iqP6mqz1TVp6vqoqo6bGKZ5yT5s0XrOG1i+hlJfmGJtivG9ZxdVX8xfr1v+x8Zb//79i1QVT9ZVbdW1c1Vdc5E+8Oq6sKq2ltVn6yqn9l3JVVV/WxVvWOi73ClUFWdN67jTeNtvWmVQ7tvXU+uqg+Mr+C6uaretMTVP2dW1Q3jcXvd5BVeVfVDVXVNVX2uqt5bVccvs7nLk3zbWuoDALYmYREAsOm6+6+T7M4oGEmSL2R029hhGQUQL6qq7xjP2xfoHNbdD+nuDySpJK9NclSSU5Icm+RnJzbxL5NcOzH9Z0m+oaoeMQ5HFpL8dpLDJtq+JV8JpiZr3bf9J4y3/9vj6UcmeViSo5P8cJJfqaqHj+f9z/G8EzO6eukFSc7JCrr7VUn+PMlLx9t66UrLLHJPkp9IckSSpyV5VpIXL+rznRnt/5OSnJXkh5JkPN6vTPJdSXaO67h4mW1dk+QJa6wPANiChEUAwP1lT5JHJEl3X97dH+3ue7v76oxCitP3t2B3X9/d7+/uL3b33iS/vKj/YUnumOh/Y5IbMwqnnpDkH8a3qP3lRNuDknxoDfV/OclruvvL3X1pkjuTPKaqDszoyqlXdPcd3f2JJL+U5AfXsO516e4ru/uD3X33eLu/lvuO4y9092fHY/L6JM8bt/9Iktd29zXjW/F+Lsk3LnN10R0ZjTMAsM3N8kGJAMADy9FJPpskVfWUjJ5l9LgkByc5JMnv7G/Bqvq6JG/MKOg5NKN/8PrcRJfPjdsn7bsV7caMrppJkr+YaPtQd39xDfV/ZtHzje5K8pCMruo5OMknJ+Z9MqP93VRVdXJGwdlCkgdn9N3uykXdblpU11Hj18cneUNV/dLkKjOqe3Jf9jk0yW1TKBsAmHOuLAIANl1VfXNGIcRfjJt+K8klSY7t7ocleUtGQUWS9BKreO24/fHd/dAkPzDRP0muTnLyomX2hUXPyFfCoj+faLvPLWjr9OmMrjqavCLnuCSfGr/+QkZBzj6PXLT8Uvu7Wm/O6MHTJ43H5ZX56nFJRrfsTda1Z/z6piQ/0t2HTfx8TXf/1X62dUqSj2ygVgBgixAWAQCbpqoeWlXfnuSdSd7R3R8dzzo0yWe7+/9V1ZOT/LuJxfYmuTej5/9kov+dGT30+ugkP71oU5fmvrdfXZHkieP2vxy3fTTJo5P8qywfFt2yaPv71d33JHlXkvOq6tDxbVz/Mcm+h1p/OMlpVXVcVT0sySvWua1DqupBEz8HZDQun09yZ1U9NsmLlljup6vq4VV1bJKXZfTspmQU0L2iqr4hGR7S/T3LbP/0JH+8ijoBgC1OWAQAbIY/qKo7Mrp65VUZ3So1+cDnFyd5zbjPf8kobEmSdPddSc5L8pfjv/L11CT/NaMHNN+e5I+SvGfR9i7M6K9+fc3Eeq5LcmuSm7v7tnHbvUn+OslDk+zvCppk9PDst4+3/72r2N8fy+gKohsyunrqt5JcMN7m+zMKaK7O6BaxP1y07BuS/NvxXyR74zLbuDPJP0/8PDPJT2UUtN2R5NfzlSBo0u+Pt/vhjMbuN8d1/W5GfyHunVX1+SQfy+ivyt1HVT0qyalJfm+Z+gCAbaK6N3LlMwDAfKiqn0tya3e/fta1bDfj5xr9Y3f/6qxrAQA2n7AIAAAAgIHb0AAAAAAYCIsAAAAAGAiLAAAAABgIiwAAAAAY7Jh1AYsdccQRfcIJJ8y6DAAAAIBt48orr/x0d+9cTd+5C4tOOOGE7Nq1a9ZlAAAAAGwbVfXJ1fZ1GxoAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAA2ERAAAAAANhEQAAAAADYREAAAAAgxXDoqq6oKpuraqP7Wd+VdUbq+r6qrq6qp60aP5Dq+pTVfWmaRUNAAAAwOZYzZVFb0vy7GXmPyfJSeOfc5O8edH8/5bkz9ZTHAAAAAD3rxXDou6+Islnl+lyVpILe+SDSQ6rqkclSVV9U5Ijk7xvGsUCAAAAsLmm8cyio5PcNDG9O8nRVXVAkl9K8tNT2AYAAAAA94NphEW1RFsneXGSS7v7piXmf/UKqs6tql1VtWvv3r1TKAkAAACA9dgxhXXsTnLsxPQxSfYkeVqSZ1TVi5M8JMnBVXVnd7988Qq6+/wk5yfJwsJCT6EmAAAAANZhGmHRJUleWlXvTPKUJLd3981Jnr+vQ1WdnWRhqaAIAAAAgPmxYlhUVRcnOSPJEVW1O8mrkxyUJN39liSXJjkzyfVJ7kpyzmYVCwAAAMDmWjEs6u7nrTC/k7xkhT5vS/K2tRQGAAAAwP1vGg+4BgAAAGCbEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADAQFgEAAAAwEBYBAAAAMBAWAQAAADBYMSyqqguq6taq+th+5ldVvbGqrq+qq6vqSeP2b6yqD1TVx8ft3zft4gEAAACYrtVcWfS2JM9eZv5zkpw0/jk3yZvH7XcleUF3f8N4+ddX1WHrLxUAAACAzbZjpQ7dfUVVnbBMl7OSXNjdneSDVXVYVT2qu6+bWMeeqro1yc4kt22wZgAAAAA2yTSeWXR0kpsmpneP2wZV9eQkByf5xylsDwAAAIBNMo2wqJZo62Fm1aOS/K8k53T3vUuuoOrcqtpVVbv27t07hZIAAAAAWI9phEW7kxw7MX1Mkj1JUlUPTfJHSX6muz+4vxV09/ndvdDdCzt37pxCSQAAAACsxzTCokuSvGD8V9GemuT27r65qg5O8rsZPc/od6awHQAAAAA22YoPuK6qi5OckeSIqtqd5NVJDkqS7n5LkkuTnJnk+oz+Ato540W/N8lpSQ6vqrPHbWd394enWD8AAAAAU7Sav4b2vBXmd5KXLNH+jiTvWH9pAAAAANzfpnEbGgAAAADbhLAIAAAAgIGwCAAAAICBsAgAAACAgbAIAAAAgIGwCAAAAICBsAgAAACAgbAIAAAAgIGwCAAAAICBsAgAAACAgbAIAAAAgIGwCAAAAICBsAgAAACAgbAIAAAAgIGwCAAAAICBsAgAAACAgbAIAAAAgIGwCAAAAICBsAgAAACAgbAIAAAAgIGwCAAAAIDBjpU6VNUFSb49ya3d/bgl5leSNyQ5M8ldSc7u7r8dz3thkp8Zd/3v3f32aRU+z2656Jbc8Kob8sUbv5hDjjskJ553Yo58/pGzLmtbGsb6k19MDkxyT3LI8dtrzBcfT4efeXg+c+ln7rPPh595eG551y255zP3JEl2HL4jJ73hpPtlHNZzzC+1TJINnTu3XHRLrnvZdfsdg2mcm6tZx7x8BkyzjuXWtdK4r2e98zKGy9kKNW4lGx3PaRxLKx3n8/J+z1Mt29W8jPG0jt+19J+Xfb+/bcb3EiC57sXXZc/5e5J7klRy4NcemHu+cM+WOaceqJ+JSVLdvXyHqtOS3Jnkwv2ERWcm+bGMwqKnJHlDdz+lqh6RZFeShSSd5Mok39Tdn1tuewsLC71r16717MtcuOWiW3Ltudfm3rvuHdoOePABecz5j3nAHFT3l6XGep/tMubL7eNq1MGVx17w2E0dh/Uc80vu10FJVaW/9JXPpLW8j7dcdEuuOeea5Mtf3b5vDJJs+Nxczb7Oy2fANOtYbl1Jlh33lf5HZan1PvKFj8w/vf2fZj6Gy5mX93m72Oh4TuNYWuk4n5f327G3+eZljNdSx1pr3sqfv5thM76XAOOg6M179jt/3s+pefl9ME1VdWV3L6yq70ph0XiFJyT5w/2ERb+W5PLuvng8fW2SM/b9dPePLNVvf7Z6WPSBEz4wuuJjkUOOPyRP+8TTZlDR9rW/sd5nO4z5Svu4Gps9Dus55teyX6utf7l1HnL8IUmy4XNzNfs6L58B06xjuXUlS4/rara13/dsfMXcWtd3f5qX93m72Oh4TuNYWs9xPov327G3+eZljNdSx1pr3sqfv5thM76XAMnlOy5f8jNl0jyfU/Py+2Ca1hIWrXgb2iocneSmiend47b9td9HVZ2b5NwkOe6446ZQ0ux88calf9Hsr531W2lMt8OYT2MfNnsc1nPMr6Wm1fZd7/amUctk+7x8BkyzjvWua93z9/OlYp7O6Xl5n7eLjY7nNI6lzf4smxbH3uablzFeSx1rrXkrf/5uhs34XgJkxaAome9zal5+H8zKNB5wXUu09TLt923sPr+7F7p7YefOnVMoaXYOOe6QNbWzfiuN6XYY82nsw2aPw3qO+bXUtNq+K21vGufmatYxL58B06xjuXVt5H3e7/wD17e++9O8vM/bxUbHcxrH0nqO81m83/NUy3Y1L2O8ljrWWvNW/vzdDJvxvQTIfj9TJs3zOTUvvw9mZRph0e4kx05MH5NkzzLt29qJ552YAx781cN6wIMPGB6Sx/QsNdb7bJcxX24fV6MOrk0fh/Uc80vu10GjeteynsXrzEH3bd83BtM4N1ezjnn5DJhmHcuta6VxX896jzr3qLkYw+XMy/u8XWx0PKdxLK10nM/L+z1PtWxX8zLGa6ljrTVv5c/fzbAZ30uA5Khzj1p2/ryfU/Py+2BWphEWXZLkBTXy1CS3d/fNSd6b5Fur6uFV9fAk3zpu29aOfP6Recz5jxk946BG9zNu5QdgzbOvGutkSK6305gvdTwd9aKjltzno150VA48/Cvx/Y7Dd2z6w633V+NK47/UMqe89ZQ89oLHrvvcOfL5R+aUt56y3zGYxrm5mnXMy2fANOtYbl0rjft61nvyr548F2O4nHl5n7eLjY7nNI6llY7zeXm/56mW7Wpexngtday15q38+bsZNuN7CZCc/Ksn56gXHfWVK4wqOfAhB26Zc2pefh/Mymr+GtrFGT2s+ogktyR5dcb/jtzdb6mqSvKmJM9OcleSc7p713jZH0ryyvGqzuvut65U0FZ/wDUAAADAvJnqA667+3krzO8kL9nPvAuSXLCaQgAAAACYvWnchgYAAADANiEsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGAgLAIAAABgICwCAAAAYCAsAgAAAGCwqrCoqp5dVddW1fVV9fIl5h9fVZdV1dVVdXlVHTMx739U1cer6pqqemNV1TR3AAAAAIDpWTEsqqoDk/xKkuckOTXJ86rq1EXdfjHJhd39+CSvSfLa8bLfkuTpSR6f5HFJvjnJ6VOrHgAAAICpWs2VRU9Ocn1339DdX0ryziRnLepzapLLxq//dGJ+J3lQkoOTHJLkoCS3bLRoAAAAADbHasKio5PcNDG9e9w26SNJvnv8+juTHFpVh3f3BzIKj24e/7y3u6/ZWMkAAAAAbJbVhEVLPWOoF03/VJLTq+qqjG4z+1SSu6vq65OckuSYjAKmZ1bVaffZQNW5VbWrqnbt3bt3TTsAAAAAwPSsJizaneTYieljkuyZ7NDde7r7u7r7iUleNW67PaOrjD7Y3Xd2951J/jjJUxdvoLvP7+6F7l7YuXPnOncFAAAAgI1aTVj0N0lOqqpHV9XBSb4/ySWTHarqiKrat65XJLlg/PrGjK442lFVB2V01ZHb0AAAAADm1IphUXffneSlSd6bUdDzru7+eFW9pqqeO+52RpJrq+q6JEcmOW/c/u4k/5jkoxk91+gj3f0H090FAAAAAKaluhc/fmi2FhYWeteuXbMuAwAAAGDbqKoru3thNX1XcxsaAAAAAA8QwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAarCouq6tlVdW1VXV9VL19i/vFVdVlVXV1Vl1fVMRPzjquq91XVNVX1d1V1wvTKBwAAAGCaVgyLqurAJL+S5DlJTk3yvKo6dVG3X0xyYXc/Pslrkrx2Yt6FSV7X3ackeXKSW6dROAAAAADTt5ori56c5PruvqG7v5TknUnOWtTn1CSXjV//6b7541BpR3e/P0m6+87uvmsqlQMAAAAwdasJi45OctPE9O5x26SPJPnu8evvTHJoVR2e5OQkt1XVe6rqqqp63fhKJQAAAADm0GrColqirRdN/1SS06vqqiSnJ/lUkruT7EjyjPH8b05yYpKz77OBqnOraldV7dq7d+/qqwcAAABgqlYTFu1OcuzE9DFJ9kx26O493f1d3f3EJK8at90+Xvaq8S1sdyf5vSRPWryB7j6/uxe6e2Hnzp3r3BUAAAAANmo1YdHfJDmpqh5dVQcn+f4kl0x2qKojqmrful6R5IKJZR9eVfsSoGcm+buNlw0AAADAZlgxLBpfEfTSJO9Nck2Sd3X3x6vqNVX13HG3M5JcW1XXJTkyyXnjZe/J6Ba0y6rqoxnd0vbrU98LAAAAAKaiuhc/fmi2FhYWeteuXbMuAwAAAGDbqKoru3thNX1XcxsaAAAAAA8QwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABsIiAAAAAAbCIgAAAAAGwiIAAAAABtXds67hq1TV3iSfnHUdsIQjknx61kXAFuYcgo1xDsHGOY9gY5xDW9vx3b1zNR3nLiyCeVVVu7p7YdZ1wFblHIKNcQ7BxjmPYGOcQw8cbkMDAAAAYCAsAgAAAGAgLILVO3/WBcAW5xyCjXEOwcY5j2BjnEMPEJ5ZBAAAAMDAlUUAAAAADIRFsIyqOraq/rSqrqmqj1fVy2ZdE2xFVXVgVV1VVX8461pgK6qqw6rq3VX19+PfSU+bdU2wlVTVT4y/y32sqi6uqgfNuiaYd1V1QVXdWlUfm2h7RFW9v6r+Yfzfh8+yRjaPsAiWd3eSn+zuU5I8NclLqurUGdcEW9HLklwz6yJgC3tDkv/T3Y9N8oQ4n2DVquroJP8hyUJ3Py7JgUm+f7ZVwZbwtiTPXtT28iSXdfdJSS4bT7MNCYtgGd19c3f/7fj1HRl9OT96tlXB1lJVxyT5tiS/MetaYCuqqocmOS3JbyZJd3+pu2+bbVWw5exI8jVVtSPJg5PsmXE9MPe6+4okn13UfFaSt49fvz3Jd9yvRXG/ERbBKlXVCUmemORDs60EtpzXJ/lPSe6ddSGwRZ2YZG+St45v5/yNqvraWRcFW0V3fyrJLya5McnNSW7v7vfNtirYso7s7puT0T+sJ/m6GdfDJhEWwSpU1UOS/O8kP97dn591PbBVVNW3J7m1u6+cdS2whe1I8qQkb+7uJyb5Qlz2D6s2fqbKWUkeneSoJF9bVT8w26oA5puwCFZQVQdlFBRd1N3vmXU9sMU8Pclzq+oTSd6Z5JlV9Y7ZlgRbzu4ku7t735Wt784oPAJW518n+b/dvbe7v5zkPUm+ZcY1wVZ1S1U9KknG/711xvWwSYRFsIyqqoyeEXFNd//yrOuBraa7X9Hdx3T3CRk9TPRPutu/5sIadPc/Jbmpqh4zbnpWkr+bYUmw1dyY5KlV9eDxd7tnxUPiYb0uSfLC8esXJvn9GdbCJtox6wJgzj09yQ8m+WhVfXjc9sruvnSGNQHwwPNjSS6qqoOT3JDknBnXA1tGd3+oqt6d5G8z+ku3VyU5f7ZVwfyrqouTnJHkiKraneTVSX4+ybuq6oczCmK/Z3YVspmqu2ddAwAAAABzwm1oAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAMhEUAAAAADIRFAAAAAAyERQAAAAAM/j+ZLtx5JTJoKwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = [1 for i in range(len(both_colours))]\n",
    "plt.plot(both_colours,y,'mo')\n",
    "plt.title(\"Data(Without Label)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STEP 1 (Guess)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimates for the mean\n",
    "red_mean_guess = 1.1\n",
    "blue_mean_guess = 9\n",
    "\n",
    "# estimates for the standard deviation\n",
    "red_std_guess = 2\n",
    "blue_std_guess = 1.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot will look like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIsAAADFCAYAAADUknJhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VNXWx/HvTgIJRaUK0kFEilQjgooiIGLFhoj9FcSCvV0Ve9drQa+iYrl2sHvtigI2igQp0kukFyO9BpLs94+VkBCSMIEkZ8rv8zznOTNnzsysmZzM7Fln77Wd9x4RERERERERERGAuKADEBERERERERGR8KFkkYiIiIiIiIiI7KRkkYiIiIiIiIiI7KRkkYiIiIiIiIiI7KRkkYiIiIiIiIiI7KRkkYiIiIiIiIiI7KRkkYiIiIiIiIiI7KRkkYiIiIiIiIiI7KRkkYiIiIiIiIiI7JQQdAD51ahRwzdq1CjoMEREREREREREosakSZP+8d7XDGXfsEsWNWrUiJSUlKDDEBERERERERGJGs65RaHuq2FoIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFIiIiIiIiIiKyU0jJIudcL+fcHOfcfOfc7QXcfpNzbqZzbppz7kfnXMM8t2U656ZkL5+XZPAiIiIiIiIiIlKyEva0g3MuHngBOAFYCkx0zn3uvZ+ZZ7fJQLL3fotz7irgCaBv9m1bvfftSjhuEREREREREREpBaH0LOoIzPfep3rvtwMjgN55d/Dej/beb8m+Oh6oV7JhioiIiIiIiIhIWQglWVQXWJLn+tLsbYXpD3yT53qScy7FOTfeOXdGQXdwzg3M3iclLS0thJBERERERERERKQ07HEYGuAK2OYL3NG5C4Fk4Lg8mxt475c755oAo5xzf3rvF+zyYN4PA4YBJCcnF/jYIiIiIiIiIiJS+kLpWbQUqJ/nej1gef6dnHM9gMHA6d779Jzt3vvl2etUYAzQfh/iFRERERERERGRUhRKsmgicIhzrrFzrjxwHrDLrGbOufbAy1ii6O8826s65xKzL9cAjgbyFsYWEREREREREZEwssdhaN77DOfcNcB3QDzwuvd+hnPuASDFe/858G+gMvChcw5gsff+dKAF8LJzLgtLTD2WbxY1EREREREREREJI8778CoRlJyc7FNSUoIOQ0REREREREQkajjnJnnvk0PZN5RhaCIiIiIiIiIiEiNCmQ1NREREREREIkxGBqxbB2vW2LJ2rS3OQWIiJCUVva5SBcqXD/pViEgQlCwSERERERGJMOnpMHcuzJwJM2bY5X/+yU0KrVkDGzbs23PExUHdutCoETRunLvOuVyvHsTHl8CLEZGwo2SRiIiIiIhImMpJCs2YYUtOcmj+fMjMtH3i4iyBU6uWJXcOOwyqVYOqVW2d93KVKtazaNs2e+y86/yX09Jg4UL46y8YNQqWLYO8JW8TEqB+fTjkEDjySDjqKFtXrRrIWyUiJUjJIhERERERkTCxZQv8+qslZ0aNgj/+yE0KxcdD06bQsiX06WPrVq2gWTMbOlbatm+HxYtzE0h//WWXZ8yAhx+GrCzbr0UL6NzZlk6dLM44VcsViSiaDU1ERERERCQg27fD779bYujHH2HcONixw3rtdOoExx4LrVtbwuXQQ62WUEkYONDWw4aVzONt3AgTJ1r8OcuaNXbb/vtbj6Ojj4aTT4bkZOvdJCJlqzizoSlZJCIiIiIiUoZmzICvv7YE0S+/wObNljzp0AG6dbPlmGOgcuXSi6FrV1uPGVM6j+89zJu3a/Jo+nTrfVS3LvTuDWecAccdpyLaImWlOMkiDUMTEREREREpZQsWwIgRtkyfbttatIBLL4Xu3S1pUq1aoCGWKOdseFyzZnDJJbZt9Wr46iv43//gjTdg6FA44AA45RRLHPXqBfvtF2jYIpJNySIREREREZFSsHw5vP++JYh+/922HX00PP88nHkm1KkTbHxlrXp1uPhiW7ZuhR9+gM8+g88/h/fesx5G3btb4ujss21/EQmGhqGJiIiIiIiUkNWr4aOPLEH00082HKtDBzjvPOjbFxo0CDpCU9rD0IojMxPGjrXE0aefWuHs8uUtYXT55RarahyJ7DsNQxMRERERESkjWVnw3Xc2rOrbbyEjA5o3h/vuswTRoYcGHeHu2rULOoJc8fHQpYstTz4JU6fC66/D22/D8OE2A9yAATZkr1atoKMViQ3qWSQiIiIiIrIX1q+32jsvvGDFnGvXtiFW/fpB27bqDbOvtm6Fjz+GV16Bn3+2GeJOP916G51wgiWZRCR0mg1NRERERESklMycaXWH3nrLZjLr3BmuvdaGTWlmr9IxZw68+qol5/75x4bz9e9vPY5irfaTyN4qTrIorrSDERERERERiXSZmTaLV48e0KqVDZPq0wdSUqzeTr9+kZUouvBCWyLFoYfCv/8Ny5bBBx/Y9XvvhUaNLGE0Z07QEYpEFyWLRERERERECrFhgyUpmja1WbrmzoVHHoElS+C//4XDDw86wr2zdKktkaZ8eUvSff89LFgAV1wB774LLVpYz66JE4OOUCQ6KFkkIiIiIiKSz/r18OCD1nPlttts/fHHkJoKd9wBNWsGHaE0aQL/+Q8sWgSDB8OoUdCxI3TvDiNH2kx0IrJ3lCwSERERERHJtm6dzWLWsCHcc4/N0JWSAqNHw1lnWZFlCS8HHmiJvcWL4amnYPZs6NkTkpPhww9tCKGIFI+SRSIiIiIiEvPWrLHkUMOGcP/90K0b/PGH1SmK1KFmsWa//eCmm6z312uvwaZNcO650Ly5zai2Y0fQEYpEjpCSRc65Xs65Oc65+c652wu4/Sbn3Ezn3DTn3I/OuYZ5brvEOTcve7mkJIMXERERERHZF6tX2xCmRo2sd8oJJ8CUKfDJJ9C+fdDRlZ7OnW2JRomJcNllNmvdxx9D1aowcCC0bAnvvw9ZWUFHKBL+nN/DQE7nXDwwFzgBWApMBPp572fm2ed4YIL3fotz7iqgq/e+r3OuGpACJAMemAQc7r1fW9jzJScn+5SUlH18WSIiIiIiIoVbvRqefBKefx42b7aiyXffDYcdFnRkUtK8h6+/tlpTf/5pScBHH7Whas4FHZ1I2XHOTfLeJ4eybygjbjsC8733qdkPPgLoDexMFnnvR+fZfzyQMwnjicBI7/2a7PuOBHoBw0MJTkQkpmRlwbZtsHWrLfkvg50qy7uUL7/7dbV6RERECrVtmyWIHnrIZjrr2xfuugtatQo6MiktzsEpp0CvXjBihCUFe/WCrl3hscfgyCODjlAk/ISSLKoLLMlzfSlQ1L9Tf+CbIu5bN/8dnHMDgYEADRo0CCEkEZEw5z1s3AgrVxa9rFplpzO3boXt20vmuStVskqPBx4ItWrlXs5/vXZtqF5dySUREYkJWVk2BOmOO2z2rJNPhscfj92eRGefbeuPPw42jrIUHw8XXGC9yIYNs2GHnTrBmWfCww9DixZBRygSPkJJFhX0K6LAsWvOuQuxIWfHFee+3vthwDCwYWghxCQiEh7Wr4c5c3ZfFiyALVt23z8hwRI2tWtDnTrWD7pyZahQIXdJStr9clKS3T89ffdl+/Zdr2/YAGlplohauBB+/92uFzQVyAEHQLNmcMghtuS9XKVKqb51IiIiZeXnn+GWW2DiRGjXzoofd+8edFTBWr066AiCU748XHMNXHopDBkCTzxhhcwvucRmwlP/BZHQkkVLgfp5rtcDluffyTnXAxgMHOe9T89z36757jtmbwIVEQnUunU2b+7UqbsmhVatyt0nPh6aNIFDD4UePSwZVLv2rku1ahAXwESUWVk2zcvff+cuy5bB/Pkwdy789hsMH249onLUrJmbQGrb1uafbdfOklsiIiIRYM4cuP12+OwzqFsX3nwTLrwwmK9iCT+VK9sQxCuvtBpGL7wA770HN99sPdDU5JFYFkqB6wSswHV3YBlW4Pp87/2MPPu0Bz4Cennv5+XZXg0rat0he9MfWIHrNYU9nwpci0jgtm6FyZOtR87EibbMm5d7e40alhDKvzRpYqeqItW2bTbX7Lx5lkCaN8+W2bNtyBxY67p5c0scJSfbXMLt2kHFisHGLiIikkdaGtx/P7z0knXQveMOuOEGfV3l1bWrrceMCTKK8LJkic2M9/bbds7viSfg/PM1Yl+iR3EKXO8xWZT9gCcDQ4B44HXv/cPOuQeAFO/95865H4DWwIrsuyz23p+efd/LgDuztz/svf9vUc+lZJGIlCnvYcYMGDfOkkK//w7Tp+cO2apbF444wpaOHS0xUqNG2ccZdItuxQqYNMl6V02aZO9VTq+quDibizY5GY46Co4/Hg4+WC0rEREpc9u327Cihx6y0eADB8K999oIcNlV0E2LcDZuHFx7rTV5jj4annsOOnTY8/1Ewl2JJ4vKkpJFIlLqUlPhxx9tGTXKTj8CVK2amxjKWerUCTbWHA8+aOu77w42jhzew/LluQmklBRLIP3zj91er54ljY4/Hrp1g4YNg41XRESi3o8/Wh2a2bPh1FOtV4gKFhcu3JoW4SYrC954w3qlpaVB//5WBPvAA4OOTGTvKVkkIpLXypUwenRugmjhQtt+0EFW3bJbN+jSRb1h9pX3Vhxi1Ch7v8eMyU0eNW6cmzw6/njrsSUiIlICli61GjMffGAjwp97zqZJFykJ69fDAw/YcVWpkhXAHjQIypULOjKR4lOySERi244dNu3JF1/ADz/YMDOw2b1yerp07261d5QcKj1ZWfbejx5ty08/wdq1dlvz5nD66XDaadC5sxUHFxERKYbt2+HZZ602UWam9QC57bbcCURFStLs2Vb36rvvrMfakCHQs2fQUYkUj5JFIhJ7NmyAb7+1eU+//tpmL0tKsh5D3bvb0r595CYlTjrJ1t98E2wc+yIzE6ZNs8TRt99az6MdO6wG1MknW/KoZ0/Yb7+gIxURkTA3apQNOZs1y847DBlivYokdNHQtChr3sOXX8KNN8KCBXD22ZawVIdpiRRKFolIbFi2DD7/3BJEo0fbKcYaNazV2Ls3nHBC9Ex7Eo1VKNevt9NzX3wBX31lvY7Kl7feXzm9jurXDzpKEREJI8uWwS23wIgRNsL52Wft60KKLxqbFmUlPR2eesrqPpUrZ7WMrr46cs9JSuxQskhEotfcufDhh/DZZ1ZUGaBpU0sO9e5ts3FF4zd1tLfoMjLgt98scfT55zBvnm1v3x769oXzzlORbBGRGJaRYYmh++6zTqm33w7/+hdUqBB0ZJEr2psWZWHBAksSff+9TQo7bJg1XUTCVXGSRXGlHYyIyD5bvhyeecZmJzv0ULjrLksIPfKI1cSZOxeefNKGnEVjoigWJCTAccfZ33HuXCsM8MQTkJhovwgaNbK5a59/HlatCjpaEREpQ3/8AR07Wo+i446zr/777lOiSIJ38ME2sn74cFiyxBJGN90EmzYFHZnIvlOySETC09q18NprVmuoXj375vXe+vwuXQrjx1sly5YtVaQ6Gh16KNx6K4wbB6mplhjcuBGuvRbq1LHaRq+/brWpREQkKm3ebAmiI46AFSusY/EXX9gPdJFw4Zx1gJ41Cy6/3M5vtmxpHaVFIpmSRSISPrZutZbgmWdC7dowYICdprnnHutpkpJiSaNYrCJ46qm2xKLGjS0xOG0aTJ8Od95pCaT+/aFWLTjjDHj/fTt+REQkKnz/PbRubeeIBgywH+LnnKPzQyUplpsWpaFqVXjpJRtVf8ABVh3hzDOtKSsSiVSzSESC5b31HnntNUsUbdwIBx1kp2jOPx8OP1wtQ9md95Y8HD7cEkXLl0OVKnbMXHYZdOig40ZEJAKlpdl5oXfegWbN4JVX4Nhjg45KpHh27ICnn4b777cKCQ8/DIMGqVqCBE8FrkUk/KWlwdtvw6uv2unCSpXg3HPhwgutIIG+TSVUmZnw0082LO3jj2HbNmjb1pJGF1wA1asHHaGIiOyB95YguvFG2LDBytXdeSckJQUdmcje++svK4D97bc2B8urr0KLFkFHJbFMBa5FJDxlZtpU6X362FCym2+23iCvvmrFCF5/Hbp1U6KoIF275k5bIruKj7fj5p137DgaOtQKZl9/vdU36tvXxjNkZgYdqYiIFCA1FU48ES6+2HoTTZ4MDzygRFFpU9Oi9DVuDF9/bedHZ8+Gdu2sDOOOHUFHJrJnShaJSOlbtMimLWnSBHr1gtGj4ZprrP7M2LFWe2a//YKOUqJBlSpw1VU2RG3KFLjySvjhB/sV0qQJ3HsvLF4cdJQiIoLl8IcMgcMOs3krnn8efv0VWrUKOjKRkuOcdZyfOdPqGA0ebLP7TZ4cdGQiRVOySERKR2YmfPklnHSSnVZ54AFo3hw++ACWLbOB3GoNSmlq2xaefdaOt/fft+PvwQfteDzjDBg5ErKygo5SRCQmzZ4NXbrYsLNu3eyH9KBBEKdfJxKlatWyZvAnn8DKlTbL3+DBNnpeJBzp41hEStaaNfDkk3DIIXDaaTB1Ktx9tw3azhmClpgYdJQSS5KSrB7Wd9/Zcfivf1mPtp49rXDAkCGwbl3QUYqIxISMDHjiCRuOM3u2Dc/54guoVy/oyETKxplnWnL0ootsSFr79tYsEQk3ShaJSMmYPNmGk9WtC7feCvXr2+mTRYtsKoiGDYOOUMSOw0cesXls334bqlWz09p168LAgZbcFBGRUjF9uhX5/de/4OST7QfzhRdq8kqJPVWrwn//a4Wvt2yBY46xUoubNwcdmUguJYtEZO9t325Tlx99tE1VPmKEVaecOtVmp+rTB8qVCzrK6HDuubZIyUhMtF8o48bBpEnQr58VyG7Xzlpsw4fb8S0iIvtsxw546CFrKvz1l40M/vhjqF076Mhim5oWwTvxREuiDhoEzz0HrVtbaU+RcOC890HHsIvk5GSfkpISdBgiUpSVK23GqWHDYNUqaNrUvuUuvdQKDItEorVr7TTf0KGwYIH9ihk0yIpk16gRdHQiIhFpyhS47DLrgNy3L/znP1CzZtBRiYSfX3+1/5V58+Dqq+Hxx6Fy5aCjkmjjnJvkvU8OZV/1LBKR0P35p32LNWxopwiTk+Gbb2DOHLjhBiWKStOWLbZI6alaFW66CebOteO6XTurt1W/viWMZs8OOkIRkYixfbtNQHnEEbB8uRX1HTFCiaJwoqZFeDnmGEuu3ngjvPgitGkDY8YEHZXEspCSRc65Xs65Oc65+c652wu4/Vjn3B/OuQzn3Dn5bst0zk3JXj4vqcBFpIx4b4WBe/a0b60RI2DAAEsQffkl9OqlqUvKwskn2yKlLy7OjutvvoEZM6wC5RtvWDHsU06BH36w/wsRESnQ5MmWJHrgARvlO3OmFfWV8KKmRfipWNEmDP75Z4iPh+OPh2uugU2bgo5MYtEef+E55+KBF4CTgJZAP+dcy3y7LQYuBd4r4CG2eu/bZS+n72O8IlJWtm2D116zwdO9etmA6kcegaVL4YUXbLYzkWjXsqUNt1yyxH71pKTACSdA27Y2ZC09PegIRUTCxo4dNqdFx47w99/w+efw1ls2l4CIhO6YY6wE6A032Oh49TKSIITSHaAjMN97n+q93w6MAHrn3cF7v9B7Pw3IKoUYRaQspaXZj+KGDa0HUUICvPkmLFwId9yhFp/Eppo1bUja4sWWJILcIZkPPgirVwcbn4hIwKZOtSTRfffBeedZx8zTTgs6KpHIVbEiPPOMzRmT08vo2ms1Y5qUnVCSRXWBJXmuL83eFqok51yKc268c+6MgnZwzg3M3iclLS2tGA8tIiVmwQK46ipo0CC3yMCPP1pf8osvhvLlg45QJHiJiVbIfepUGDkSDj8c7rnH/m+uvx4WLQo6QhGRMrVjh+XMc2oTffopvP22zi2JlJQuXazZcf311rm/TRtLIImUtlCSRa6AbcUp1tAgu9r2+cAQ59zBuz2Y98O898ne++SaqnonUrYmT7ZTgM2aweuv23TiM2daPaJu3cAV9BEgEuOcgx494KuvbIhmnz7WT/zgg+1/aNq0oCMUESl106dDp06WMz/7bOtNdEaBp4ZFZF9UrAhDhthQNOega1e47jr1MpLSFUqyaClQP8/1esDyUJ/Ae788e50KjAHaFyM+ESkN3luvoZ49oUMHK+R766021OyVV6yQr4SXSy+1RcJPq1ZWADs11U77/e9/VtPopJOsVadi2CISZTIy4NFHrXPlkiXw0UcwfDjUqBF0ZFIcalpEnmOPtV5G114L//mPTdz6669BRyXRyvk9NGKdcwnAXKA7sAyYCJzvvZ9RwL5vAF967z/Kvl4V2OK9T3fO1QDGAb299zMLe77k5GSfkpKyly9HRIqUmWlz1z7+OEyaBLVrW+W8K6+EAw4IOjqR6LB2rc15++yzVuH1iCPgX/+y0+3x8UFHJyKyT2bOtATDxInWqfKFF6ysm4iUrTFjrHziwoVw443w0ENQoULQUUm4c85Nyh75tUd77Fnkvc8ArgG+A2YBH3jvZzjnHnDOnZ79hEc455YCfYCXnXM5iaQWQIpzbiowGnisqESRiJSSbdtsRqfmzeHcc2HDBrv+11/2I1aJovD3zz+2SPirWhXuvNNaby+9BGvWwDnnWI+9116D7duDjlBEpNgyM+GJJ6xDcmoqvP8+fPCBEkWRTE2LyNa1q416v/JKePpp62U0fnzQUUk02WPPorKmnkUiJWjTJvux+tRTsHKlejhEsq5dba15UyNPZqZVfH3sMevRV68e3HKLzTZYqVLQ0YmI7NGcOdabaPx4OPNM6zxZq1bQUcm+UtMievzwA/TvD0uXWhPj/vshKSnoqCQclWjPIhGJQGvW2LdEgwZWi6hVK6tRNGGCVaBUokik7MTHW8+iiRPhu++sCPYNN0CjRvDww7BuXdARiogUKDMzt8fC3Lnw3nvw8cdKFImEmx494M8/LWH0xBNWT2zixKCjkkinZJFINFmxAm67DRo2hPvusyp448fb6QbNbCYSLOesqPyYMVaN8sgj4a677P/1zjutvpGISJiYNw+OOw5uvtk+umbMgH791JQQCVf7729VJr79Ftavh86drZmRnh50ZBKplCwSiQYLF8LVV0Pjxjbk7PTTbRDzZ5/ZD1IRCS9HHw1ffgmTJ9usaY89Zkmj666DxYuDjk5EYlhWltXnb9vWEkRvv23Nidq1g45MREJx4okwfTpcdJF1YE5OtlHwIsWlZJFIJJs9Gy65BJo2hVdfhYsvtsIC774LrVsHHZ2I7Em7djBihP0vn3++FQI5+GCb3mTevKCjE5EYs2ABHH+8jZTt1s2SRRdeqN5EIpGmShX473/tvNSaNXbu+O671ctIikfJIpFINHWqzWrWsiV8+CFce61NTTJsmCWOJPpcdZUtEp2aNbOZ0hYssL/z8OE2e+H559vpQRGRUpSVBS+8AG3awJQp9iPziy+gTp2gI5PSpKZF9DvlFGtGXHghPPSQehlJ8Wg2NJFI8vvv9kn/xRew336WJLrhBs1bKxJtVq2yqrJDh9qshmeeCYMHW8VKEZEStGCBFcX96ScbvvLqqzZpo4hEl6++goEDrYlxxx1WzygxMeiopKxpNjSRaPPzz1Zd8sgj4bff4IEHYNEiG4isRFFsWLLEFokNtWrB449bPbJ77oHRo+104Mknw9ixQUcnIlEgpzZRTm+i116Db75RoiiWqGkRW9TLSIpLySKRcOU9jBxpM5odd5wNPXviCfvxePfdULVq0BFKWbroIlsktlSvDvffb//3jzxi8+AefbQVExk1yj4nRESKKWemsxtusBpFM2ZYqTTVJootalrEnqpV4Y03VMtIQqNkkUi48d6GmXXqZL2JUlPhuefsx+Ktt9rwMxGJLQccYH3GFy604WmzZ0P37pY4+vprJY1EJCSZmTZpaps21sPgzTetyVG3btCRiUhZUi8jCYWSRSLhIivLilW3bw+nnw5paVawesECq01UoULQEYpI0CpVghtvtCTyCy/AsmXW4ktOhk8/tc8REZECzJ4NxxwDt9xi56JmzrRJVNWbSCQ2FdTL6PbbYevWoCOTcKFkkUjQMjLgnXfgsMNshrNt2+yTe84cuPxyVZ4Tkd0lJcHVV9tYktdeg/Xr4ayzoG1bGDHCug+IiGDNjMcfh3btYO5cePdd+OwzOOigoCMTkXBwyik2FPXSS3M/K379NeioJBwoWSQSlO3b7Ude8+Y2YDwhAd5/3z6tL7kEypULOkIRCXfly1uhkdmzLemcmQn9+kHLlja+ZMeOoCMUkQD9+SccdZT1FjjlFOtNdP756k0kIruqUsVmQhw50n6idOkC11wDGzcGHZkESckikbK2bZsNH2naFAYMsD6gn31mU5Gcey7ExwcdoYSjm2+2RaQgCQlwwQVWgOCjj6BiRTtF2KwZvPyyKleKxJj0dCta26ED/PWXnYv66CObaFEkh5oWkl+PHpZkvv56GDrUBj58913QUUlQnA+zopjJyck+JSUl6DBESt7mzfaj7d//hpUr7VTf3XfDiSfqFJ+IlCzv4auv4MEH4fffrXrtrbfa0NaKFYOOTkRK0W+/2bmo2bOt4/LTT0ONGkFHJSKRZuxY6N/fPksuucQ+S6pVCzoq2VfOuUne++RQ9lXPIpHStm4dPPwwNGpkp29atYLRo20wcK9eShRJaObMsUUkFM7BqafC+PHw/fdw8ME2R3ajRvDYY7BhQ9ARikgJ27ABBg2yItZbt8K338JbbylRJIVT00KKctRRMHkyDB5sI91btoSPPw46KilLShaJlJa0NPt0bdgQ7rrLphgYOxZ++AG6dlWSSIrniitsESkO5+CEE+Cnn+CXX+Dww+GOO+xz6Z57YPXqoCMUkRLw5Zd2LurFFy0vPH26dVwWKYqaFrInSUnw0EOQkgJ16sA558DZZ8Py5UFHJmVBySKRkrZsmU1t3agRPPq4K9j6AAAgAElEQVSotdYmT7aWXOfOQUcnIrHqmGPgm29g4kQ4/ngbotawoQ1PW7ky6OhEZC/8/bfVtD/tNCtQO24cPPMMVK4cdGQiEk3atbNR7Y8+Cl9/DS1aWE2jrKygI5PSpGSRSElJTbXTM02awH/+Y6n3mTPhgw/sE1ZEJBwkJ8Mnn1gFy969rQhBo0Y27cnixUFHJyIh8N4mPGzRwv6dH3gAJk2yTswiIqUhIcFmVvzzT+jY0Ya9Hn20XZfopGSRyL6aOdMqSDZrBm+8YdNYz5tnrbjmzYOOTkSkYIcdBu++awUrLroIhg2z2kb/939WzVJEwtLcuTa69NJLLVk0ZYrNl1G+fNCRiUgsaNrUyiG+/TbMn2+zLt5+O2zZEnRkUtJCShY553o55+Y45+Y7524v4PZjnXN/OOcynHPn5LvtEufcvOzlkpIKXCRwEybAGWdYkYBPPrE5Jv/6ywoGNG4cdHQiIqFp2hReeQUWLICrrrI5tlu2tKIEmp1UJGxs2wb33gutW9u/5tCh8PPPljASESlLzsGFF9q5pYsvhscft3NQ330XdGRSkpz3vugdnIsH5gInAEuBiUA/7/3MPPs0AvYHbgE+995/lL29GpACJAMemAQc7r1fW9jzJScn+xQ1TiVceQ8jR9psQqNHQ9WqNnTjuus03YiUrh9+sHWPHsHGIdEvLQ2efRaefx7Wr7cuDHfcocL8IgH6/nu4+mrL6V5wATz5JNSuHXRUEunUtJCS8tNPVo1jzhyro/bMM1CrVtBRSUGcc5O898mh7BtKz6KOwHzvfar3fjswAuiddwfv/ULv/TQgf4mrE4GR3vs12QmikUCvUAITCSuZmfDRR3DEEVawes4ca6ktWmSFApQoktLWo4dac1I2ata0qU8WL7ZThdOmQbduVqD/f/9TNUuRMrR8OZx3njU94uPtx/077yhRJCVDTQspKccdB1Onwn33wccfWyWOV15RkyHShZIsqgssyXN9afa2UIR0X+fcQOdcinMuJS0tLcSHFikD27fDa6/ZkIw+fWDDBvvkS02Fm2+G/fYLOkKJFVOm2CJSVvbfH267DRYutOG1f/9tQ29bt7ZCBTt2BB2hSNTKzLS5Mpo3h88+s/NS06ZB9+5BRybRRE0LKUmJiTZUdupUaNsWBg6Eo46y4vsSmUJJFhXU57zosWvFvK/3fpj3Ptl7n1yzZs0QH1qkFG3caP0nmzSBAQOgUiWb1WzWLLuemBh0hBJrbrjBFpGylpQEV15pVXXffRfi4qxAQdOmNlxt06agIxSJKikpNqvZdddZh77p062AtZoeUtLUtJDS0Ly5Vet4800733TEEVYScc2aoCOT4golWbQUqJ/nej1geYiPvy/3FSl7K1ZYbY769eGmm+CQQ6xS26RJ1rMoPj7oCEVEgpGQAOefb90bvvgCGja0XxkNGsBdd8GqVUFHKBLRVq+2qag7drThZ++/D99+a3lZEZFI4pydV5ozx+YAeuUVmzh62DDrOSmRIZRk0UTgEOdcY+dceeA84PMQH/87oKdzrqpzrirQM3ubSHiZNQv694dGjeCJJ6BnT/j9d0uL9+ypoq4iIjmcg1NPtWmYxo2D44+HRx6x5NEVV1gPJBEJWUYGvPCCnZ96+WW49lqbYejcc9X8EJHIdsABNlhj8mSbQPqKK6BTJ/uZJeFvj8ki730GcA2W5JkFfOC9n+Gce8A5dzqAc+4I59xSoA/wsnNuRvZ91wAPYgmnicAD2dtEguc9/PILnHaa1SQaPtyGmM2da0POjjgi6AhFRMJbp05WyXL2bLj0Uutz3rw5nHUWjB8fdHQiYW/0aOjQwSZWbd/e6sc8+6yVDBMRiRatW8OYMfDee7BsmTUfLr/cJmCV8OW8D7X8UNlITk72KSkpQYch0Swz06pF/vvfMGGCzWR2zTXW91uzmkm46trV1mPGBBmFSNFWrbKqvEOHwtq1cMwxcOut1hMpLpTOzCKxYdEiuOUWm2i1USN46ik480z1JJKypaaFBGHjRivaP2SIzRX00EPW40jVPsqGc26S9z45pH2VLJKYsWkTvPGGnbKbP9+KV998s50Nr1gx6OhEijZ2rK2POirYOERCsWmTzST59NOweLEVXbnhBrjkEqhcOejoRAKzZQs8/riNeI+LszKJN98MFSoEHZnEIjUtJEgzZ9qw21GjrOfRU0/BCScEHVX0U7JIJK9Fi+D5562y2vr1NsXILbfYKTylsEVESs+OHfDJJ1awYMIEqFLF5tK95hqbSEAkRngPH35ozY8lS+C88yxhpH8DEYll3lsz4bbbIDUVTj7ZBn+0bBl0ZNGrOMki9QmX6DVunFWHPPhg+6HSq5dtGz8ezjlHiSKJLGPH5p4CFIkU5cpB3772uTt2LPToAU8+CY0bQ79+qnApMeH33224T9++UL261YYfPlyJIgmemhYSNOfg7LOtl9GTT8Jvv0GbNnD11fD330FHJ+pZJNFlxw4rtjpkyK5nsQcNsumdRSKVCgtItFi40OoavfKKFS44+mi48UY44wwl8SWqLFgAd95pc2YceKDV6BgwQIe5hA81LSTc/PMP3H8/vPiiVQkZPBiuvx6SkoKOLHqoZ5HEnjVrrAhAkyZ2tnrNGht6tmSJbVeiSEQkPORU81261Hp9LltmvT2bNrW+56tXBx2hyD5JS4PrrrOJAb/8Eu65x0olqoCriEjRatSw80nTp1sy8/bb7bN0xAgbsiZlS8kiiWx//AGXXQZ169qnSbNm8MUXNo3zoEEqpCoiEq7239+KXs+fb1NC1a9vRQvq1bPuF5MnBx2hSLFs2QIPP2yj34cOhf797fC+/36b8UdERELTvDl8/jn8+CNUrWp9ATp3hl9/DTqy2KJkkUSe9HR45x37xDj8cHj/fZthZ9o0+0TRFM0iIpEjPt4KFvz8M0ydChdfbAVdOnSAY46x04nbtwcdpUihMjLg1VfhkEPgrruge3c7K/7SS3DQQUFHJyISubp1g5QUeP11m1y1Sxc45RSdTyor+kUtkWPJEhu4Wr8+XHSRDTV79llYvtxaZK1bBx2hiIjsizZt4OWXbYja00/DypV2OrFhQ7jvPvu8FwkT3ltn5rZt4fLLbcT7L7/Ap5/aWXEREdl38fHwf/8H8+bBY4/ZfEUdOkCfPjBrVtDRRTcVuJbw5j2MGgUvvAD/+59tO+00G2LWvbt6EEnsmDLF1u3aBRuHSFnKyoLvvrMadN98k9sL6cor4bjjbBoVkTLmvR2W991nc2kccgg8+iicdZYOSYksalpIJFq3zs4nPfOMDf+96CK4916baFX2rDgFrpUskvCUlgZvvmmz5cyda9XOBgywHwgNGwYdnYiIlLX58216lNdft5Zis2Y22+Ull9h3hEgp8x5GjrQk0bhx1pNo8GA7412uXNDRiYjElrQ0m8fo+eft3NLll9tQYA3/LZpmQ5PIlJVlrbBzz7WC1bfeCjVrWtJoyRI7badEkcSqH36wRSRWNW1qs6gtW2bfCzVrwi232PdFv34werSmSpFS4b19/HbpAieeaE2SF1+0c1kDBypRJJFLTQuJZDVrwpNPwoIFNqHAsGE2wcBtt2li1ZKinkUSvBUr4L//teqQf/0F1arZmeIBA6Bly6CjEwkPXbvaesyYIKMQCS8zZlgP1LfegrVrbTzQ5Zfbd8iBBwYdnUSB0aNteMMvv1he8s477UdJYmLQkYnsOzUtJJqkplrPz3fegUqV4Kqr4KaboHbtoCMLL+pZJOEvMxO+/hrOPNMKVg8eDI0awXvv2Vnjp59WokhERIrWqhUMGWLfG2+/DbVq2SnFevWgb1/rrZqZGXSUEoF++sl+SHfrZmet//MfGwl59dVKFImIhKMmTezc0fTpVuL2qafs5+WgQbBoUdDRRSYli6RszZljp+UaNbJ5D3/7DW6+2fpyjxplQwmSkoKOUkREIkmFCnDhhdb9Y8YMaxmOHAk9e9r3zZ132vePSBGysmx2s2OPtUTRnDk26eqCBXDNNWqeiIhEgpYtrf/BnDlw8cXWAblpU7j0Upg9O+joIouSRVL61q61wf2dOtlcso8/btMjf/CBTY/8+OM2dEBERGRftWxpU6QsX27fM23a2PdM8+Zw1FFW1GDduqCjlDCyfTu88Qa0bg2nn25noIcMsSEN112nJJGISCRq2tS+8lNTLeH/wQfWROjTByZPDjq6yKBkkZSOHTvgyy/tv7F2beu3vWWLVSFbtgy++spuK18+6EhFRCQaJSXZ98xXX9mJiX//G9avhyuusO+lfv1s/nMNU4tZGzZYs6RxY5vRLCHBal3Mnw/XX28d1kREJLLVq2fnkBYtso7G338PHTrAySfDr78GHV14U4FrKTnew7RpNkvNu+/C33/bdMYXXGDFRtu1A+eCjlIkMuUMoTn00GDjEIlk3sOkSdaN5L33rOdrnTo2hK1fP2jbVt9TMWD5cnjuOev0vGGD1SW67TYbtag/v8QSNS0kFq1fD0OHWgIpLQ06drQTBOecExv9GIpT4FrJItl38+bBiBEwfDjMmmVzyJ52miWITjpJc8qKiEj4SU+3AjVvvAHffms9jFq0sKRRv37Wf12iysyZNn/G229DRob9MLj1VkgOqcksIiLRZMsWm5D7ueesfO5BB1nJw4EDoWbNoKMrPSWeLHLO9QKeBeKBV733j+W7PRF4CzgcWA309d4vdM41AmYBOVUlx3vvryzquZQsihCLF9vAz+HD4Y8/7FRcly5w3nlw7rlQvXrQEYpEly++sPVppwUbh0g0+ucf+Ogj+077+WfblpxsSaO+fW3OdIlIO3bA//4HL7xg04MnJcFll9l0ygcfHHR0IsFS00LEJjf47jurVff99zbj5YUXWm+j1q2Djq7klWiyyDkXD8wFTgCWAhOBft77mXn2uRpo472/0jl3HnCm975vdrLoS+/9YaEGr2RRGFu1Cj780HoR/fabbevY0RJEffrYgFARKR1du9p6zJggoxCJfkuXwvvv2zC1nJMhxx5riaNzztHJkAixYoXNgPPyyzbsrGFDuPJK6N8/us8YixSHmhYiu5o1y3oavfWW9Tw6/nhLGp16KsTHBx1dyShOsiiUAtcdgfne+1Tv/XZgBNA73z69gTezL38EdHdOo76jwsqV1tI64QSr63DttTbA/+GHrQLkhAlw441KFImISHSoVw9uvtlqG82ZA/fea9+FV15phbF79oSXXrJtEla8t45hfftCgwb2p2vdGj7/HBYsgNtvV6JIREQK16KF1bNbssQmUp0/H844A5o1sz4TsSaUZFFdYEme60uztxW4j/c+A1gP5Jx6a+ycm+yc+8k516WgJ3DODXTOpTjnUtLS0or1AqQUpKbCU0/B0UdbgujKK3PLx0+fbkWs77xT/bdFRCS6NWtmGYdZs6yX0S23wMKFcNVV9v14zDFWIXPhwqAjjWmbNlnjvk0bOO44G0Zw3XVWUvHbb22ITbScERYRkdJXrZpNfJCaapVXDjooNidPTQhhn4J6COUfu1bYPiuABt771c65w4HPnHOtvPcbdtnR+2HAMLBhaCHEJCXJe/jzT/j0U1umTrXt7dvD/ffDWWdBy5aaIkRERGKTc/ad2L49PPKIVUr++GP45BMrfnPTTTYP79ln23dm8+ZBRxz1srLgl19sAtYPP7SEUfv28OqrNmKwYsWgIxQRkUiXkGDVVvr0sZ/MsSaUZNFSoH6e6/WA5YXss9Q5lwAcAKzxVhApHcB7P8k5twBoBqgoUdAyMmD8eKv6+Omn1j/bOetN9PTT1t+uceOgoxQREQkvzkGrVrbcc499f37yiS2DB9vSooV1Zzn1VOjc2VqbUiJSU62WxFtvwV9/wX772bCzAQPgyCN1XktEREpHLH6/hFLgOgErcN0dWIYVuD7fez8jzz6DgNZ5Clyf5b0/1zlXE0saZTrnmgC/ZO+3prDnU4HrUpSWZv2xv/7aSr6vXWvT2nfvDmeeCb17Q61aQUcpIgVZkj0auH79ovcTkeAsW5bbS/fnn+3ETNWqcOKJljjq1UsFsvfCxo02Wd0bb9jb6pw1XS691Jov6kUksnfUtBCJPSU6G1r2A54MDAHigde99w875x4AUrz3nzvnkoC3gfbAGuA8732qc+5s4AEgA8gE7vXef1HUcylZVIKysmDyZPjqK0sQ/f679Z+rVQtOPtmWE06AAw4IOlIREZHosn49jByZ+x38998QFwedOsEpp1jyqHXr2DxVGYLMTJuh6c03bcTfli1WQuqSS+Cii/TjVkREZG+UeLKoLClZtI/++QdGj4ZvvrFl5UpriHbsaMmhU06xQf1xodQ2F5Gw8f77tu7bN9g4RKT4srIgJcUSR199ZTOtgc281qsX9OhhXWVq1Ag2zoBlZFjPoQ8/tFF9f/9t57P69rVeRJ06KbcmUpLUtBCJPUoWxZItW+DXX+GHH2yZMsV6D1WpYt3eTznF1gceGHSkIrIvuna19ZgxQUYhIiVhxQrrbfTVVzBqlPVCAjuZc8IJljw65hioUCHYOMtARoZ9rH34oY3eS0uzYWWnngrnnGPrGHgbRAKhpoVI7ClOskgVFyNNRoadnfzxR0sOjR0L27db7aGjjoIHHrCzk0ccoYKaIiIi4eigg6B/f1syMqyn0ciR9r3+zDPwxBOQmGgJox49LIHUrl3UzP++Y4flyD76yBJEq1dDpUpWE/ycc+Ckk1SHSEREJGjqWRTuduyAP/6w+WF/+cVS/xs22G3t21tiKOcMZKVKgYYqIqVIp/9EYsPmzTYWK6fH8LRptr1KFfuuP/ZYWzp0sBNFEWLNGptb4+uvbVmzBipXhtNPtwRRr17qQSRS1tS0EIk96lkUyTZvtintc5JD48fbUDOApk3hvPMsOXT88TFf20BERCTqVKpkXWtOOsmur1xpvYl/+smSSF9+adsrVrQexV26WPLoyCPDKtviPUydmjvabvx4K91Uo4aNkD/7bBsln5QUdKQiIiJSECWLgrZq1a7JoT/+sC7pcXHQti0MGGBnErt0gdq1g45WREREylLt2nDBBbaAtRt++cUSR7/8AvfdZ5mZcuVsMotjjrFK0J06lXm7YeNG6wyV03to+XLbfvjhMHiwJYmSk6NmNJ2IiEhU0zC0srRli9UlmDDBprGfMAEWL7bbEhOtkdeliy2dO2tKexHJ9c8/tlaPQhHJa906+O233OTRxIl20gmgQQPrcdSpk607dCjR3kfp6Xa+a/RoG8YydqyNnt9/f+jZ05JDvXrpXJdIuFLTQiT2aDa0cJCZCbNm7ZoYmj7dtgM0amTJoSOPtHVysvpii4iIyL7ZuhUmT7Z2x4QJls1ZtMhuS0iwXstHHmnL4YfDoYeGPCFGero1acaMsQTRuHGwbZtNZ9+hA3TrZgmio46KqHJKIiIiMUPJonBw+eXw6qt2uUoVSwjlJIeOOAJq1Qo2PhGJLG+8YetLLw0yChGJRKtW5SaOJkyw3kcbN9ptSUnQpo1NmtG+vWV9WreGpCS2bLHR8T/9ZMmhsWMtF+WcTc7WtauVUOzSxZo6IhJZ1LQQiT1KFoWDX36BhQstOdS0qdUgEhHZW5qyRERKSmYmzJ5tPZD++AMmTybrjynM23AgEziS8a4zExKPY1p6MzK89Tpq22oHx/dIoOvxji5doFq1gF+DiOwzNS1EYo9mQwsHObWHRERERMJJfDyra7fi9+qtmLD/hYxPhN/jPWtxAOxXbhsdK83mXwnPc+SmHzmKsVSfsQZWVofJreD7w6BVKzgse129esAvSEREREqakkUiIiIiUSorCxYssGnsc5YpU2DJErs9Ls5yPuec43ZOota8eRJxce2AdrCqn9VcnD4dZsyw9TvvwIYNuU9Sq5Y9SIsWVgOpWTNb6tfX1GciIiIRSskiERERkSiwfj3MnLlrYmjaNNi82W6Pj4fmza3jc9u2ufNrVK5cxIPWqmVL9+6527yHZctyE0g5SaQ338ythQQ202vTprnJo5xE0iGHQM2aVvxIREREwpKSRSIiIiIRIivLegXNnr3rMmcOrFiRu98BB1hCqH9/W7dtayPGSmTiVeegXj1bevXK3e69FdOeO3fXZfZs+PJL2LEjd99KlWxm2MaNbWnSJPdy48aw334lEKiIiIjsLRW4FhGJBFu22LpixWDjEJFSl5kJS5dCair89Zct8+ZZzmXuXJuRLEeVKjb6q3nz3KVtW2jQIMw67mRkwKJF9gLmzct9YTlL3h5JYHWQGje2F1Kvng1py1nq1YM6dSBB5zxF9oWaFiKxR7OhiYiIiISpjAxYudJGci1cuGtSKDUVFi+2fXLExUHDhrsmhHKWqBjN5T2sWZP7BuRNIi1ZYsumTbveJy4ODjooN5FUp45dP+ggqF0793L16pqRVkREJJuSRSIi0WboUFtffXWwcYhIkTZtspFYy5ZZ76ClS3e/vGKFDSfLq2bNXUdh5R2V1aABlCsXzOsJG+vX5yaOli7NvZyzrFixe+8ksN5HtWrlJpBq1bI3u7ClQoWyf20iAVHTQiT2KFkkIhJtuna19ZgxQUYhEnO2bbNOL6tX27JqFfz9t60LupwzrCOv/feHunWtE0zedd26VranUSOV6CkRmzdbl62VKy15tGJF7uWc9apVkJa2a9etvCpVyk0cVasGVavaOu+Sd1vVqjYWMCkpCrp4SaxR00Ik9hQnWaTB3iIiZWjVu6tIHZxK+uJ0Ehsk0uThJtS6oFbQYZW4na9zUTrEA5mQ2DDY15v3vY+vFk/Wtiz85uwTJuWBDCALiIc6A+vQbGizEnmugv7O+WNxODLWZBR5TBT2nlY/uTqrv14d0jFVUFxAqR2TxT3e8+9f/eTqrPpgFZmrMwFIqJ7AIc8esttjFHS/ne9J/URq392E+J61WL+encuGDbmX162zRNCaNbmJoZx13vpAecXHQ+9Kq7hkaypVd6SzqVIi87o2IeO4WtSqlZsMqlvXkkWRpDh/tyLf+3z3LfXjr1IlOPhgOPhge66nCnlc7+0Pn5ZW9LJ2LSxYYAfD2rV2v8KUK2dVxYta9t/fpp7bb7/cdd7LlSvbUoxaTKX1nVLU4xbntuJ8PgUh0uKV4My9ei7Lhy2HzNxtZdWuiZW2o4SfkHoWOed6Ac9izdNXvfeP5bs9EXgLOBxYDfT13i/Mvu0OoD/2r3Wd9/67op5LPYtEJFqtencVcwbOIWtL7viTuIpxHDrs0D1/6UfQ6b+CXmeOkF9vGcZUmDpX7V3CaE9/5z3FUtB7VJz4C3uPC3yMcuCcw2/3e7x/cYV6vHsP27fD8rdXsfi6OfitRb/GrHjHgrObs+iQWmzaBDWmrOLIX+dQLjP3fh7I28djG3E8yaH8SMGvqVw5K22T01kk53L+dbVquTPJ7/h2FfOu2Mv/5zBWnM+pUI7LnPsCZXb87dNnbWGysizDmJNRzLvkzUIWtmzYEPpzJSVZxeFKlWyd93Ke9arFhzDnh/ZkZcTnvs7yWRx62Vpq9ch+nAoVbElMLHwpX36Xuk5FvX+w+9+xyL9xPuH0P1Kc4zcc4i0NEdS0CNTcq+ey/MXlBd5W2sdIqXyeSUwr0WFozrl4YC5wArAUmAj0897PzLPP1UAb7/2VzrnzgDO9932dcy2B4UBHoA7wA9DMe5+Z/3lyKFkkItFqXKNx1iskn8SGiXRe2LnoO0dQi66w15kjpNdbwvYUU4HioWtG1xJ7rpzXHUos5esn0mZmZ7KyLJky/bBx7FgaevxxByVywDedycxk57LtjHGwKrTH2FE9kb8e6UxGhs12Xth6+3ZbCrrcf+Q4Dti2+/P9E5/IoOqd2bYN0tNtARjOOGoTWnwrSeSCuM5Urgyvbh5Hzcw9329HtUTWvdS5wA4fFSoUfwTRPv0/h7HivK5Q/68SGyYChPw/uK/vYVj+bbKybJjcpk1WW6mo9aZNNp5x8+Zd1/m2jVv1LOn+wN2eKpGVdKZf8eIrV25n4mjcupdIz6q5++MmrgUc6elVdr+tkhUgT99ceY9PlVg1nc73zbEeVAkJ1k0vZx3qEheXuy7qchHLuI4LSF9ayHDEvPHWL0/n+cn2IZGzxMVFxbDDCGpaBGpMwphdehTlV5qfLWH5eSYRraSHoXUE5nvvU7MffATQG5iZZ5/ewH3Zlz8CnnfOueztI7z36cBfzrn52Y83LpTgItmjj8Lw4UFHIRKdwqzUWsieW5ROQU3LbYvSad266Pu+nmrry/awH4T+/pTWfkMLeZ05ti1Kp3nz0J4n77b8txd2W87lvNte30NMBfGZVhPX+10fs7DrOctnG9MpaO6lrYvSKVcOvsso+Pa8ti1J36WGzY/s+T55ZaxIp127XbcV5zHiV6dzxRVF3B5vnRFylnLldr++fwGJIoDqmemcdVZup4akJFvXGhx6Mqy2Sycjw36rjYkL7X7l1qbTp0/IT7FH6YsLft7CtkeK4ryuUF9rcd+TfX0Pw/JvExeXO+zsoINK5CHT48YUvN3VhilTreDW1q225GRm97Rs3076CzUKftwCkkQ7b9tcKfS415aD668Pef/Sks6PEMKnYvqSbfYhVZiCkkg5S/7b8y/5b897vaDLJb1e8Z5dbnL+rsmv/ImwghJjJb2tIMVJyJXGY+bI/A8U0YpIX7SNPTbk9lL6oucKfO7SfM5ARELy9brrYMCAoKMoU6Eki+oCS/JcXwocWdg+3vsM59x6oHr29vH57ls3/xM45wYCAwEaNGgQauxhrWZNaNo06ChEolckfKfkt3FBIvtv3f3HysYKiTTbw2inx5qNASDUQVEl3WYqzn6bFiey35bCf5RtqphI27Z716Ysqv1a0OWc9ea0RCpvLt4PRe+gd+/Q2/U5y7bXEqm4sYAf1vsncts1kP58IhU2FB3LjiqJPHFn7knwjAcTKb829GWLtQwAAAh0SURBVPh9jUQ+GbbrSXb3f4nwd4iJlbqJLJ1gSZ+EhF3X8fGhHQ/jGiUWeDY0qWEiL75YwP7DCt6/IIkNEnfGkNggtPslNijix95eKOx5S/p5ylpxXldx3/vi/H33RbT+bfIr8nW2abP3j/tlYT0ZkoCC/45F3bbbvvUTYcpq66KYd8nbFbKoJSsrd72ny97nXs+3JN6ZQfqa8nuOt8p2uPVhdnb1zFnyX8+/HQq+PZQzEAVdLoX1GLKnQ/NH577gos7MlNa2ghTnzGBpPGZeM701CgqRWGEje2zI7aXEBRtJ37p74bvSfM4yFylngatWDTqCMhdKsqig/4z8f9HC9gnlvnjvhwHDwIahhRBT2BswIOYSjyKyB6vebVLguPMjhzXh9AsCDKyEFfQ6c8RVjKPjsCacVsavd1XPwmMqTL0r6/Dy0L14ruSC/84dhjbhpAtgVcuiY4mrGEe755twYp73aNWBoccfVzGOFkOaUOvMfHE9XcBjFFIzptnjTai126md4mnycMHvQ05R41D2L4gr73Z5jFDuV9Tz7q3ivr5IUZzXVdz3PtTjb1/fw2j92+RXWq9zT49bnNvyi6sYR5NHm1oBsIA1qRxazaImz7eHC3qVYWQSbursoWZRk2FHwgWnl8pzNymkZlFpPqdIjlB6pC8F6ue5Xg/I/9+ycx/nXAJwALAmxPuKiMSEWhfU4tBhh1r9DmfjzaOxQOEurxNsagSCfb353/v46vG4SnnOZ5Qn9xsxfu+LWxf0XPlfd0GxJFRPKPKYKOo9rXNVnZCOqYLiavHfFjR/vXmpHJPFPd4L2r/OVXWIr55bvDehegLNX2++y2MUdr/S/j+L1v/n4ryu4rz3ZXn8RevfJr/Sep1FPW5xbyuL/8W9FWnxSnCaDW1Gnavq7PzuzVEWx0isfJ5JeAqlwHUCVuC6O7AMK3B9vvd+Rp59BgGt8xS4Pst7f65zrhXwHrkFrn8EDlGBaxERERERERGRslOiBa6zaxBdA3yH5VNf997PcM49AKR47z8HXgPezi5gvQY4L/u+M5xzH2DFsDOAQUUlikREREREREREJFh77FlU1tSzSERERERERESkZBWnZ1FxZuIVEREREREREZEop2SRiIiIiIiIiIjspGSRiIiIiIiIiIjspGSRiIiIiIiIiIjspGSRiIiIiIiIiIjsFHazoTnn0oBFQcchpaIG8E/QQUhY0rEhhdGxIUXR8SGF0bEhhdGxIYXRsSFFiZbjo6H3vmYoO4Zdskiil3MuJdRp+iS26NiQwujYkKLo+JDC6NiQwujYkMLo2JCixOLxoWFoIiIiIiIiIiKyk5JFIiIiIiIiIiKyk5JFUpaGBR2AhC0dG1IYHRtSFB0fUhgdG1IYHRtSGB0bUpSYOz5Us0hERERERERERHZSzyIREREREREREdlJySIREREREREREdlJySIpU865fzvnZjvnpjnnPnXOVQk6JgmWc66Xc26Oc26+c+72oOOR8OCcq++cG+2cm+Wcm+Gcuz7omCS8OOfinXOTnXNfBh2LhA/nXBXn3EfZbY1ZzrnOQcck4cM5d2P2d8p059xw51xS0DFJMJxzrzvn/nbOTc+zrZpzbqRzbl72umqQMUowCjk2YvI3rJJFUtZGAod579sAc4E7Ao5HAuSciwdeAE4CWgL9nHMtg41KwkQGcLP3vgXQCRikY0PyuR6YFXQQEnaeBb713jcH2qJjRLI55+oC1wHJ3vvDgHjgvGCjkgC9AfTKt+124Efv/SHAj9nXJfa8we7HRkz+hlWySMqU9/57731G9tXxQL0g45HAdQTme+9TvffbgRFA74BjkjDgvV/hvf8j+/JG7Adf3WCjknDhnKsHnAK8GnQsEj6cc/sDxwKvAXjvt3vv1wUblYSZBKCCcy4BqAgsDzgeCYj3/mdgTb7NvYE3sy+/CZxRpkFJWCjo2IjV37BKFkmQLgO+CToICVRdYEme60tRQkDycc41AtoDE4KNRMLIEOA2ICvoQCSsNAHSgP9mD1F81TlXKeigJDx475cBTwKLgRXAeu/998FGJWGmlvd+BdhJK+DAgOOR8BQzv2GVLJIS55z7IXsseP6ld559BmPDTN4NLlIJA66Abb7Mo5Cw5ZyrDHwM3OC93xB0PBI859ypwN/e+0lBxyJhJwHoALzovW8PbEbDSCRbdv2Z3kBjoA5QyTl3YbBRiUgkibXfsAlBByDRx3vfo6jbnXOXAKcC3b33SgzEtqX8f3t3yOJFFEZh/DnBVQwmy8IqGsRqFDUsrEFEtolxEQx+AEHWIgar+AlMK4LIohuMFotgMAhaZAVd0FXEJphew4yzKli9F//PL81MOmGY4R7uvAMHfjlfwC3hGiXZxVAU3a2q9dZ51I2TwHKSs8AeYF+Stapy0actYKuqfu5CfIBlkXacBt5W1WeAJOvACWCtaSr1ZDvJfFV9SDIPfGodSP2YxTWsO4v0TyU5A1wFlqvqW+s8au45cCTJ4SRzDIMmNxpnUgeShGHuyOuqutU6j/pRVatVtVBVhxieGU8sigRQVR+B90mOjpeWgFcNI6kv74DjSfaO75glHICu320AK+PxCvCoYRZ1ZFbXsJmRUkydSPIG2A18GS89q6rLDSOpsXF3wG2Gv5LcqaqbjSOpA0lOAU+Bl+zMpblWVY/bpVJvkiwCV6rqXOss6kOSYwyDz+eATeBiVX1tm0q9SHIDuMDwGckL4FJVfW+bSi0kuQcsAvuBbeA68BC4DxxkKBfPV9WfQ7D1n/vLvbHKDK5hLYskSZIkSZI08TM0SZIkSZIkTSyLJEmSJEmSNLEskiRJkiRJ0sSySJIkSZIkSRPLIkmSJEmSJE0siyRJkiRJkjSxLJIkSZIkSdLkBymMV7R9rxGDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-3, 12, 100)\n",
    "plt.plot(x, mlab.normpdf(x, red_mean_guess, red_std_guess),'r',x, mlab.normpdf(x, blue_mean_guess, blue_std_guess),'b',alpha=1)\n",
    "plt.ylim(-0.02,0.28)\n",
    "y = [0.001 for i in range(len(both_colours))]\n",
    "plt.plot(both_colours,y,'mo')\n",
    "ymaxr =max(mlab.normpdf(x, red_mean_guess, red_std_guess))/0.28\n",
    "ymaxb =max(mlab.normpdf(x, blue_mean_guess, blue_std_guess))/0.28\n",
    "plt.axvline(red_mean_guess, color='r', linestyle='--', ymin=0.067,ymax=ymaxr);\n",
    "plt.axvline(blue_mean_guess, color='b', linestyle='--',ymin=0.067, ymax=ymaxb);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2.\n",
    "### Likelihood of each point under current guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)\n",
    "likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "likelihood of point 1.7617657967268245 being red 0.18884523534521547\n",
      "likelihood of point 1.7617657967268245 being blue 2.7155455454919998e-05\n"
     ]
    }
   ],
   "source": [
    "print(\"likelihood of point {} being red {}\".format(both_colours[1],likelihood_of_red[1]))\n",
    "print(\"likelihood of point {} being blue {}\".format(both_colours[1],likelihood_of_blue[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3\n",
    "### Convert Likelihood into Weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "likelihood_total = likelihood_of_red + likelihood_of_blue\n",
    "\n",
    "red_weight = likelihood_of_red / likelihood_total\n",
    "blue_weight = likelihood_of_blue / likelihood_total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4\n",
    "\n",
    "### With current estimates and new weights we can compute new estimates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key bit of intuition is that the greater the weight of a colour on a data point, the more the data point influences the next estimates for that colour's parameters. This has the effect of \"pulling\" the parameters in the right direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_mean(data, weight):\n",
    "    return np.sum(data * weight) / np.sum(weight)\n",
    "\n",
    "def estimate_std(data, weight, mean):\n",
    "    variance = np.sum(weight * (data - mean)**2) / np.sum(weight)\n",
    "    return np.sqrt(variance)\n",
    "\n",
    "# new estimates for standard deviation\n",
    "new_blue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)\n",
    "new_red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)\n",
    "\n",
    "# new estimates for mean\n",
    "new_red_mean_guess = estimate_mean(both_colours, red_weight)\n",
    "new_blue_mean_guess = estimate_mean(both_colours, blue_weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Red Mean New: 3.1844279003578717 and Blue Mean New: 7.651045781252283\n"
     ]
    }
   ],
   "source": [
    "print(\"Red Mean New: {} and Blue Mean New: {}\".format(new_red_mean_guess,new_blue_mean_guess))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5 (Repeat)\n",
    "### Lets Run for 5 iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.1844279003578717 7.651045781252283\n",
      "3.4456051738982914 7.053630693427869\n",
      "3.2933810421494556 6.884581733005299\n",
      "3.1699583349617333 6.848715840027546\n",
      "3.0789827515638435 6.852124794685631\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIQAAADGCAYAAACjDVtCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VGXaBvD7pHcCqRASekcUjCIIiAqIXbHXVVRsKK5Y1rJr2c+GCjZUEDtrWzs2EEURpAWU3nsIkEAa6ZnM+f64PXtmwiSZhAlp9++6znVg6kmZyXnved7nNUzThIiIiIiIiIiItBx+DX0AIiIiIiIiIiJydCkQEhERERERERFpYRQIiYiIiIiIiIi0MAqERERERERERERaGAVCIiIiIiIiIiItjAIhEREREREREZEWxqtAyDCM0YZhbDQMY4thGP+o5nYXG4ZhGoaR6rtDFBERERERERERX6oxEDIMwx/AVABnAugN4ArDMHp7uF0kgDsBLPH1QYqIiIiIiIiIiO94UyF0IoAtpmluM02zDMBHAM73cLt/A5gEoMSHxyciIiIiIiIiIj7mTSCUBGC3y//T/7rsfwzD6A8g2TTNb3x4bCIiIiIiIiIiUg8CvLiN4eEy839XGoYfgCkArqvxgQxjHIBxABAeHn58z549vTtKERERERERERGp0fLlyw+YphlX0+28CYTSASS7/L89gAyX/0cC6AvgF8MwACARwNeGYZxnmmaa6wOZpjkdwHQASE1NNdPS3K4WEREREREREZEjYBjGTm9u582UsWUAuhmG0ckwjCAAlwP42rrSNM080zRjTdPsaJpmRwCLARwWBomIiIiIiIiISONQYyBkmqYDwHgAswGsB/CJaZprDcN43DCM8+r7AEVERFq8GTO4iYiIiIj4iDdTxmCa5ncAvqt02b+quO3wIz8sERER+Z+lS7m/8caGPQ4RERERaTa8CoRERETkKHE6gV27gI0buW3ZAixZAoSFAR98ACQm2lvr1oDhae0HEREREZHqKRASERFpKGVlDHys8GfjRmDzZqCkhNcHBgKdOgH5+cD+/cDkye73DwoCEhLsgOiYY4DRoxkeiYiIiIhUwzBNs+Zb1QOtMiYiIi1SWRnw1VfAF18wDHI6eXl4ONCjh/vWqRMQEACMGweYJvDccwyG9u07fMvIAA4c4OOccw5w8cW8v4iIiIi0KIZhLDdNM7Wm26lCSERE5GgoLmYI9N57DG769AGuu84Of9q1A/yqWOshOpr7Vq24de9++G1ME1i9Gvjvf4HPPwc+/hg44QTgkkuAU04B/P3r7UsTERERkaZHFUIiIiL1qbCQIc1//gPk5ACpqcANN3BfX/1/srNZhfTZZ6weio8HLroIuOACICamfp5TRERERBoFbyuEFAiJiIjUh/x84KOPuOXnA4MHMwg69tijdwxOJ7BgAfDJJ8DixZx+dtppwJVXAn37Hr3jEBEREZGjRlPGREREGkJODquBPvkEKCridK0bbgB69677Y77yCvfjx9fufn5+wLBh3HbtAj79FJg1C5gzB7jiCj5ecHDdj0tEREREmiwFQiIiIr7gdAIzZwLTprFx9IgRwNixQLduR/7Yq1Yd+WOkpAB33w3ceiswdSrw4YfAokXA448fWVglIiIiIk1SFd0rRURExGu5ucDf/w689BJw0kmsDnrqKd+EQb4WGgrccw/w6qusYLr+euCNNwCHo6GPTERERESOIgVCIiIiR2LlSvbkWboUuP9+Lg3fFJZ7P/FE9jcaNYpVTWPHAjt3NvRRiYiIiMhRokBIRESkLpxOLiF/001AYCDw9ttc4r2+Vg6rD1FRwL//DTz9NJCezmDrk0/4tYmIiIhIs6ZASEREpLby8tiP56WXgOHD2US6Z8/6e76EBG6VvfUWtyM1YgTw8cfA8ccDkyYBd9wBZGYe+eOKiIiISKOlZedFRERqY9Uq4IEHgOxs9g1qalVB1TFN4PPPgSlTWPV0//3A6NENfVQiIiIiUgtadl5ERMSXnE5WAr3yCqt13nzz6KzOdfAgkJ8PlJTUvIWFAUlJQPv23CIiavdchgFcdBFwwgnAI48ADz8MbNkC3H578wm9RERERASAAiEREZGa5eUBjz4K/PYbcNppwD//CURG1s9zlZUBGzcCa9dyy8y0l53v18/9ttu2AUFBwLHHAsHBQEEB8Pvv9vXR0QyGkpKA5GTuExMBvxpmjKekADNmAM88A7zzDr/+Bx6o+X4iIiIi0mQoEBIREanO/v3ArbcCGRnAvfcCl17q22oZ0wT27WP4s2YNsHkzl4APDGRfotNOY5WQvz8rdkJDGf6EhAAvvshjmTjRfrz8fGDPHjaJtrb164GKCl4fEAB06AAMGQKkpjJQ8sTfnyFQdDT7FOXnswF1VbcXERERkSZFgZCIiEhVMjKAW25hhcy0aazE8QWnE1i9mgHQ2rUMfACgbVs2qe7bF+jalaEQAHz4IffJye6P4ymYiori1quXfZnDwWDLCohWrwbefZcrig0aBAwbxueuzDCA224DWrViX6FDh4DnnuPUNBERERFp0hQIiYiIeLJ7N8Og4mLgtdd80y/I6QQWLwa++w7IymKlT69ebNzcpw8QE3Pkz+FJQACniyUlAQMHAmPGsDfQ/Pncfv4Z6NYNOOUUoH9/3t7VVVcxFHr8cVZLvfQS/y8iIiIiTZYCIRERkcp27GAY5HAAr78OdO9+ZI9XUQEsWQJ8+y1w4AB79Nx6KyuBKocvR4NhMADq1o1T4BYtYjA0YwYbUQ8eDAwdCsTH2/c55xxWHv3jH8CNNwJTp7pfLyIiIiJNigIhERERV1u3MqwBOE2sS5e6P5anIOj224FjjqldH6IOHTxfnpBQ92OzREYCo0YBI0cCGzYAv/4KzJ0LzJnDqqgxY+ypasOGcZW1u+4Cxo5lKFTVsYmIiIhIo2aYptkgT5yammqmpaU1yHOLiIh4tHEje+YEBbEyqK5hR+UgqEMHVtjUNghqKLm5wMKFwLx5XLls+HDg/PPZ0Brg92n8eP775ZfZ/FpEREREGgXDMJabppla4+0UCImIiIDNncePB8LDGQa1b1/7x2jqQVBlRUXAV1+xaigyErjkEuCEE/i17NrFaqe8PDacPv74hj5aEREREYECIREREe+tXAnccQeXWJ82zfOKWzXZu5fLs+/axSDo3HPZI8gXQdATT3D/0EPul8+cyf3VVx/5c1Rn507ggw/YW6lHD+CKK/g9ysxkiJaeztXHBg+u3+MQERERkRp5Gwj5HY2DERERabSWL2eoERvLpsq1DYNMkz13/u//gOxsYNw44IEHfFsVtHMnt8r27+dW3zp0AO6/n6uN7d4N/PvfwJdfMkB74w2gc2fg3nu5nL2IiIiINAlqKi0iIi3XkiXA3XdzOfZXX2UoVBvZ2cA777CnTr9+wDXXcCWu5sjPj02l+/cHPvsM+P57YOlS4LLLuAz9DTcAEyYAb74JdOrU0EcrIiIiIjVQhZCIiLRMK1cCf/87V/56/fXahUGmyaXaH3uM06iuvZbNqJtrGOQqMhK47jrgnnuA4GAGaR9+yGltgYHsK3Q0qpZERERE5IgoEBIRkZZn506GQYmJwGuvAW3aeH/fQ4fYZ+idd7gc+7/+BZx8ctNsGn0kunUDHn4YuPhiVki9+SaXoy8sZD+m/PyGPkIRERERqYamjImISMuSnc3Awt+fU52io72/76pVwHvvAcXFwEUXASNGcCpVfevRw/Plycn1/9zV8fcHRo5kv6TXX+eKZJdcAvznPwyHXn0VCAlp2GMUEREREY+0ypiIiLQcxcXAzTcDW7cC06cDffp4d7+SEuCTT4CFC7kc/dix7DvU2Did3Ft/2133rv/29wcCfPyZUGkp8O67bNIdGgr89hswdChXH/P1c4mIiIhIlbxdZUxnaCIi0jJUVHD1rw0bGFJ4GwZlZgKvvML9mWcC55zT8AFHYSGQk2Nv2dlAbi5QVub9Y0REsDqq8hYWVrdjCg4GbrqJK4599hn38+axt9C//tXyptSJiIiINHJendEahjEawIsA/AHMME3z6UrX3wLgdgAVAAoAjDNNc52Pj1VERKRuTBN45hlgwQKGQsOGeXe/TZvYY8jPD5g4kX1zjiankw2aH32UVUpjxjAAsoKfhQvZyPmCC3hsoaEMXqzwZdky7gcO5N66rqwMyMtjiLRxI1Bebj9nYKAdDrVuzabb3vZYMgxOo0tJYQXWwYPAxx/z/nfc4ZNviYiIiIj4Ro2BkGEY/gCmAhgJIB3AMsMwvq4U+Hxgmubrf93+PACTAYyuh+MVERGpvXfeAT7/nKtjXXSRd/f5/Xdg5kwgLo5hRm2XpK+rigpgzx5g+3auYFZaCqxbx2lefn5A164Malq3ZmVQYCCrljzZvZv7fv2qf87CQoZDrltGBrB5M5eWj4wEOnbklpBQc9+k7t3ZcHraNFYLvfACj/fqq2v5zRARERGR+uJNhdCJALaYprkNAAzD+AjA+QD+FwiZpum6lEg4gIZpTCQiIlLZd98BU6cCo0dzafiamCbwxRfA7NlAr17AuHF1n0blLYeD4c327VwBrbwcCApipU2nTsDixZymdu657vcLDPTN84eHc6vcF6m4mMezYwdDqdWr2SQ6JYXhUPv2VU+fi45mVVVyMquzHniAFUzeBnIiIiIiUq+8CYSSAOx2+X86gIGVb2QYxu0A7gYQBOA0nxydiIjIkVi2DHj8cSA1lX1saqpsKS0F3n4b+OMPTiu7/HJW5tSHsjJg1y5g2zYgPZ2hUEgI0KULQ6B27eznbqieRaGhQM+e3MrLeZw7dnDbtInH1b49w6FOnQ4PqAICgCuvZIB0882stDJNLlUvIiIiIg3KmzNMT10gD6sAMk1zKoCphmFcCeBhAH877IEMYxyAcQCQkpJSuyMVERGpjS1bgHvuYRjx7LOsuKlObi6XSd+1C7j0UuC00+qnEXJuLpev37SJPYLCwrisfKdOQGLi0VnGvi4CA3mMnTrxuPfutcOhHTs4xa5XLy5BX7miasgQTtm76CJgwgRef9ZZDfBFiIiIiIjFm0AoHUCyy//bA8io5vYfAXjN0xWmaU4HMB3gsvNeHqOIiEjtZGYCd97J4OHll9kDpzq7d3MlseJiTiurqedOXY/pzz8Znvj7s+qmWzcgPr7m4Kmq4+ncufr7JSTU6VBr5OfH6WVJScDJJ7Px9erVDLpWr2afo3793JtR9+gBfPklV2q7/XY2nR45sn6OT0RERERqZJhm9bmMYRgBADYBOB3AHgDLAFxpmuZal9t0M01z81//PhfAIzWteZ+ammqmpaUd4eGLiIhUUlgI3HADmyLPmMEGx9VZuZK3Cw8Hxo/nFChf2r2bQdDevVyavU8fbqGhvn2exiA/H1izBtiwgVPgkpOBY4/l9DfL0qWswAoPB557jgGRiIiIiPiMYRjLa8pkAC8qhEzTdBiGMR7AbHDZ+bdM01xrGMbjANJM0/wawHjDMEYAKAeQAw/TxUREROqd0wn885/sy/PyyzWHQXPnAp9+CnTowMqgVq18dxxbtzJsys5m+DFoEKuCfNUIujGKigIGDwYGDADWr2c49M03XKGtXz9WNJ14IjBlCnDXXezvVFwMXHhh/UzPExEREZEq1VghVF9UISQiIj43bRrwxhvAffexCqUqpgl89RXw/ffA8ccD11/vm6DG4WB1zKpVQEEBl1o/9lhOoTqS3kD33cf9pEnul7/+Ove33OL5fj/+yH1DTc2qqODS9atWsXdSRARw3HEMxiZPZs+mbt3YePqaaxpv/yQRERGRJsRnFUIiIiJNws8/Mww67zzgkkuqvp1pAv/9L/DTT8DQocBVV/mmOmXbNi4PX1DA3j0nn8yG1r547Nxcz5cXFlZ/v5KSI3/uI2H1SurRg826V64EFizgEvaXXcYqqnnzgFmz+LXcdFPzrqASERERaUQUCImISNO3dSvwyCNA377AP/5RdQhjmsAHHwDz53MVsUsvPfLAJicHWLiQPYvatAHOOce9Z47we9yhA7ft24FFi1iddcYZ7LG0YwfDtKIiNpxujv2VRERERBoZ1WaLiEjTlp8PTJzIFcWqW17e6QTefZdh0OjRRx4GlZVxqfVPPwUOHmRF0JgxCoNq0qkTv/epqeyvNGQIfzb79wObNrHRdH5+Qx+liIiISLOnCiEREWm6nE7gwQeBffu4jHlcnOfbVVQAb70FpKVxStlZZ9U9DDJNYONGrpZVUgL06gWccAIQElL3r6OlCQhg4+kePYAlSxioffYZEB3NQG/SJDadjo1t6CMVERERabYUCImISNP1yiucavTww1zFyhOHg2HRypXARRcBo0bV/fkyMzk9LCsLSEzkilpHI7Q48UTPl/fsWf39kpJ8fyy+FB7OqXu9e3O1sVmzgOHDGeA98wxwzz3sxyQiIiIiPqdVxkREpGmaPRt46CE2kL7/fs+3KSvjSlxr1wKXXw6cemrdnqu4mJUsmzZxatrAgVwdS3zH6QRuuw2YM4e9hXJyGGjdd59CIREREZFa8HaVMfUQEhGRpmfjRuDxx4H+/YG77/Z8m9JS4OWXuaLVtdfWPQzatg345BNgyxYuIX/ZZQqD6oOfH/Dii8CgQaz66tgR2LABeOwxVmaJiIiIiE8pEBIRkaYlJ4dNpFu14rQiT8uUFxcDL7zAEGfsWDZ8rq3SUi5lP3cun+vii1kZ1BDLot95J7fKXnqJW1W+/55bUxEcDEyZAsTHA3/+CYwYAWzeDPz978DevQ19dCIiIiLNigIhERFpOhwOTg/Lzgaef57LvFdWWMhQYedO4Kabqu6/U530dK4etm0bV8M67zw2PG4oJSXcKisv51YVh4NbUxIfz9XicnP5/b/1Vq5AdsstrPYSEREREZ9QU2kREWk6pkwBVqzgdLFevQ6/vriY04727GGAUFWj6ao4HOwVtHYtA6AzzmjaK12ZJhs0FxXZ4VF5+eEhkad+gqbJldgCAlgVZe2tf/v7199x9+vHKrBnnuG/n3oKeOQRXnb33WxEXZ/PLyIiItICKBASEZGmYdYs4OOPgauu4rLxlVk9g3bvZlVJbcOgzExg3jwgLw845hguJR/QBP5MOp1Afj5Dn9JS9+Bn61beZv169/tUDlMMgxVVANChg32ZaVZdYeTn5x4UhYQAoaFsuh0czPsfiYsvBlav5gpxL77In+0DDwCTJ7OC67zzgJiYI3sOERERkRasCZzpiohIi7dhA6tETjjBcy+d8nJg6lROMbrpptqFQU4nq47++IPLoJ9zDtCune+O3ZfKyhj8WFtGBo9/82Zeb4UzQUEMZlq35mVdutjVPYGBnsOa0FDuu3d3v9wKhRwO9woj1+CppIRBmlVp5OfHgCgszA6JQkNrV9VjGMCDD7IP1MMPA++/z6lk//438MUXDMGGDWOjbz/NgBcRERGpLQVCIiLSuOXnc+nx6GjgyScPDxUcDuC117gk/PXXA8cf7/1j5+SwKujAAQYhgwczTGkMTBMoKODX37Urw6DVq+3rQ0IYhgQH89jDwg7/3gwYwL03/Y9atfJ8uWHYQZIVGlV1vCUlDKqKi7nPzeX31hIczNAtKgqIjKz5ex0SAkyaBFxzDXDvvcDbbwP//CeDoSVLGIbt3MkpZFFRNX+NIiIiIvI/humpb8BRkJqaaqalpTXIc4uISBPhdAJ33QUsXQrMmAH07et+fUUFMG0asHIlQ4MhQ7x/7I0bgQULGHQMG8ZlzhtaaSkrbfLzgUOH+PUbhl1lY22hoU2nKqa83D0kKiiwG2GHhDDIsQKiqr6mBQv4e3DOOewltHMn+0k5HEDv3gyZhg9vHD9DERERkQZmGMZy0zRTa7qdKoRERKTxmjED+P139o6pHAY5ncBbbzEMuvxy78MghwNYuJCBUFISq0uqq3ypTxUVDH7y87mVlvLy4GA2s64pKGkKAgNZfeRagVRcbH/NBw6wf5NhABERdkAUFmbffsgQTgV84w32d7roImDCBPYW2rmTl82Zw6mCJ57YtL9fIiIiIkeJAiEREWmcfv+dAcDZZwNjxrhfZ5rAu+8CaWm87tRTvXvMvDzgxx+5bP2AAZxedqTNj2vL4eBUtZwcVsuYJgOMqCggIYH74GD3+4wbx/306e6XP/889xMnen6uWbO4P/fcmo9r0ybuK/cQqg+hodwSEhjsWVPj8vO5QtyePQySWrcG2rRhBdBNN3H1t2efBXr0YEA4YQLwwgvsMzRkCLBqFcOl00/nfURERESkSgqERESk8cnIYCPhbt1YHeQa2pgm8MEHwOLFDDrOOMO7x9y2Dfj1V4YvZ54JJCfXz7F74nQyjDp4kKGHaXK6lBUARUQc/WCqsbDCMKsHUHk5q6Zyc4GsLAY8wcEMhh5+GLjxRvaUmjkT6NwZuO02rkC2bBlwwQX8vfjsM1Z+tW/fsF+biIiISCOmmmoREWlcSkvZQNg02VA4JMS+zjSB//4XmD+fQdDZZ9f8eE4nsGgRMHcuK04uuujohEGmyWBjxw5Wrmzbxh468fHse9OnD6esRUa23DDIk8BAhj+dO7NpdseODIT27gV27wZuuYVB0f33c8pdz57AzTfzuu++Y5+hsDD+Oy3NXvlMRERERNyoQkhERBoP0wSeeYb9faZMca/wME3gyy+Bn35i9ceFF9YcpBQU8Pb793OK0Ukn1X9/meJiVgJlZ7Paxd/fnvrUkiuB6sLfH4iJ4VZezml24eHAVVdx+twjj7DZdJ8+wNixwJtvAh9/zOllixcDK1bwZ9+QfaJEREREGikFQiIi0nh8+SXw9decFjR0qPt1338P/PADL7/00pqDlfR04OefWUUyYgQrTuqLaTIAysxkFZBhsIlymzbcq8nxkQsMZHVVfDzQqRO/1198wdDwpJOAxETg4otZQfbuuwyF2rblCmWffcbfgcTEhv4qRERERBoNBUIiItI4rFvHKWKDBtlNlC3z5gFffQUMHMjqkOrCINNkZcjy5QxkRowAoqPr55jLyzl9KSuLzaJDQoCUFFYEBfjwT+zIkZ4vP/746u/XpYv3z9G6tfe3bWjBwcCjj7LX1Pvvcwqe08mvYdAg9ooKDgb+9jeu1jZ3LhtsDxzIlchEREREBIbZQHPrU1NTzbS0tAZ5bhERaWRyc4Grr2bQM3Om+xLlixcDb7/NfjI338xpRFUpK2NV0K5dXC1ryBDfBjOWggJWqOTmMoCKjmblSmSk759LqpaZyYCwVStOFysqYjj3009cpW7ECE4lM02GRNu3s1H50KH183shIiIi0ggYhrHcNM3Umm6nsyEREWlYTifw0EOccjVjhnsYtHIlp//06MEpQNWFQXl5wOzZXMXr5JPZV8bXx5mTY08L8/dnCBQXd/gy8b5WUsK9a4NtgAEYAAQFeb6fw8G9N+GH08l9U5reFh8PPPUUVxp78kng6ac5Lax9e4aLs2ezp9OFFwKDB7MXUVoag7xRo7Q0vYiIiLRoTeisT0REmqXXXweWLOFS4r1725dv2MDGwSkpHPAHBlb9GLt3s59MSQlw1lm+DYMcDk5NWr2aK4aZJo+pXz8GD/UdBgHAnXdyq+zll7lV5fvvuXljyxZuTU1qKjB+PKuCZs5kENSmDTBhAn8Xli0Dvv2Wv08REQyGcnOBzz9nw2kRERGRFkoVQiIi0nB+/RV46y3g/POBCy6wL9++HXj1VVaA3Hnn4ZUxrlatYqDUujWXovfVtK3ycgYGWVmsntG0sMbrmmuANWsYjvXqxZDIMDhdzOlkVVBSEptRl5czeFy7ln2Fhgzh0vUiIiIiLYwqhEREpGHs2gX8618cnN9/v315Rgbw0ksMXiZMqHpaj8PBZtOLFwMdOzJU8kVYU1bGiqM1axgIRUfzGLt0URjUWBkGm0ynpAAPPMBpfQCnv40dCxxzDCuliot5m7AwBkcOBy9fsMCeMiciIiLSQigQEhGRo6+oCLjnHk4DmzTJ7oGTlQW88AIv//vfq14drLCQ1R2bN7MaZOTI6qeUeaO0FNi5k0FQVhYrjvr2ZVVJaOiRPbbUv7Aw4Nln+XO87z67v1JAAHDLLWwm/c47wL59/Ll2785eU61bcxWyjz5iYCQiIiLSQigQEhGRo8s0gccfZz+eJ59kE2CAfV1eeIFVG3fdxeXCPdm/n/2CrMbAAwYc2fGUlPBY1q4FDh7k8/bty6qjo9EfSHynUyfgkUcY6k2ebF8eGAjcfjuQnAxMm8YgMSaGP+fzz+dy9OvXAy++CGzcqGohERERaREUCImIyNE1cyYrMsaPB048kZcVFjIMOnSIPYPatfN8340bWRkUEMCBfMeOdT+OkhL2Klq7lquHxcVxalFKStWrdjWUc8/lVtmgQdyq0qMHN2/ExHBr6k4/Hbj2WuDTT4FvvrEvDwnh71ZsLDB1KqvBDIMVQmeeCVx3HVeO+89/uDrZgQMML0VERESaKcNsoJOd1NRUMy0trUGeW0REGsiyZazUOPVULhFuGAxmpkwB0tM5YPcUYDidbBy9ejWbA48YUffqnfJy9ik6eJDPHx8PJCR4tzS7NA0VFfw9W7WK08S6d7evy8nh1LKSEuDee4G2be3rioqAr78GNm1iz6jjjmM42br1Uf8SREREROrKMIzlpmmm1ng7bwIhwzBGA3gRgD+AGaZpPl3p+rsB3AjAASALwFjTNHdW95gKhEREWpj9+4GrruLg+t132fOlvJwrQ23eDNx6K5dyr6ysjEuK797NKT4nncRmwbXlcPAYMjNZ+REXx+lqR9p76GjIzeW+ck+lggLuIyI836+khPvqVmmzOBzcN5dgLDsbuPpqfj0zZwJRUfZ1mZkMhfz82G/ItTKqogKYPx9YuZL36dOH+6Qk98cQERERaaS8DYRqPKM2DMMfwFQAZwLoDeAKwzB6V7rZHwBSTdPsB+BTAJNqf8giItJslZWxGqOsDHjuOYZBFRXA9OmsxrjuOs9hUH4+8NVXwJ49wNChwODBtQ+DnE42El6zhvvWrTnIT05uGmEQwNDivvsOv3zaNG5V+fFHbt7Yto1bc9GmDfDMMwx/Hn7YvS8Imj6+AAAgAElEQVRQfDxXsCsrY3VaXp59nb8/K9iGD2dguWkTpzJu3sx/FxYe9S9FREREpD54c1Z9IoAtpmluM02zDMBHAM53vYFpmvNM0yz667+LAbT37WGKiEiTNmkSsG4dm0l36MDB+dtvc0rPlVeyqW9l+/YBX37JAfhZZ3GZ8NowTa4WtmYNA6WICC4fr2bRLccxxzCI/P13YMYM9+vatwfuuINh0IsvHh70HHcccMYZdigUFsaKqw0bgC1btCKZiIiINHneBEJJAHa7/D/9r8uqcgOA7z1dYRjGOMMw0gzDSMvKyvL+KEVEpOn64gsGO2PHsurCNIEPPmA/oTFjgGHDDr/Ppk1sCBwcDFx4YdVNpquSnc1m0bt28TF69AC6dtXy8S3RmDHAOeewGu2339yv69wZuO02TiV85RUuWe+qY0fgvPPYa2rRIoaKSUmcqrduHZuSW8vbi4iIiDQx3gRChofLPDYeMgzjagCpAJ71dL1pmtNN00w1TTM1Li7O+6MUEZGmae1aVgcNGgTccgvDoM8+48D8zDNZgeHKNIGlS4FffmGz3wsuAFq18v75Cgq4fPj27Zxa1rUrw6CqeuxI82cYwAMPAD17curY9u3u1/fqBdx4Iy9/7TW7l5IlNpahZJs2XB1v3z72skpMZIPqtWtZgVZRcfS+JhEREREf8CYQSgeQ7PL/9gAyKt/IMIwRAB4CcJ5pmqWVrxcRkRYmO5vTdeLigP/7PwY0333Hnjannspl412Vl/O6P//k1K4zz/R+aldpKfvfbNzIAX3HjnyM2oRJ0nwFBwPPP8/9xInsTeWqf38uVb9+PaeWufYbAjhd7JxzuPLY0qXAggUMhPr2ZaPvffsYDGmpehEREWlCvFlKZBmAboZhdAKwB8DlAK50vYFhGP0BTAMw2jTNTJ8fpYiINC0OB6sycnPZK6hVK64U9vXXXCXssstYuWEpKABmz2aINHgwB9reqKgA9u5l42DD4NSyhIS6rULWmF18sefLTzml+vv1rrwGRDWae+VuQgKr1W65BXjoIfYNcv09GTyYfYE++YSr4F13nfvvaEAAcPrpbEqelsZQaeRIoFMnNqlOTwd27uTvYvv2WpFMREREGj1vl50/C8AL4LLzb5mm+YRhGI8DSDNN82vDMOYCOAbA3r/usss0zfOqe0wtOy8i0ow98wzw3/+yifRZZ7Gp77vvshJj3Dj3gXhmJsOgigoOuJOTq35ci2myGiMjg+FTbCzDoKayapg0nC++AJ54ArjmGq40Vtk33wCzZrG31ZVXuodClm3bOK0xJAQYPZrTyQBOIduzhxVrrVqx35D6VomIiMhR5u2y895UCME0ze8AfFfpsn+5/HtErY9QRESap08/ZRh07bUMg1asAN57z+7V4hoGbdkC/PorEB7OKTmtW9f8+Hl5rMYoKQEiI1mNERZWf19PY7B/P/cJCe6X5+RwX9X3raCAe296KFnNkYOCan98TcmFF3IJ+fffB7p14++oq7PP5vdi9mx+Ly6++PBQqHNn/u7Nng189RWDzJQU/hyioxly7t3LxtNxcQwrA7w65RIRERE5anR2IiIivrNsGaflDBkCjB/PviozZnAAfeut9qDYNIHlyxkWtW3LqTchIdU/dnExg6D8fPaC6dKFg++W4J//5H76dPfL33qL+4kTPd9v3jzuzz235ufYsYP77t1rfXhNzt13A1u3sreV1W/KYhgMjcrK2EQ6OJgrjVUWF8fbzZ4N/PADp0L268f7JyQAMTEMhbKyOBWybVtOLfNUcSQiIiLSAJpZkwUREWkwu3cD99/PAfYTT3DA/dprrI4YP95uEO1wsJ/QihVcAezss6sPgxwOPva6dUBhIaeU9enTcsIg8b2AAE5rjI1lmHbggPv1hsE+VyefDHz7LQMfT8LDGRZ16gQsXgzMn283pA4I4O9q796s0EpP5+9wXl79fm0iIiIiXlIgJCIiR66ggFUXhgFMmcKqiFdeYW+VCRPsKV2FhWwsvW0bKypOOaXqBtCmycdZu5ZTcOLi2GxaVRbiC9HRXHmsoICr4VlT5iyGAVx9NXDCCew79PPPnh8nIAAYMQIYMADYsIEBUkmJfX1ICNC1KzeA0yS3bHG/jYiIiEgDUCAkIiJHxukEHnwQ2LWL08UcDuCFF1gV8fe/s9cKwHDniy9YITF6NKfXVOXQIS4BvmsXm/L27s0eLerDIr7UrRvw2GPA6tXAU08dvmS8nx9w/fXAcccBH3/M5eY9MQwgNRU47TSGl19+yRX2XLVqxd/j9u0ZQq1bx6qhior6+dpEREREaqAzaxEROTIvvcRVxB58kNU7kyezIujuu+1mx66rMp1/vr0qU2VlZRwk5+SwoW/nzt41mhapq9NOA266CXjjDfZPuuIK9+v9/Xn9q68CM2dyJbuBAz0/VteuDEDnzGEoNGIEAyCL1V+oTRuukLd/P/sLJSWx55C0GE4ns8CKCv7b4bD/XVHB/1v/BphN+vtzq+rfgYFVF1yKiIh4okBIRETqbtYsDpIvvZSD5OefZ6+gu++2B7grVgBpaRwIjxrleRlup5OD4337+P927Xh7jW7o6qs9Xz5yZPX3q64Kq7LKK5i1JDfdxGlcU6awWfmJJ7pfHxDApugvvwy88w5H3gMGeH6shAQ2m/7hB+D774HBg9nzylVgINChA6dB7trFht5ZWew5FB5eH1+h+IBpsrd9fr7n7dAh7ouLOSOwpMT+t+u+uPjwGYq+EhjI3D00lPvKW2got8hIblFR7pt1WWSk3n5FRFoCw6xcHn2UpKammmlpaQ3y3CIi4gMrVwK33AL078/qoBde4MfUEyeyUsjh4JLyW7dyas6wYby+spwcVgWVlbEaqH375r/0uTQ+RUWcHpaVBbz3nntlj6W0lL/nO3cCt93GnlZVKS9n36GdOzlVbPDgqkfY2dl8DZSXM0hNSuLIXuqdaTLEOXiQvcVd99Z24ABnAOblVT/Dz9+fM2VDQvjjCwzkW5n178BA3iYggMVifn7cDMPeAM9719N107S3yv+3qo4qKvjr5HBwb21lZQymCgpqDqUiIljMFhNj72Njubcui43l27Z+XUVEGhfDMJabppla4+0UCImISK3t3Qv87W8cMUyaBEybxsvvuYcVEkVFnDaTmclqi+OOO/wxiou5etihQ/zIOjnZ7jck7nbu5L5DB/fL9+/nvqrqHquPjTcrsllNjqtb8a2527MHuOYajnTfeoulEpUVFbGSKCMDuOMOoGfPqh/PNIGlSxmetmvHKWRVfX+dTr6u9u9nSqBl6n2iqIiFh/v28dtb+d8HDjA0qSwggG2fIiJYtBUczNAjIMCepmWxgpjycu9aQpmmvRida5hTOehxvb11H+t+VXENl1zDpsqssMqadmZlla7PYwVIxcUMkA4d4roAnrRpAyQmVr1FR6viSETkaFIgJCIi9aOoCLjxRg6In3sO+O9/OQqaOJGD2MxMhkFlZezP0rGj+/0rKnjfrCyOEJKS+DGzBr5VGzeO++nT3S9//nnuJ070fL9Zs7g/99yan2PTJu67d6/98TUnK1YAt9/O6p+pUz1XqxUU8Ht/4AAwfjzQo0f1j7lpE/Dbb+ytdcYZVffQAliFtHs3y1FCQhiUegqmBACDi337WGC1e7e9WYFPfr59W+uUNzqa39KwMDvocQ1EHA47THG9r9PJy4KDGQxZ1T8BAXblj+vzWPexAiPXIAhwr/xxPR13fd7vvuP+rLMOf+zKmxXkuFYJufYhso7Ful3liiKHw64osoIkKzCyvkarV5F1ufVYDoddeZSXZ38Pra8nKIi5ddu2LL5LTuaWksI/AcHBtf/Zi4hI1bwNhNRDSEREvOd0Ao88wn4rjz0GfPopz/ytMGjTJmD+fH6kfsEFhw98Dx7kyM3hYP+Udu20cpg0LgMG8Hf7wQf5u/7EE4eXNlgr6E2Zwr5Ct98O9OpV9WN2784Uwmo27SkotQQHszl1Xh6Tjc2bed/27VvsqNk0mTNv386WS67hz549DDAAhhp+fpzCFBXFypTkZF5mBSEhIXbgYgU8ISH81lrhUOWpXdZ9raqgymGLtQG8j+tmhSjW/4OCuFmhknVd5bDFuv2mTTyGa67h47uGRZX/Xfl4PG2uoU9pKXN7aysttZtZOxz29DLXnkhlZbyutNSuICor4/coIIDHHx/PYwoKcv/+lZRwcb3Fi3kff387eIqPZzhkbcnJfIkkJ2s6mohIfdJZuIiIeO+ll4B589iEd/58ntXffTdHXosWcfluT1NjCgs5eissZFjUrRs/nhdpjEaN4tStF19kWcNddx1+m6go/u5PmcJKoltvPbx5tKv4eGDMGIZCc+Zwmfr+/auujGvVis+xfz/LXdat47EkJjbbuTdOJ4sHt21jn+1t2+x/FxW5V++0bs1crn17BgumyZDFCnwqKhgkWE2UXcMeq/rHCl9cV/myQiXXQMeaXuXaDygoiG9hoaHch4XZgUjloMfa1+XHtns399XNTPQl1zDIColKSux95a20lD+bggJ7Kyy0N+sxSkv52EFBfClYlVKBgfzZlZXxV3zpUj6uNTUvIIAVRJ07A5068U9Ht26cPduSZ7eKiPiKAiEREfHOzJnczj2XPW1KSlglERfHFZXS0znN5qST7JGPw8GP8A8c4Jl/x45aXluahquvZhAzcyar3y677PDbREbaodCrrzIUqq7RdFgYXz/z53PlvYMHgeHDqy6BMAwGQG3a8HW0dy/v0749E5EmyjQ5Y3TzZlbAbN1qBz9WBUp5ObPjVq24tWnDtxWrabPTaf/fWlErKMiu6rECICtEcjgOb/JshT2uFUKulUJWCOR6H9fQ52iqzexPoO4zQP387K+5NioqWElUVGRvhYXu/87L458CKziyVmUrKeHPJyyMLykrILKCu7w8YMECLtpn9W/y82Ow1LEjFwbs0YOhWbduLbaQTkSkThQIiYhIzb7/nqsrDR3KkVphIcOgyEjgiy94dn/KKXYvFWvEl5HBkZvVPMLTKmMijZFhsEl6ZiZ7ZcXHA6eeevjtIiIYCr3wAvDaa8DNNwP9+lX9uP7+fJyYGGDJEuCrr9hXqLqG6kFBLI+IjWXJyLZtvH1yMtOQRqysjIe7ebMdAG3axEG+tfJVeLi93Lmfnx0MWAP/yEheZq3eZYUWwcH2lCQrBHJd4cuqzLECIytgsO4bHHx46GNNc5LasVZYi4io/nZOp3s1kRUMHTrEPxn79vF3w7osL48/t9BQfvZgrZxmGPz3ypXMV62gKCiIeWnnzgyKundnUBQXZ/9uiIiITU2lRUSkeosXAxMm8Kw6MZEf506YwLPvn3/mqGvUKHulq0OHOGgtLuYILzlZtf1HaulS7k880f3y9eu5r6p/zZ493Ccl1fwcVvddNTB2V1LCyp+NG4HXX6867CkqYiiUns5Q6Nhja37s9HRg7lyOZkeM4HTLmlQOW61eXI0gbC0qAjZs4LZ+Pb9l27czFCov56G3asXBfUAAf9UiIzm4t1bzCg21Qx+r6bO1YpYV7FhBjlWtY22WgIDDK36sECgwsGkFPh9/zP2QIdx781IGmu7L2TT5krOqh/LzGQplZPDlkpVlX25VF/n7202t/fz4u1ZQwGlqVm+mNm04zaxTJ/4p69mT30vrd8q1t5SISHOgVcZEROTIrVvHwW1MDM+kDYNh0MGDnPISF8cwKDyco770dCAnhyOv9u29W+5cpLHLyQHGjuUI9O232fXWk6Ii9tnauZMrw/XvX/Nj5+UBs2fzsQcPBnr39u6YHA57tT6r0UpMzFEb1VYOf9avZyWQ1Ww4OJghj2vwExTEvTXNy6rIiYy0Z5m6BkOuFTyuK3IBvL1VEWSFRNa/G0E2JvWkvNw9KMrOtpuNZ2YCubm8vLDQriIqK+N9S0r4e2uFRK1a8c9ax47c9+jBbNXqO+VaWSYi0tQoEBIRkSOzezcHwabJxgwREVxNacsWjvy6deMUMj8/1vnv38/7JSayWkhn0b5TVUMQq+NscrLn+x08yL03fZuKirhXs2/Pdu8Grr+er4O33qp66fiSEoZC27cDN94IHH98zY9dVsZqu127WLpw8snepxpFRTy2ggL+7FJSmKj4UHk5fwXXrAHWruW2dasd/lgNlsPCOMgOC2MWHBlpBz/BwQyGrMF45aXfg4LsSg+Hw37ugAAOyl0rOax+QS3Br79yb7Wm8rYFW0t8OZeUMBDKzbVXpdu5k7lpTg63wkK7R5XTyZdecbE9xTA6mpVEXbpw69yZMzWt30HXJuWqKBKRxkyBkIiI1F12Nge/WVmcjhQXx3Bo5UqebQ8cyKkzOTmsCior4wA5KanljNSOpnHjuJ8+3f3y55/nfuJEz/erTSfaunahbUlWrwZuuYXLwk+bVvVUyJISLke/bRtfNyecUPNjmyar7v74gyPQkSOr7ytUWXY2X4vl5UwNkpLq1DDFNDnTcM0abqtXc2+tGBUQwJDBqvoJD+fbQ3i4PT0rOtoOeFq35uZ6fUAAB+FW5Qbg3uvHtTqjpfd8Oekk7h96iPv6birdHFVUsGooJ4fZ6ebNfGnu2WMHRSUl9tTGigqGnVZ1WlISwyErKEpJsZuRW0GotdqcqtNEpLHwNhBSU2kREXFXVATceScHlz17shn0hRcCCxfybPessxj+bNrERg+hoay1r6mbqEhTd8wxwBNPAPfeyxH6s896roQLCeFr6JVXgDffZCnCwIHVP7ZhMDiKjwfmzQM+/xw47bSqq78qa9OGSczevazWy8nh/Jf4+GpLGYqKGPisWsXwZ8UKrgRlLRNurfQVE8N9fDxf6oGBdiNof3/urUMID+e3wAp+iovt5cwdDl5nTRuzBtItPfiR+mP1EGrThoHO8OG83OHgy2TfPoZEW7awomjvXr4GDh1iULRtG/thmSZ/d8PDWSDbuTOnmnXubBcMWtVy1u91WJg+IxGRxk2BkIiI2MrLOdhduZJnvF27spvp8uWcBjZ8OD9qXbeOo72UFFYzqHZeWorhw/kamTSJgdB993n+/Q8OBsaP53L0b7/NhGXYsJofv0MHYMwY4McfubrfgAGcdubNa8zPj+UM1mpkVhfe5GSgVSuYJge7q1bxJb50KXsAWdURISEMdhISuG/XjsFNUBAHvFYfoJgYO/yxpnw5HKwiKimxN6vqJzbWHiCrJ4s0FgEBrG6Li2PWayku5pSzjRv5+ti6lbM59++3Vz9bsgRYsICPERrKlnl9+rAiq0MHvoasaiF/fzscCg+3XzMiIo2BAiERESGnE3jsMVYndOzIM+TevdmIoU8ffrS6fTvr6a2VjQL0Z0RaoEsvZVnBe+9xNHjHHVWHQrffDrzxBvCf/7AT7tln1xzuREUB55/PEeeKFQx1Tj3V+9X6goOBrl3hyM7Hxt8ysfKrg0jb4sSyTVHYl+n/v+ofa9pXu3Z2D/iQEDv8CQtjmGMFQOHhHNwWF7OyqLiYORdgTyWzeghp0CtNVWgoQ50OHbhmAsDA88ABu4H6li38c2i179q8mb21AL6GoqOB446z/3R26MDbWZ06FBKJSGOhM3kREaGXXgI++4wfbQ4YwOqfwkIudR4YyIYLUVEcOYaGNvTRijSsO+5gIvLee0xDbr3Vc9ATFMTr3n+fPZ3y8oArrqi5TCYggNVICQmcrvn55+wrFBdX5V2Ki+1pXwsXAn/+GYVDhyJRVuJEoFGBVhGliI30Q2LfILRt6/e/wCc6mpVAcXEMf2Ji+H/D4FtAUZG9DDjAtwOrcbSmxUhLEBDA9RISE5nNAgx3srL4mlu7lrOoN2xgJVF+Pgv8Zs1i0BMaynZ8J53EGdZdu/KzlcxM95DICofCw7lpKqWI1Dc1lRYREU5pefppjvCGDOFHmlFR7CFkGDyj/WvaiTSAVau479fP/fKtW7nv0sXz/ayV3xISan6OggLu1QvKe04n8NRTwBdfsPG31fzbE9MEvvqKo8T+/YEbbvB+tJeVxSlkRUV8ffbsCYD93f/8k9NXFi7kYNRq1hwezpdw69bMcNtEV6BN4CHEBOYjKsJEXJcoxHWLRly8gagoHp4V/hQV8UsD7EGq60BVg9Sj59tvuU/9qy2oNy9lQC/nhmKarBpauZIzqzds4D47m6+v4mL3aWZDhvBn2707//xaPbes4ZnVk8j1NajG1SLiDa0yJiIi3nnzTYZB4eFc7rpXL7v7ZnAw55PExalPkIgnTifwf/8HfP01cNttXFWsOj/9BHzyCXt03Xab9+uCl5Qg64sFWL6oDL8e7IvFe9pjx04/lJZyOos1/SshgQNNqzdKdDQbQVv/bxNWAv+96SjILEShMwyFkYlwhHI1Mz8/92ks1spgIlJ3FRXsR/Tnn6wm+uMPYMcOhnaFhfZnLrGx7Ct/2mlA376sRiou5m2sqZmA3dg6IsJu4K4/zyJSmQIhERGp2bRpDIMiIzk1rGdP9g/q3Jlnp+oT1DioQqhxczqBRx8FvvuOq4tde231t1+2jFV5iYm8fXS0x5vt28fGzz//zH16uonSAgfgKEdUeAVaJYSiXUoAkpP5UDEx/FHHx9shUEgIK34KC+0KBQBAwSGE5O5HuF8xwuPDEd61LULbhGpg2cioQqh5ys9nMGQ1d1+3jj2KCgsZ8AYF2X2IRo4EBg9m0GsFRNbtAAa5rgFReLj+bIuIAiEREanJK68AzzzDMGjwYIYKxx/PdXTVJ6hxsaYiTZ/ufvnzz3M/caLn+82axf2559b8HJs2cd+9e+2PT1gG8M9/AnPmAHffDVx5ZfW3X78eeO01juImTIAZn4C9e4HffmNf9+XLuSJYaSkHfJGRzGiTkoDOsflIyNuExMgCJJzUGfGpKYiPZ2FfaaldeVBYyMMCOEC0BovW5u/3VxOUjAze0AqBNSes0TjpJO4feoh7b17KgF7OTY3TySbVK1YAv/8OpKXZDatLS/mSjIpik+pTTmGz6+7d3QMi16lmwcHuIVFoqKqIRFoabwMh5cciIi2NaQKTJwPPPcczxUGD2NNk0CCGQuoTJFJ7/v7Av//Nj+0nT2YCc+mlVd++Vy/sueo+/PLYr/j59B1Y7oxGZm4wSkv5UFFR7OueksIGtG3bso1XfDyQkBCF+IjeCFz4C8q2/YDC1Z2Q3XcYMjLs+V1hYfbKYBERVU39MvC/JGnvXoZD2dksQ0lM1PrwIkeJnx///HbpAlxyCS/LyWFAtGABsGgRVzZbsgSYPx948km+tnv2BIYNA0aP5loQroHwoUN8OVuPb70XWCGRehGJCKBASESkZTFN4IkngBdf5IhxxAh+1Dh0KKsD9BGiSN35+3Okdv/9wKRJDIXGjPnf1Xv3ArNnA7/+yoHe/v3tUV56BYLKDiE2cB96do5C8jGt0aULZ262a8dsJiGBmU1JCQd7BQXAlowwONqfidCS1YhYtxSt93+KNqefhrAubWvfeDYgwE6b9uzhgR44wBRK7wsiDaJ1a+D007kBfP1v3Mgqwl9+4cpmq1ZxytnkyfyT3q0bWwGeeaZdXWa9ZxQW8qVtCQ11D4nUL0ykZVIgJCLSUpgmB6rTprE5wdVXA1ddxTVw9VGhiG8EBHDlsfvuw/5HX8Pcee0wO+8k/PEHWzo5HPy0PjqaUz66dfNDnxQ/jN7wHtqUZqDi7KsRffbJCAriAK6ggJ/yp6fb00FCQljIFxFhIKJvP4SMaMdm1cu+ARz9WSqAOlT3BAezf1hhIZ9w1y6ui52UVGWfIxE5OkJCgGOP5TZ+PGd57trFgHnuXPYj2rCBvYmmTmXg06ULMHAgA6LBgxkaWe8rhYWsQjpwgI8fEGCHQxERvK2yYJHmT4GQiEhLUFEB/O1vwJdfsgrg0Uc5nSUkpKGPTKTZ2LePTYB//TUIf66YjJKt6QheV4CCsAPwj49F794MgY45hp/kJyZyi42NhOPQPXBMnQbzy/eQvXM3Mk+5BPD3h2HwU/yEBPvT/MMaxobEshJp4UKWHu3Zw6WKIiPr9oWEhzMozs3lY23dyidu357XiUiD8/dny79OnYDrrmNgvG8f3wZ++IFvBVu3cmWzN99kQNSxI9ePGD2aFURdutjTzKwtN5ePb733uE4zU7NqkeZHTaVFRJoz0+Qn/JdfzuYDnTsDH37IEak0HVV1iN29m/vkZM/3O3iQ+5iYmp+jqIh7b5dBF+zdy4XFrClge/cye/Xz43SPTsnlGJf1BIYXfY8DY++D/6UXISGB1T3Wp/TWJ/UOBwCnEzG/fIbWK+bCr2cPGDePQ3hCRO0+pd+6lXNKAE4FrWoFOm+ZJksIMjJ4kK1bs2JI80uOil9/5b5vX+69eSkDejkL5eSw/9A333Bxw1277ObTYWH803H88QyITjyRwbNpur8/FRW5Vye6VhHpbUCk8fLpKmOGYYwG8CIAfwAzTNN8utL1wwC8AKAfgMtN0/y0psdUICQiUs9ycthw4JZbgM2bedb37bd1rxoQaeEyMhgAzZ/PACgjg6sDGQZ7/HToAPTqxSavnTr91f+nVQmC/vUPOH9bgIJLxmL/mFtRWGR4HGCFh/9VtLd4MfD++0yObruNlTm1cegQp5BlZjJEHDyY61gfiYoKznnbv5+jw5gY9hg60scVkaPm0CEGRN9+y7eZHTvcA6L27bnGxBlnAKmprGAMDrZXMrNCImv1wsBA94BIq5mJNB4+C4QMw/AHsAnASADpAJYBuMI0zXUut+kIIArAPQC+ViAkItKA8vI4Ul22jKse7d0LnHUW8NFHWjWoqVq6lPsTT3S/fP167nv18ny/PXu4T0qq+Tny87mPiqr98TVTGRnA99+zgatrBZBhMA9JSWE/j1NP5VQMTv/ibVynYBQXVCDxnacR/csXKD7tHJTc8zAiogM8T/+y7NjBZemLijgf5Pjja3fwTicP+o8/OI6NoMgAACAASURBVNIbOpQHfKTKyzkvJSuL/4+L4xeuperrxccfcz9kCPfevJQBvZzFO0VF9hSzRYtYYFhSwreP8HD+vvXrx7UnBgxgyB0dzRDJ9T2urIyP5+fnHhCFh+u0Q6Sh+DIQGgTgUdM0z/jr/w8AgGmaT3m47TsAvlEgJCLSAPLzOYI9eBCYNQuYOZNnatdey1XF9LFd0zVuHPfTp7tf/vzz3E+c6Pl+s2Zxf+65NT9HVdPSWgjTPDwA2rfP/iQ8NpYVQAMGcHE+KwBq1cp99S/XwZG/v0sPjnAT4R/OgN8b04BBg4Bnnql5Pk9+PvD66xylnXkmcP75tX8dZ2Vx3lF2NhsXDR7sm3keZWVMyA4e5DHFx3O0qCYjPmWtFPXQQ9x781IGWvzLWeqotJQVkHPmAL//zqXurYAoIoIBUa9efA/s398OwR2OSiF4MR/PMFg1FBlpB0TKjkWODm8DIW/+aicB2O3y/3QAA+t4UOMAjAOAFF98SiUiIjz7yshgLfiePfxI+ccfeSY2eTIwdmxDH6FIo2Oa7Kfx44+cXfXnn/YqYFYF0HHHsShr1Ch7ClhICD9VLyhga50dO9ynT1gDn4gI3tbObwzg5puAxHjgiScY8r34YvVNYaKigLvvZt+v77/n63vsWI6wvBUXB1x4ISuF/vyTq4cNHcpE60gEBTEhS0zk+49VNZSQwHBIKxeKNDnBwcDIkdwA5r5WQLRoEYPGLVuAr76yK4i6dweGD+f7Zdu2QNeufPm7TjHLyuL7q/UcrlVEWttCpGF5Ewh5+iiqTp2oTdOcDmA6wAqhujyGiIj8paiIA8T8fE7jWLeO08JWrGDj17fe4sBPROB0MryZO9deojkz030KWP/+LKAZNYq9mOPjeV9rULNr1+ENVtu0sQc2XrXTOf98fqR+//3A9dcDL7/MYKUqAQHA1Vez++vHHwNPP82+QgkJ3n/x/v5sCNKpE8uf5szhF3jyyUc+GgsO5uNawVBGBr+xiYkMozRfRKTJCgpiNdCIEfx/eTkLDn/80T0g+uYbOyCy3lr692dA1K4dr7OC9IICzmy31jzQcvciDcubQCgdgOvyJe0BZNTP4YiISI2KijhVIzeXZ1JFRfwI76uvgO3buVz0u+/yYzqRFsrh4Eyrn35iALRmDSt6ysuZUbRpw4xkyBAGQN268bLycnvQsnEjp0sA9hLMCQk+WIL55JM5/e+uu1jxM2UKG3VUxTD4EXy7dsC0acBTT/F+1d3Hk5gYVgutXAksX85A+eSTj3wlMoBVS1262EF1ejqrhhISGAypYkikyQsMPDwg+uUXBkRLlvA9c+tW9iQKD2cg1Lkz32sHDOBbWEICL6tquXs/v8OXu9fbh0j98eZUZhmAboZhdAKwB8DlAK6s16MSEZHDHTrEAVZ+Ps+OIiNZFfTTTwyEMjOBU07hgNHbtYlFmonSUi6mN2cO8PPPwIYN/ATa4eDLpXVrYOBArgA2ahSnOURE2J9aWy24ysv5eP7+vD4mxv7U2qfFLr17s4rvjju4EuBTT/H1W53u3YEHH2Sz6alTGRJdfHHtmnL4+fGj+44dOZL76SeO4IYM8c0a5WFhTNcKChhc79njHgypx5BIsxEY6D7FrLQU+O03vg8vWcL35O3bGRiFhdlh0LHH8m0oKYmX9ejB+7sGRHv32s8TFuYeEGlxQxHf8XbZ+bPAZeX9AbxlmuYThmE8DiDNNM2vDcM4AcAXAFoDKAGwzzTNPtU9pppKi4h4KS+PA6qCAg6m4uIY/ixaBCxYAKxezdHsmDHAk09qQn5ztHMn95WnFllNGaqaPmR95BodXfNzWKUwTeT3p7CQoc8PPzAP3byZfZMrKhjmxMQwcxk+nIOVHj348nHta1FYyKlkAGc+uX4qXZs2PUckJ4eVQuvXcxrZRRfVfJ/ycuDLL1n61LYtcMMNnFJWW04n3z/S0vjNOeEEdoz15XyNwkK+f+XmMoxS8+lasU6VrYJPb17KQJN7OUszVVbGVcysgGjjRvdl7mNjGRD16sWCx5QUvj0kJPA9uKr366Ag94BIy92LHM5nq4zVFwVCIiI1yMnhQKqoiGc/CQn8+G3xYk7aX76cg0iHgz1F7rhD/TqkWTJNhj1r1gDffceXwI4dzEorKpgtxMYCxxwDnH46pzN068bcxHVAUXnlG9e+FQ268k1xMSt/fvsNuOwyYMIE7z4CX78eePttfnEXXMDkqy6jotxchssZGUzSTj6ZPYB8qbiY72fZ2Xyfio3le5o+6hdpMUpLgaVL7SbVGzfafdlCQvj206EDA6IePexm/omJrPKsPM3MtaKz8jQznQ5JS6dASESkKbJGvvv28SPe4GCeCQUF8eO1LVtYf71tG6eLBQYCjz7KviDSfM2fz/2wYe6Xr1rFfVW9ZKqqLPIkL4/7Vq1qf3w+ZpqcabR8OQcOVrub/HxeFxDAQcJxxwGnnsopYB06uDctLSxsYoOFigquOvbBByxtevppNtyoSWEh8P77XEWsRw82qm7dum7HsG0bR2mFhSxJGTiQ3yxfKimxgyGAI8DERL7XyWHefpv7007j3puXMtCoXs4iVSotZRXc3Ll869mwge/fTidD++hoFj92785wqGtXFkVaixkahufl7gFWILm+7yt7lpZGgZCISFNSUcGGJ/v3s8Y6NJRnPRERbAC7ciUrhrZt48h4yxYO+iZP5rrY0ryNG8f99Onulz//PPcTJ3q+36xZ3J97bs3PsWkT99271/74jlBpKdvYLF4MzJvHWUyZmTzBNwzmnklJbEw6YgQHx23aMLewKoAKC+3VvypP/3Jf/r2R++UXhrwA8MgjTLxqYprA779zFTJ/f65KdvzxdXt+h4PL069cyW/agAEsvfJ1V9eyMgZDBw7w+KOj7Y7d8j8nncT9Qw9x781LGWjQl7NInZWV8a3np5841WzdOrZPdDh4WtSqFd8munZlG7TOnfm3IT6eubJrXzjrb4M1zSww0P5AQKuZSUvgbSCkCdwiIg2ptBTIyuKgqKKCZykpKTzr2bqV82Py8xkWbd3KeTKZmUCfPqwg6Ny5ob8CkVrLy2Pos3AhZypt3sy8s6iIuUNwMCshBg7kLKihQ3mZdZKflcVFrACe0IeFcUBgnew36PSvIzV8OKuE/vEP4N57gSuuAO68s/ovyjA4zatbN+DNNxkcDhoEXH557ZvIBAQweevenQnd0qWc1zFoEN+bfCUoiI/Xti3f07KyOHUtPJw/zNatNVoTaWGCgtjK7IQT+P/ycoZC8+bZLRM3bGBxbFAQ19aIi7PDoeRkblYfopQUnlq5fnCQk8PHtlaOtLYGnzos0kAUCImINISCAlYD5ebyrKR1aw6CwsMZDn39NT89dzpZFZSezqliZWXA3/7GHiOaYiFNgNPJmWsrVrCIZflyYNcu5pylpTwBDwtjxjl0KAtiBgzg/awKoB073D/lDQ/nICA8vB5W/2oM2rVjsPPSS8CHH/Ijc2+mkMXHA/fdB3z7LcPkzZu5PH1dlpWPiuJcvN27+YP74QeOrgYP5nW+YpV/tW3L4Dszk+91e/bwh6wl60VarMBArkh27LHsve908m1t4ULg119ZzLhzJzPrgAB7Zch27fi2l5DADxcSE/n2mJTEvxuuVUSZmXZ1aVCQHQ5Zf1+US0tzp0BIRORoMU1+NLV/P89GAgJ4lhIXx7OQoiI2lV2/nteVlrKLbmYmR8QxMcC//nV4HxmRRqSwkPmFFQCtX89xfkEBy/6DglgAN2QIi2FOOYWf7loNoAsLecIP2NU/sbH2SXqL6QMRGMipgAMGAI89Blx1FaeQDR9e/f38/YHzzmPC9uabwLPP8pt94YV16weUnAxccgnfi5YvBz75hI993HG+XYrNz88OgPLy+D65Zw/Xno6N5WhOIbhIi+bnx1ZpPXow6zZNZtZLl3K2bVoae+Pv2MHQKCKC+XVCAnsQtW3LsKhdO7sPUUqK/QGEpyoia8l7q5KoxfwNkhZDgZCISH1zODgdIiuLo96QEJ6BxMTw7KakhFMz1q1jbXObNhxRp6dzULRvHz+Vf+wxDoxEGgmrgO3PPxkApaXxRPzQIeabfn52X/SRIxn+DBjA/xcX88S7qIiPAdifzjbr6p/aOvVUTt964AHgnnuAK6/kioI1zW3o0oUB8qxZwM8/8wc0ZgynltX2I28/PzYu79oVWLaM4dD69UDfvvzo3tdBTatW3IqL+R6YlcVgPDqavxy+rFASkSbLMHg6lZICXHwxL8vMZI/9hQvttTgyM/kWaIU7cXEMhZKTeVqVlGQHRHFxfPuxVjSzpinv38/HDwhwn2oWHq4iRmna1FRaRKS+HDpk98UwTftjKmswU1LCifBr1jAI6tCBZxwLF3IgZE0RGz+eg8AWPzJuwawz0YQE98utjzGrWlWqoIB7bxr1lpVxX83HnwcO8ER71SqeXK9fz8Wiiop48hwUxKKRLl2YOwwaxE9yw8Ls8Mfh4GP5+R1+Uq3+DdUoK+MUso8+qt0qZADD5Q8/5MioUye+nxxJP6DcXP4CbNnCH3q/fgyH6uuj8/JyO1S3ysxiY7k181+ajRu5T0ri3tue2168nEVahKIinmYtXcrTq1WreHpWUmI3mo6K4gcVVjBkhURxcXZIFBzMx7IqWUtK7OcICbH/joWFaaqZNA5aZUxEpCGUl3N+zIEDHCEHBLASKDbWbu5aVsYzktWrefvOnTnI+flnDvDLynj20qED8OSTQM+eDfs1SYtUXGz/mq5ezSqgvXvtYMcweIIcHc0ikSFDOJPIWha7sNAelAIMiqyT5Sa38ldj8vPPwOOP8z1j3DiGOwFeFHybJj8u//RTBoWnnAKcfz5/IHWVnW2XhYWE8BehTx/vjqcuTJNhVFYWR3SGwY/yVTUkIl6qqGCW/eefdl+7jAyespWX29PDYmI4xaxtW57CxcUxLLJmtrZpw8og15DI+sDDMOy/edbfPf3Nk6NNgZCIyNGUn89BSl4eBy2RkTyDiI62K3vKyzmyXrWKI+VOnXjmMXcuP8GPjORKYrt3AxdcwP4hvuzRIU3XnDncjxrlfrn1dzS1ir/3W7dyX0NT4bIyYO2ifKxe54cNuyOwciUbdVrl8uXlDH9CQvhQgwZxJfLOnXnSXFLi/mlpcDBPgF1PhlXg5kP79gGTJgHz5/MH8uCDDGO8UVTEpvW//ML3nzFj+AM9kpFKVhZ/F3fv5g+7f38G2fU5j8JaofHgQbtqKC6Ov5DNqGrolVe4P/NM7r3tD15T8aCI2A4e5OnZihWcwb92Lf/+WZ/rRUbaCyBaoVBMzOEhUatWfCu1PjgpKmIABfBvoPV30dqrLZrUJwVCIiL1rfz/2zvz6LiuOs9/b1WpqrTLkm1ZsiQvQd4mJnFCcEIgCYSwOCRhaMLS0ITAjJumT3egw5npTM6ZpjnkwOkm3TD00D1OBpJApucMTMgCYWi2LIYAIQuJsWN5t2VLsqzFsl2Sanl3/vjWz/dVqapUkiVXWfX7nHPPq3r1lvuqbt3le3+/303QEujECY6oc1kDyXE7dzIu0MQETSiWL+c6qjt2sKcxOcnZ+/p64O67gbe9rXTPpZQfW7dyu21b5v577+X2zjtzn/fEE9zedNPZXckkl+199VW6fO3YwRg+p4bjiI0HMJEMIRRiR7W5Gbj8csb9uegiehkZkyn+yCphfvFnvgxElCyefprC0MAAReS/+AuOSIrhyBEub79/P3/cP/5joKPj3PLT388YQ319FJsuvZTxj+azQEiw/hMnnNVQUxPr4oaGC35K/sorub37bm59f+WC9PRwu2bN3OdJURY6yST/Q+Jq9vzzjHE/OcnunrR5dXVuBbPFi9lmijgkXq11dZmWROPjbtXMYNC5mEmKRC74akspE4oVhLTLpiiKMhNSKQ4+hoc5+AAo4nR0cBDib8VPn2Zv4rXX2IPo7GRAleeeA77/fbb6ra1cWezMGUZE/NSn1PVBmTPicaCnN4qXv8OVu3btotGQBH0eH2eRjUaBujBw1aZxXHZNPVavpvjT1MQOsOB57LA2N7vO6wIyxrjwuPZa4IorgPvuAx5+mFY/n/kMcOON048oOju5RP2vfgU88gjwxS8Cb3wjz82OVVUsy5ZRsTh6lCOo7du5Xb+eMYbOxT0tH8awQDY3U6k8cYLT/SMjFKIWLaI4NJsV1hRFqUhCIYZq27AB+MAHuG90lJZDO3awWvv972mk2NvLrqG4mtXXuzhEUv2INdHixa46qqpiGxyL8ToiEoklkT+pu5kyn6ggpCiKMh2eR1ew4WHnEhaJMKBrc/NUm9+BAZpfHDjA96tWcUD0yivAV79Ki6GVK2mb/Mtf0gTjc58DurvP+6MpCwfRH3ftovizezfQ//LlGBsPYyRAYUeKbjTKju6mTXT7WrEC6Bg/CFhgckU9AHrgRKPsvKr4U8bU1AB33AFs2QJ86UvA5z9Pl7C77mLdUwhjGP370kuBH/2IgtJvfwts3kxhaOnS2eVp+XKm/n4XgOqVV2iJtHHj/K2WGI1SnF++nG68Q0NMg4PO5K25OdOCU1EUpQiamlhdXn018Kd/yvb06FGKRK++SpFo507uO3CAXUcRierqKBK1tFAkkqpIBKLmZh4XjbKLGItR2xaRSGIS1dRkbnV1M2UuUEFIURQlF9ZyhD00xGmhVIqjYYkkmD3b7HnsAbz6Ktc3lZV31qzhvq9/nTPWF13E6IXf+x57B1/+MnD99Tr1oxSN57HDuXMnULudRfT+p+gBJDELxKqn2jYiUuVh3ToaaKxaRcOQtrbM1YeiUaChnx3McLe6fV2QdHcD998PPPoo65sPfxj42MeAT35y+kAVtbW0UHzHO4Af/5iuaH5haMmS2eVp2TKmsTGqlbt3U61sb6cwJD6Ic40Em25sZN09OkpBv6+PqabGjcxU5VQUZRYYQ/25owN45zu5L5mkF+5rr7HKe+klvhaRyFonEtXUUHMX79bmZopOfpGovp7tsjG0JhodpVAkSLw+v0ikKwsqM0VjCCmKogjW0nVrZIQpkeD0S1OTa5mzBy+TkzTJ+MMfeG5jI0fey5fTFePnP6d/TkcHr/XkkxzRf+xjwG23adBopSCxGMfPdX+1FSfHgG9csg1791IEisWAjhMvYSi1CIO1KxEOs3O4ciWtf1auBNYMPIPORaeRvGELAJqiV1e7zqN0IAMBaNCRhcTwMK0Rn3ySgswnPkFXrmLFj7ExJwylUgxks2XL7IUhIR53I6XTp1lfbtw4/3GGhESC383wsFsqr66OwlBTU1mKQxpDSFEubEQk2rWLEzkvvshqMBZjlZhKOYGoutq5l9XXs2qSlC0ShcM81z8JBLCrKe28P6k1UeWhQaUVRVGKwfM4+BkdpTtYMulml5ubuc21PFJ/P3vce/fynOXLObCpqaEItH07W+gNGzjgeOIJOpq/9a3AZz/LGXJFSTM5CRw+TAOKffu4wteBAywyZ84AodOjiMeBYa8JwSCFn+pqoHtlAq+7yEPH6gja2qg7SpiWSASoCUwgGgWqF0VRU8MOZF6DDFkvV02DFg6/+x2XqdqxgxaJt98O3Hxz8VPIJ086YcjzuBrZli3n7vIlFpWvvEJ3rnCY1pNr1sw+ftFMmZhwFqASJb22lsJQU1PZuJX19nIrX3mx2dK/s6KUL4kEcPAgu5GyyMMf/sDqKB7n/zcSYZUky9fLKmYNDayi/K5njY1OUKqq4jXGx90KZwCrWTkmGnVbNVBfuKggpCiKko9kkgOd0VGKQZ7HqZPGRrayDQ25p1JOnaK5Rk8PzwuFXEyMWIxLgz//PI/dtIkWRz/8IV3EVq/mSlCbN5/fZ1XKislJeqzs2cMZw4MHKf4cPswiOT7OsWkqxRQOMzU1cay8ciW1RAmTEgrpbKAyDdZyHeX77qMAs3Qp8PGPc1WyYoWh0VEKQ888w/rysssY0Lq7+9xHE/39nDo/cIB1sxT27u7zFwh6YoLPODLCuhzgSEnEIQ1IrSjKPGMt+wc9PZwckthEx45RQIrHWd36RZ2mJlbpomWLp2xzMwWj+noeF41SKJqcdPEEAV5P4gqqULTwUEFIURTFz8SEswQ6fZqtYVWV6/DncgcD2Arv388Wuq+P+9rb3eh8/34OlHbuZKt6+eW8z+OPc3CxcSMHX295S25LI2VBEosxZsDevSwihw9zpr+31xkkTEyweKVSFG8kiPPKlQzy3N7O1NkJXLTzCQQCQOJdN50VfKJRoOblX9EC4E1vyp2RmfiNDA1x29IyF1+BUm5YS8F62zYGeV68mG6r73vf9DGGhNFR4Cc/oTtsLMZgVNdeS6H7XFcQi8cpCu3eTZHIGKqea9fyD3G+TF3icT5nrraioYFtxXlUW//u77h973u5LdYFTP/OirIwGBtjX2LvXk4m7dxJq6JTp1hdJRKsosSaKBJhdSULQohIJPuamlwfQoQiz3OGkkCmUJSddLLpwkEFIUVRKptkkq3l2BhTPM79xcz6ytIRPT004Ugm2ZrKrLXnAb/5DZePP3qUg4TLLmNU3x/8gAOlq6/mYGvTJp1mWaAkkxwzinvXoUOcyevt5fbkSXawJied8BMIUPipquKs3qpVbkGmzk4KQP4ZOkl1f7WVnbBt2zIzce+93N55Z+5MPvEEt8UEHtGgI5WBtXQlu+8+BrNobmZd9Ud/VLw/UjzOazz9NOvIcBi44gqKQytWnHsex8ZYHnt6KMqEw8DrXseyOduVz2ZDLmtSY9h2NDQ4P415RGMIKYqSjedRNxeRqKeHQtGBA67PkUg4N7FolAJPYyPnAmpqGM1AXNBaWlx1JseGwxR/4nFnUQSw/5JLKNJg1uVHsYKQehYrirIwkFXBRAASs/9gkDO6y5ax1cs3E55KcRR/6BAHOLEYW7fubvaom5uB3/8e+Na32Op6Hk05rr+edr1f+xr3veMdDBitvfAFgedRVzx+3Ll4ieBz7BhX+5CAjuL3LxY/VVUsQp2dXEyprc2JP+3tLkxJdsoZ11Zn5JS5whiKN1dcQUFo2zbgH/8ReOABqg0338y6rRDhMK3S3vQm1pmyKtkvf8lzr70WeMMbZj9CaGjg+ZdfTsvM3bvdiKe2lqKT+E/Op+VlKMSRUkvL1DZGKoFQiG2MjKzKMDC1oigLi0DAWRFfc43bn0hwbvLAAcYjlGDWBw64akssiqTPISuVNTezChO9W2IVNTfz80jEpWSS8Q09LzNP8rlcV7ZaLZY3KggpinJhYi1H4qdPc8R+6lTm7G17O1uzmpr8FjoSyffQIbagiQQ7952djA3U1cXPf/pTzobHYmwd3/pW3uu554BHHmFL9773AR/9qAaLvgARY7KBAXaajhxhp6m/n2lwkJ+L771Y+4gnSSjEDs9FF7kOmgg/HR304/d3juS1ehAqJeeyy4B/+ReujfztbwPf+Q7w0EPAJZcAt9wCvP3t01vArFhBEfz976fl5NNPAw8+CHz3u07UWbNmdgXeGPenuvpqKrISiXXnTgpOXV0Uhzo753fUYQyFn/p6/rmTSScOjY3RRRjgH7y+ntPvdXU6ba4oynmjqoohK1ev5nylkEzSoH3fPvZzDh6kZdG+fTSAHBzka8/L7KeEw04Yqq1llSbV4KJF1MHFTU2En6oqTor5u96yGEaupFVk6VFBSFGUC4NUiuKPpFjMTU1EIpzBLSa+w6lTzgqor4+j+poauiOsWMGO/ugoBzYPPECVoKoKuPRSXn/XLq7aE4vRdeGTnwQ+8AFOoShlibWcyTp5ktY9R48y9ffT8mdggJ2h06czLX2SSRaxUIipqoqm1h0dHJ+2tTmrn87OTFNrEX90hR/lgmDTJqYTJ7hU/aOPAl/4AvCVrwA33EBxaOPGwu6vNTUUy6+7jiONZ55hMOtnnuEo4tJLKUCtWze7IBThMIWlNWvc6Eaisu/dS8Fp+XIXhGueXbkQCrklfgBGhBdxaHiYlYrkW8ShujqOnhRFUc4joRCrxWyPXmsZb+zQIdc1PnCAmntvL6uz/n63aqH0bWTBi7o6CkMiFvmTWBaFw5likSRB4hVli0R+tzVlftGuqqIo5cnkJEfxIgCNj3O/LLGwZIlrgQrNCk9OsjU7dowDiOFh7l+0iLPgK1dSTDp4kO4TDz7o1vnt7gYuvpgmI488QvWgpoaz5lu2cHCjZh4lx/Ooz506xZ+6t5fbvj4n9gwOUucTwUf861MpJ/qItU9LC8eVbW1cAbu1lQJQVxc/y+64qOijLBgWL6a1z5/8CVcke+wxrp742GMMeHXzzcCNNxYWwI1xwk08TpfaF15gQOvt21mHXnIJ688NG2b3B/KPbqzlH12sh559lqmlxVkXtbXN/zS0RHtvbWWexsdd+3XqlGt7QiHXdtXW8vvQEY+iKCXAGFb7ixfTmNNPPM5+1JEjLh06REui3l5WbQMDFIvEYtovFuWyLqqry4xfJMdLXyoUchNrQijkrpctFqlgNDdoUGlFUUrP5CRH9GfOcBuLcaQOUHDxTznU1hYWYaQF6+tzQV4Athitrc69IBiky8GrrwI7dvDegQD9fjo72Xl/9llOkwQCjJWxZQtjYxS7Io8yJ0xM8OcZHeVP2tfnrHtOnOB2eJgWQLkEH4npI52NSITGXcuWcSuCT0cHi0dT09SOR8njgsvyH9lBfyVYer7BrkzrFTPoFos7FTkVP7EYVxZ77DGKRMEgRw5vfjNTV1dx10kkaGH5wguMxzY+zvL8+tfTemjtWtbx58rwMEctR49ytJJKuVGPCETLlp3/oBaTk5lWrv4lfaJRjpJqatxSQelRztgYDxGDp2L1M/07K4oyXyST7IeJUCSrqB48yPdjY87SOpl0C2r4LYXEwkiM+6WLH41yK6uhyUIcMnEnMY+kLpRVWv1Jri+p5H24EqGrv9KlrgAAGa9JREFUjCmKUn5Y68QffxLxxxjW/tkd40I1+cQEFQEJ8Dk0xPuIANTWxgHAkiUcHLz6KtO+fTyuro4DkXCYZiQvvkh3B4Cz1zfeSJcJdQmbc8SyJxZj+A0RegYHKfQMDXE7PEwxKBaj/iEuXYkEryE/t4g+VVU0DhALHxF/RPRpa8vsZGg8H0Upkv37gR/+kGL5/v3c19npxKFNm4qzxEkmGSj6hReAl1+m4gvwD7p2LVN397m7faVSrPelfTh+nJVGIMA2QcShJUuKX2FtrpCorP62UARewIlE0g5WV6s5oqIoZY+1NIoU93wx0JeQnceOcT5AxCKx1PaLOBKHyG9MKQGvxYJIRKOaGicYiTWRPw6SMZkiVL7XC9HSSAUhRVFKhwg/ExOs9cfH+Xpiwq1dKa5f/jSd+DM56dQC8QM6fZqfBQLO1KO9nbaoR45Q+JGlFmTlsY4ODgBiMX720kvsiIdCnKnevJmxMKZbaUfJSSLBnzwWo5AzMOCseYaHmYaGmMbG2HEYH3eWPamUE3s8z4k9IvjU1zvLniVLnOAj2l9ra2ZAxFBoAcwOffe73N56a+b+p57i9rrrcp+3cye3GzZMfw+JebJkyUxzp1Qix47RBWz7dgbdj8dZj7/xjcBb3sIg0IsXT38dz+O08muvUSTat48VgDEUm9ato0D0utedu2gj09oiEA0Oujapvp5lX9Lixec/2mkikSkQnTmDz/9TCwDgAzeMAsEgNlwc4PcgIlGeYGX6d1YUpRyxlhbdsnCHJBGNpGpOpZxg5HdLEwEnFGK/0G9g2dhIAclv4S1BsMWqSMSjXCuhSVq9emGsjKaCkKIo808y6ZZempzMLfwArsaVzqtE353O8mdoKFMAOnXKfd7QwA67dNxDIbYkIgD19jqb+WXLGDMoleJ1Xn7ZuZKtWgVceSXTZZdpwM8cWOt+VhF5xIJnaIjWPSMjbOAHB7mVsBmyKpcIPdK4W0sNLxjkVmZ2mprcT9rayq2IP62t/CnFpSs7MOGCZutWbrdty9x/773c3nln7vOeeILbm26a/h49PdyuWTPz/CmVzfi4ixG0fTsVYIC96osvZkDqjRv5fjpzvGSSUU1376ZIdOCA8zkQl9+uLqa2tnOzmkkk3OSCJH8709TkKqHFi/n+PFsSXbmZyvjdf3kKiCdw09XDrIzFshZwwTt8y/30HGZwjjVrL3Q1XFGUSsNfNR8/7rZibdTXx31SFUoSL3m/uCOTiWI5VF2dGeNfBKbqameFdPvtc+PBXGqKFYTU9lRRlMLE45mijz/5O6SAE34kUpwIQPkGALL808gIVYbRUffaH19BZm7Xr6ewk0xSkTh2jIOGw4d5DsBavbmZysHEBI975hnXyW9ocALQ5s1UGSoIz+NPKgLPqVNOzPF//WNjTvgZGXFxvc+cce5asgqX3+Q3EHBCj6RwmOOplpbMsZUksfJpbXWN84Kw6lGUSqG6GrjmGiZruerX9u0U3596Cnj8cR5XU0NrtY0bnVCU7Y4bCtFdrLsbeM97WNns28e6fs8erlwmlnGhECPAd3XRmqirixagxSrFVVXOqlSQdkNGI319fB4hGmU7tGgRBSLZ1tbO9tsrjAkAwQDQtIjv1y3lNh53kzCyHRpy7fKhdKy7BPKv9bwQfSQURbngyVU1Z2Nt5gSlbAcHaZkuC4ucOMG+65kzmcKRhJfzu5tJv/Xd76ZRaqWggpCiVDLJpAvKIv462e/9lj7GOBvM2tqpHcxcwo/nOXORU6f4WsSf0VEn5wPsaDc1cQa4ocGZpgwNcUDw7LMuBoQgcRWiUSoXe/a4KJwysLjhBopJGzbw/QUeMEZEHdHqxGpHVjweG8v8ykdHnavW6dNsFEUQEmHHP8PieZlWPH5LHtmK73ZzsxN2ZAXm5maKPy0tFHqamxeQ65aiKIUxxgk6ACuSI0cYvH/HDsZwe+ghJ1y0t7N+XrWKdf+qVRR2xFqzqopuY+vWuesNDnIiQNKLL7J9AFhBiQ+pKM1LlzLV1U1fCUWjFJU6Oty+WIztkF8137ePFbAQDrNSbGpi+1VX53wVamvnvvKTYBuNjZn7k0lW7sl0e14bYz5HRjLbW0CjsSqKcsFijNPmpbnJx/h4pkW79ImHh90iJYODLmblokXn5xnKBRWEFGWh4XlulC++Ov7X/vd+YQVwgk9VFTuz/nUeZUYxu3MormIjI5mij7w+cyZTVAJ47aYmzuxay4HB5CSP37PH1doSRCYWc9J9KsX3J0/yGLl2MEg5//rrObhYv54rhp3vGBAFyP4Z/KsSy+zF6dNuwbWxMTf2kK9VLHUmJjJX1PILOfJa3hvDMZIxTtCRFApx0r6+nmMYEXQWL3YNraTmZreVRXB0rKAoSkGMcS5eW7Zw3+Qk3cFklcfXXgN+8YvMNmnZMicSiVC0YgUrIRF43pC2hLeWPfkjRygQydI3L72Uec2aGneuiEUtLbxmY2N+ixkJUNHZmblf2j4RikZG6K4s8eoE/2qZfrFIlp2vqZm7tioU4vVlQLPK95nfzTt7EkjMP3NdT3wv/K+z32vAa0VRypTq6qk6fz48r/L6tlp7K0q54h/di32jP7pavvfZIo/gj6RWWzt1FlA6dZ7nOoyy3veJEy7IpUQLltfZbmP+JZ+M4f3EryiR4HWPHGHnXXyXROEQlUJELVFF/MpDXR0HFmvWODeBlSspBp1jh1qyKV+r37VKHtefXfl6Tp7M1ML8MUFziTeySpb8tKKJiYgjyVo+tiT5evyvA+n4ojK2aGjguKahwU1Wy4R1UxM/k21zM89TYUdRlPNKJAJccgmTEI+zbTh40KUDByjq+F2Iq6roZ7psmRN25HVrK1c7a2pi5ZhKcYLh+HEX3X5ggNY9zz8/1QK2vn6qK5i8lmiltbXOylQCO2f7NaRSbmLk9OlMs83Dh6cKRgArYv8CC/J6fD0QMMDIOJ99PDj7pRHFLyKfe5uYn/qFIv8shgSFK9TP8Pte+N9n+2T4kzZAiqKUCRe4E8GsKCqotDHmXQC+BiAI4H5r7ZezPo8AeAjA5QCGAHzQWnuw0DU1qLSyoLB26mjeb6JR6LVf9PGn6f6buTpdwaA7z9rM19Kh83f2RKkQ8Uc6e9Y6ZcTfMRRlwn99eQZRR8bGcruhiamKqB/pZaVsMAgbZP4taMLitSxBatESJFqWYXJpJyaWdmF8SRfONHdiPNKEybg5m/1YzG3PnHEuVP5HEjHGH/7I/1VkG035xRoRZgpt/X1ZYwAkPdi4B1j244PVAQSjgbNCTijkjK5kPCGrC4ulTmOjG5vU17txiN8DQYSfUoaBGHh4AD139CA1RGEw1BJC99e60fqR0sVmGnh4APvv3o/JQ5NstVJwWwMg/bc417yevc/hSUS6Ilh9z+op1/IfE2wOwsAgOZzMe7yck+s7BTDt/Qrlq5j8ztf3MN3xAIoqR9nntmxpwdCTQ3P+TOf6jBcKM3muYstyvt931t+f59GW/+BBrlk8MODEnf5+xHceRfxIDDbuwYQDiCwPo2pptRN0RAXPfl1b6yp8/wSIP55dLuEGYIUtFj+S/Osj+1cAk+SPqZdMTp09yJ5wkYYMwNhvxjD46CCSw0mEmkNY8t4laHjLYteQyDI5YsnrX1M53zYUwsD/Op73dylYNlIpDDzUi4N/sxeJIzFEOgJoflsDRn92Aomj44i0B9F5RzsWv7up8ESVkB14zp/8vsqFXmenWZZfRen5dA+ObTvGPkOayIrzU04WalujlIY5W2XMGBME0APgBgC9AJ4H8GFr7U7fMZ8G8Hpr7aeMMR8C8O+ttR8sdF0VhJQM8okY0+3L9T57n4zcp0vZI/1cr0UhSCZhk3xtE1wL0XpM8Dx41sCmPNiUBy9leRnPwkt6SCUtP4dBKmXh2QCSnkEyYeGlLFIpi0QSfJ0AUp4HL2WQSnhIxvl5MuHBxlNITKZgE0kk4yl48SRswkMqkUQqYWETKdikh1TSg02lYFO8v5fy4HkW1gOs58Gm0nlPWcBLpV97sACsF4CFBSyPhQfAivkKfJ9ZWASQRAgpE0QKQSQR4nuEkEQQE4EaJlQjFqjDGVOLmKnDuKlJp2pMmGrEbRieNTl/1lw/cS6yJxv9Opbf4ibb+sY/iek3oJK+tn+RNH8Skaa+Hki8MIrh/3EE4XgSEXiIwEM0YvHvvrISqz62FLW1CyeO58DDA9h1+y4GLfVhwgbrvrmuJJ2YgYcHsHvrbnixaQYgaWab11z3CdQEsHbb2owBVaG8ZB8v5+T6ThEETNDAxu205+fK17LblqH/wf6C+Z0NxXwP0x2PKrBuyTI2zP5tivlt5+KZisnzfNznfDOT5yq2LAPI+fsaM33ZnfUz/MddCIyPIYwRpvAoOj5UhYaLUhR2JGL+yZNMhfq9EiBNFHpZ0l1MM/39Ar/1rswsiIVsPnPOQMA1JrLapog42aJOOujz6NMjOPqNXiCeAtXsFIJhoO3jrWjcXJvbglieJZ95qTEY++0p9D98Aqm4gUUQFkGYcBXa/7wLNhBE7z/1IzVpAARhEYCJhNB51yq03LgEQz8excEvHoE3gfS5BhYBAAFO7CAAEw3ior/vxpJb07+xv49VaIYl+7X0vWaK77lPPDmMg397GN6ETefPpPMKAAYmGsTqL1+EJe9rndoxyNdpyJXkvoX2y+tC+/xb5bzS8+keHPvnYzk/m+96f6G2NUrpmEtB6CoAn7fWvjP9/i4AsNZ+yXfMj9PHPGeMCQHoB7DEFrj4QhGEbo08iqQtzvOuCGOsGcFGbe72+T4sYpfJ8ZnJd3rGPW2uc232PpN1ztTryD1t1v7M62e9tlPvL50De9Z0wF3H5riuPXts7vvKfbyz7wMZ14PvfNmX/9uc+o1kP3nm/hl2IEzuO/oPsDDwTAAWAW5NAB6C3Kb3W2PS74PwTBAmEIANuBk+GwwiEDAZfWKZ/JO+tt+S3C/MZMe4rKqaGstaJl+lfy3JH5pB+tj5JlTnykT0uZXP0TIli8iKCK46eNXc3KRMyPesQOmet1Ce8jGbvBbzOxeTl+x7zzT/RZ8vFlLTnD9TZlrez+X5ij13rsveQv1Pz+S5ii3LAIr+fefi+5vxb+N5dNsSkWh0NDOIm5iaytb/enw808w0V7ydbCEjn8iR/Tp7tsN3vfhgEjaV2RMAABMMILI8klt08F8n1wyKtZg4Ogmk3D55FQiyL2Sn1BcWJmgQaQtjsi8Ob8rnU/sQJmgQ7Yxm7SzQT8n12XTCiPhWZ+/zMX5oAjZVqOfLZ4t2RfMekZfzKdxMd69yEpHKKS8zILZnvODnJmRQvWoW5aQIxg9MwCZzjLjm8Z5KHn70I8YhvcCZy2XnlwM44nvfC2BzvmOstUljzEkALQBOZGVqK4CtANDV1VXErcuflBdE0s69s6HB1AqhmOOMT3aYzfV4bHH5cXW9xVRBYqrM4BcuTLawke1y4zsme+uulpZyTLYk5AjAAiZbhsm8b4ZsY/x55LWN73hj3Pbs4cZJSjDm7DnGyP25D4buQ+yvUQYygfS+ADMSCBgYSQbpz3liIGBgggaBIGACAQSCAZhQACZgEAwZmED6fTiEQFUIgUh6Gw4hGA4hGA7ChIKoigRgQkEEI0FUVVchlE7B6ipUVYcQrA6jKhpCsCaMcDSYYTHjt+LODgUgok4l+t76mTycezCUb/+FTKFnKtXzzua+c3mOf38x180+ZqZ5Kfr8PBPs5/o7zbS8n8vzFXvuXJe9hfqfnslzzaYsz/XxM7lG3msHAnQba2xkDLpzwR9vLzuJ1ZB/m73PvxqAuFVlL/eYTKL3c68BSMEgBQPPbVMpdN7SlrmGslwnnxu7WDZZi9OHB339Kg9n+1ApgP0gL6vfZWFSQGR1I8Z7R+GfNuP5cNeQ9ykg2lw31bS3kGCVb5t9fL7zs7EWNiUTc+5ZphyWgjPfzXWtfPtyPUuefMya+bruud57ARLANPVSEsD4/HwnJjmZW7Scx3sqecgl+C9gihGEcpXN7FJZzDGw1m4DsA2ghVAR9y57HkncVOosKIpSRkS6IrlnrLsiJcjN/JLvWeWzUlAoT4XOmav7+K9VTF6y7z3T/Bd9fj4LoXP8nWZa3s/l+Yo9d67L3kL9T8/kuWZSlou2EJqD76+kv00g4GIFzSMDX89vBdX532ZvYbWngHUVkPt3jKyI4Kqnr8KumVjrvVB6K7rfF2utuaf0eVVKy69DT+WdQAHm1zL0pQVqjaqUP8XM5fcC8K+z2QEg27ny7DFpl7FGAMNzkUFFUZQLidX3rEagJrNqDdQEzgZWXUisvmc1479kYcKmZM+b6/svxGzzWszvPF1ecpWLfN8pgsxrMefnylf71vZ5KZczLe85v5MqULDKIvu3Kea3nY//2kL9T8/kuYoty/l+32LK7mxYqL+Nn/l6xkLXne6epfovzpbZ1MVKZdK+tT3vZ/NdTiqhPlPKk2J6zc8D6DbGrDLGhAF8CMDjWcc8DuC29Ov3A/h5ofhBiqIoC5XWj7Ri7ba1nGU1nNlZqAEBWz/SivXfWo9gixvNh1pCJQsoLXk6+/0DTmiQrW9cei55LeZ3zj4m2BJEqCVUsFzk+07XP7ge6765btpylS9fa76xZl7K5UzLe67j139rPdY/OH05ynVu+5+1z/t/baH+p2fyXMWW5Xy/bzFld76f4UJlvp6x0HWnu2ep/ouzZTZ1sVKZrPnGGrT/WfuUSYrzUU4qoT5TypNil53fAuCr4N/jm9bae4wxXwDwO2vt48aYKIBvA9gEWgZ9yFq7v9A1F0pQaUVRFEVRFEVRFEVRlHJhLoNKw1r7JIAns/b9V9/rCQC3zjSTiqIoiqIoiqIoiqIoyvmnwtcDUhRFURRFURRFURRFqTxUEFIURVEURVEURVEURakwioohNC83NmYQwKGS3FyZTxYDOFHqTChli5YPJR9aNpR8aNlQCqHlQ8mHlg0lH1o2lHwspLKxwlq7ZLqDSiYIKQsTY8zviglepVQmWj6UfGjZUPKhZUMphJYPJR9aNpR8aNlQ8lGJZUNdxhRFURRFURRFURRFUSoMFYQURVEURVEURVEURVEqDBWElLlmW6kzoJQ1Wj6UfGjZUPKhZUMphJYPJR9aNpR8aNlQ8lFxZUNjCCmKoiiKoiiKoiiKolQYaiGkKIqiKIqiKIqiKIpSYaggpMw5xpi/N8a8Zox5xRjzfWNMU6nzpJQWY8y7jDG7jTF7jTF/Xer8KOWDMabTGPMLY8wuY8wfjDF3lDpPSnlhjAkaY14yxvyg1HlRygdjTJMx5nvp/sYuY8xVpc6TUh4YYz6bbk92GGP+1RgTLXWelNJhjPmmMea4MWaHb1+zMeYnxpg96e2iUuZRKQ15ykbFjWNVEFLmg58AuNha+3oAPQDuKnF+lBJijAkC+O8A3g1gA4APG2M2lDZXShmRBHCntXY9gCsB/LmWDyWLOwDsKnUmlLLjawD+n7V2HYBLoGVEAWCMWQ7gLwG8wVp7MYAggA+VNldKiXkAwLuy9v01gJ9Za7sB/Cz9Xqk8HsDUslFx41gVhJQ5x1r7b9baZPrtrwF0lDI/Ssl5I4C91tr91to4gP8N4JYS50kpE6y1fdbaF9OvT4GDuuWlzZVSLhhjOgDcCOD+UudFKR+MMQ0ArgHwPwHAWhu31o6WNldKGRECUG2MCQGoAXCsxPlRSoi19hkAw1m7bwHwYPr1gwDee14zpZQFucpGJY5jVRBS5ptPAPhRqTOhlJTlAI743vdCB/xKDowxKwFsAvCb0uZEKSO+CuA/AfBKnRGlrFgNYBDAt9LuhPcbY2pLnSml9FhrjwL4CoDDAPoAnLTW/ltpc6WUIa3W2j6AE1MAlpY4P0p5UhHjWBWElFlhjPlp2jc7O93iO+Zu0B3k4dLlVCkDTI59uryhkoExpg7A/wXwGWvtWKnzo5QeY8x7ABy31r5Q6rwoZUcIwGUA/tlauwnAGajLhwIgHQvmFgCrALQDqDXGfLS0uVIU5UKjksaxoVJnQLkwsda+vdDnxpjbALwHwPXWWh38Vza9ADp97zug5tuKD2NMFSgGPWytfaTU+VHKhqsB3GyM2QIgCqDBGPMda60O7pReAL3WWrEm/B5UEFLI2wEcsNYOAoAx5hEAbwLwnZLmSik3BowxbdbaPmNMG4Djpc6QUj5U2jhWLYSUOccY8y4A/xnAzdbaWKnzo5Sc5wF0G2NWGWPCYHDHx0ucJ6VMMMYYMA7ILmvtP5Q6P0r5YK29y1rbYa1dCdYbP1cxSAEAa20/gCPGmLXpXdcD2FnCLCnlw2EAVxpjatLty/XQgOPKVB4HcFv69W0AHithXpQyohLHsaYCRC/lPGOM2QsgAmAovevX1tpPlTBLSolJz/B/FVzt45vW2ntKnCWlTDDGvBnAswBehYsT81+stU+WLldKuWGMuQ7A56y17yl1XpTywBhzKRhsPAxgP4DbrbUjpc2VUg4YY/4WwAdBd4+XAPwHa+1kaXOllApjzL8CuA7AYgADAP4GwKMA/g+ALlBEvNVamx14Wlng5Ckbd6HCxrEqCCmKoiiKoiiKoiiKolQY6jKmKIqiKIqiKIqiKIpSYaggpCiKoiiKoiiKoiiKUmGoIKQoiqIoiqIoiqIoilJhqCCkKIqiKIqiKIqiKIpSYaggpCiKoiiKoiiKoiiKUmGoIKQoiqIoiqIoiqIoilJhqCCkKIqiKIqiKIqiKIpSYaggpCiKoiiKoiiKoiiKUmH8fxcedM4F3UXIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# estimates for the mean\n",
    "red_mean_guess = 1.1\n",
    "blue_mean_guess = 9\n",
    "\n",
    "# estimates for the standard deviation\n",
    "red_std_guess = 2\n",
    "blue_std_guess = 1.7\n",
    "\n",
    "\n",
    "y = [0.001 for i in range(len(both_colours))]\n",
    "plt.plot(both_colours,y,'mo')\n",
    "plt.ylim(-0.02,0.40)\n",
    "\n",
    "for i in range(5):\n",
    "    likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)\n",
    "    likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)\n",
    "    \n",
    "    likelihood_total = likelihood_of_red + likelihood_of_blue\n",
    "\n",
    "    red_weight = likelihood_of_red / likelihood_total\n",
    "    blue_weight = likelihood_of_blue / likelihood_total\n",
    "    \n",
    "    blue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)\n",
    "    red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)\n",
    "\n",
    "    red_mean_guess = estimate_mean(both_colours, red_weight)\n",
    "    blue_mean_guess = estimate_mean(both_colours, blue_weight)\n",
    "    \n",
    "    print(red_mean_guess,blue_mean_guess)\n",
    "    \n",
    "    x = np.linspace(-3, 12, 100)\n",
    "    ymaxr =max(mlab.normpdf(x, red_mean_guess, red_std_guess))/0.39\n",
    "    ymaxb =max(mlab.normpdf(x, blue_mean_guess, blue_std_guess))/0.39\n",
    "\n",
    "    plt.plot(x, mlab.normpdf(x, red_mean_guess, red_std_guess),'r',x, mlab.normpdf(x, blue_mean_guess, blue_std_guess),'b',alpha=0.2*i)\n",
    "    plt.axvline(red_mean_guess, color='r', linestyle='--', ymin=0.067,ymax=ymaxr,alpha=0.2*i)\n",
    "    plt.axvline(blue_mean_guess, color='b', linestyle='--',ymin=0.067, ymax=ymaxb,alpha=0.2*i)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see estimated mean is going near the real mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x20fa5732978>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIQAAADGCAYAAACjDVtCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecFPX9x/H39/YKRRHpwtGUrhTxpGiKRmNs0RiNYtSIJVijRqOxJTHGxJZYYid2jaCxJKgxxfxEEwXlQFCKoGLh6B2lHFe+vz8+N+7esXu3d+zdbHk9H4/vY3ZnZmc+uze7N/OZb3HeewEAAAAAACB35IUdAAAAAAAAAFoWCSEAAAAAAIAcQ0IIAAAAAAAgx5AQAgAAAAAAyDEkhAAAAAAAAHIMCSEAAAAAAIAck1RCyDl3uHNuoXPuI+fclXGWj3fOrXbOza4pZ6c+VAAAAAAAAKRCfkMrOOciku6R9G1JZZJmOOemeO/n11n1ae/9hc0QIwAAAAAAAFIomRpCoyR95L1f7L3fLmmypGObNywAAAAAAAA0l2QSQj0kLYl5XlYzr67jnXPvOeeedc71TEl0AAAAAAAASLkGm4xJcnHm+TrPX5Q0yXtf7pw7V9Jjkr61w4acmyBpgiS1bdt2v0GDBjUyXAAAAAAAACQyc+bMNd77zg2tl0xCqExSbI2fYknLYlfw3q+NefonSTfH25D3fqKkiZJUUlLiS0tLk9g9AAAAAAAAkuGc+yyZ9ZJpMjZDUn/nXF/nXKGkcZKm1NnZHjFPj5G0INlAAQAAAAAA0LIarCHkva90zl0o6Z+SIpIe9t7Pc85dL6nUez9F0kXOuWMkVUpaJ2l8M8YMAAAAAACAneC8r9sdUMugyRgAAAAAAEBqOedmeu9LGlovmSZjAAAAAAAAyCIkhAAAAAAAAHIMCSEAAAAAAIAcQ0IIAAAAAAAgxzQ4yhgAAEiRTZukt9+Wpk2TFiyQhgyRxo6VRo2S2rULOzoAAADkEBJCAAA0B++lDz+05M9bb9l07lyb75zUo4c0ebKt65y0zz6WHBozxqYDBkh5VOQFAABA8yAhBABAqixZIj35pCWApk+X1qyx+bvtZome44+XDjjAagTttpu0YYPVGJo+3RJGTz8tTZxor9l992hyaPx4qWfP0N4WAAAAso/z3oey45KSEl9aWhrKvgEASKnqaum++6Qrr5S+/FIaNMgSOQccYNPBg3es7XPVVTa98cba21m40JJDQZk/X2rbVrr5Zuncc6k1BAAAgHo552Z670saWo8aQgAA7IyFC6Wzz5b+9z/psMOk+++X+vZt+HXTpu04Ly/PkkeDB0tnnmnzPvlEOucc6YILpEmTpAcflAYOTO17AAAAQM7hNiMAAE1RUSHddJM0fLg0b5706KPSP/6RXDKoMfr2lf75T9v+vHm2vxtvtP0DAAAATURCCACAxnr3XWn0aGv2ddRR1qzr9NOtc+jm4Jxtf/586bvfla6+2vohmjWrefYHAACArEdCCACAZG3bJl1zjbT//tKyZdKzz0rPPSd169Yy++/WTfrLX6Tnn5dWrLCk0JVXSlu3tsz+AQAAkDVICAEAkIw335RGjJB+9zvp1FOtts7xxzd9e8XFVpriuONs/+PHW2fTw4dLb7zR9FgAAACQc0gIAQBQH++tidbXv241cf7xD+vPp0OHndvuk09aaardd7cOpl99VaqslL75TastFNLooQAAAMgsJIQAAEjEe+lnP7NOnM88U5o7V/rOd8KOqrZDDpHef1+aMMFqC112GUkhAAAANIhh5wEAiMd7q3Fz223ST34i3XlnajuNvuQSm95xx85vq21bG+6+qEi6/XYpP9+SQ83VyTUAAAAyHgkhAADq8l669lrplluk885LfTJIkmbPTu32nLM4q6qkW2+VIhHr74ikEAAAAOIgIQQAQF2//rUlUyZMkO6+O3OSKs5Jd91lSaGbbrKaQtdfnznxAwAAoMWQEAIAINZvfmMJoTPPlO67T8rLsO728vKke++1pNANN1hS6Fe/CjsqAAAApBkSQgAABG68UfrlL6XTT5f+9KfMSwYF8vKkBx6w0ceuu86aj117bdhRAQAAII2QEAIAQLJ+d66+WjrlFOmhh5o/GTRgQPNuPy/PhqWvqpJ+8QtLCl11VfPuEwAAABmDhBAAALffLl1xhXTSSdKjj1rypLlNnNj8+4hEpEcesaTQ1Vdb87HLL2/+/QIAACDtkRACAOS2u+6SLr1UOuEE6cknLWmSTSIR6bHHLCl0xRX2/NJLw44KAAAAIUuqPrxz7nDn3ELn3EfOuSvrWe8E55x3zpWkLkQAAJrJffdJF10kHXec9NRTLZsMmjDBSkvIz7dk1wknSJddJv3xjy2zXwAAAKStBs98nXMRSfdI+rakMkkznHNTvPfz66y3q6SLJL3dHIECAJBSr7wiXXCB9N3vSpMnSwUFLbv/RYtadn/5+Zb0qqqSLrlE6t9fOuKIlo0BAAAAaSOZGkKjJH3kvV/svd8uabKkY+Os9xtJt0jalsL4AABIvU8/lU49VRo61JJBhYVhR9QyCgqsptCwYfb+P/007IgAAAAQkmQSQj0kLYl5XlYz7yvOuX0l9fTev1TfhpxzE5xzpc650tWrVzc6WAAAdlp5ufSDH9iQ7M89J7VpE3ZELatNG+nZZ62m0AknSNu4jwMAAJCLkkkIuTjz/FcLncuTdLukyxrakPd+ove+xHtf0rlz5+SjBAAgVS65RCotlR5/XOrXL+xowtGvn3U0PXOmdPHFYUcDAACAECSTECqT1DPmebGkZTHPd5W0j6SpzrlPJY2RNIWOpQEAaefxx6X775d+/nPp2Hitn1vQiBFWwnLssdKVV0oTJ0qPPhpeHAAAAAiF897Xv4Jz+ZIWSTpE0lJJMyT90Hs/L8H6UyX9zHtfWt92S0pKfGlpvasAAJA6770njRkjjR4t/fvf2Te8fFNUVkqHHSZNmyZNny4NHx52RAAAANhJzrmZ3vsGK+k0WEPIe18p6UJJ/5S0QNIz3vt5zrnrnXPH7HyoAAA0s40bpeOPl9q3lyZNIhkUyM+3z6NDB/t8NmwIOyIAAAC0kKTOiL33f5f09zrzfplg3YN2PiwAAFLEe+mMM6RPPpGmTpW6dQs7InPqqTZ98slw4+jaVXrmGemgg6TTT5deeEHKS6ZFOQAAADIZZ3wAgOz2hz9YkuPWW6WvfS3saKLKyqykgwMPlH7/e2nKFPucAAAAkPVICAEAstfrr1vHySecYKOLIbGLLpJOPFG6+mrptdfCjgYAAADNjIQQACA7LV8unXSStNde0kMPSc6FHVF6c0568EFpwABp3Dhp6dKwIwIAAEAzIiEEAMg+FRWWDPriC+m556R27cKOKDPsuqt9Xps3W22hioqwIwIAAEAzISEEAMg+V18t/fe/0sSJ0j77hB1NfGPHWkk3Q4ZYTaG33pKuuCLsaAAAANBMGHcXAJBdXn7ZOki+4ALplFPCjiaxG28MO4LExo2Tpk2T7rhDOvhg6Zhjwo4IAAAAKea896HsuKSkxJeWloaybwBAllq71moEdekivfOOVFQUdkSZa/t2adQoacUKae5cqVOnsCMCAABAEpxzM733JQ2tR5MxAED2uOACSwo9/nj6J4OOP95KuiostM9x3TrpvPOkkG4gAQAAoHmQEAIAZIenn7Zy3XXS8OFhR9OwtWutpLNhw6Rf/1p69llp8uSwowEAAEAKkRACAGS+5cul88+XRo+mI+RUu/xyacwYq321bFnY0QAAACBFSAgBADKb99KECdKWLdJjj0n5jJeQUvn59rlu2yadfTZNxwAAALIECSEAQGZ75BHppZekm26SBg4MO5rsNGCAdPPN0iuvSA89FHY0AAAASAESQgCAzPXZZ9Ill0gHHST95CdhR9M4hxxiJVNccIENQf/Tn0qffBJ2NAAAANhJDDsPAMhM1dXSoYdKM2ZI778v9ekTdkTZ77PPpKFDpZEjpf/7PymP+0oAAADphmHnAQDZ7Z57pNdek26/nWRQS+ndW7rjDun116U//jHsaAAAALATSAgBADLPokXSz38uHXmkdNZZYUfTNEccYSXTnHGGdPTR0lVXSR98EHY0AAAAaCISQgCAzFJZKZ1+utSqlfSnP0nOhR1R02zdaiXTOGefe5s20o9+ZH8PAAAAZBwSQgCAzPL730vTp0v33it17x52NLmpWzfpvvus/6abbgo7GgAAADQBCSEAQOZ47z3pl7+UfvAD6aSTwo4mt514ojRunPTrX0uzZ4cdDQAAABqJhBAAIDNs325NlDp0sNpBmdpULJvcfbfUqZN02mlSeXnY0QAAAKARSAgBADLDDTdIc+ZIEydaEiLTHX20lUzWsaP04IPS3LlWUwgAAAAZw3nvQ9lxSUmJLy0tDWXfAIAMM2eOVFIinXyy9PjjYUeDus44Q3riCemdd6SRI8OOBgAAIKc552Z670saWi+pGkLOucOdcwudcx85566Ms/xc59z7zrnZzrn/OeeGNCVoAAB2UFkpnXmmNRW7/fawo0E8t90mde4snXWWVFERdjQAAABIQoMJIedcRNI9ko6QNETSyXESPk9574d670dIukXSbSmPFACQm37/e2nWLOs3qGPHsKNJnYMOspINdt/d/j6zZ0u33BJ2NAAAAEhCMjWERkn6yHu/2Hu/XdJkScfGruC93xTztK2kcNqhAQCyy8KF0nXXSccfbwXp67jjbPS366+X5s8POxoAAAA0IJmEUA9JS2Kel9XMq8U5d4Fz7mNZDaGL4m3IOTfBOVfqnCtdvXp1U+IFAOSK6mprgtSmjY1mhfR3113SLrvY362qKuxoAAAAUI9kEkLxxvXdoQaQ9/4e7/1ekn4u6dp4G/LeT/Tel3jvSzp37ty4SAEAueWee6Q335TuuEPq1i3saJCMrl2lO++Upk+35BAAAADSVjIJoTJJPWOeF0taVs/6kyV9b2eCAgDkuE8/la66Sjr8cOm008KOBo1xyinSkUdK11wjLV4cdjQAAABIIJmE0AxJ/Z1zfZ1zhZLGSZoSu4Jzrn/M06MkfZi6EAEAOcV76cc/lpyTHnjAptnoxBOtZBvnpPvvlyIR+zt6uhUEAABIR/kNreC9r3TOXSjpn5Iikh723s9zzl0vqdR7P0XShc65QyVVSFov6fTmDBoAkMUeeUR69VUbtapXr7CjaT7nnx92BM2nZ0/p1lulc8+VHnzQEkMAAABIK86HdOeupKTEl5aWhrJvAECaWrZMGjJEGj5ceu01KS+ZiqwZassWm7ZpE24czaW6WjrkEGnWLGnePKm4OOyIAAAAcoJzbqb3vqSh9bL4TBsAkFG8t1oz5eVWqySbk0GS9bNz5JFhR9F88vLs71hRIZ13Hk3HAAAA0kyWn20DADLGM89If/ub9JvfSP37N7w+0t9ee0k33CC99JI0aVLY0QAAACAGCSEAQPjWrJF+8hNp//2lSy4JOxqk0sUXS6NHSxddJK1aFXY0AAAAqEFCCAAQvosvljZskB5+WMpvcLwDZJJIRHroIWnTJksKAQAAIC2QEAIAhOull6SnnpKuuUbaZ5+wo0Fz2Htv6Re/kJ5+WvrrX8OOBgAAAGKUMQBAmNavtyRQhw7SzJlSYWHYEbWcRx+16fjxYUbRcioqrEngypU26liHDmFHBAAAkJUYZQwAkP4uvdQSBI8+mlvJIMkSQbmSDJKkggLpkUesv6iLLw47GgAAgJxHQggAEI6XX7ZE0JVXSvvtF3Y0LW/NGiu5ZN99pauvlp58UpoyJexoAAAAchpNxgAALS+2qVhpqVRUFHZELe+gg2w6dWqYUbS87dut6diqVTQdAwAAaAY0GQMApK/YpmK5mAzKZYWF9nen6RgAAECoSAgBAFpWrjcVA03HAAAA0gAJIQBAy1m/XpowwZqL/eIXYUeDMF1zjTRsmHTOOdK6dWFHAwAAkHNICAEAWg5NxRCIbTp2ySVhRwMAAJBzSAgBAFoGTcVqO+88K7ksaDr2xBPSiy+GHQ0AAEBOYZQxAEDzY1QxJMKoYwAAACnFKGMAgPRBU7EdLVliJdcFTcdWr6bpGAAAQAsiIQQAaF40FYvvtNOsgKZjAAAAISAhBABoPowqhmRde62NOjZhAqOOAQAAtAASQgCA5kNTMSSLpmMAAAAtioQQAKB50FQMjUXTMQAAgBZDQggAkHpr19JUDE0T23RszZqwowEAAMha+WEHAADIMt5LP/6xNf158UWaiiVy2WVhR5CeCgulxx+XRo2SzjxT+tvfJOfCjgoAACDrJFVDyDl3uHNuoXPuI+fclXGWX+qcm++ce8859x/nXO/UhwoAyAh/+pP0wgvSjTdKI0eGHU36+u53rWBHw4dLt9xiCcX77gs7GgAAgKzUYELIOReRdI+kIyQNkXSyc25IndXelVTivR8m6VlJt6Q6UABABliwwDoEPuww6ac/DTua9LZwoRXEd9FF0hFHWE2quXPDjgYAACDrJFNDaJSkj7z3i7332yVNlnRs7Are+9e891tqnk6XVJzaMAEAaW/bNunkk6VddpEee0zKo5u6ep1zjhXE55z0yCNSu3Z2XG3dGnZEAAAAWSWZs/UekpbEPC+rmZfIWZJe2ZmgAAAZ6KqrpDlz7CK+W7ewo0E26NrVkotz50o//3nY0QAAAGSVZBJC8Xpy9HFXdO5USSWSbk2wfIJzrtQ5V7p69erkowQApLdXXpHuuMOa+Rx1VNjRIJscfrg1Q7zrLunll8OOBgAAIGskkxAqk9Qz5nmxpGV1V3LOHSrpGknHeO/L423Iez/Re1/ivS/p3LlzU+IFAKSblSul8eOloUOlm28OOxpko5tuso6mx4+Xli8POxoAAICskExCaIak/s65vs65QknjJE2JXcE5t6+kB2TJoFWpDxMAkJaqq+0ifdMmadIkqVWrsCNCNioqsuNr82Y73qqrw44IAAAg4zWYEPLeV0q6UNI/JS2Q9Iz3fp5z7nrn3DE1q90qaRdJf3HOzXbOTUmwOQBANvnjH6V//EO67TZp773DjiazXHutFSRn8GBrlvivf0m33x52NAAAABnPeR+3O6BmV1JS4ktLS0PZNwAgBWbPlkaPtqHBX3jBRoUCmpP30vHHSy+9JE2fLo0cGXZEAAAAacc5N9N7X9LgeiSEAACNtmWLtN9+1lRszhypU6ewI8o8s2fbdMSIcOPINGvXWn9CbdtKs2bZFEgz3kvl5dLWrdGyZUvt55WVUlVV7VJdveO8vDwpEpHy860kelxYKLVps2Np3dq2AQDIHckmhPJbIhgAQJa59FJp4ULp1VdJBjXVJZfYdOrUUMPIOB07Sk88IR1yiH2Gf/pT2BEhC1VXSxs2SKtWSatXS2vWSBs32rxgGvs4mG7aFE38hHTPNa7WraMJorZtpd12k3bfXWrfPjqN97hzZ6lLF3s9ACD7kBACADTO889LDzwgXXml9K1vhR0NctHBB9vxd+ON0ne+I51wQtgRIQN4b4mbsjJp6VKbrlhhCZ9Vq6LJn1WrLAFUWZl4W7vuakmV9u1tusce0qBB9jiolRMkYeI9bt06WsOnbglqBAXF+2htosrKHR8Hz4MaSVu21C6bN9d+/uWX9jmsXSt9/LElstavb/j9dukide0af9q9u1RcbNOCgtT/7QAAzYMmYwCA5H3+uTVx6tdPevNNzvx3xkEH2ZQaQk1TUSF97WvSokXSu+9KffqEHRFC9sUX0iefSJ9+Ki1ZEk36xCaAtmzZ8XVBsiOoDVP3cZcuVjEtqDnTrp0larKJ9/bZrF8fTRCtXx9NkK1cueN0zZoda0E5J3XrZsmh4mKpZ8/o4+JiqVcvm2bb5wcA6YYmYwCA1NqyRfre96wtxVNPkQxCuAoK7DgsKZGOPdYSlLvsEnZUaEbbtkmffWZJnyDxEzz+5BOr8RIrPz9ac2X4cOmoo6QePex5jx5W9thDatUqlLeTVpyzpmRt29rnk4yqKksKrVwpLV9uSbggAVdWZq2K//Mfa0YXKz9f6t1b2nNPqW/f2tM997SmaoxRAAAtg4QQAKBh3ktnnmkdIb/8stUQAsK2117S00/bSHfjx0vPPEPvuRmustKSPosWWfnww+jjzz+vXSOlsNAqhvXta33c9+0bLT17Ws2eTD4c0r0SYSRiTca6dpWGDUu83qZNVkNryZJoQm/xYps+/7wllWK1a2eJoQEDpIEDrQSP27Vr3vcEALmGhBAAoGE33mgX3jffbBff2Hm/+13YEWSHww6Tbr1Vuuwy6YYbpF/+MuyIkISNG6UFC6R586QPPogmfT7+2FoDBnbbzZIBX/ua1L9/tDZJ375WuyeTEz4NyZafiHbtrAweHH/5pk3RWl6LF1v5+GOptFR69lmrlBro2nXHJNHgwXY80AwNABqPPoQAAPV78UVrkvPDH9roTtTlR7rxXjr9dDs+n39eOu64sCNCjQ0bpPnzrcybF50uXRpdp6jIkj0DBuxYOnXiJyeXlZdbcmjhQiuLFkWnsTWLWre2xNDee9cuvXtnd9IQABJJtg8hEkIAgMTmz5fGjLHbsG+8wdjDqfTWWzY94IBw48gW27ZJ3/ymZRumTZOGDg07opyyfbvV9JkzR3rvPStz50rLlkXXib1oHzIkOu3Th9od8fATUb916yw5FCQZgxKbbGzbNnrMDR1qfUkNH26dhgNANiMhBADYOevWSaNG2RjFpaXJ9zSK5KR7ByGZaNky62S6VSvpnXeseglSbuXKaOInmC5YEG3qVVhoF+D77BOtqREkfqitkTx+IpomqJUWmySaN886vg7ssUc0OTRsmE0HDrQOrwEgGzDKGACg6SorpZNOsl5Ap04lGYTM0L279MIL0je+IZ14ovTPfzIa3k7w3kbymjXLyrvv2nTlyug6PXrYBfWRR0YvrPv352NHeNq3t1pVdWtWrVljCczY8p//RBOZRUWWvBwxQho50srw4VKbNi3/HgCgpZAQAgDs6PLLpVdflR5+WBo7NuxogOSNHi1NnGijjl12mfTHP4YdUUaoqrIRveomfzZssOX5+VbL5/DDpX33teTP0KFUwkLm6NRJOuQQK4GgqWNQ223OHGnKFPvXJ1mNtsGDowmi/fazhNGuu4bzHgAg1UgIAQBqe/RR6Y47pIsvls44I+xogMY7/XS7srv9drvFf9ZZYUeUVry3jnpnzLDWoDNmWPJn82ZbXlRkCZ+TTopeCO+zj7XEA7JJYaEd68OGSaeeavO8l8rKosnRmTPt/sgTT9hy56wW3MiR0v77Wxk50vorAoBMQ0IIABA1bZp0zjnSoYdKv/992NEATXfLLdar8XnnSYMGSQceGHZEoQgubmOTP6Wl0Zo/rVpZjZ8zz7TaDyNH2sdFky/kKueknj2tHHtsdP7y5VZzbuZMSxS9+aY0ebIty8uzGnRBgmj//S3JVFgYznsAgGTRqTQAwJSV2Vls27bWIW+HDmFHlN1mz7bpiBHhxpHN1q+3jtE3bbIsSM+eYUfU7DZssKTP229bmTEj2udPfr4184q9aB0yhORPuuInIv2tXGnfsdiyZo0tKyy0Cor7728tWUePtppFdKwOoCUwyhgAIHlbt0pf/7qN4Tt9uvWsCWSD+fOlMWPsSuy//82qHmIrKqzvkyD588471h9KYNAgy4cFyZ/hw2n2BTQn76XPPosmh0pLrXzxhS1v396+k0GCaPRo+uEC0DxICAEAklNRYSMy/e1v0l//Kh1zTNgR5YZXX7XpoYeGG0cuePFFa/tx5JHS889nbDuOJUssXzttmiWAZs2Stm2zZV26RC8wgyRQ+/bhxoudw09EdqiqskRtkLh9+23p/fel6mpbvuee0e/umDHWhDNDf6IApBESQgCAhlVWSqecIj3zjI3G9JOfhB1R7jjoIJtOnRpmFLnj/vutP6Hvf196+mlrP5XGtm2zvkqCBND06dLSpbasVSvr7yc2AdS7t/V9guzBT0T22rw5+v0OkkTB97uoyL7fY8bYIJ9jxkjFxeHGCyDzJJsQSu+zIQBA86mqslHEnnlGuvVWkkHIbueea1mWn/5U+tGPbMigSCTsqCRFm5kEiZ9p06z/mIoKW963r/TNb0YvEOmsFshsbdtK3/iGlUBZmSWGpk2zcs890m232bLi4mhyaOxYq0VE808AqUBCCAByUXW1XSA/+aR0ww3Sz34WdkRA87vkEksKXXWV3YZ/6KFQenjdts2ae731VvTib/lyW9amjTX3uuwyu/gbM0bq2rXFQwTQwoqLrRx/vD3fvt0Sw0GSeNo06S9/sWWFhTYi4AEHWILogAOk7t3Dix1A5iIhBAC5xnvpooukBx+Urr1WuuaasCMCWs6VV1pG5te/tlvs997b7G2tli6NJn/eesuSQbG1f771LbuoC2r/pHlrNgAtoLDQmoOOGmX/siVLHMcmiGJrEfXqVTtBNHw4IwgCaBinHACQS7y32kD33GPT668POyKg5f3qV5YUuvlmqyl0++0pSwpVVEhz5ljiJ0gCff65LWvVSiopsYpKBxxgtX+6dUvJbgHkgD32kI47zooUrUUU/Nb873/S5Mm2rHVrq214wAHRRBEjmgGoi06lASCXXHut9NvfShdeaJ1I0wtteBYutOnAgeHGkau8t/6E7rzTag397ndN+j6sWWN37IME0DvvSFu32rLi4tp37EeMoO8fJI+fCDTFkiXRGkRvvSW9+260RuKAAbUTREOGhNJqFkALSOkoY865wyXdKSki6UHv/U11ln9D0h2Shkka571/tqFtkhACgBZ2ww3SL34h/fjHNuISZ4HIdd7byGMPPGBNyH75y3pXr66WFiyoXfsnuGjPz7eOXmMvtnr2bIH3AAD12LrVRjQLfrfeektavdqW7bab1VQMfrdGjZLatQs3XgCpkbJRxpxzEUn3SPq2pDJJM5xzU7z382NW+1zSeEn0SgoA6ejWWy0ZdNppJIPSxYsv2vS73w03jlzmnPUhVF5uzciKiqSf//yrxZs2WY2f4CJq+nRp40Zb1qmTJX3Gj7cLqZIS6xAaSBV+IpAKrVtLX/uaFcny4B+phxu6AAAdzElEQVR/XDuxfd11Nt85aejQ2jUb99qLysRANmuwhpBzbqyk67z336l5fpUkee9vjLPuo5JeooYQAKSRu+6yHilPOslGFaPH2vRw0EE2nTo1zCggSVVV8qecqo+eLtW0H96tae2+o7fekt5/P3qRtM8+0QukAw6Q+vXjIgnNi58ItJSNG3dMfm/aZMs6d67921dSYkkmAOktZTWEJPWQtCTmeZmk0U0MaoKkCZLUq1evpmwCANAY995ryaDvfU964gmSQUCNzZulGTOCvjYimjbtKa2Rk56S2rUq15hvFOm44+wCaPRoa1oBANlot92kb3/biiRVVe3YPHbKFFuWn2/9oQUjI44dK/XuTYIcyFTJXBnE+3o3qSdq7/1ESRMlqyHUlG0AAJJQVSVdfrmNnnT00TbsCOPPIkd5Ly1eHO1oddo06b337GsiWae9Rx/tNHb/So199jINee1uRQZdKF37B5KoAHJOJGK1IvfZR5owweatXl17yPuHHrIKyJKNlhibINpvP2oRAZkimbOcMkmx3SIWS1rWPOEAAHbaxo3SySdLr7wi/eQn0m23cVGLnPLFF1b7Z/r0aAk6Ud1lF6vxc9VVduEyerTUsWPwynzpnNukyyOWTF240JKp7duH9VYAIC107mz9WQV9WlVWWrPa2ET7Cy/YsoICq0U0Zky09O1LLSIgHSVzhTBDUn/nXF9JSyWNk/TDZo0KANA0H39sZ2uLFkn33Sede27YEQHNqrpa+uCD2smfuXOtVpAkDRokHXWUJX7GjrU73pFIPRuMRCyJuvfe9v0ZO9Z69+3Xr0XeDwBkgmBkxX33lc4/3+atWhWtRTR9uvTww9FaRJ061U4Q7b8/I5oB6SDZYeePlA0rH5H0sPf+t8656yWVeu+nOOf2l/SCpN0lbZO0wnu/d33bpFNpAEix11+Xjj/erpCffVb61rfCjgj1WVLTPR9jkzfKihXW+enbb9t0xozoyF/t29e+4Bg1Stp9953Y2euvS9//vj1+7rloL79AC+AnApmuslKaN8+SQ2+/bdMFC2yZc9KQIfZbPXq0/V7vvTcVmoFUSbZT6aQSQs2BhBAApNBDD1lthr32stoM/fuHHRGw0zZvlmbNiiZ/3n5b+vxzWxaJSMOG2UVEkAAaMEDKy0txEEGtuw8/tE7af/zjFO8AAHLHhg32ex6bJFq3zpa1bm39DwUJolGj6LAaaCoSQgCQC2I7jz7sMOnpp+nvJFM8/bRNTzop3DjSREWFNfWaMSNa5s6Ndvzcp0/0ImH0aGum0KZNCwW3caM0bpz0j39Il1wi3Xort7HR7PiJQC4IOv2Prfk5a5ZUXm7Lu3SJJof239+Gve/UKdyYgUxAQggAst2mTdZ59N//TufRmShofjR1aphRhKKqyvprDhI/paXS7NnRC4Ddd7eT/tGjo0mgLl3CjVmVlZZ8veMO6YgjpEmTGIsezSqHfyKQ47Zvtw6r33knmij64INo33B9+tj/iJISSxKNHMm9MKCuZBNCXDkAQCZavDjaefT990vnnBN2REBcVVXW2mrmTLvrW1pq0y+/tOW77GJNBC68MHr3d88907CJQH6+1cQbPFi64IJoZ9N77RV2ZACQVQoL7f/CfvtJ551n8zZtsv8dwU2EGTOsu8RA//7R/yH77We1SHfdNZz4gUxCQggAMon30hNPWLMVSfrXv6SDDw43JqBGZaXdxZ01K5oAevdd6wtIklq1koYPl8aPtxP3/fe3fn/qHfUr3UyYYFceJ5xgVx5/+IN0xhlpmMECgOzRrp3Vmovt23/tWvtfEySJ3nhDeuqp6PIBA6z2UGzZqYEGgCxEQggAMsWnn1pNoH/9SzrwQOnRRxkKG6HZts36+Jk925I+s2ZJc+ZIW7fa8jZt7A7tmWfa3dqRI61yTVa0ajz4YGvHcOaZ0llnWfOxBx6wqk0AgBbRsaN1n3jYYdF5K1bUvinx1lvS5MnR5X37RpND++4rjRghdetGTh+5KxtOywAgu1VVSXffLV1zjZ2x3H231aFO+XBKQHxr1liy5913LQE0e7bVBAo6fN51VzuxPvdcO8neb78MrPnTWHvtJb32mjRxonTFFdLQodINN0gXXZTlbxwA0le3btKRR1oJrFlj/7+CJNGsWdJzz0WXd+liiaERI6wW64gR9j8sK25gAA2gU2kASGfz51sNhOnTrSPb+++XevUKOyqkwpo1Nk2j4VIqK61bqvfesw4933vPkj9lZdF1ioujJ84jRlgiqE+fHM9PLlliSdqXX7ZesB98UNpnn7CjQoZLw58IIGts2GA3OoKbHHPmWK3Xigpb3qqV5fmDBNHQoVZocoZMwShjAJDJtm+XbrxR+u1vreH8nXdKP/whdZqREt5Ly5dHkz7BdMECO/QkuzM6cGC0Sn1w55SL0wS8t3YJF11kw9RffbWVwsKwIwMAJGH7dqv9GiSIgmTRunXRdXr0iCaHhg2z6aBBUlFReHED8ZAQAoBM9fbb0tln262qk0+2ZFDnzmFHhVR79FGbjh/frLtZtUqaN88qm82bFy1r10bXCU5wg5PbYcPsBJdcRhOsXm2dvj/1lLT33tJDD1mtIaCRWugnAkA9vJeWLrUbJ7El9gZKJGI3UIYOtZ/9vfeWhgyxlsUFBeHGj9xFQggAMs26ddJvfmMJoB49rHnYUUeFHRWaSzBUytSpO70p7y3xM39+7cTP/PnRZieSVTYLTlaHDYsmgDp02OkQUNfLL1unSkuXWq2hX/zCekAFkpTCnwgAKVZRIX344Y6Jok8+ia5TUGB9EQUJoiFD7HG/ftxwQfMjIQQAmWLlSun226V77pG+/NL6IrnpJrt6R/ZqwtVeRYW0eLFVaa9bNmyIrheb+AlOQPfeW+renVaHLWrTJunKK6X77pPatpXOP1+69FLr9RRoAAkhIPNs3mz/k4ObM8GNmsWL7eaNZE2y+/WzWkWDBtWecoMGqUJCCADSXVmZdOutNkpRebl04onW58iwYWFHhpaQ4GrPexs298MPo2XhQjvB/Ogj6/g5sMcetU8mBw8m8ZOW5s61PsEmT7bbwmefLV1+OR3Eo14khIDssWWL/S8PavDG/l8POrKWrIeA4H/6oEFWw6h/f6lvX/opQuOQEAKAdPXxx1YD6LHH7Or/tNOsFsGAAWFHhhbivVRx4EHaulV6/qKpXyV+PvrIypdfRtctKLCTweDkMPZO4m67hfce0AQffmjf/ccft4zdj35k3/1+/cKODGmIhBCQ/SorrZnZwoXRJFEwXb06ul5ent1D6NfPzgn6948+JlmEeEgIAUC6mTfPaglMmmRX+WedJV1xhdS7d9iRoRls2yZ9+qlVE49XXtp8kCTpYE1VJCLtuWf0JC+29OplHVYii3z2mdUOfPBBuzV88slWO3DIkLAjQxohIQTktvXrpUWLojeMYmsOxzYTz8uTeva084ig9O0bfdypE7WGcxEJIQBIB1VV0uuvW/9Azz9v/Yice6502WXW3gcZq7xcWrLEkj6ffRYtQRJo6dJofwGS1Lp17ZO1/j22aK+9pL2GtlGfPoxEkpOWL5duu836GNq8Wfr+960PsYMPJgsIbdli0zZtwo0DQPpZu7Z2ouijj6ym0eLF1jVlrLZtayeK+vSxm029e1vp0IGEUTYiIQQAYfHeho6fNEl65hnrEGa33WykoYsvZqShDOC9nWwtWWJdPS1ZIn3+ee2kz4oVtRM+eXnWd0+fPjbUbGzyZ889pa5dOeFCAmvW2OiCd90lbdxoB8uJJ1rNoTFjOHAAAEnbvLl2DeUgURRMg2RzoG3baHKod+9osqhXL6m42M5tGBUt85AQAoCW5L2NNzppknUc++mn1qD7yCOlceOko4/mNm+aqKy0IdqXLbNSVhZN+sROy8trv66gwKpk9+5tSZ/gxCl4XFzcyFo+995r0/PPT9E7Q8bbutWGq588WXrpJTsIe/e235CTT7YO50kO5Qx+IgCkWnDDK7Zmc92ybt2Or+va1c6Biot3LD16WKX3tm1b/v0gMRJCANASPvrIkkCTJkkLFlgzj0MPtYu3732PXn9bUHm5VZMOSpDwWb48+njZMksGVVfXfm1+vp3QFBfXPuGJfdy1a4pb8dBBCOqzaZP017/ab8u//23NTwcPtt+WceOsgylkNX4iAIThyy8tMRTcMItXYvswCrRrZ4mh7t1rT2Mfd+1q63Fvo/mREAKA5lBWJv33v1beeMM6ipakr3/dLtSOP17q0iXcGLOE99Z6ZvVqa1GzerUlc2KTPitWRB/HOzlxzoZw7d699glK7PPiYvuTtXiXLVztIVmrV0vPPWfJoTfesHmDB0vf+IaVr3/dspfIKvxEAEhXX34ZTQ4tXWo334IbcLHTbdt2fG1RkZ13de26Ywnmd+5spWNH+lhsKhJCALCzvLee+t54I5oE+uQTW7bLLtKBB0rf/rb19cHFWL2qqy25s3atVUUOpsHjNWuiSZ+grFljzbvi2W03O2Ho1q32iUTs8+7dbZq2JxJc7aEpysqsb7J//1t6803piy9sfp8+lhgKEkQDBnALNsPxEwEgkwU39oLk0PLldgOv7s29YF6ic7727W2ktE6dLEkUO+3Y0TrF7tCh9uOiopZ9r+mIhBAANEZ1tXUes2CBNH++9NZblgBatcqWd+5sF1lBGT7c2hnlkPJya8WycaPVxtmwwYZEjZ3WnRckfdav37GZVqz27aN3g4IS/MOvW7p0kVq1arn33Wy42sPOqqqS5syJ1lj8738tmyrZF+XrX5fGjrXh7AcPth5C8/LCjRlJ4ycCQK6orrZzxyBBFNwkTDRdvVravj3x9tq0qZ0k2n13O9eMLfHmtW9vfSFlw/2UZBNCuXU1AwDbt1u/PwsWSB98YNPgceywC717S9/5TjQBNHBgxv138N76qP3ySytffBGdxj6OnRckfDZt2vFx3U6W64pEdvzn2ru3/TMO/iHHm7Zvn3O5NSA1IhFp5EgrF19sX/pFi6LJoTfesKZmgTZt7Lds8ODapV8/hpABAIQmLy9au2fw4IbX997OX2NrmweP65a1a+00P7ihuXlz/dtessS6E8gVSdUQcs4dLulOSRFJD3rvb6qzvEjS45L2k7RW0kne+0/r2yY1hACkXEWFdSqzdGm0B+Hg8dKl9gu/eHHtOqk9e+54cTR4sFVFaaKVf16pxdcsVvnn5SrqVaQ9f7unup7S9avlVVWWqNm2rXYJ5m3ZYo+3bo0+jjfdvLn+smVL7WHR65OfL+26qzXFatcu8TR4HCR8qt5co7U3f6w2G7arlapU0DFf/e/sX+v9trSvPv/Pyu2/VpWiUyep5jPJ38lYG/o7110n0iEiJ6fKdZUJ1w9es+jiRapaW1UrTkkN7q++uJKJt7k+h4bWlxT3Pdf3eRb1KlLHIztq7d/Xpvw97ex7TCtr1kQT37Hl88+/WqVaEZXn91DewN4qGtnHeljv3j067d5dK1/L0+JfLWnwWE70922uzy+j/zZJCuO729A+w/ouNlVTfouRmxadv0jLJi6zc4YaRb1b5jjJhd+zdFBRUbu2e92a7hdeKLVuHXaUOy9lTcaccxFJiyR9W1KZpBmSTvbez49Z53xJw7z35zrnxkk6znt/Un3bJSEEYAfe26/0li1WvvgiflukmMd+/QZVrtmgymWrVLlqnSqUr8qaUqECVUZaqbJLd1V23kMVnburoueeVor7qmKPXqrIb62KCqs4VFFRu2zf3nApL69dvliyXZsWbdP2aqcK5VlxearapUDbfZ62bUvcRjoZRUX2T6p1a6vSmkzZZRdL9tSdxj4uLGx8BaiVf16pBWcskCpqz3eFToMeHhTKSczKP6/UwgkLVb2lnvZpMZoaa7z95LXJ08CJA2tdUNUXS931g9fE+0wVkVzEyW/3Db4+XlzdTu+mFY+tqDfepkjmc2hofRVIqlatk29px79NMn/bVLynZGJujv20tFUPLVbZhf9R622fqo0+UxstUWHeWu2y+wZFNq7c4YfKy6lC7VWuTqpQO1Vql69KVUE7dTxlL1W3ba8lD27S9vK2Nl9tVZVfJO9aqboiIsvIpu7zy9a/Tazmeo/1bVdSvfsM67vYVE35LUZuWnT+Ii27b1ncZc19nOTC7xlaVioTQmMlXee9/07N86skyXt/Y8w6/6xZZ5pzLl/SCkmdfT0bz5aE0PSH5qm6KrkLj2R5n9xVWXN0/xRvm3HnqekxJr+PRNt0MY8TbyN4HhvrV/PivM5/9dzFXTfetO4+vnpe7eXlEq4bb70dSk0sted7mxe8rtrLy9rd2rpevlrRedXBsh2n9tiruqrmcZVXtY++prra18yPrrfjtKZUeVVV27pVVd62VWU1YaqrvaqqnE2r3Vfzq6qkqupg6qJFkVolSO7Ee1yliKrVckNDFRRY4iSYBqWoKFrKZ29SfnmlClWtAlWrQF4FqlarXZz6/LibWre2/m9iS915rVtby46601atQhgJqx7T+kyzWjhxFPUu0thPx7ZwRPXHlEhTYk20n9htJRNL3X03Nv6kXx/UkGrg9Y2VzOeQzPqJNPbzrG/fTdXY95gp6n1fi0dbzaKa2pUf/+i/iqxboUKtUZHWKl9f1JQvVaAvlLdDBnNHXnmqUitVq0hVaiVf0EptRnSKZriDH9S6P7CxzwsKrDpjJCJFIvr0xiWqWF8tr4ikiLzy5JWn/I5F6nfHQGuD4JxNgxLvuWTTRCV2efA4mWmsZOfVMff4udq+csfPt6BroYY+v0+TtilJ73//fVWsiLPdbtYLf6JlQ58fmvC1idYPWzLxpkusCNesr82yGxQJNOdxUt93kmOzhe23X1Y0o05lH0I9JC2JeV4maXSidbz3lc65jZI6SlqTXLiZ69tn99KX2jXsMIAmyVOVnLzyVK08Vcd97OQVUdVX8+srEVUpz3lFXLXynJSXF30cyZPy8qqVl2ePC1tJkYhTJN8pki9FCmSPC7wiBdWKFDhFCp0K2nhFWkv5bSKKtC5Uftsi5bfKVyRi1wXBtUFwnZCfX/tx3ecFBYlLcL1R93okmJ/MefbUvFnxF2yWDrqtW2r/gCEr/zzxxXl9y5pTU/abytfEzk9mu3XXaWwsSb8+TjKoKftL9vWNnZ/M9pN9baqPvVS9l3RT7/vKy7NOqbt0kUaM0JL1bevdVp62K19ffpUkChJGEW1VRNuVp22KqLz2tKJcbTq1tdqg69fvWC0z3vOq2gdyn0QBrZV0WqM/krQUJ+VjVko6sOnbTXh5uaKeF62QdEA9r02wftiSijdNYkW4Rja0QjMeJ/V+Jzk2W9ayZdIee4QdRYtJJiEU7xKobuWNZNaRc26CpAmS1KtXryR2nf7+evMiVVWmvqpOsk03XMJ6NKnddxNvbNl6cWJMfh/x35+rtU7ibQT7jrdOrXl11vtqWnd+ntthfsLX5DlbnmDdeOslvEEpH10vKDXPgxucsfO/ugkacdFpMD/PffVaRSLx75jWnVdzN9ayMAV1nkfv1jJ6jVTUqyj+Xfde2Tf+ZaL3GiwLQ30x1feaVO0ndlvJxFJ3342NP+nXJ6ohtJN/p8Ye7zvz/pJ9baqPvWz9TjfmfTX02VerUK73HqrSHtqS5N+3qHeROv+9kTWsqmuqlVZWSlVVemfING1fslVO1XKqkqtpe1jUPV8jXx8RrcYaVI+t+zxIMMWtphtbtbfO42SmsZKdF8f8U+epYlWcWgNdCjTkiSFN2qYkzf/R/ITblZR4n48PSfjaROuHLZl40yVWhGvOkXPqryHUjMdJfd9Jjs0W1qFD2BG0qGQSQmWSesY8L5ZUt3FlsE5ZTZOx3SStq7sh7/1ESRMlazLWlIDTzSFX7Bd2CADSyJ6/3TNuG/CgY9Vssudv90zYh1BY7zfe51+fpsaazN+5oVjiHReJPtNEfQjFe31j+hDa2b9TY4/3uJ9JPX0INebzbGjfTZWt3+nGvK/GHMvx/r7ONXzsJiW4SVFgSYveNw6L/x5uGSj1y44+NzreFv89DrxtoHRY099jx9tGJN6uEvQhdNtA6Ttd4762rtj1w9ZQvOkUK8LV+py+9fch1IzHSb3fSY5NNKNkbufPkNTfOdfXOVcoaZykKXXWmSLp9JrHJ0j6v/r6DwKAbNX1lK4aOHGginoXSc7ugmdrh4BdT+mqwY8MVqRjtGOj/I75oXUoHcT01ecv6asupoJpTM3AnYk1mb9z3XUiHSPK75hf73GR6DMd/NhgDXp4UIPHVaK4Btw7oFmOy8Ye7/HWH/zIYA1+rOHjKN5ru5/Xvdm/a9n6nW7M+0r2WE70903m2G3u95Cpmus91rfdhvYZ1nexqZryW4zcNODeAep+XnfV7Z6yJY6TXPg9Q3pKdtj5IyXdIft6POy9/61z7npJpd77Kc65VpKekLSvrGbQOO/94vq2mS2dSgMAAAAAAKSLVHYqLe/93yX9vc68X8Y83ibpB40NEgAAAAAAAC2PHmABAAAAAAByTFJNxpplx86tlvRZKDtHc+okaU3YQSBtcXwgEY4NJMKxgfpwfCARjg0kwrGBRLLp2Ojtve/c0EqhJYSQnZxzpcm0VURu4vhAIhwbSIRjA/Xh+EAiHBtIhGMDieTisUGTMQAAAAAAgBxDQggAAAAAACDHkBBCqk0MOwCkNY4PJMKxgUQ4NlAfjg8kwrGBRDg2kEjOHRv0IQQAAAAAAJBjqCEEAAAAAACQY0gIIeWcc7c65z5wzr3nnHvBOdc+7JgQLufc4c65hc65j5xzV4YdD9KHc66nc+4159wC59w859zFYceE9OKcizjn3nXOvRR2LEgfzrn2zrlna843FjjnxoYdE9KDc+6nNf9P5jrnJjnnWoUdE8LjnHvYObfKOTc3Zl4H59y/nXMf1kx3DzNGhCPBsZFz17EkhNAc/i1pH+/9MEmLJF0VcjwIkXMuIukeSUdIGiLpZOfckHCjQhqplHSZ936wpDGSLuD4QB0XS1oQdhBIO3dK+of3fpCk4eIYgSTnXA9JF0kq8d7vIykiaVy4USFkj0o6vM68KyX9x3vfX9J/ap4j9zyqHY+NnLuOJSGElPPe/8t7X1nzdLqk4jDjQehGSfrIe7/Ye79d0mRJx4YcE9KE9365935WzeMvZBd1PcKNCunCOVcs6ShJD4YdC9KHc66dpG9IekiSvPfbvfcbwo0KaSRfUmvnXL6kNpKWhRwPQuS9f0PSujqzj5X0WM3jxyR9r0WDQlqId2zk4nUsCSE0tzMlvRJ2EAhVD0lLYp6XiQt+xOGc6yNpX0lvhxsJ0sgdkq6QVB12IEgre0paLemRmuaEDzrn2oYdFMLnvV8q6feSPpe0XNJG7/2/wo0Kaair9365ZDemJHUJOR6kp5y4jiUhhCZxzr1a0za7bjk2Zp1rZM1B/hxepEgDLs48hjdELc65XSQ9J+kS7/2msONB+JxzR0ta5b2fGXYsSDv5kkZKus97v6+kzaLJByTV9AVzrKS+krpLauucOzXcqABkmly6js0POwBkJu/9ofUtd86dLuloSYd477n4z21lknrGPC8W1bcRwzlXIEsG/dl7/3zY8SBtHCjpGOfckZJaSWrnnHvSe8/FHcoklXnvg9qEz4qEEMyhkj7x3q+WJOfc85IOkPRkqFEh3ax0zu3hvV/unNtD0qqwA0L6yLXrWGoIIeWcc4dL+rmkY7z3W8KOB6GbIam/c66vc65Q1rnjlJBjQppwzjlZPyALvPe3hR0P0of3/irvfbH3vo/sd+P/SAZBkrz3KyQtcc4NrJl1iKT5IYaE9PG5pDHOuTY1/18OER2OY0dTJJ1e8/h0SX8LMRakkVy8jnU5kPRCC3POfSSpSNLamlnTvffnhhgSQlZzh/8O2WgfD3vvfxtySEgTzrmvSfqvpPcV7Sfmau/938OLCunGOXeQpJ95748OOxakB+fcCFln44WSFks6w3u/PtyokA6cc7+WdJKsuce7ks723peHGxXC4pybJOkgSZ0krZT0K0l/lfSMpF6yJOIPvPd1O55GlktwbFylHLuOJSEEAAAAAACQY2gyBgAAAAAAkGNICAEAAAAAAOQYEkIAAAAAAAA5hoQQAAAAAABAjiEhBAAAAAAAkGNICAEAAAAAAOQYEkIAAAAAAAA5hoQQAAAAAABAjvl/NZQuXim8968AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# estimates for the mean\n",
    "red_mean_guess = 1.1\n",
    "blue_mean_guess = 9\n",
    "\n",
    "# estimates for the standard deviation\n",
    "red_std_guess = 2\n",
    "blue_std_guess = 1.7\n",
    "\n",
    "for i in range(20):\n",
    "    likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)\n",
    "    likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)\n",
    "    \n",
    "    likelihood_total = likelihood_of_red + likelihood_of_blue\n",
    "\n",
    "    red_weight = likelihood_of_red / likelihood_total\n",
    "    blue_weight = likelihood_of_blue / likelihood_total\n",
    "    \n",
    "    blue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)\n",
    "    red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)\n",
    "\n",
    "    red_mean_guess = estimate_mean(both_colours, red_weight)\n",
    "    blue_mean_guess = estimate_mean(both_colours, blue_weight)\n",
    "    \n",
    "y = [0.001 for i in range(len(both_colours))]\n",
    "plt.plot(both_colours,y,'mo')\n",
    "plt.ylim(-0.02,0.50)\n",
    "\n",
    "x = np.linspace(-3, 12, 100)\n",
    "ymaxr =max(mlab.normpdf(x, red_mean_guess, red_std_guess))/0.5\n",
    "ymaxb =max(mlab.normpdf(x, blue_mean_guess, blue_std_guess))/0.5\n",
    "\n",
    "plt.plot(x, mlab.normpdf(x, red_mean_guess, red_std_guess),'r',x, mlab.normpdf(x, blue_mean_guess, blue_std_guess),'b',alpha=1)\n",
    "plt.axvline(red_mean_guess, color='r', linestyle='--', ymin=0.06,ymax=ymaxr,alpha=1)\n",
    "plt.axvline(blue_mean_guess, color='b', linestyle='--',ymin=0.06, ymax=ymaxb,alpha=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ! pip install prettytable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+--------------------+--------------------+----------------------+\n",
      "|           |      EM Guess      |       Actual       |        Delta         |\n",
      "+-----------+--------------------+--------------------+----------------------+\n",
      "|  Red mean | 2.9096206771965387 | 2.8132116984626867 | 0.09640897873385201  |\n",
      "|  Red std  | 0.8542016076490087 | 0.8001119672360739 | 0.054089640412934736 |\n",
      "| Blue mean | 6.838227247495596  | 6.972553898467638  | -0.13432665097204133 |\n",
      "|  Blue std | 2.2271904083245206 | 2.0388437926931493 | 0.18834661563137134  |\n",
      "+-----------+--------------------+--------------------+----------------------+\n"
     ]
    }
   ],
   "source": [
    "from prettytable import PrettyTable\n",
    "\n",
    "t = PrettyTable(['', 'EM Guess', 'Actual', 'Delta'])\n",
    "t.add_row(['Red mean', red_mean_guess, np.mean(red), red_mean_guess-np.mean(red)])\n",
    "t.add_row(['Red std', red_std_guess, np.std(red), red_std_guess-np.std(red)])\n",
    "t.add_row(['Blue mean', blue_mean_guess, np.mean(blue), blue_mean_guess-np.mean(blue)])\n",
    "t.add_row(['Blue std', blue_std_guess, np.std(blue), blue_std_guess-np.std(blue)])\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expecation maximization algorithm should always we there in data scientist toolbox. This tutorial has tried to explain the intuitvely how EM works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Intuitive Explanation of EM](https://stackoverflow.com/questions/11808074/what-is-an-intuitive-explanation-of-the-expectation-maximization-technique)\n",
    "\n",
    "[Expectation Maximization Primer](https://www.cmi.ac.in/~madhavan/courses/dmml2018/literature/EM_algorithm_2coin_example.pdf)\n",
    "\n",
    "[How EM works](https://math.stackexchange.com/questions/25111/how-does-expectation-maximization-work)\n"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}