{
 "metadata": {
  "name": "",
  "signature": "sha256:d043207a7d27c4f6f2236ec7066634b384138daeba1531a9665cfa1069f37475"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Expectation Maximization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Expectation Maximization (EM) is a powerful technique for creating maximum likelihood estimators when the variables are difficult to separate. Here, we set up a Gaussian mixture experiment with two Gaussians and derive the corresponding estimators of their means using EM."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Experiment: Measuring from Unseen Groups "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's investigate the following experiment:\n",
      "\n",
      "Suppose we have a population with two distinct groups of individuals with different heights. If we randomly pick an individual from the population, assume we don't know which group the individual is from. So, we measure that individual's height and choose another individual. The goal is to estimate the mean heights of the two distinct groups when we have an unlabeled distribution of heights sampled from both groups.\n",
      "\n",
      "Group **a** is normally distributed as\n",
      "\n",
      "$$ \\mathcal{N}_a(x) =\\mathcal{N}(x; \\mu_a,\\sigma) $$\n",
      "\n",
      "and likewise for group **b**\n",
      "\n",
      "$$ \\mathcal{N}_b(x) =\\mathcal{N}(x; \\mu_b,\\sigma) $$\n",
      "\n",
      "Note that we fix the standard deviation $\\sigma$ to be the same for both groups, but the means ($\\mu_a,\\mu_b$) are different. The problem is to estimate the means given that you can't directly know which group you are picking from.\n",
      "\n",
      "Then we can write the joint density for this experiment as the following:\n",
      "\n",
      "$$ f_{\\mu_a,\\mu_b}(x,z)=  \\frac{1}{2} \\mathcal{N}_a(x) ^z \\mathcal{N}_b(x) ^{1-z} $$\n",
      "\n",
      "where $z=1$ if we pick from group **a** and $z=0$ for group **b**. Note that the $1/2$ comes from the 50/50 chance of picking either group.  Unfortunately, since we do not measure the $z$ variable, we have to integrate it out of our density function to account for this handicap. Thus,\n",
      "\n",
      "$$ f_{\\mu_a,\\mu_b}(x)=  \\frac{1}{2}  \\mathcal{N}_a(x)+\\frac{1}{2}  \\mathcal{N}_b(x)$$\n",
      "\n",
      "Now, since $n$ trials are independent, we can write out the likelihood:\n",
      "\n",
      "$$ \\mathcal{L}(\\mu_a,\\mu_b|\\mathbf{x})= \\prod_{i=1}^n f_{\\mu_a,\\mu_b}(x_i)$$\n",
      "\n",
      "This is basically notation. We have just substituted everything into $ f_{\\mu_a,\\mu_b}(x)$ under the independent-trials assumption. Recal that the independent trials assumptions means that the joint probability is just the product of the individual probabilities. The idea of *maximum likelihood* is to maximize this as the function of $\\mu_a$ and $\\mu_b$ after plugging in all of the $x_i$ data.  The problem is we don't know which group we are measuring at each trial so this is trickier than just estimating the parameters for each group separately.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Simulating the Experiment"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We need the following code to setup the experiment of randomly a group and then picking an individual from that group."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "from numpy import array, linspace, random\n",
      "from scipy.stats import bernoulli, norm\n",
      "from matplotlib import cm\n",
      "from matplotlib.pylab import figure, subplots\n",
      "#random.seed(101) # set random seed for reproducibility\n",
      "mua_true=4 # we are trying to estimate this from the data\n",
      "mub_true=7 # we are trying to estimate this from the data\n",
      "fa=norm(mua_true,1) # distribution for group A\n",
      "fb=norm(mub_true,1) # distribution for group B\n",
      "fz=bernoulli(0.25) # each group equally likely \n",
      "\n",
      "def sample(n=10):\n",
      "    'simulate picking from each group n times'\n",
      "    tmp=fz.rvs(n) # choose n of the coins, A or B\n",
      "    return tmp*(fb.rvs(n))+(1-tmp)*fa.rvs(n) # flip it n times\n",
      "\n",
      "xs = sample(1000) # generate some samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here's a quick look at the density functions of each group and a histogram of the samples"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f,ax = subplots()\n",
      "x = linspace(mua_true-2,mub_true+2,100)\n",
      "ax.plot(x,fa.pdf(x),label='group A')\n",
      "ax.plot(x,fb.pdf(x),label='group B')\n",
      "ax.hist(xs,bins=50,normed=1,label='Samples');\n",
      "ax.legend(loc=0);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VFXe+PHPNyH03kkjdAEJoSShSqiCoujyiKKwdnEt\nqz7Psz93V1dgdXd17a67rro+diwoYkWQEkJLAtIhlAiBFAi9E0g5vz/uJExCkplJ7syduXPer1de\nZO7cc+93QvKdM9977jmilELTNE2ztxCrA9A0TdO8Tyd7TdO0IKCTvaZpWhDQyV7TNC0I6GSvaZoW\nBHSy1zRNCwIuk72IjBeRHSKyW0Qer2a/eBEpEpHJnrbVNE3TvKvaZC8iocDrwHigFzBVRHpWsd9z\nwI+ettU0TdO8z1XPPgHIVEplKaUKgU+BSZXs9zDwBXC4Bm01TdM0L3OV7COAbKfHOY5tZUQkAiOJ\nv+HYVHpLrsu2mqZpmm+4SvbuzKXwCvB7Zcy7II4vd9tqmqZpPlDHxfO5QJTT4yiMHrqzAcCnIgLQ\nGpggIoVutkVE9JuCpmlaDSilxPVel3au8gvjzeAXIAaoC2wEelaz/7vArzxpa4RgXzNnzrQ6BNOc\nOKFUYqJS06crtXevsW3mzJlqwQKlIiKU+tOfLA3PK7z5//fTLz+pNn9vo55Z/owqLC5USil1+sJp\n9fKal1Xb59uqtblrvXZupez1u1kZu78+R+6sNoc7f1Xbs1dKFYnIQ8BCIBR4RymVISIzHM+/6Wlb\nt9+FNL9y+jRcfTXEx8Nrr4E49SfGj4dNm2DwYOjUCe6807o4A8Wuo7u4bd5tzL1pLiNiRpRtb1y3\nMY8OepTOLTpz7ZxrWXDbAvp36G9hpJpduCrjoJRaACyosK3SJK+UurPC48vaaoHp97+Hbt0uT/Sl\nWrWCr7+GESOgZ08YNMj3MQaKUxdOccOnN/D0yKfLJXpn1/e4nnOF57jli1vYdP8mGoQ18HGUmt3o\nO2i9LCkpyeoQam3FCpg/H/7xj8sTvfPr69kT/v1v+PWv4eJF38boLd74/3ty6ZMMjhzMfQPuq3a/\nW668hb7t+/Ln5X82PQawx+9mdez++jwlyuLFS0REWR2DVrWCAoiNhb//HW64wb0248fDhAnwyCPe\njS0QZRzO4Kr3riLjwQxaN2ztcv/8M/nE/jtWl3O0y4iIRxdodbLXqvWPf8CiRfDtt+632bYNkpJg\nxw6jvKNdcu2caxndaTT/Pfi/3W7z1s9v8WXGlyycttCLkblHKqvhaV5XWY7UyV4zzYUL0LUrfPUV\nDBzoWdsHHoD69eGll7wTWyBauncp9317H9se2Ea9OvXcbnex+CJdX+vKF1O+ICEiwYsRuuZIMJbG\nEGyq+pl7mux1zV6r0vvvw5VXep7oAZ54At57D06cMD2sgPX86uf54/A/epToAeqG1uXxoY/zlxV/\n8VJkWjDQyV6rVGEhPPssPPlkzdpHRBh1+3feMTeuQJVxOIP1B9Zza59ba9T+rn53sTZ3LRsPbjQ5\nMi1Y6GSvVerbbyE8HIYOrfkxHn3UqPkXFZkXV6B6Ne1V7h9wP/Xr1K9R+wZhDXgk8RFeS3vN5Mi0\nYKGTvVap//wHZsyo3THi4yEy0hi2GcyOnjvKZ9s+44H4B2p1nDvi7mBexjxOXThlUmRaMNHJXrvM\n/v2QlgaTJ7ve15WHHzbG3gezDzd/yMTuE2nXuF2tjtOucTtGdx7NJ1s+MSkyzZfOnDlD48aNueaa\nayw5v0722mXefRemToWGDWt/rEmTYP16yLlsCrzg8eHmD7m97+2mHOve/vfynw3/MeVYwazIgtri\nl19+SXR0NMnJyeTn5/v8/DrZa+UUFxsXVe+5x5zj1a9vfEKYM8ec4wWabYe2kX8mn5ExI0053tjO\nY8k/k68v1FZi/fr19OvXj6ZNmzJlyhRuvvlm/vSnPwGQnJxMZGQkf//73+nQoQN33303Fy9e5NFH\nHyUiIoKIiAgee+wxLjpu/X7vvfcYPnx4ueOHhISwZ88eAO644w7uv/9+xo0bR9OmTUlKSmL//v3V\nxvf+++9zzz33MHToUD766CMv/ASqp5O9Vk5KinEjVFycececPh0++ACCcXj2h5s/5LY+txEaEmrK\n8UJDQrkj7g4+2PSBKcezi4sXL3LjjTdy1113cfz4caZOncr8+fPL3QSWn5/P8ePH2b9/P2+++SbP\nPPMM6enpbNq0iU2bNpGens4zzzzj9jnnzJnDU089xZEjR4iLi+O2226rct99+/aRkpLClClTmDJl\nCh98YMH/nydTZHrjC5tPcRxo7r9fqb/9zdxjFhcr1bGjUuvXm3tcf1dcUqwiX4pUW/K3mHrcLflb\nVNRLUaq4pNjU47rD1d+r8ZZe+y9PLV++XEVERJTbNmzYMPUnx7zby5YtU3Xr1lUXLlwoe75Lly5q\nwYIFZY8XLlyoYmJilFJKvfvuu2rYsGHljici6pdfflFKKXX77berqVOnlj135swZFRoaqnJyciqN\n7+mnn1aDBw9WSil15MgRVadOHbVhwwa3XltVP3M8nOJY9+y1MsXFMG8e3HSTuccNCYFp08CCT66W\nStmXQuuGrbmy7ZWmHrd3m940qtuI9Nx0U49rBrPSvafy8vKIiCi/6mlUVFS5x23atKFu3brl2nTs\n2LHscXR0NHl5eW6dT0SIjIwse9yoUSNatmxZZfsPPviAmxx/WK1atSIpKYn333/frXOZRSd7rUxK\ninEzVJcu5h/7ppuMaReCqZQzL2MeN/Uy+Z0TI9FM6TWFz7d9bvqxA1WHDh3Izc0tt61iDb3ivD7h\n4eFkZWWV2z88PBwwkve5c+fKnjt48GC5tkopsrMvLbF95swZjh07Vtbe2erVq8nMzOSZZ56hQ4cO\ndOjQgTVr1jBnzhyKi4s9e6G1oJO9Vubzz2HKFO8cOzbWSPRbtnjn+P5GKcX8HfO54Qo3pwr10E29\nb+KL7V9Qokq8cvxAM2TIEEJDQ3n99dcpKiri66+/Zu3atdW2mTp1Ks888wxHjhzhyJEj/PnPf2b6\n9OkA9O3bl23btrFp0yYKCgqYNWvWZe1/+OEHVq1axcWLF/nTn/7E4MGDL/t0AcaF2XHjxpGRkVF2\nfWDr1q2cP3+eBQt8t9yHTvYa4L0STikRuPFGo3cfDNYfWE/9OvXp2bqnV47vz6UcK4SFhTFv3jze\neecdWrRowccff8zEiRPLlW0q9uyffPJJBg4cSGxsLLGxsQwcOJAnHfODdO/enaeeeooxY8bQo0cP\nhg8fXq69iHDrrbcye/ZsWrVqxYYNGyodYVNQUMDcuXN5+OGHadu2bdlXTEwM06dP9+mFWj3rpQbA\nypXw4IPG8oLekpJizHG/YYP3zuEvnlz6JIXFhTw39jmvnqOopIhnxzzrtXNUFEizXiYmJvLAAw9w\n++3m3OPg7M477yQyMpKnn37a9GNXpGe91Ez1/fcwcaJ3zzF0KOTmglOZ1La8WcIpdW23a/l+9/de\nPUcgSUlJ4eDBgxQVFfH++++zdetWxo8f75VzBcobnjOXyV5ExovIDhHZLSKPV/L8JBHZJCIbRORn\nERnl9FyWiGx2PKc/b/qx77+Ha6/17jlCQ+G66+w/V87uo7s5ev4oiZGJXj1PQkQCB88cZP/J6m/m\nCRY7d+4kLi6OFi1a8PLLL/PFF1/Qrl3tpqioiogE3EIu1ZZxRCQU2AmMAXKBtcBUpVSG0z6NlFJn\nHd/3Ab5SSnV1PN4LDFBKHavmHLqMY7H9+6F/f8jPNxKyN82fD//6l7H6lV29mvoqm/M3884k78/v\nPP2r6QyJHMJv4n/j9XNBYJVx7MJXZZwEIFMplaWUKgQ+BSY571Ca6B0aA0cqxuRuMJo1fvjBWDfW\n24keYNQoWLMGnEa12c7CXxYyodsEn5xLl3I0d7lK9hFAttPjHMe2ckTkBhHJABYAv3V6SgGLRWSd\niNxb22A17/juO++XcEo1bWp8ili+3Dfn87XzhedZsX8FozuN9sn5ru5yNcv3LedcoY3fPTVTuEr2\nbn1eU0rNV0r1BK4DPnR6aqhSqh8wAXhQRIZXegDNMufPG6Nkrr7ad+ccPx5+/NF35/OlFftXENsu\nlhYNWvjkfC0atKBf+34s27vMJ+fTAlcdF8/nAs73HEdh9O4rpZRaISJ1RKSVUuqoUuqAY/thEfkK\noyy0omI75xsWkpKSSEpKcvsFaLWzapWxzmzLlr475/jxxhTKdrQwcyHju3hnBEhVxncdz097fuLa\n7j76eKZZIjk5meTk5Bq3d3WBtg7GBdrRQB6QzuUXaLsAe5RSSkT6A3OVUl1EpCEQqpQ6LSKNgEXA\nbKXUogrn0BdoLfSHP0BYGPz5z747Z0kJdOgAqanQqZPvzusLvf/Vm3cnvUtCRILPzpmem87d39zN\nlt94//ZkfYHW93xygVYpVQQ8BCwEtgOfKaUyRGSGiJQuWjcZ2CIiG4BXgVsc29sDK0RkI5AGfFcx\n0WvWW7wYxozx7TlDQoyy0cKFvj2vt2WfzCb/TD4DOgzw6XkHdBhAzqkcDp456HpnLWi5HGevlFqg\nlOqhlOqqlPqbY9ubSqk3Hd//XSl1pVKqn1JquFJqrWP7HqVUnOPrytK2mv84dgx27oRBg3x/7rFj\njTcaO1mydwmjO482be56d4WGhDIyZiRL9izx6Xk192RlZRESEkKTJk1o0qQJ7du358EHH/T5aln6\nDtogtmwZDBsGTtOH+MzIkZCcbJR07GLp3qWMihnlekcvGN1pNIv32uzd04usWJbw5MmTnD59mi1b\ntrBmzRr++c9/+vT8OtkHMStKOKUiI42Lwlu3WnN+symlWJa1jJGdzFl+0FNjOo9hyZ4lQV1P9/dl\nCUu1adOGsWPHsn37dhNfvWs62QcxK5M9GL37pUutO7+ZMo9lopSiW8tulpy/e6vuKBS7j+225PxW\n8/dlCeHSfDp5eXksXLiQwYMH1+zF1pCroZeaTWVnw8mTxrBLq4wcCZ98Ao8+al0MZint1Vs1X4qI\nMLrTaJbsWUL3Vt0tiQFAZpvz+tVMzz6hpKamUlxczMMPPwzAjTfeSEJC+RFRISEhzJ49m7CwMMLC\nwpgzZw6vv/46rVu3BmDmzJnMmDGDP7s5NG3ixIkMGzYMgL/85S80a9aM3NzcSue0B8rOc/LkSYYM\nGcLkyZM9eo21pZN9kFq+HK66yhgZY5WRI+E3vzHm0vfFVA3etCxrGeM6j7M0hhEdR7BozyKfzZNT\nGU+TtFn8aVnCqpL90aNHCQkJoaCggKeeeoqrr76a1atXu3U+M+gyTpBavhxGjLA2hnbtIDw88Oe3\nV0qxbK919fpSI2JGsDxreVDW7f15WcKK6tevz+23305qairHjlU5R6TpdLIPUikpRs/eanao2+84\nsoMGYQ2IaR5jaRydmnciRELIPJZpaRxW8OdlCUuVvglfuHCBDz/8kA4dOtDSh7eu62QfhA4cgMOH\noU8fqyMxPl2suGwCjcCSsi+FER0t/piE0XMdETOClH0pVofic/66LKGz5s2bl42zT0tL45tvvjHx\nJ+CaXpYwCH3+OXz0Efj4d61SBw9Cr15w5Ii11w9qY9q8aSTFJHFP/3usDoW3fn6LlftX8sGN3lnb\nNJCmS9DLEpYXoH9eWm34Q72+VPv20KoVbNtmdSQ1t2L/CoZH+8eEriM6jmD5PpvOH+2CXpawejrZ\nByFvJvvS5doqflVn+PDALeXsP7mfgqICS4c7OuveqjsXii6QdSLL6lB8Ti9LWD1dxgkyx45BTIzx\nbx0vDLwVkcsWQRCq7wm9+66xTOEnn5gfj7d9vPlj5u2Yx5dTvrQ6lDJT5k7huu7XMb3vdNOPHUhl\nHLvQZRytRlavhsTEyhN9Vb1yb/dgSnv2gZhD/KmEU2po1FBWZa+yOgzNz+hkH2RWrYKhQ6t+XlXy\n5W1duhg3Vu3b54OTmcwvk320Tvba5XSyDzKukr0VRAKzbn/03FFyTuXQt31fq0Mpp2+7vmSdyOJE\nwQmrQ9H8iE72QeTCBVi/3pr5610ZNsx4Iwokq7NXkxCRQJ0Q/5p1JCw0jPjweNZkr7E6FM2P6GQf\nRNavh27doEkTqyO53ODBxvWEQLImZw1DIodYHUaldN1eq0gn+yDijyWcUnFxsGcPnDpldSTuW529\nmiFRfprsdd3eErNmzSqbcsHf6GQfRPw52YeFwYABkJZmdSTuKSwuZF3eOhIjE60OpVKDIwezLm8d\nhcWFXj1PdSO4zPpy18qVKxkyZAjNmzenVatWDBs2jHXr1nnx1V/On8feu0z2IjJeRHaIyG4RebyS\n5yeJyCYR2SAiP4vIKHfbar6jlFEm8ddkD4FVytmcv5mY5jE0r9/c6lAq1ax+Mzo178TGgxu9fq7K\nRnCZ9eWuU6dOMXHiRB555BGOHz9Obm4uM2fOpF69eia8Qvf58z0I1SZ7EQkFXgfGA72AqSLSs8Ju\ni5VSfZVS/YA7gLc8aKv5yN69xtj66GirI6nakCGBk+z9uYRTakjUEFJzUq0Owyd27dqFiHDzzTcj\nItSvX5+xY8fSp08ffvnlF0aNGkXr1q1p06YN06ZN4+TJk2VtY2JieOGFF4iNjaVJkybcfffd5Ofn\nM2HCBJo1a8bYsWM5ccIY2VS6ePjbb79NREQE4eHhvPjii1XGlZqaypAhQ2jRogVxcXEsX35pKov3\n3nuPLl260LRpUzp37sycOXO89wPCdc8+AchUSmUppQqBT4FJzjsopc46PWwMHHG3reY7a9YYPWer\nuPMRffBgo4wTCIuQr87x/2Q/KHIQa3KCY0ROjx49CA0N5Y477uDHH3/k+PHj5Z5/4oknOHDgABkZ\nGWRnZ5ebslhEmDdvHkuWLGHnzp189913TJgwgWeffZZDhw5RUlLCa6+9Vu54ycnJZGZmsmjRIp57\n7jmWLFlyWUy5ublMnDiRp556iuPHj/PCCy8wefJkjh49ytmzZ3nkkUf48ccfOXXqFGvWrCEuLs4r\nP5tSrpJ9BJDt9DjHsa0cEblBRDKABcBvPWmr+UZqqrVDLt35iN6mDbRtCz5eh7lG1mSvCYhkHyw9\n+yZNmrBy5UpEhHvvvZe2bdsyadIkDh06RJcuXRg9ejRhYWG0bt2axx57rFwPG+Dhhx+mTZs2hIeH\nM3z4cAYPHkzfvn2pV68eN954IxsqrLAzc+ZMGjRowJVXXsmdd97JJ5XM9fHRRx9xzTXXlE3GNmbM\nGAYOHMj333+PiBASEsKWLVs4f/487dq1o1evXt77AeF6WUK3ClBKqfnAfBEZDnwoIld4EoTzu2xS\nUhJJSUmeNNfckJoKt9xidRSuDR5sXEi2cm1cV/JO53Hm4hnLFhd3V/dW3TlRcIL8M/m0a+ydCcH8\nyRVXXMG7774LGJOiTZs2jUcffZRXXnmF3/72t6xcuZLTp09TUlJy2aIhzhOmNWjQoNzj+vXrc+bM\nmXL7Oy95GB0dzZYtWy6LZ9++fcydO5dvv/22bFtRURGjRo2iYcOGfPbZZ7zwwgvcfffdDB06lBdf\nfJEePXpU+fqSk5NJTk5274dRCVfJPhdwXsgxCqOHXiml1AoRqQO0dOznVtvKVoHRzHP+vNFb7t/f\n6khcGzTIKOXMmGF1JFVLzUklMTLRr0deAIRICImRiaTmpDLpiuCqoPbo0YPbb7+dt956iz/+8Y+E\nhISwdetWmjdvzvz588sWJq+Kqwut+/fvL0vM+/fvr3SFqujoaKZPn85bb71V6THGjRvHuHHjuHDh\nAk888QT33nsvKSlVLzxTsSM8e/bsamOsyFUZZx3QTURiRKQucDNQbskLEekijt96EekPoJQ66k5b\nzTd+/tlYIKRBA6sjca002fuztJw0BkX44W3IlRgUERylnJ07d/LSSy+VrUObnZ3NJ598wuDBgzl9\n+jSNGzemadOm5Obm8vzzz9f6fM888wznz59n27ZtvPfee9x8882X7TNt2jS+/fZbFi1aRHFxMQUF\nBSQnJ5Obm8uhQ4f4+uuvOXv2LGFhYTRq1IjQ0NBax1WdapO9UqoIeAhYCGwHPlNKZYjIDBEp7XtN\nBraIyAbgVeCW6tp652Vo1UlNtfbirCf69DEmRHMaLOF3UnNT/XZ8fUW+uEgrXvxyV5MmTUhLSyMx\nMZHGjRszePBgYmNjefHFF5k5cybr16+nWbNmXHfddUyePNnlp7KKSxBW3H/EiBF07dqVMWPG8Lvf\n/Y4xY8Zctm9kZCRff/01f/3rX2nbti3R0dG8+OKLKKUoKSnh5ZdfJiIiglatWrFixQreeOMND16x\n5/R89kHgv/4LfvUruPXW6verbC56cD0fvatjCJVf/KnquFddBU89BY6/H79SVFJEi+dakP1Ytt+O\nsXd2/Pxxol+J5vjjx02ZwyfY57PPysqic+fOFBUVEeKjdTT1fPaa29as8c/Jz6qSmGh8GvFHWw9t\nJbJpZEAkeoAWDVoQ1TSKLfmXX0DUgotO9jaXkwOFhdCpk9WRuM+f6/ZpOWkMigygd04gMTKRtFw/\n/YEGIH+/MF8VnextLj0dEhKMOeNrw5erVyUmGsneH6sFablpJEYERr2+VEJ4Aum56VaHYQsxMTEU\nFxf7rIRjpsCLWPNIWpqRPGvLl6tXRUZC3brGFA/+JjUnVffstYCkk73NlfbsA40/1u1PFpxk/8n9\nXNnWj+/4qkSftn3IOpHFqQsBNH+0Zjqd7G2suNgYYx8fb3UknktMhLVrrY6ivHV56+jfob/frUzl\nSlhoGHHt41iX59vpfjX/opO9jWVkQPv2UOHO8ICQkGB8KvEn6bnpJEQE4MckzK3b+2IOe/1l/rWx\nwOqiaB5JSwvMEg4YC5ls3GiMJAoLszoaQ3peOlOvnGp1GDWSGJnIZ9s+q/VxgnmMfaDTPXsbS083\n5+KsFZo0gZgY2LrV6kguCeiefYQekRPsdLK3sUDu2YN/lXJyT+VSWFxIx2YdrQ6lRjo178TF4ovk\nnKpyHkPN5nSyt6lz52DXLujb1+pIas6fkn1prz5Qb6gREd27D3I62dvUhg3QuzfUr291JDXnj8k+\nkCWEJ7A218+GOGk+o5O9TaWnB+aQS2d9+sCePXD6tNWRGBdnAz3Zx0fEk57nJ++ems/p0Tg2Ur7E\n8DHwE2+88R4QmKMo6taF2FhYvx5GjLAujhJVwrq8dcSHB/a7Z3x4PD/n/UyJKiFEdD8v2Oj/cZsp\nnc6gK/FsJd3rUxt4W3y89aWcnUd20rpha1o1bGVtILXUplEbWjRowe6ju60ORbOATvY2dIwWHKQ9\nV7DD6lBqLT7e+jtp1+atDfhefan48Hh9kTZI6WRvQ+sYyAB+JpQSq0OpNb9I9rn2SfYJEQmszdMX\naYORTvY2lE4C8djjD7p7dzh+HA4fti6GtXlrA/7ibCndsw9eOtnb0FriScAef9AhIcbUCessmsPr\nYvFFthzaQv8O/a0JwGQDwgew5dAWLhZftDoUzcdcJnsRGS8iO0Rkt4g8Xsnzt4nIJhHZLCKrRCTW\n6bksx/YNImKP7OPnFL7r2ftqQRMrL9Juyd9C5xadaVS3kTUBmKxx3cZ0at5JL1MYhKpN9iISCrwO\njAd6AVNFpGeF3fYAVymlYoGngbecnlNAklKqn1LKHp+D/VwuERRRh47s8/q5Ki5o4q2RP1bW7e10\ncbZUfES8rtsHIVc9+wQgUymVpZQqBD4FJjnvoJRao5Q66XiYBkRWOEZg3l8eoNYSTzxrbfVDL032\nVtwqYKeLs6Xiw+P1nbRByFWyjwCynR7nOLZV5W7gB6fHClgsIutE5N6ahah5ojTZ20lUlPFvdnb1\n+3nD2ry1xEfYMNnrnn3QcXUHrdt9KREZCdwFDHXaPFQpdUBE2gA/icgOpdSKim1nzZpV9n1SUhJJ\nSUnunlarYC3xPMKrVodhKpFLvfvoaN+d9+zFs/xy/Bdi28W63jmAxLaLJfNYJmcvnrXNtYhgkJyc\nTHJyco3bu0r2uUCU0+MojN59OY6Lsm8D45VSx0u3K6UOOP49LCJfYZSFqk32Wu2sY6DtevZwKdlP\nnuy7c244uIHebXpTN7Su707qA/Xq1KN3295sPLiRodFDXTfQ/ELFjvDs2bM9au+qjLMO6CYiMSJS\nF7gZ+MZ5BxGJBuYB05RSmU7bG4pIE8f3jYBxgB4C4FVdaMJp2nHI6kDc5u7SbFZcpLVjvb6ULuUE\nn2p79kqpIhF5CFgIhALvKKUyRGSG4/k3gaeAFsAbjqF3hY6RN+2BeY5tdYCPlVKLvPZKNAjAen3F\nOqFUsS0+3lg8vaTEGHvvC2vz1jKuyzjfnMzH4sPjWbx3sdVhaD7kctZLpdQCYEGFbW86fX8PcE8l\n7fYAcSbEqLkt8JK9u9q0gebNYfdu6NHDN+dcm7eWJ4Y/4ZuT+Vh8RDx/W/k3q8PQfEjfQWsr9k32\n4Nvx9sfPHyf/TD5XtL7CNyf0sZ6te3LgzAFOFJywOhTNR3Syt4miIoA4BvCz1aF4jS+T/bq8dfTr\n0I/QkFDfnNDHQkNCiWsfx7o8i+ah0HxOJ3ubyMgAyKM5J13tGrB8mezteOdsRXqZwuCik71NGEnQ\n/T9cX81rY6YBA2DTJigs9P65giHZ62kTgotO9jbhabIH38xrY6amTY2bqrZv9/651uWts92dsxXp\n4ZfBRSd7m6hJsg9EvpgB88DpA5wrPEen5p28eyKLdW7RmXOF5zh45qDVoWg+oJO9DVy4UNrb3Wh1\nKF7ni7p9aQknEEpbtSEiDAwfqOv2QUInexvYtMlY0QnOWx2K1yUk+CDZ2/jO2Yp0KSd46GRvA2vX\nGj3eYNC3L+zcCee9+L5mp2UIXdHJPnjoZG8DwZTs69eHnj1ho5cqVkopW05rXJX4CGNue2XFYgGa\nT+lkbwPBlOzBKOV46yLtnuN7aBjWkPaN23vnBH4mvEk49erUY++JvVaHonmZTvYB7vRpyMqCK6+0\nOhLf8eZF2mAYX1+RXrkqOOhkH+B+/tmoY4eFWR2J73g12QfRxdlSCREJum4fBHSyD3Dp6UZZI5j0\n6gV5eXBw4CgwAAAekElEQVTCC3N4BdPF2VIJEQmk53r55gXNcjrZB7hgTPahodCvH6wzeQ6v4pJi\nNhzcwIDwAeYe2M8NDB/I+gPrKSopsjoUzYt0sg9wwXZxtpQ3xttvP7ydiCYRNK/f3NwD+7nm9ZsT\n0TSCjMMZVoeieZFO9gEsP9+4QNu1q9WR+J43RuSk56YHXQmnVHx4vC7l2JxO9gGstFdv87v6K6WT\nvbn0RVr708k+gAVjvb5Ux47GVMe5ueYdMz0vuJO97tnbm8tkLyLjRWSHiOwWkccref42EdkkIptF\nZJWIxLrbVqud9PTgrNeD8WnGzN79ucJz7Dq6i77t+ppzwAAT1z6OHUd2cL7Q/vMrBatqk72IhAKv\nA+OBXsBUEelZYbc9wFVKqVjgaeAtD9pqNaSUUcYJ1p49mJvsNxzYQO82valXp545Bwww9evUp2eb\nnmw4uMHqUDQvcdWzTwAylVJZSqlC4FNgkvMOSqk1SqnStfDSgEh322o1l5kJjRtD++C4q79SZib7\nYK7Xl0qMSCQtJ83qMDQvcZXsI4Bsp8c5jm1VuRv4oYZtNQ8Ec72+VHy8Mda+pKT2xwrmen2pxIhE\n0vN03d6u6rh43u2p8ERkJHAXMNTTtrNmzSr7PikpiaSkJHebBq20NEhMdH9/Oy7E0aoVtGkDO3YY\nd9XWRnpuOjNHzDQnsACVEJHA7OWzrQ5Dq0JycjLJyck1bu8q2ecCUU6PozB66OU4Lsq+DYxXSh33\npC2UT/aae9LS4Kab3N+/4juvXVJ/aSmnNsn+8NnDHD13lO6tupsXWADq0boHx84f4/DZw7Rp1Mbq\ncLQKKnaEZ8/27I3ZVRlnHdBNRGJEpC5wM/CN8w4iEg3MA6YppTI9aavVzIULsHUr9O9vdSTWS0yE\n1NTaHSMtN434iHhCJLhHIodICPER8aTl6rq9HVX7262UKgIeAhYC24HPlFIZIjJDRGY4dnsKaAG8\nISIbRCS9urZeeh1BZdMm6NYNGjWyOhLrDRpkfMqpjbScNAZFDDInoACXGJGox9vblKsyDkqpBcCC\nCtvedPr+HuAed9tqtedpvd7O4uJg1y44e7bmb35puWn8NvG35gYWoBIjEnl97etWh6F5QXB/bg1Q\nOtlfUq+esXDLzz/XrH2JKiE9N53ECP0DhUt30pYoE4Y4aX5FJ/sAlJamh106S0yseSln55GdtGzQ\nUl+QdGjXuB3N6zdn99HdVoeimUwn+wBz5AgcOlT7oYZ2Upu6fVpuGoMidb3e2aDIQaTm1PKqt+Z3\ndLIPMKW9+hD9P1emNj37tJw0XcKpYFCETvZ2pFNGgFmzxujJapd07gwFBTWbATM1N5XESJ3snQ2K\nHMSanDVWh6GZTCf7AJOaqpN9RSI1G29/9uJZdh3dRVz7OO8EFqDi2sex+9huzlw8Y3Uomol0sg8g\nxcXGTJc62V9u8GDjU48n1uatJbZdLPXr1PdOUAGqXp169G3Xl3V5Ji/yq1lKJ/sAkpEB7doZc8Jo\n5dUk2a/JXsPgyMHeCSjA6Yu09qOTfQDR9fqqJSTAxo3GVBLuWpOjk31VdN3efnSyDyC6Xl+1xo2N\nKSQ2uLn2hlLKSPZROtlXprRnr5Tbk9dqfk4n+wCie/bV86SUk3ksk/p16hPZNNL1zkEoqmkUdULq\nsOf4HqtD0Uyik32AOHYMcnIgNtb1vsHKk2SvSzjVExGGRA3RpRwb0ck+QKSmGnXpOi6nrgteHiV7\nfXHWpaFRQ1mdvdrqMDST6GQfIFatgiFDrI7Cv3XtatxclVPpEjnlrclZw5Ao/QOtzpCoIazKXmV1\nGJpJdLIPEKtX62TviojxM1rtojN66sIpMo9l6pupXIhrH8cvx37h1IVTVoeimUAn+wBQWGgsrO18\ncVZELvvSYOhQWLmy+n1Sc1Lp36E/9erU801QAapuaF0GhA8gLUevXGUHOtkHgE2bICYGmjcvv11V\n+NKMZL/KReVh1f5VDIse5puAAtyQSF3KsQud7AOALuG4b+BA2LkTTp+uep+V2SsZGjXUd0EFsCFR\nQ/RFWpvQyT4ArF5t9Fg11+rVM5YqrDjlcVm5K1RYumMpE/tO1KUvNwyOGkxabhrFJcVWh6LVkstk\nLyLjRWSHiOwWkccref4KEVkjIgUi8j8VnssSkc3OC5FrnlEKVqzQyd4Tw4ZVXspRwNr20PskqAKf\nhxWQWjdsTXiTcDbnb7Y6FK2Wqk32IhIKvA6MB3oBU0WkZ4XdjgIPAy9UcggFJCml+iml9EJ6NbB3\nr5HwO3e2OpLAUd1F2lVRMGy/b+MJdMOjh7Ni/wqrw9BqyVXPPgHIVEplKaUKgU+BSc47KKUOK6XW\nAYVVHEN/Vq6FFStg+HBjWKHmniFDjDJOUdHlz62MhqE62XtEJ3t7cJXsI4Bsp8c5jm3uUsBiEVkn\nIvd6Gpx2Kdlr7mvVCiIjYXOFyoMCVkXD0OxKm2lVGN5xOCv2rdCTogU4Vzff1/Z/d6hS6oCItAF+\nEpEdSqnLugizZs0q+z4pKYmkpKRantY+VqyAhx+2OorAcemi6xsMGLADeLXsucyWEFoCnY5bElrA\n6tisI2GhYWQey6Rbq25WhxO0kpOTSU5OrnF7V8k+F4hyehyF0bt3i1LqgOPfwyLyFUZZqNpkr12S\nnw+HDsGVV1odSWBRwBxS+IL/Yp4j2QuQ0hGu2qfrip4SkbJSjk721qnYEZ49e7ZH7V2VcdYB3UQk\nRkTqAjcD31Sxb7m/IRFpKCJNHN83AsYBWzyKLsitXGlcbAwNtTqSwHMVKaRwFSVOv5YpHWHEPguD\nCmC6bh/4qk32Sqki4CFgIbAd+EwplSEiM0RkBoCItBeRbOAx4EkR2S8ijYH2wAoR2QikAd8ppRZ5\n88XYja7X11wkuTTjJBlcGjxW2rPXPFdat9cCl8sJc5VSC4AFFba96fT9QcqXekqdAfRMU7WQkgKv\nv251FIFrBMtJ4Sp6sx2awdm60OOI1VEFpl5tenG84Di5p3KJaOrJGA3NX+g7aP3UsWOQmQnx8VZH\nEriuIoXljDAe6Hp9rYRICCM6jiA5K9nqULQa0sneT61YYSzGERZmdST+obJZPl3N/Fnas1cAHWFE\nlhWR20dSTJJO9gFMJ3s/lZwMegTqJRVn+FRVbHcWQxahFJNJ17KevVZzI2NGkrwv2eowtBrSyd5P\nLVsGI0daHUVgE2Aky5jXJAEaQp9DVkcU2Hq37c2JghPknHJ79LXmR3Sy90PHjsGePTBggNWRBL5R\nLOWrTs0hC0L0DaC1ouv2gU0nez+UkmLM76Lr9bU3kmVs6nQM9lodiT3oun3g0sneDy1bpuv1Zolm\nP0WdUmBvV6tDsYWkmCSWZS2zOgytBnSy90NLl+p6vVn2toCw0DNwZLzVodhC7za9OXPxDFknsqwO\nRfOQTvZ+5uBByMkxltfTam9pJ+i3twUw2upQbEFEGN1pNEv2LLE6FM1DOtn7mSVLjF69ng/HHEs7\nweS9x4ERFOtfd1OM6TyGxXsXWx2G5iH92+9nfvoJxoyxOgp7KBEj2d+49xSQx89cPrzJ1Y1Z2uXG\ndB7Dkj1LKFElVoeieUAnez+iFCxeDGPHWh2JPWxpC00uQKcTAIv4ict/sNXdlKVVLrpZNC0atNDr\n0gYYnez9yM6dRvmmqx44YopFXWDcL2WPWMQ4K8OxlTGdxrB4jy7lBBKd7P3I4sVGCUdXEsyxsCtc\nXZbsl7Oe/pyiict27sy7E+zGdNbJPtDoZO9HdL3ePOfCIDUSkrJKt5wnkTSSSXLZtqp5eLRLkmKS\nWJ29moKiAqtD0dykk72fuHDBmPxM1+vNkdIR+h+AphcubRunSzmmadGgBX3a9SFlX4rVoWhu0sne\nT6xcCT17QuvWVkdiD+Xr9Qad7M01oesEFuxe4HpHzS/oZO8nFiyACROsjsI+FlaS7GPZzCmasodO\n1gRlMxO6TmBBpk72gUInez+hk715sprDkYYwMK/89hAUE1jAD1xjTWA2069DP44XHGfvcT3LXCBw\nmexFZLyI7BCR3SLyeCXPXyEia0SkQET+x5O2mmH/fjh0SE+RYJbvu8GEzMqnNL6W7/mea30flA2F\nSAjju47XvfsAUW2yF5FQ4HVgPNALmCoiPSvsdhR4GHihBm01jF791VdDiP6cZYrvusO1uyp/biw/\nsYqhnKWhb4OyqQldJ/Bj5o9Wh6G5wVV6SQAylVJZSqlC4FNgkvMOSqnDSql1QKGnbTXDDz/oEo5p\nwmBl9OX1+lLNOMVA1rGUUT4Ny911cwPNuC7jWL5vOecLz1sdiuaCq2QfAWQ7Pc5xbHNHbdoGjXPn\njPnrdbI3SSeIz4NmF6reZSLf8R0TfReTgx3H77ds0JK49nEs3bvU6lA0F+q4eL42v49ut501a1bZ\n90lJSSQF0codixcbtfqWLa2OxCaqKeGUupbveZnHfBNPEJjUYxJf7/yaa7vrayHelJycTHJyco3b\nu0r2uUCU0+MojB66O9xu65zsg83XX8P111sdhT2UCNAdJr5f/X7d2UUDzgP9gA1uH9+TkotS5vbb\nqzu32efy1PU9ruf51c9TokoIEX3hyVsqdoRnz57tUXtX/zPrgG4iEiMidYGbgW+q2Lfib6MnbYNS\ncTF8951O9mZJjwAuQI+j1e8nwA3MB2706PiVlWB8WZrx1zJQ15ZdadmgJWtz11odilaNapO9UqoI\neAhYCGwHPlNKZYjIDBGZASAi7UUkG3gMeFJE9otI46raevPFBJq0NGjbFjp3tjoSe/jqCsDN37Ab\n+QpPk71Wteu7X883O3Vfzp+J1R8BRURZHYNVfv97Y0rjv/yl6n2q/fhecd9KtlW13d1t/rCvO+0V\n0ONh2P0lqLzq9wUoQQglh12MoBuZpsVatr3C77SIuL1vZWrb3ttSc1K5+5u72fbANqtDCRoiglLK\n7dqiLrBZRCn44gv41a/c2LeSL6287W2goA5QSaKvTAgK+JqvdO/eFAkRCZwoOMGOIzusDkWrgk72\nFtm4EUpKoH9/qyOxh696wo0eFwm/Yh5uvNtqLoVICJN7TmbutrlWh6JVQSd7i/Tv/1f27n2OkBD7\n3GBjpS97wo0edyqT2UV3sok0PZ7a3Djl6Y1X/nKj1k29bmLudp3s/ZVO9hYwSqw3sZa5ujRjgp2t\nIL8xDN/nactCbuQr5nKT6THVtuzmSXt/KfENjR7KkXNH2Hlkp4VRaFXRyd4CmzYB1GEAP1sdii18\neiVM2QahNch0t/Apn3KL+UEFobJSju7d+yWd7C3w+ecAcy+7MQFq9/E/GCmMZH/L1pq1H8ky9tGR\nX9DjX80wpfcUPt/2udVhaJXQyd7HSkpgzhyAOZU+7y8fyQPF5nZwPgwS3b2vu4I6FHMTc/mMm80N\nLEgNjR7KiYITbMnfYnUoWgU62fvYypXQpAnAJqtDsYXSXn1tPgPdwqd8wlTTYgpmIRLCbX1u46PN\nH1kdilaBTvY+9uGHMG2a1VHYhMDHsXBrLTuRQ1jNKZqykb7mxOWnfDXF8rTYaXy85WOKS4pNP7ZW\nczrZ+1BBAXz5Jdx2m9WR2EQnaH0OYvNrd5gQFLfzPu9xhylh+TNflAl7t+1N20ZtSc5K9tIZtJrQ\nyd6HvvsO+vWDSPOHdQenOLjT/Ukrq3U77zOHW4Ewcw4Y5KbHTuejLbqU4090svehd9+F6dOtjsIe\nTtYDusPUGo7CqagLe+hJBuj1aU0xtc9U5u+Yz+kLp60ORXPQyd5H9u+H1FSYMsXqSOzh897AHqOM\nY5Y7eA+CoJTjC+0bt2dkzEg+2fqJ1aFoDjrZe1npRbCOHWdx7NjrNGqkx8+b4Z3+wEZzj3kTc4Gr\nyKODuQd2wa73Vtw34D7e+vktq8PQHHSy94FCQonkbjbyth4/b4KfO8CBxsBuc4/bmLPAp/yHe8w9\nsAt2vbdibOexHDl3hPUH1lsdioZO9j7xI+OJIJe+bLY6FFv4Vzzcvw4vZcY3eIv7KHS5YqfmSmhI\nKPf0v0f37v2ETvY+8E8e5H7+bXUYtnC8PszrCXebNArnclvoxF6+5TpvncAn/GUmzLv63cXn2z7n\nRMEJS86vXaKTvdddwQb6cQufWh2ILbwXB9fuhrZnvXeOB/kn/+RB753AB/xlwZvwJuFc0+0a3ln/\njkURaKVcJnsRGS8iO0Rkt4g8XsU+rzme3yQi/Zy2Z4nIZhHZICLpZgYeOB7hfv5NfS5YHUjAKwqB\n1xLhIS//Jv2KeezgCjbTx7snChKPDnqU19Jfo6ikyOpQglq1yV5EQoHXgfFAL2CqiPSssM81QFel\nVDfgPuANp6cVkKSU6qeUSjA18gBw9CjAzfym3I9Eq6kve0LkKRhUw0nP3FWXQn7La7zA/3r3REFi\nYPhAoppGMX/HfKtDCWquevYJQKZSKkspVQh8CkyqsM/1wPsASqk0oLmItHN63j5jyTz09tsA82nH\nIatDCXgK+PtQ+N1q35xvBm/yPdeynyjfnNDmHhv0GK+mvWp1GEHNVbKPALKdHuc4trm7jwIWi8g6\nEbm3NoEGogceAPij1WHYwrJOcC4MJu7yzfmac5I7eZdXeNQ3J7S5SVdM4t1J71odRlBzlezdva5T\nVe99mFKqHzABeFBEhrsdmQ00bQpw0OowbOFvw4xefYgPrzQ+yiu8z+0coZXvTmpTdULq0LVlV6vD\nCGquBhPnQrnPsVEYPffq9ol0bEMplef497CIfIVRFlpR8SSzZs0q+z4pKYmkpCS3gteCREfY0wKm\n+3gJgEhyuZnP+Dv/D6h0bIKm+UxycjLJyck1bi9KVd1VEpE6wE5gNJAHpANTlVIZTvtcAzyklLpG\nRAYBryilBolIQyBUKXVaRBoBi4DZSqlFFc6hqosh0IlIpR+PhMs/Nrm7zVv7+vJc7u6rgJA74f31\n8OtN1e9rRlwVt+USTiybOUZvFPnV7uvLuEzd18Z/f3YmIiil3L4mWm0ZRylVBDwELAS2A58ppTJE\nZIaIzHDs8wOwR0QygTeBBxzN2wMrRGQjkAZ8VzHR242/3MhiJ4s7A43gNotuPo4gj1/zAfAHawLQ\nNJNU27P3SQA26tlX1ov36x6dhedyZ98Sgfh7Yf0qUNt8E1dl2/JpS3u2kUkiXdhj2rlqG5dp+9rk\n7y/YmNqz1zQrvRcH9YuAShK9LxlDZ1/idzxvbSCaVgs62Wt+6VQ9eGIUvPKj1ZGUeokN9GMZSVYH\nomk1opO95peevgqu/gXi86yOpNQFnud3PMKrekZMLSDpZK/5nXXh8EFfeO4nqyMpbzJfEkEuz/M7\nq0PRNI/pC7Qm0hdoTdg3BGLvg/9dDdM3+z4uV+33Ec0AfmYlw+jJTr+Jqzb7Vqayv8mqRpfZ5e83\n0OgLtFpgGw7hp2Gan67z0pH9zGYmd/F/QKjV4ZjCk6mQ/WHaZK1mdLLX/MbKaCAe3v7Wv2fP+w1v\n0JgzwFNWh6JpbtPJXvMLxxrAbb8CvjGmMfZnISg+ZDpwD0sZaXU4muYWnew1yxWFwLRfweQMwEez\nWtaWMfb+dqbxkZ4GWQsIOtnXQGXTIuipEWruf8dBYYj/jb5xbTH/w4tczzecoZHVwWhatXSyryF/\nWeMz0P0zHn7sCnPnQliJ1dF47r95iYGsYyqf6PH3ml/TyV6zTj94bhj88DE0L7A6mJoR4F+Ouf+m\n8RFFNhmho9mPTvaaJd7vC4yExR9A5+NWR1M7dSlkLjdxjJbcwXu4XibC/+kypf3oZK/5lAKeGwpP\njQQ+hO5HrY7IHPW5wNdM4jgtgK8Dvoavy5T2o5O95jMFdWDGdfBRLKz6P+Cw1RGZqyHnmc8NwEFG\nsJwsOlodkqaVsX2y1yNn/MOeFjDsLjhe30j0/j6WvqbCKALuZhofkUga33Gt1SF5XVV/Y/rvzr/Y\nfm6cKuerqcU5a7vUoCf7+tO8MDXZt0QgNB5aj4AnU+C3aZfujrX7z2AVQ5jGRySRzHv8N4oTfhGX\npefS8+iYRs+No/mNFdGQeA9wpdGbf8Qp0QeDoaxmC31oxFlgO/9mhh6to1kmaHv2tWWLXpaXzpUS\nDX8dDtvbwLOL4dZtUNl/sZ1/BpfvG0cSL5FHOL/nWW7jY+pS6AdxBU7PvroykNV5zAqm9+xFZLyI\n7BCR3SLyeBX7vOZ4fpOI9POkrVWqGm3g7jatvBP14a0BwH1w5w3wqwzY9Q+YuhX9QwNgI0sZxRv8\nhjncSgxZPMEzQBerAwsoepRQzVU7IFhEQoHXgTFALrBWRL5RSmU47XMN0FUp1U1EEoE3gEHutK2J\njz/8kNSUlMu2R3XuzP/7wx9qc2jNQ/uawWvtYNcAWB4DY/YAS2HnLxBqk7/CZDBtIUIBRrGMUSxj\nOz15i/uAVQwghxv5iqtZyAB+xlcpLNknZ7FOcnIySUlJVofhN6ot44jIYGCmUmq84/HvAZRSzzrt\n829gmVLqM8fjHRh/H51ctXVs96iMc9eUKRTMnctgp237gReqaRNQH3P99FwHGxllmY3tjZWk5kRD\nmzBo/y08HgbX7YKmF+z3M5jl+PJeXKEs4Sq+5ToWcjUH6MAJUplFKnFs5Eq2EkMWdSgx/WcwC5jt\ntdflnTJOleeqcNxgKPl4WsZxdatfBJDt9DgHSHRjnwgg3I22NTIGuMvp8TqMZF/VL4JWiVA4Hgan\n6xmLe5+oDzSAdxrCoUZwoAnkNDV677SEXiXQ6zDEHYTRe2HOcsg/aiSL2yx+KYGtuKy3D5BPW9oz\nhIsM5E1msJ1e5NMOyGICWUSRTTh5tCOfNhwGjrKBEzTjJE047ZhnP0DnnjDZTIw3NGfBnA9cJXt3\n3wJ99zMMDeXVBg34KiysbNOJ4mI4e9ZnIbhrTh9AwcTYyp+f6PjX+Yd8LVD6Xl26fbw46pNO/wKM\nEmNoY+kXAgkhUCxQHGJMHVwYAoRCZChcCDVubCpw/K93KoQmF4weeYsCoABWnoU256DLMRi+D2JO\nQMIJOHqufOx3E9x/ON5iTJ08n78wv2zbWRrSmBgepBO5RJBLBBuJ4zBtgBbcSQtO0oxTNOUsjYAw\nGlFAfQqoxwXqcpG6XAQK6UsRoRQTSjEHeBO4k6soJoQSBEUIJYBiLAqp8AVwjeNf5z72tRUel5p4\n2ZZM035OmudclXEGAbOcSjF/AEqUUs857fNvIFkp9anj8Q5gBEYZp9q2ju32+EylaZrmY2aWcdYB\n3UQkBsgDbgamVtjnG+Ah4FPHm8MJpVS+iBx1o61HwWqapmk1U22yV0oVichDwEKM1ZXfUUpliMgM\nx/NvKqV+EJFrRCQTOAvcWV1bb74YTdM0rXKW31SlaZqmeZ+l0yX4801XtSUiUSKyTES2ichWEfmt\n1TGZTURCRWSDiHxrdSxmE5HmIvKFiGSIyHZHidI2ROQPjt/NLSIyR0TqWR1TbYjI/4lIvohscdrW\nUkR+EpFdIrJIRJpbGWNtVPH6nnf8fm4SkXki0qy6Y1iW7J1uuhoP9AKmikhPq+LxgkLgMaVUb2AQ\n8KDNXh/AI8B27Hkj46vAD0qpnkAsYJsSpOM62r1Af6VUH4wy6y1WxmSCdzFyibPfAz8ppboDSxyP\nA1Vlr28R0Fsp1RfYBVR7V6mVPfsEIFMplaWUKgQ+BSZZGI+plFIHlVIbHd+fwUgW4dZGZR4RiQSu\nAf6DzUZhOnpIw5VS/wfG9Sel1EmLwzLTKYzOSEMRqQM0xLjLPWAppVYAFdc8ux543/H9+8ANPg3K\nRJW9PqXUT0qp0pWb04DI6o5hZbKv6mYs23H0pPph/IfYxcvA74AAXCbcpU7AYRF5V0TWi8jbItLQ\n6qDMopQ6BryIcfN5HsYIusXWRuUV7ZRS+Y7v84F2VgbjZXcBP1S3g5XJ3o4f/S8jIo2BL4BHHD38\ngCciE4FDSqkN2KxX71AH6A/8SynVH2OUWSCXAMoRkS7Ao0AMxqfNxiJi6xuhHXOy2DLniMgTwEWl\n1Jzq9rMy2ecCUU6PozB697YhImHAl8BHSqn5rvYPIEOA60VkL/AJMEpEPrA4JjPlADlKqbWOx19g\nJH+7GAisVkodVUoVAfMw/k/tJl9E2gOISAfgkMXxmE5E7sAop7p8s7Yy2ZfdsCUidTFuuvrGwnhM\nJcZMTO8A25VSr1gdj5mUUn9USkUppTphXNhbqpT6tdVxmUUpdRDIFpHujk1jgG0WhmS2HRgz0zZw\n/J6OwbjQbjffALc7vr8dsFOHCxEZj1FKnaSUcjkhkmXJ3tGjKL3pajvwmc1uuhoKTANGOoYnbnD8\n59iRHT8ePwx8LCKbMEbj/NXieEyjlNoEfIDR4drs2PyWdRHVnoh8AqwGeohItojcCTwLjBWRXcAo\nx+OAVMnruwv4B9AY+MmRX/5V7TH0TVWapmn2p9eg1TRNCwI62WuapgUBnew1TdOCgE72mqZpQUAn\ne03TtCCgk72maVoQ0Mle0zQtCOhkr2maFgT+P/NsekoUh5eRAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xedb3fd0>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Just from looking at this plot, we suspect that we will have to reconcile the samples in the overlap region since these could have come from either group. This is where the *Expectation Maximization* algorithm enters."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Expectation maximization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The key idea of expectation maximization is that we can somehow pretend we know the unobservable $z$ value and the proceed with the usual maximum likelihood estimation process."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The idea behind expectation-maximization is that we want to use a maximum likelihood estimate (this is the *maximization* part of the algorithm) after computing the expectation over the missing variable (in this case, $z$). \n",
      "\n",
      "The following code uses `sympy` to setup the functions symbolically and convert them to `numpy` functions that we can quickly evaluate. Because it's easier and more stable to evaluate, we will work with the `log` of the likelihood function. It is useful to keep track of the *incomplete log-likelihood* ($\\log\\mathcal{L}$) since it can be proved that it is monotone increasing and good way to identify coding errors. Recall that this was the likelihood in the case where we integrated out the $z$ variable to reconcile as its absence. \n",
      "  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sympy\n",
      "from sympy.abc import x, z\n",
      "from sympy import stats\n",
      "\n",
      "mu_a,mu_b = sympy.symbols('mu_a,mu_b')\n",
      "na=stats.Normal( 'x', mu_a,1)\n",
      "nb=stats.Normal( 'x', mu_b,1)\n",
      "\n",
      "L=(stats.density(na)(x)+stats.density(nb)(x))/2 # incomplete likelihood function "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next, we need to compute the expectation step. To avoid notational overload, we will just use $\\Theta$ to denote the $\\mu_b$ and $\\mu_a$ parameters and the data $x_i$. This means that the density function of $z$ and $\\Theta$ can be written as the following:\n",
      "\n",
      "$$ \\mathbb{P}(z,\\Theta) = \\frac{1}{2} \\mathcal{N}_a(\\Theta) ^ z \\mathcal{N}_b(\\Theta) ^ {(1-z)} $$\n",
      "\n",
      "For the expectation part we have to compute $\\mathbb{E}(z|\\Theta)$ but since $z\\in \\lbrace 0,1 \\rbrace$, this simplifies easily\n",
      "\n",
      "$$ \\mathbb{E}(z|\\Theta) = 1 \\cdot \\mathbb{P}(z=1|\\Theta) + 0 \\cdot \\mathbb{P}(z=0|\\Theta) =  \\mathbb{P}(z=1|\\Theta)  $$\n",
      "\n",
      "Now, the only thing left is to find $  \\mathbb{P}(z=1|\\Theta) $ which we can do using Bayes rule:\n",
      "\n",
      "$$  \\mathbb{P}(z=1|\\Theta)  = \\frac{ \\mathbb{P}(\\Theta|z=1)\\mathbb{P}(z=1)}{\\mathbb{P}(\\Theta)} $$\n",
      "\n",
      "The term in the denominator comes from summing (integrating) out the $z$ items in the full joint density $ \\mathbb{P}(z,\\Theta) $\n",
      "\n",
      "$$ \\mathbb{P}(\\Theta) = (\\mathcal{N}_a(\\Theta) + \\mathcal{N}_b(\\Theta))\\frac{1}{2} $$\n",
      "\n",
      "and since $\\mathbb{P}(z=1)=1/2$, we finally obtain\n",
      "\n",
      "$$  \\mathbb{E}(z|\\Theta) =\\mathbb{P}(z=1|\\Theta)  = \\frac{\\mathcal{N}_a(\\Theta)}{\\mathcal{N}_a(\\Theta) + \\mathcal{N}_b(\\Theta)} $$\n",
      "\n",
      "and which is coded below.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ez(x,mu_a,mu_b): # expected value of hidden variable\n",
      "  return norm(mu_a).pdf(x) / ( norm(mu_a).pdf(x) + norm(mu_b).pdf(x) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, given we we have this estimate for $z_i$,  $\\hat{z}_i=\\mathbb{E(z|\\Theta_i)}$, we can go back and compute the log likelihood estimate of\n",
      "\n",
      "$$ J= \\log\\prod_{i=1}^n \\mathbb{P}(\\hat{z}_i,\\Theta_i) = \\sum_{i=1}^n \\hat{z}_i\\log \\mathcal{N}_a(\\Theta_i) +(1-\\hat{z}_i)\\log \\mathcal{N}_b(\\Theta_i) +\\log(1/2)  $$\n",
      "\n",
      "by maximizing it using basic calculus. The trick is to remember that $\\hat{z}_i$ is *fixed*, so we only have to maximize the $\\log$ parts. This leads to\n",
      "\n",
      "$$ \\hat{\\mu_a} = \\frac{\\sum_{i=1}^n \\hat{z}_i x_i}{\\sum_{i=1}^n  \\hat{z}_i } $$\n",
      "\n",
      "and for $\\mu_b$ \n",
      "\n",
      "$$ \\hat{\\mu_b} = \\frac{\\sum_{i=1}^n (1-\\hat{z}_i) x_i}{\\sum_{i=1}^n  1-\\hat{z}_i } $$\n",
      "\n",
      "Now, we finally have the *maximization* step ( above ) and the *expectation* step ($\\hat{z}_i$) from earlier. We're ready to simulate the algorithm and plot its performance!\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Lf=sympy.lambdify((x,mu_a,mu_b), sympy.log(abs(L)),'numpy') # convert to numpy function from sympy\n",
      "\n",
      "def run():\n",
      "    out, lout = [], []\n",
      "    mu_a_n=random.random() * 10 # initial guess\n",
      "    mu_b_n=random.random() * 10 # initial guess\n",
      "    for i in range(20): # iterations of expectation and maximization\n",
      "        tau=ez(xs,mu_a_n,mu_b_n)                 # expected value of z-variable\n",
      "        lout.append( sum(Lf(xs,mu_a_n,mu_b_n)) ) # save incomplete likelihood value (should be monotone)\n",
      "        out.append((mu_a_n, mu_b_n))             # save of (pa,pb) steps\n",
      "        mu_a_n=( sum(tau*xs) / sum(tau) )        # new maximum likelihood estimate of pa\n",
      "        mu_b_n=( sum((1-tau) * xs) / sum(1-tau) )\n",
      "    return out, lout\n",
      "\n",
      "out, lout = run()\n",
      "\n",
      "fig=figure()\n",
      "fig.set_figwidth(12)\n",
      "ax=fig.add_subplot(121)\n",
      "ax.plot(array(out),'o-')\n",
      "ax.legend(('mu_a','mu_b'),loc=0)\n",
      "ax.hlines([mua_true,mub_true],0,len(out),['r','g'])\n",
      "ax.set_xlabel('iteration',fontsize=18)\n",
      "ax.set_ylabel('$\\mu_a,\\mu_b$ values',fontsize=24)\n",
      "ax=fig.add_subplot(122)\n",
      "ax.plot(array(lout),'o-')\n",
      "ax.set_xlabel('iteration',fontsize=18)\n",
      "ax.set_title('Incomplete likelihood',fontsize=16)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "<matplotlib.text.Text at 0xf16bab0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAEjCAYAAADuR70GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VNX9//HXhwAJIFBRRFAUjWDd910xVQERtVptRa17\nrdUa0m9bW2Up/Kq2bq0SqlbcsYLaClqNIogGsIq44L6AURQQAVFBloQsn98f9wYmmck2mclkJu/n\n43EfM3Pmnns+M9GbDydnMXdHRERERESSp12qAxARERERyXRKukVEREREkkxJt4iIiIhIkinpFhER\nERFJMiXdIiIiIiJJpqRbRERERCTJlHSLiIi0AmZ2oZlVmdmuqY4l2cysX/hZL4ij7m/M7PQExzPO\nzKpqlVWZ2Z9qn2Nmzc6dzGyxmd0f8br6Z79TrXMeam5biRArPmm69qkOQERERNqseDYL+Q0wB5iW\n5FgOB5Y2cE5z2oq81tNhe1/Vc46kOSXdIiIikm4s2dd09/kt1C7u/jXwdTKuLa2HhpeIiIi0UmZW\nbGZzzewEM3vTzNab2btmdlqMc/czs2lm9rWZbTCzj8zs6lrn/J+ZfWxmZWb2pZlNMLOutc6pMrNr\nzewqM/vCzNaZ2dNm1tPMepvZ42a2xsw+N7M/1KpbPQzhGDN7wsy+D+P5h5nlNOLzHmtms8xsbdju\ndDPbK+L9xcBOwLlhO1Vmdl+t7+C/ZvZN+B28ZGZHN/oLj/4exjZwzolhnIVmZmHZT8xsXviz+tbM\nHjOzvg1cp67hG2Zmw83sw7Cd18zsqBj1f25mb5vZRjNbZWaTzGz7Wud0MLPrwmErZWb2Wfhzbl/r\nvF3NrCiMf6WZ3QZk1xe/NI6SbhERkdbLgVzgNuAW4CfAcuDfZpZbfZKZHQq8AuxCMPziJODvwA4R\n5/wF+BvwHHAycBNwIVBUnTBGOB8YCFwG5APHAP8CngTeAE4HngVuMLOhMeL+F7AwPO9W4FLgzvo+\nqJkNA2YBa4FzgXOArsBcM9sxPO00giEY0wmGYxwOXBvWPxB4GfgB8AvgDGA18Hz4XjzqHN5hZucT\nfB9/cfcR7u5m9ivgP8B7YfuXAXsDs81sqya2bQTf+/8Bo4CzgCzgaTPrHhHHL4FJwPsE3/fVwJCw\nzS4R13sQ+CPwADAsfPxjWF59rY7ATGA/4AqC/z52AUY3MXaJxd116NChQ4cOHSk+CBKcKmDXiLJi\noAzIjSjrCVQA10SUzQE+B3LquHaP8Dr31So/N2zzlIiyKuAjoF1E2d/C8pERZVnAishrRnyGO2q1\nMzKMuX/4ul943vkR53wCzKxVryuwCrg1ouwzYFKMzziLIPFsH1HWDvgAmNbAdz8OqKpVVgX8qfY5\n4ef+A7AJuDji/a2ANcA9ta7TL/zuC2p9hljf204RZYsJ/tHQPaLsoPC8s2v9DGbVavOo8Lz88PXe\ntT9PWD4qLN8nfH1p+PrQiHOM4B8RlZHx6Wj6oZ5uERGR1m2Ru5dUv3D3VcBKoC+AmXUGjgQedvfS\nOq5xONCBoAc60qMEyfDAWuUz3T1yNY+Pw8fnIuKoJEiUdyTaYzHaaQccEis4M+sP7ApMNrP21Qew\nEZgXI77a9TuF5/w7fF1dvx1BMl5v/Sa6jSABP8Pd74soP4LgHwm1P8NSgu8vnhhecfc1Ea/fCx+r\nh6vsTvCPsIcjK7n7/wj+EVbdZvVj7Z//v2q9fwTwhUeMZ/cg8/43SRrP3pZoIqWIiEjr9k2MsjKg\neoz01gTJZe2VNiL1CB+XRxa6e4WZrY54v9q3tV5vqqO8PCKOSCvqeL1D7RND24WP94ZHbZ/XUa9a\nD4Je3z+FR22JXAVkOPAuQTIfqfozPF9HvdVNbMep9bN397JwJFD1dx7z5xpaEfF+XeetqPV+b6J/\ndtRRJk2kpFtERCS9fUswJCBWj3O16uStN/BhdWHYE7sNsRP75tg+sh2gV/i4rI7zqxPSq4mdtG6K\nURbpO4Lv4B8E45uT6TiCcc/PmtlJ7r4+LK/+DBcQDHOp7fskxBL5c61te+C1GOd9WuucyPeXA3vG\nuFavGGXSRBpeIiIiksbcfQPwEvDzelYIeYUgcR1eq/wsgg644gSH9bNar4cTJMWv1nH+xwRjmPd2\n9zdjHO9FnFsGdI6sHCa+c4H9gQWxrpGIDxV6H8gD+hMk3tWTFf9HkFj3r+MzLEpgDNU+JuiFrvFz\nNbMjCVZ5KQ6L5oSPtX/+54aP1ee9DPQ1s8MirtWO4OepNcObST3dIiIirVussbS1y34PzAZeMbO/\nEfQo7wrs58HKGt+G5deY2XqClUf2IFj5Y667FyU4vqFmdhNBj/ChBEM+Howcmx7J3d3Mfg08Ga6g\n8W+Cdat7EYxX/9zdbw1P/wA4JlztZAWwyt0/B35LkFw+Z2b3Eqxysi1wIMGk0Gua8Rlrx/uRmeUB\nL4btneju35vZVcDtZtaTYIWVNQRDao4FXnT3KeElGjM+usFz3L3Sgl0z77Jg98qHw/auJ1g95r7w\nvPfMbAowLvzrxisE47dHA5Pdvbpn/kGCvzZMNbORBJNYf0UwVl1juptJPd0iIiKtR+3exLp2JaxR\n5u6vE6xYsQSYABQBvwtfV58ziiAxHQo8RbACx4MEy8fFE1t98f0cGABMJVjybiLBEnR1X9z9WYIJ\nfV2AuwmS1hsJxkq/HHHqNQQ9vI8B84GxYf0FBBM1VwOFBJM+bwP2IvgHSUOfraGe3BrnuPtCgmR6\nZ4LEu6u7TwROJZjgOIng5zCWIN9aUOtasa5f3+vYQbnfDZwH7AM8QfCdPQcc6+4bI069MHzv4jCu\ni4AbCIbDVF+rHBgEvAXcQbCsYAlwXWPjkbpZMCm19TOz3YFHIop2Bca4e2GKQhIRyRhm9lOCFRl+\nCBxS/ef4sNfxLrYsVVbg7rPD9w4i+KWcAzzj7gVheTZBwnEgQQJ0VtgTKRnOzC4k6F3dzd0/beB0\nkTYlbXq63f1jdz/A3Q8guPlvAKalOCwRkUzxLsHGGnNqlV9KsH7xvgQ9YH+LeO9O4BJ37w/0N7MT\nw/JLgNVh+a0EvWsiIm1a2iTdtZwAlLj7kgbPFBGRBrn7R+Gfy2vbg2DcavX60N+Z2SFm1hvoGrGe\n7ySC3QIh+PN69S53jwPHJy9yaYXS40/oIi0sXZPu4cDkVAchItIGvA2camZZZrYLwV8adySYrBW5\nLvQytqzBvAPhWGJ3rwDWmFntdaAlA7n7A+6epaElItHSbvWScHzhKcAfUx2LiEg6MbOZbFmXN9JI\nd3+qjmr3EfR2v06wQcnLBNtBqzdTRKQJ0i7pJph1/Ub4Z84azEy/BEQkbbl7UpfkcvdBcdSpJFjx\nAgAz+x/BUmRrqLkZy45s6fleRrBG8Jfh8mTd3T1q8xXds0UknTX1np2Ow0vOBqbU9aa762jkMXbs\n2JTHkE6Hvi99X8k8WpnNv0jMrFP15h9mNggo92D893JgrZkdZsG+1OcBT4bV/suWZcjOJHq77M1S\n/b2n06H/p/R96ftqPUc80qqnO7zxn0Awm15ERBLEzE4nWNt4W6DIzBa4+1CCzUmmm1kVQU/2eRHV\nriBYMrATwZKB08Pye4GHzGwRwZKBtXfBkxZSVDSHwsIZlJW1Jzu7ghEjBjNs2MCk101l29V1P/74\nJV55ZXSLt6vvK30+c7xxxy3V/1JI8L86XBpv7NixqQ4hrej7ahp9X00T3r9Sfh9tyUP37MZ5+unZ\nPnjwKN9552N98OBR/vTTsxtdLzd3pINvPnJzRzaqfnPqprLtmnXHpqhdfV/p85mbHne1eO7ZabM5\nTmOYmWfS50m24uJi8vLyUh1G2tD31TT6vprGzPAkj+lubdrSPTvenrWiojkUFDxHScn1QDGQR27u\nKMaPH8LQoQMpK4OyMigtDY7I5/n5o3ntteuirrnPPmP43e+upbISqqqgsnLLUf36rrtGs3BhdN3c\n3DGce+61EelK7OOxx0azeHF0/Z13HsMZZwT1q9V+Pm3aaL74IrruTjuN4bTTro1Zr9oTT4xmyZLq\nusH3VV331FOvja4Q4cknI+vWbDeZdVPZds26xbTU95XYuBNTd8iQMUyf3nDc1eK5Z6fV8BJJLCVE\nTaPvq2n0fYkEaibOgZKSUQCcdNJAvv0WVq2ClSuDI/L5f/4zgxUrquvlhXWv55RTxgAD6dgRcnKC\nIzu75vOPP479K/6rr7KYNQuysoKjXbstz6tfb9gQu255eRZmwTlm9R2x62dlZdGnT/DcItKVyOcd\nOsSu26FDFrvuWrPMaqU8HTtG1s2rUXfAgJiXraNuzXaTWTeVbafq+2pu/WTULS3Nqr9iAijpFhER\nSaK//31GjYQbgsT5Jz8Zg/tAunSB7baDnj2Dx+rnu+8OPXq0Z8WK6GsefXQWs2dHJ52RhgypYMaM\n6PIDD6xk0qT6Y/7ggwqWLo0u32OPSsaNq78uwIsvVvDZZ9Hl/ftX8rvf1V93+vQKSkqiy3fbrZKC\ngvrrFhXVXTc/v/66Tz+dmrqpbFufeYucnMr6KyZCU8ejtOYDjQ8UkTSFxnRnjG+/dX/qKfff/c79\n4IPd27UbG3MQxuGHj/WysvqvNXjwqJh1hwwZ3WAcsceuXtOMca+Nq5vKttOxbrrG3RY/c6R47tnq\n6RYREWmEusZlf/stzJ0LxcUwezYsXAiHHQZ5efD3v8Of/1zB889HX69790o6dqy/zREjBlNSMqpG\nT3lu7kjy809sMN7qMeMTJoyhtDSLnJxK8vNPbNRY8ubUTWXb6Vg3XeNui5+5uTSRUkSkFdBEytYt\n1rjs7t1H0aPHEFatGsgRR8CxxwaJ9iGHUCOZjlU3N3ck48c37hd9UdEcJkyYGZEgDGqZ5c1EpE7x\n3LOVdIuItAJKulu3IUNGM2NG9IoHhx02hjlzrm2wx1qJs0hm0eolIiIiSfD557F/XebkZDWYcEPw\nJ20l2SJtWzpuAy8iItIiysrg8sthyZKKmO+3yIoHIpIRlHSLiIjEsGQJDBwIK1bA/fcPJjd3VI33\ngwmNg1IUnYikGw0vERERqWXWLDj3XPjtb+Gqq8AsWE87FSseiEhm0ERKEZFWQBMpWwd3uPFGGD8e\nHn4Yjjsu1RGJSGukiZQiIiJxWrMGLrwQli+H+fOhb99URyQimURjukVEpM17771gfe0+fYINbpRw\ni0iiKekWEZE2bcoU+NGPYMwYuP12yM5OdUQikok0vERERNqkTZuCSZJFRfD887DffqmOSEQymZJu\nERFpE4qK5lBYOIOysvZABStXDiY3dyCvvQZbb53q6EQk0ynpFhGRjFdUNIeCgucoKbl+c1mPHqO4\n8UbYemst+yciyacx3SIikvEKC2fUSLgBvvnmem6/fWaKIhKRtkZJt4iIZLxgSEm00tKsFo5ERNoq\nJd0iIpLxsrMrYpbn5FS2cCQi0lYp6RYRkYw3YsRgevUaVaMsN3ck+fmDUhSRiLQ12gZeRKQV0Dbw\nyXfEEXNYt24m22yTRU5OJfn5gxg2TJMoRaTp4rlnK+kWEWkFlHQn19dfw267wRdfQLduLdKkiGSw\neO7ZGl4iIiIZ75FHYNgwJdwikjpKukVEJONNmgTnn5/qKESkLVPSLSIiGe3DD2HpUjjhhFRHIiJt\nmZJuERHJaA89BD//OWRpSW4RSSFNpBQRaQU0kTI5qqpg553h2Wdh772T2pSItCGaSCkiIhKhuBh6\n9lTCLSKpF3tfXBHZrGhmEYWTCynzMrItmxHnjGDYoGEZWzdd407XzyzJpQmUItJaKOmWtJCqhKpo\nZhEFtxdQckDJ5rKS24PnDdVPx7rpGne6fmZJrvXr4ckn4cYbUx2JiIjGdEsLSmTim7sgl/G/Hp/0\n+kMuGsKMfjOiyo/77Dgeuf0RKqoqqPRKKqsqox4v/b9LeWXAK1F1D/34UG7+681UeVXMw90ZM2YM\nC/ZcEFV33/f35Zo/XYO743jUY5VXcct1t/D+Pu9H1d3z3T0puKZg87lAzOe333A7H+33UVT9H779\nQy77w2WbzwVq1AWYePNEFu6/MKrugLcHcMnvLtn8uvb/p45z39/uY9H+i6Lq9n+rPxf99qIa59Z2\n/9/u55MDPolZ94L/uyCqvLYH/v5AzPq7vbVbg/Uf/PuDMesO+XwI0++b3mDb1TSmO/EeeggefRSe\nfjppTYhIGxXPPVs93dIkLd1jXFlVyd//9fca9QBKDihh9D2j+a7Xd5RWlLKxYmPwWL6x5uuKjUy/\nazpfHfZVVP2zbz6bfp/0o7yqnPLKcjZVbtr8vPpx42cboV90XLO/mM2ed+xJlmWR1S4r5uPnKz+H\nAdF1P1r9EaNfGE07axd1mBntrB2fr/085vfx5bovmfbRNAzDzGo8Vtf/uvTrmHW/Kf2G1798fXMd\nIObzteVrY9b/vuJ7Fn+3GGPLPSayLsD6yvUx626o2MDXG2rGFXkdgNKq0ph1S6tKWVtWM6bqdquV\neVmddTeUb4j5XmPqb6raRGlF7Lg2n+Ob6mxbUmvSJLj00lRHISISUNItjRZP4uzufL/pe/764F9j\nJs5XTLiC4zYcx9qytVHHmtI1bKzYCEuAXaOvveT7JRQtKqJT+07ktM+hU4fgsXOHzvTo1GPz63md\n5/EVX0XV779Nf+47/T46ZHWgQ7sOdMzquPl59eNpn5zG8zwfVfeEficw/ar6ezGHvD6EGUT3kh+x\nwxFMv6iBujNj1z1o+4N49MxH66/7ZOy6+223HxNPmVhvXYBF/17El3wZVb73tntz24m31Vv3vUfe\nYxnLosr32nYvbhp0U71135r8FktYElW+5zZ78tcT/lpv3Tf+9Uadda8//vp66wK8/q/XY9bfY5s9\nuO646+qt+9pDr/EFX0SV57TLabBdSZ6lS+HNN+HUU1MdiYhIQEm3NFrh5MKYifOou0fxWffPWLFu\nBV+t+4oV61cEx7rgsX279lR8VQH9o6/ZoX0HBu40kG7Z3WIeXTp2YehnQ2MmkQdvfzCTz5jcYNyP\ndn2Uj4geLtGzU0/2236/euv+5pzf8Nntn9UcmvJmLvlX5jfY7ohzRlBye0la1U3XuNP1M0vyPPww\nnHkm5OjfPiLSSqRV0m1mPwDuAfYCHLjY3eelNqr009ghIt+Xfc/C1Qs3HwtWLog51GLpuqV8uOpD\nem3Vi4P7HEyvrXrRq0svtt9qe3pt1YvOHTozZGHs3tfdfrAbFx1wUfRFI6Qyoar+XiZMmUBpVSk5\n7XLIvzK/UUNq0rFuusadrp+5NTGzm4GTgU1ACXCRu68J37sGuBioBEa4+4yw/CDgASAHeMbdC8Ly\nbGAScCCwGjjL3WOPmUoC92BoycSG/7AjItJiEj6R0sx6AgcDHYG57v5NAq/9IDDb3e8zs/ZAl+pf\nCuH7mkjZgFhDRHZ6fScuPONCug3oFiTY3yzk468/Zk3ZGvr36M+AbQYwYJsBPHnHk7y3z3tR12zM\nhLGYkxnfzGX8lY2fDFkjqTm7aUlNc+uLJFuqJ1Ka2SBglrtXmdkNAO5+tZntCUwGDgF2AJ4H+ru7\nm9l84Ep3n29mzwCF7j7dzK4A9nb3K8zsLOB0dx8eo82k3LPfeAN+9jP45BOwNjU1VURaSjz37CYn\n3WZ2OFAAvOXuN9Z67zzgDqAzYMAG4DJ3f7hJjcRutzuwwN1jjO7dfI6S7npUVlVyzM+P4ZXdo1fU\n2GbeNvy84OebE+zdt9mdHbrtQDvbsn9SqhNnkUyW6qQ7kpmdDpzh7j8Pe7mrqu/3ZjYdGAd8Drzg\n7nuE5cOBPHf/VXjOWHd/NewgWe7uPWO0k5R7dkEB9OgBY8cm/NIiIkDLrV7yc+AsYG6txncD7g2v\nWUHwZ8jOwP1m9ra7R3eRNs0uwCozux/YD3gDKHD3hpcmaKMqqyp5e8XbFC8upnhxMS998RKlq0ph\n9+hz9+7V8CS5RPwJX0m2SFq4GJgSPu8DRA7jW0rQ410ePq+2LCwnfFwC4O4VZrbGzHok8i+fdSkv\nhylTYJ4GHopIKxNP0n1M+PhUrfLLwuvNIRgXWA48CPwM+A3wizhjrNaeYHzgle7+mpndBlwN/Cny\npHHjxm1+npeXR15eXjObbX3qGpNdUVXBW1+9RfHiYmZ/PpuXvniJ3lv1Jq9fHufucy4TT5nIBR9d\nEHNsdWNXWlDiLJIYxcXFFBcXt2ibZjYT2D7GWyPd/anwnFHAJndveJZyAiT6nv3ss7D77rBrnX8T\nFRFpukTcs+MZXrIC2BrIjvy7oJl9SNCH+iN3nx2W9QM+BT5x9xgrFjep3e2BV9x9l/D10cDV7n5y\nxDkZP7wk1hCPbV7ehl0O3oWFXRfSt1tf8vrlcezOxzJw54H02qpXg/WbMkRERJKjNQwvMbMLgUuB\n4929NCy7GsDdbwhfTwfGEgwveTFieMnZwEB3v7x6CIq7z2vp4SVnnglDhmh9bhFJrpYaXtIDWFsr\n4e5BkHCvIWLYibsvNrONwI5xtFODu39lZkvMbIC7LwROAKK33ctwsZbtW33kavp+0JdP7vuEnl2i\nfq/VkCkrLYhIYpnZicBVwLHVCXfov8BkM/s7wbCR/sD8cCLlWjM7DJgPnAcURtS5gGBYypnArJb4\nDN98A88/D/fc0xKtiYg0TTxJ93qgu5l1dN+8Fdux4eM8d6+qdf6mONuJJR942Mw6Ei5plaDrtnqb\nKjcx9cOpzPtyXsxl+7p36t5gwl1NQ0REJIYJBKtOzQx3/HzF3a9w9w/M7DHgA4L5OldEdLpcQbBk\nYCeCJQOrlzG6F3jIzBYRLBkYtXJJMjz2GJx4IvzgBy3RmohI08STDH8AHE7Qe1E95u/C8LE48kQz\n6wp0Bz6JL7ya3P1tgmWr2ozF3y1m4hsTuW/Bfey13V7s0n0X3ubtqPO0+52INIe7x9i+avN7fwH+\nEqP8DWCfGOVlBPN5WtSkSTB6dEu3KiLSOO0aPiVK9R7Ud5nZ7WY2DTiFoAek9v7UR4SPi+KMr02q\nrKqkaGERJ08+mYMnHszG8o3MvnA2s86fxfWXXk/ugtwa5+e+mUv+2dr9TkTarkWL4NNPYfDgVEci\nIhJbPD3ddwKnAwOByyPK/+zui2udW/0nxRfjaCejxVqB5JAjD+G+Bfdx1xt30bNzTy4/+HIe++lj\ndO7QeXM9jckWEYn20ENwzjnQPq32WRaRtiSuHSnD2ejnEAwzWQM86+5zap3TAXgIyCZYZeTj5ofb\nYFxpsXpJrBVEtpq7FVW5VZw97GwuP/hyDupzUAojFJGW1hpWL2lpibpnV1VBbi5Mmwb775+AwERE\nGtAiO1K2ZumSdA+5aAgz+kWvlX3cZ8cx64EWmeQvIq2Mku74zZ4N+fnw9tva9l1EWkY89+x4xnRL\nM62vXB+zvJLKFo5ERCT9TZoE55+vhFtEWrdmjX4LN6zJA/oCndz9z4kIKpPNLJnJG0vfgNzo97QC\niYhI02zYAFOnwgcfpDoSEZH6xdXTbWadzOyfwBcEywbeSLBDWeQ5W5vZt2ZWaWa7NT/U9PZ92fdc\n9tRl/OKpXzDywpFagUREJAGefBIOPxx69051JCIi9WtyT3c4ibKIoId7A8EOlEcQTJjczN2/NbP7\ngP8jWMXkuuYGm65mfTqLS/57CYN2HcQ7v3qH7jndObD3gVqBRESkmaqHloiItHZNnkhpZpcRLBu4\nCBjq7p+a2VdAT3fPqnXu4cDLwIvufnyCYq4vtlY1kXLdpnX8YeYfeGrhU9x9yt2cuNuJqQ5JRFop\nTaRsuuXLYc89Ydky6Ny54fNFRBKlpSZSnhc+5rv7pw2c+yZQBewZRztp7cXPXmTfO/dlY8VG3r38\nXSXcIiIJNnky/OQnSrhFJD3EM5Fyb4LdJxvc8MbdN5nZGqBHHO2kpfWb1nP181cz7aNp3HXyXQwb\noCEjIiLJMGkSFBamOgoRkcaJJ+nOATa6e3kjz+8ElMbRTqtXe1fJ4044jomrJ3L0Tkfz7uXvsnWn\nrVMdoohIRnrrLVizBo45JtWRiIg0TjxJ93JgJzPr4e7f1Heime1HkKR/Ek9wrVmsXSVnTZzFyAtH\n8ufTtHKiiEgyFBXNobBwBu+/357s7AqefXYww4YNTHVYIiINiifpfhG4ELgI+FsD544LH5+Po51W\nrXByYY2EG6DyuErmz50ffDMiIpJQRUVzKCh4jpKS6zeXFRSMAlDiLSKtXjwTKf8OODDGzAbFOsHM\n+pjZw8CPgU3A+PhDbJ3KvCxmeWlVRo6kERFJucLCGTUSboCSkuuZMGFmiiISEWm8Jifd7v4eUAB0\nA6ab2dtAd8DMbKqZvQEsBs4mSM5/5e6fJy7k1qGioiJmuXaVFBFJjrKy2H+cLS3NilkuItKaxLUj\npbv/A/gJsATYhy0b45wGHEAwbGUJcJq7P9D8MFuXL9Z8wcIfLKTnKz1rlGtXSRGR5MnOrqOzI6ey\nhSMREWm6eMZ0A+DuT5jZUwQ7Ux4J9CZI4r8CXgFmuXvsO2QaW7FuBSdMOoFrzr2GAesGaFdJEZEW\nMmLEYEpKRtUYYpKbO5L8fO2DICKtX5N3pGzNkr0j5TcbvyHvgTzO2OMMxuaNTVo7ItL2aEfKxikq\nmsPZZ89kwIAstt22kvz8QZpEKSItLp57tpLuRvq+7HsGPTSIo/oexS2Db8GsTf1uFJEkU9LdeNtu\nCx9+CD17NnyuiEgyKOlOUtK9sXwjwyYPY7ceu3HXyXcp4RaRhFPS3TgbN8LWWwePuhWLSKrEc89u\n8phuM6siWJWkSdw9LaeXl1eW89N//5Ttt9qeO4fdqYRbRCSFli2DHXZQwi0i6SfeiZRt4nZXWVXJ\nedPOo52148HTHiSrXVr+u0FEJGMsXQo77pjqKEREmi6epHvXBt7vDhwC/IZgRZOLgLfjaCelqryK\nXz71S1ZtWEXROUV0yOqQ6pBERNq8ZcuUdItIempy0u3uixtx2ttm9hDwDHAvcFBT20kld+e3z/2W\nD77+gJnnzSSnvTa8ERFpDZYuDYaXiIikm7g2x2kMdy8DRgDbAGm1vt644nEULy7mmXOeYauOW6U6\nHBERCWnyOA0FAAAgAElEQVR4iYikq6Ql3QDu/j7wPTAkme0k0i0v38Kj7z/KjPNmsHWnrVMdjoiI\nRFDSLSLpKu4dKRvDzDoCnYGOyWynOYpmFlE4uZAyL2PV96tY1XsVb1z/Btt12S7VoYmISC1KukUk\nXSU16QbOAbKAL5LcTlyKZhZRcHsBJQeUbC7b6fWdeGf+O/Qd1DeFkYmISCxKukUkXTV5cxwz26mB\nU3KAvsCPgUuBbOBmd/9jXBE2LbYmbbQw5KIhzOg3I7r88yFMv296IkMTEamXNsdpWHk5bLUVbNgA\nWVrBVURSqEU2xwEW0/DmOJFBvApcG0c7SVfmZTHLS6tKWzgSERFpyPLlsN12SrhFJD0lY3OcKuBb\n4B3gMeAed6+Is52kyrbsmOU57bREoIhIa6OhJSKSzpq8eom7t2vgaO/u27r7ce7+z9aacAP86Pgf\nkfVizS6T3DdzyT87P0URiYhIXZR0i0g6S/ZEylarrKKMB9c8yNUXXM3rc1+ntKqUnHY55F+Zz7BB\nw1IdnoiI1KKkW0TSWZtNum/6300M2GYA1w2/LtioXkREWjUl3SKSzpK6OU6imdliM3vHzBaY2fx4\nr7No9SLGvzqeCUMnJDI8EZG0ZWY3m9mHZva2mU01s+5heQ8ze9HMvjezCbXqHGRm75rZIjMbH1Ge\nbWaPhuXzzGznRMS4bJm2gBeR9FVvT7eZXUDDK5U0irtPSsRlgDx3/6YZcXDFM1dwzdHXsFP3hlY/\nFBFpM2YAf3T3KjO7AbgGuBooBUYDe4dHpDuBS9x9vpk9Y2Ynuvt04BJgtbv3N7OzgBuB4c0NUD3d\nIpLOGhpecn+C2nEgEUk31L9ySoOmvDeFVetXUXB4QYLCERFJf+4+M+Llq8AZYfkG4H9m1j/yfDPr\nDXR19+q/Ok4CTgOmA6cCY8Pyx4F/JCJGJd0iks4aSrrnJKidhPSWh9d53swqgbvc/e6mVP5247f8\nfsbvmXbWNNq3a7PD2UVEGnIxMKVWWe37+A7A0ojXy8Ky6veWALh7hZmtMbMezfkrZWVlsE53nz7x\nXkFEJLXqzTzdPa+F4miso9x9uZn1BGaa2UfuPrexla+ZdQ2n//B0DtvxsCSGKCLSOpnZTGD7GG+N\ndPenwnNGAZvcfXKLBteAlSuhRw/o2DHVkYiIxCetunvdfXn4uMrMpgGHAjWS7nHjxm1+npeXR15e\nHgAvL3mZpxY+xQdXfNBS4YqI1Km4uJji4uIWbdPdB9X3vpldCJwEHN+Iyy0DIgd77MiWnu9lwE7A\nl2bWHuheVy93XffsqMaWaWiJiKROIu7Z5p6okR/JZWadgSx3/97MuhBM+vl/7j4j4hyP9XnKK8s5\ncOKBjD5mNGftfVbLBS0i0khmhrs3a85KM9s/EfgbcKy7fx3j/QuBg9w9P6LsVWAEMB8oAgrdfbqZ\nXQHs4+6Xm9lw4DR3j5pIWdc9O5YnnoD774cnn4zjw4mIJFg89+x06unuBUwzMwjifjgy4a7PrfNu\nZYeuO/CzvX6WzPhERNLZBKAjwdA9gFfc/QoIlmsFugIdzew0YJC7fwRcATwAdAKeCVcuAbgXeMjM\nFgGr0colIiLNS7rNbCBwJNAH6EI9K4u4+8XNacvdPwP2b2q9xd8t5qb/3cT8S+cT/iIREZFa3L1/\nPe/1q6P8DWCfGOVlQEJ7OZR0i0i6iyvpNrN9gMnAXo2s4gSz4VuUu3PlM1fyuyN+x65b79rSzYuI\nSIIsXQp7NfY3johIK9TkpDtcm/V5oCfwQfh8BLAeuJVgZvxxwK4Ef1b8J1CRoHibZOqHU/nsu8+Y\netbUVDQvIiIJop5uEUl38fR0/54g4X4O+LG7bzKzEcD37v4nAAvGcVxKsCHCQcDJCYq30daWraVg\negGPnPkIHbO0xpSISDrTFvAiku7axVFnaPg4yt03xTrBAxOBUcCJwK/jjC9uo18YzYm7ncjROx3d\n0k2LiEgCuQc93Uq6RSSdNXnJQDNbD3QAsqvXejKzKuBbd9+m1rndgG+BN9z90MSEXG9s7u68/uXr\nnDz5ZN6/4n226bxNwxVFRFIs1UsGpkJjlwxcvRr694dv4t7PUkQkseK5Z8fT010FrKl1p1wHdDOz\nrMgT3X0tsBYYEEc7camoquCypy/jpkE3KeEWEckAGs8tIpkgnqR7GUGCHVl3MZAF7Bt5opl1B7oT\nrP3aIm6ffzvds7tz3r7ntVSTIiKSREq6RSQTxJN0f0wwvGSPiLI54eNVtc69Nnz8MI524nLd3Ou4\nc9idWpNbRCRDKOkWkUwQT9L9XPgYuSLJPwiWBRxuZu+a2cNm9g5wZfj+nc2IsUm6/a8bnyz4pKWa\nExGRJNPKJSKSCeJZMvAxgjW4N1QXuPtHZnY+cDfBhjnVWxg4cJu739PcQBvr0wM+peD2AgCGDRrW\nUs2KiEiSLF0KRx2V6ihERJqnyUm3u39NsFZ37fJHzGwWwZKCOwLfAc+7+8JmR9lEJQeUMGHKBCXd\nIiIZQMNLRCQTxLUNfF3cfRUwKZHXjFdpVWmqQxARkQRQ0i0imaDJY7rN7BQzS2iyngw57XJSHYKI\niCSAkm4RyQTxTKR8EvjKzO4ys2MTHVAi5L6ZS/7Z+akOQ0REmmnt2mBHym7dUh2JiEjzxJN0rwV6\nAJcCL5jZF2Z2s5kdkNjQ4jPk8yGMv3K8xnOLiGSA6pVLtAqsiKS7eLaBzwFOAs4GhgHV4zicYA3v\nKcBkdy9JYJyNja1RWwqLiLQ22gY+tpkz4cYb4fnnWygoEZFGaJFt4N291N2nuvtPgV7ABQRrd1cB\nPwT+H7DQzF41swIz69XUNkREREDjuUUkc8QzvGQzd//e3R9y96FAb+DXwEuAAYcAtwJLzWxmsyMV\nEZE2R0m3iGSKZiXdkdz9a3e/090HAjsDfyTY/j0LOC5R7YiISNuhpFtEMkXCku5qZtYBOJigp3uX\n6uJEtyMiIplPW8CLSKZIyHrbZtaOoDf7HOB0oHvE218AjySiHRERaVvU0y0imaJZSbeZHU6wisnP\nCCZVVlsN/JtgFZOXmtOGiIi0XUq6RSRTNDnpNrO9CHq0hwP92DJ0ZAPwBDAZmOHuFQmKUURE2qCN\nG2HdOth221RHIiLSfPH0dL8b8bwCmE6QaD/p7hsTEpWIiLR52hhHRDJJvMNL5hAk2v92928TGI+I\niAigoSUiklniSbp3cvelCY9EREQkglYuEZFMEs+OlEq4RUQk6dTTLSKZJOHrdIuIiCSCkm4RySRK\nukVEpFVS0i0imURJt4iItEpKukUkkyjpFhGRVmnpUk2kFJHMYe6e6hgSxsw8kz6PiLQdZoa7t6kV\nqeu7Z5eXQ5cuwQY5WVktHJiISAPiuWerp1tERFqd5cuhVy8l3CKSOZR0i4hIq6Px3CKSaZR0i4hI\nq6OkW0QyjZJuERFpdZR0i0imSauk28yyzGyBmT2V6lhERDKJmd1sZh+a2dtmNtXMuoflg8zsdTN7\nJ3z8UUSdg8zsXTNbZGbjI8qzzezRsHyeme3c1Hi0BbyIZJqkJ91mtruZzTSzDWa20sz+bWYD47xc\nAfABoCVKREQSawawl7vvBywErgnLVwEnu/u+wAXAQxF17gQucff+QH8zOzEsvwRYHZbfCtzY1GDU\n0y0imaYlerpvACYT3KwLga7AdDO708wavdSKme0InATcA7SpZbVERJLN3We6e1X48lVgx7D8LXf/\nKiz/AOhkZh3MrDfQ1d3nh+9NAk4Ln58KPBg+fxw4vqnxKOkWkUzTPpEXM7NOQGmthVfnuvv9tc7r\nAfwNGAlc38jL3wpcBXRLRKwiIlKni4EpMcrPAN5w93Iz2wFYGvHeMqB6QMgOwBIAd68wszVm1sPd\nv2lsAEq6RSTTJLqn+xXgWzMrMrM/mNlhQJmZ7Rp5krt/4+4XAds15qJmdjKw0t0XoF5uEZG4hEP9\n3o1xnBJxzihgk7tPrlV3L4K/XF6W7DgrK4N1uvv0SXZLIiItJ6E93QTDR0YAcwn+nDg6bONPZlYI\n/A94zd3Xm1lHIKeR1z0SONXMTgrrdDOzSe5+fu0Tx0WMWMkLDxGR1qY4PFqSuw+q730zu5BgGN/x\ntcp3BKYC57n7Z2HxMsIhKKEd2dLzvQzYCfjSzNoD3evq5R43btzm53l5eeTl5bFyJfToAR07NvKD\niYgkWXFxMcXFxc26RkK3gQ8T6dPc/bHwdXvgYOBYgvz3SGArYAVQBRS4++NNbONY4PfufkqM97QN\nvIikpVRvAx9OgvwbcKy7fx1R/gNgNjDW3Z+oVedVgo6W+UARUOju083sCmAfd7/czIYT/F4YHqPN\nmPfs11+Hyy6DN95I4AcUEUmgeO7ZCe3pdvdNwGMRryuAeeFxo5llAf0JEu8P3X19vE01N1YREalh\nAtARmBnOcX/F3a8ArgRygbFmNjY8d1CYmF8BPAB0Ap5x9+nh+/cCD5nZImA1EJVw10fjuUUkEyW0\np3vzRc2ygYnATGBaM5Lrprarnm4RSUup7ulOhbru2f/4B3z4Idx+ewqCEhFphHju2c2aSGlmPzOz\n081sm8hydy8DfgNsAO41s0ub046IiLQd6ukWkUzU3NVLbiZYg3Wlmb0frr093Mz6uPu37j41HMd3\nevNDFRGRtkBJt4hkouaO6T6EYDOEo4FjCJaSugzAzD4F3gC+A7ap6wIiIiKRli7VFvAiknmalXS7\n+0qCsdsTAcLNEo5hSxJ+JlACRC3tJyIiEsuyZerpFpHMk5SJlJsvbnYcwcz3c9y9NGkNbWlPEylF\nJC1pImXAHTp3hq+/hi5dUhSYiEgDWnwiZUPc/QXgr8C4ZLYjIiKZ4ZtvoFMnJdwiknmSmnQDuPtr\nwK4NnigiIm2eJlGKSKZq1phuM5sLbGLLjsbz3L08xqnZzWlHpKhoDoWFMygra092dgUjRgxm2LCB\nLVK/rdVN17jT9TNLTUq6RSRjuXvcBzCeYC3uqvBYD8wC/gT8GBgM3ArMaE47TYjHJbmefnq2Dx48\nyo89dqwPHjzKn356dtLrPv30bM/NHenBaM/gyM0d2SL121rddI07XT9zpPD+lfT7ZGs6Yt2z77rL\n/Re/aNJXJyLS4uK5Zzd7IqWZdQQOBY4FBgJHArVH4/2bYIWTl919Y7MarD8Wb+7naQvi7ZUrKppD\nQcFzlJRcv7ksN3cU48cPabB+Y+q6Q3l59DF8+Gjmzr0u6pqHHz6GW265lspK6j3+/OfRvPNOdP29\n9hpDQcG1uENVVezjrrtGs3BhdN3ddhvD+edfuznNqqqKTLmCY8qU0Xz2WXTdfv3GcMYZW+pWf/bI\n59OmjWbJkui6ffuO4ZRTrq1xbrXq50VFo1m6NLrujjuO4cQTr61RFqv+c8+NZtmy6Pp9+oxhyJBr\no86PNGPGaL78MnbdQYOuja4QYebMuuuecELDdZcvj67bu3fDdQGefz7++nXVHTJkDNOnN9x2NU2k\nDPzpT5CVBWPH1lFJRKQViOee3dx1unH3TcBL4XG9mbUHDmJLEn408NPwKDez14HZwBR3f7e57bdF\nzf0zeO3kt6RkFEDUNcrKYO3aLce4cTNq1AvqXk9BwRjmzh1IaSls3AilpUQ9f/vtGaxdG133xz8e\nQ4cOAykvDxLkDh2ij9WrY/9n+t57Wfz+99C+ffBLuq7jyy9j1//66yxefRXatYs+zILH9etj1y0r\ny6KsbMt5kY/Vh3vsumZZ9OpV89ygfMvzjh1j183OzmLPPbecH/lY/fyFF2LX7dQpi0MOqXl+rPov\nvRS7/lZbZXH00XXXBXj55dh1u3bNIi8v5lubzZtXd93jjqu/7quvtmf58ujybt2yOOGE+usCzJ8f\nf/266paWZjXcsERZuhSOOirVUYiIJF6zk+7a3L0CeDU8bjKzdsD+bEnCjwGOAC4E+iS6/UzXlKS5\n2rp1sHw5fPUVjBwZO3E+77wx7LzzwBpJtjt07w7dugXHp5/G/s+lvDyLH/wAcnKCo1Onmo85OfCb\n37RnwYLouocfnsXMmUFynZUVncQBDBlSwYwZ0eVHHVXJ9Ol1fFGNqL///pXcc0/9dT/4oIJly6LL\n99yzkuuiOzdrmDu3gsWLo8sHDKjkqqvqr/vssxWUlESX5+ZW8utf11932rQKFi2KLt9110p++cv6\n6wI89lgFCxdGl++ySyUXX1x/3UceiV23X79KLryw/rpTplTw8cex615wQf11J0+uu+75jdgl4OGH\n469fV92cnMqGG5YoGtMtIhmrqeNRmnsABuwDHJ2Ea8c/OCdNDB48yqMHMrjvvfdo/8tf3PPz3X/6\nU/ejj3bfbTf3Ll3cO3Vy33VX96OOct9227Ex6++331h/4w33RYvcV650Ly1tfNtDhoyOO+7G1I09\nZvaaZo7XbVz9tlY3XeNO188cCY3pdnf3PfZwf++9Jn11IiItLp57dsJ7uhuR5DvQ5oeVNGWIyMaN\n8P778M478P77sX9kK1dmsWYN7Lpr8KfZ7beH3r2Dx65dt/Qg19Xru/32lRx4YP0xjxgxmJKSUbXG\nZY8kP//EBj9vc+pWfy8TJoyhtDSLnJxK8vNPbPSQmubUb2t10zXudP3MEk1bwItIpkrqjpQtLV0m\nUtY1qfC224aw114DeffdIMGuPj7/HAYMgH33hfnzY0/sa+ykrdhtj2T8+MYlCUVFc5gwYWZEcjGo\nSePJ460rkuk0kTIY1tanD3z/feyhZiIirUU89+xkbwN/NHAjsA3wCtAB+BD4h7uvSUJ7aZF0Dxky\nmhkzohPnrKwxbL/9tey7LzWOAQOgY8fgnOYmzdXXUPIr0roo6YYPP4TTT4ePPkphUCIijZCS1Usa\ncCPBpMnv3f0iADPbH7jLzCa7+3+T3H6rU14Oy5bF/toPOyyL//2v/vqJ+FP2sGEDlWSLSKujSZQi\nksmSnXRPAnYGHqoucPe3zOxs4Lwkt92qfPIJ3HsvPPAAbNxYEfOcrl0bt9qBkmYRyURKukUkk7VL\n5sXd/S5339Hdr6lV7u4+KZlttwalpTBlChx3HBx5ZNDL/cIL8PDDg8nNHVXj3GBS4aAURSoiknpK\nukUkk7X46iWZpK4VSN57D+65Bx5+GA44AC6/HE49FbKzg3p77KHVDkREalu6FPbfP9VRiIgkh5Lu\nOMWa0LhgwSh69IB16wZy0UUwfz7sskvs+hoiIiJS07JlcPLJqY5CRCQ5lHTHqbAwemfHVauup2/f\nMbz33kDa65sVEWkSDS8RkUyW1DHdmaysLHZW3bVrlhJuEZE4KOkWkUympDtO5eWxVyDJyWncCiQi\nIrLFxo2wbh1su22qIxERSY6UJN1mdpKZ9UhF24nwzjvwwQeD6dVLK5CIiCTCsmXB9u/aiVJEMlXS\nB0KY2aXAYmCeu38fFi8EzjGz9u5+W7JjSKQFC2DoUJg4cSCdO2sFEhGRRKhOukVEMlVLjD7eC/gn\nUGVmHwL/A14GXgcub4H2E+b112HYMLjzTvjJTwC0AomISCJoPLeIZLqWGF6ynmBXyt7An4ANwAjg\neeCNFmg/IV59NUi47767OuEWEZFEUdItIpmuJXq6N7j70vD5E+GBmV0CLG+B9pvt5ZfhtNOCLdxP\nOinV0YiIZJ6lS6F//1RHISKSPC3R093PzKJmF7r7vcBhLdB+s8ydGyTcDz2khFtEJFnU0y0ima4l\nerrHAsVm9j7wODDb3ZeYWTbQqwXaj9uLL8JZZ8GUKXD88amORkQkcy1dqomUIpLZkt7T7e5fAkcC\n3wL3Ap+b2Xrga+DVZLcfr+efDxLuxx5Twi0ikmzLlqmnW0Qym7l7yzVmtg1wOJADvBox1jtR1/dE\nfJ7p0+H882HqVDj66AQEJiLSADPD3dvUKtXV9+zycujSJdggJysr1VGJiDQsnnt2UoaXhENHJgIz\ngWnuvh7A3VcDRcloM1GefhouvhiefBKOOCLV0YiIZL7ly6FXLyXcIpLZmpV0m9nPgHJgTphQA+Du\nZWb2G+BHwL1mNsvd725eqMlRVDSHwsIZlJW1Z82aCj77bDAzZgzk0ENTHZmISNugSZQi0hY0t6f7\nZqAv4Gb2ETAHmE2QhH8JTAWmmtkzQKtLuouK5lBQ8BwlJddvLttxx1GsWgWgTW9ERFqCkm4RaQua\nO5HyEOBXwMNAZ+AyYDKw1Mw+MbNHzewuYJtmtoOZ5ZjZq2b2lpl9YGZ/be41Cwtn1Ei4AZYuvZ4J\nE2Y299IiImnFzG42sw/N7G0zm2pm3cPyQ81sQXi8Y2ZnRdQ5yMzeNbNFZjY+ojw7vP8vMrN5ZrZz\nfW1r5RIRaQualXS7+0p3n+ju57v7LgS93ucAdxDsRHkmwRCTguYG6u6lwI/cfX9gX+BHZtasaY5l\nZbE7+ktLNbBQpKnMTEcjj1ZqBrCXu+8HLASuCcvfBQ5y9wOAwcDtZlZ9k7wTuMTd+wP9zezEsPwS\nYHVYfitwY30Na+USEWkLErpkoLsvc/dH3P3K8MY9CHgPeCtB198QPu0IZAHfNOd62dkVMctzciqb\nc1mRNsvddTRwtFbuPtPdq8KXrwI7huUbI8o7AWvcvdLMegNd3X1++N4k4LTw+anAg+Hzx4F6F17V\n8BIRaQuSuk63u78A/BUYl4jrmVk7M3sLWAG86O4fNOd6I0YMJidnVI2y3NyR5OdHbaApItKWXAw8\nU/0iHGLyPvA+8NuweAcgctnXZWFZ9XtLANy9AlhjZj3qakxJt4i0BUnfkdLdXzOzqxJ0rSpg/3Cs\n4XNmlufuxZHnjBs3bvPzvLw88vLy6rxe794D2WorGDhwDGVlWeTkVJKffyLDhmkSpYgkV3FxMcXF\nxS3appnNBLaP8dZId38qPGcUsMndJ1e/GfZm72VmPwSmm1lxomIaN24c770H//kPVFTUf88WEUmV\nRNyzm7U5jpnNBTYBxeExz93LY5z3pLv/OO6GYrc9Btjo7rdElDVpc5yLLoIf/hD++MdERibSNlmw\nUUCqw2j16vqerBVsjmNmFwKXAseH82hinTML+ANBz/aL7r5HWH42MNDdLzez6cA4d59nZu2B5e7e\nM8a1vLLSycmBdeugY8ckfTARkQSL557d3J7uNwlu0D8KX280s3kEywa+DWwEhhKMA2wWM9sWqHD3\n78ysE8F48f8X7/VWroQnnoBPPmluZCIi6S+cBHkVcGxkwm1m/YCl7l4RrkLSH1jk7mvNbK2ZHQbM\nB84DCsNq/wUuAOYRTKifVVe7K1fC1lsr4RaRzNespNvdC8KhI4cCxxIsbn0kW5Lwav82s+OBl919\nY5zN9QYeNLN2BGPRH3L3Om/kDbn7bjjzTNim2YsZiohkhAkEk9RnhiusvOLuVwBHA1ebWTnBZmi/\ndPe1YZ0rgAcIOlaecffpYfm9wENmtghYDQyvq1GN5xaRtqJZw0tiXjD4U+JBbEnCjwa6hW+XA68T\n9IRPcfd3E9x2o4aXlJdDv37w7LOw776JjECk7dLwksZpzcNLWpqZ+bRpzv33w5NPpjoaEZHGS8Xw\nkijhTPVXw+OmsGd6f7Yk4ccARwAXAn0S3X5jPP44DBighFukJRQVzaGwcAZlZe3Jzq5gxIjBTZ6s\nnIhrSOuknm4RaStaYvWSKoKx328Ct1rwd8u9ge7JbrsuhYVwVULWUxGR+hQVzaGg4LkaO7+WlATL\ndDY2aU7ENaT1UtItIm1FUtfpjsUD77r7Sy3dNsBrr8GXX8Ipp6SidZG2pbBwRo1kGaCk5HomTJjZ\nYtfo168ft9xyC/vuuy9du3blkksuYcWKFQwdOpTu3bszaNAgvvvuO4qLi+nbt29U3Vmz6p86Mn/+\nfI444gi23npr+vTpQ35+PuXlUYs4SR20BbyItBUtnnSn2oQJ8OtfQ/uk9/GLSFlZ7P/RnnsuCzMa\ndcyYEfsapaVZMctrMzOmTp3KrFmz+Pjjj3n66acZOnQoN9xwAytXrqSqqorCwsKY27M3Ztv29u3b\nM378eFavXs0rr7zCrFmzuOOOOxoVm2gLeBFpO9pU0v3VV/DUU3DJJamORKRtyM6uiFk+ZEgl7jTq\nGDw49jVyciobHUd+fj49e/akT58+HHPMMRxxxBHst99+ZGdnc/rpp7NgwYK4Ph/AgQceyKGHHkq7\ndu3Yeeed+eUvf8ns2bPjvl5bo+ElItJWtKmke+JEOOss6FHnZsQikkgjRgwmN3dUjbLc3JHk5w9q\n0Wv06tVr8/NOnTrVeJ2Tk8O6desafa3aFi5cyMknn0zv3r3p3r07o0aNYvXq1XFfr63R8BIRaSva\nzCCLTZvgn/+EmY0fSioizVQ90XHChDGUlmaRk1NJfv6JTZoAmYhr1BZryb4uXbqwYcOGza8rKytZ\ntWpVg9e6/PLLOeigg3j00Ufp0qULt912G48//njcsbU1nTpBly6pjkJEJPnaTNL9n//AnnvCXnul\nOhKRtmXYsIHNXmUkEddoyIABAygtLeWZZ55h0KBB/OUvf6GsrKzBeuvWraNr16507tyZjz76iDvv\nvJPtttsuqbFmEg0tEZG2os0MLykshBEjUh2FiLQGkZMjqydLduvWjTvuuINf/OIX7Ljjjmy11VZR\nq5nEcssttzB58mS6devGL3/5S4YPH97g5EvZQkNLRKStSPiOlKlU146Ur74KZ58NixZBVuMWPBCR\nJtKOlI2jHSm3MDP/xS+cu+9OdSQiIk0Tzz27TfR0T5gAV16phFtEpLXR8BIRaSsyPulevhyKiuDi\ni1MdiYikq6FDh9K1a9eo44Ybbkh1aGlPSbeItBUZP5Hyn/8Mhpb84AepjkRE0tWzzz6b6hAylpJu\nEWkrMjrpLiuDu+6CF15IdSQiIhKLJlKKSFuR0cNLHnsM9t03WCpQRERaH/V0i0hbkbFJtzuMH69l\nAkVEWrPu3VMdgYhIy8jYpHvePPjuOzjppFRHIiIiddGS5iLSVmRs0l1YGCwT2C5jP6GIiIiIpIuM\nTKAT7ycAAA4ISURBVEmXLYPnnoOLLkp1JCKS6caNG8d5552X6jDS1pAhoykqmpPqMEREki4jVy/5\n5z/h3HM1VlCkNSiaWUTh5ELKvIxsy2bEOSMYNmhYi18jWbTle/PMmHEdJSWjABg2bGCKoxERSZ6M\nS7pLS2HiRJijjhORlCuaWUTB7QWUHFCyuazk9uB5Y5PmRFwjmWJt6S5NU1JyPRMmjFHSLSIZLeOG\nlzz6KBx4IOy+e6ojEZHCyYU1kmWAkgNKmDBlQotdo1+/ftxyyy3su+++dO3alUsuuYQVK1YwdOhQ\nunfvzqBBg/juu+8oLi6mb9++UXVnzZpV7/XNjNLSUoYPH063bt046KCDeOeddxr9+SRQWpqV6hBE\nRJIq43q6x4+H669PdRQiAlDmZTHLn/v0Oez/NXJYxmdAv+ji0qr/3979B1lV3nccf3+E1DWAMQjd\nSFaDJRJs/iCkjrUjCGmGLGBNF9MktTWGCdpWm9XE/oixNTD9QUw7JogR86OEmAxGmwrRRGWgBMfM\npMRKSLQaA6k6gZQiTZsAIavIfvvHeZa9XM7d3bu7l3vPvZ/XzBnufZ7nnPu9h3O/+8w5zzlPz5BW\nl8T69evZsmULR44cYdasWezYsYO1a9cyY8YMFi1axKpVq5g7d27uuoMNH4kIHnjgAe69917WrVvH\nypUr6erqYufOnYwd23Qptmba2o7WOwQzs5pqur8Ihw5BZ2e9ozAzgFN1am555691snHZxiFto/OF\nTjax6YTytlPahhxHd3c3kydPBmDOnDm0t7czc+ZMABYvXsyWLVtyO91DdcEFF3D55ZcDcOONN3Lb\nbbexbds2Zs+ePexttpJp026mu3tBvcMwM6upphteMmbMX/PIIx7QbdYIrv+D65m2Y9pxZdO+O43u\nK7pP6jba29uPvT7ttNOOe9/W1sahQ4eGvK08HSXTKkqio6ODvXv3jmibraKz8xZuv32Bx3ObWdNr\nujPdzz77d9xwg++EN2sEfTc63vGVO+jp7aHtlDa6P9hd1Q2Qo7GNcnk3P44bN47Dhw8fe3/06FH2\n798/pO3t3r372Ove3l727NnDlClThh1fK9m48W/rHYKZ2UnRdJ1u8J3wZo3k0vmXjvgpI6OxjcFM\nnz6dnp4eHn74YebPn8+KFSt46aX8Menltm/fzoYNG7jssstYtWoVbW1tXHTRRTWN18zMiqXphpf0\n8Z3wZlZJ6c2RfTdLnn766axevZqrr76ajo4Oxo8ff8LTTCptq6uri/vuu4+JEyeybt061q9fz5gx\nzkFmZtZPzfSMWUkB2ffp7LzFly3NTiJJfmb1EFTaT6m8pWbakRQ+ZsysiIaTs5vyTHd2J/z8eodh\nZmZmZgY0Yafbd8Kb2WhbuHAhEyZMOGG59dZb6x2amZkVRNMNL2mm72NWJB5eMjQeXtLPOdvMisrD\nS8zMzMzMGpA73WZmZmZmNdaUz+k2s/oofRSfmZmZ9SvMmW5JZ0vaKulpSf8h6fp6x2Rm/SLCyxCX\nRiTpHyX9QNL3Ja2X9Jqy+nMkHZL0ZyVlvyHpKUm7JN1eUn6qpPtS+TZJbziZ38XMrBEVptMNHAE+\nHBFvBi4C/lTS+XWOqdAeffTReodQKN5f1fH+KpxNwJsjYiawE/hoWf0ngYfKyu4ClkbEecB5khak\n8qXAT1P5p4BP1C7s1uHfVHW8v6rj/VV7hel0R8R/R8T30utDwA+AKfWNqtj8A6uO91d1vL+KJSI2\nR0RvevsdoKOvTlIX8BzwTEnZWcCEiHg8FX0J6Eqv3wncnV7fD7y9hqG3DP+mquP9VR3vr9orTKe7\nlKSpwCyyPwxmZja6PgA8DCBpPPCXwPKyNq8H9pS8/0kq66vbDRARrwA/lzSxhvGamTW8wt1Imf4A\n/AtwQzrjbWZmQyBpM/C6nKqbI+Lrqc1fAS9HxD2pbjnwqYg4LN8pa2Y2bIWaHEfSq4BvAI9ExMqc\n+uJ8GTOzMvWeHEfSEuAa4O0R0ZPKHgPOTk3OAHqBW4D1wNaIOD+1uwK4JCKulbQRWB4R2ySNBfZG\nxOScz3PONrPCqjZnF+ZMdzrDsgZ4Jq/DDfX/g2VmVlTpJsi/AOb2dbgBIuKSkjbLgIMRsTq9PyDp\nN4HHgfcBq1LTB4H3A9uA3wO25H2mc7aZtZIijem+GLgSeJukHWlZMNhKZmY2JHcA44HNKb+uHsI6\n1wH/BOwCfhQRG1P5GuBMSbuADwE31SJgM7MiKdTwEjMzMzOzIirSme6KJC2Q9GyaiOEj9Y6n0Ul6\nQdKT6WzW44Ov0VokfUHSPklPlZRNlLRZ0k5JmySdUc8YG02FfbZc0h5fmTpepYm+WukYc86unvP2\nwJy3q+OcXZ3RytuF73RLGgN8GlgA/DpwhSfNGVQA8yJiVkRcWO9gGtBasuOp1E3A5oiYTjY+1ZfL\nj5e3zwL4ZDrOZpUMPWh1lSb6aoljzDl72Jy3B+a8XR3n7OqMSt4ufKcbuJBsLOELEXEEuBf43TrH\nVAS+gamCiPgW8H9lxaWTfdxN/yQgRsV9Bj7OTlBhoq/X0zrHmHP28Pn3VIHzdnWcs6szWnm7GTrd\nxyZhSPbQP0GD5QvgXyU9IemaegdTEO0RsS+93ge01zOYAumW9H1Ja3xp90RlE321yjHmnD08ztvV\na5Xf1Ghyzh7ESPJ2M3S6fSdo9S6OiFnAQrJLJHPqHVCRRHb3sY+7wd0FnAu8BdgL3FbfcBpLmujr\nfrKJvg6W1jX5Mdas36vWnLdHoMl/U6PFOXsQI83bzdDp/gn9EzeQXu+p0NaAiNib/t0PbCC73GsD\n2yfpdQCSzgJerHM8DS8iXoyE7LFyPs6SNNHX/cCXI+JrqbhVjjHn7GFw3h6WVvlNjQrn7IGNRt5u\nhk73E8B5kqZK+hXgvWQTM1gOSa+WNCG9Hge8A3hq4LWM/sk+SP9+bYC2xrEE1GcxPs6AASf6apVj\nzDm7Ss7bw9Yqv6lR4Zxd2Wjl7aZ4TrekhcBKYAywJiI+XueQGpakc8nOkkA2I+k676/jSfoKMBeY\nRDZG62PAA8A/A+cALwDviYif1SvGRpOzz5YB88guUwbwPPDHJWPfWpak2cBjwJP0X4r8KNmsji1x\njDlnV8d5e3DO29Vxzq7OaOXtpuh0m5mZmZk1smYYXmJmZmZm1tDc6TYzMzMzqzF3us3MzMzMasyd\nbjMzMzOzGnOn28zMzMysxtzpNjMzMzOrMXe6rSFJmiepV9L7B2/deCQtSfHPrXcsZma15pxtNjh3\nuq2RRVpIs9ctlzSzzjEdk/7ILJP0mpzqKFnMzFqBc7bZADw5jjWkNOXqq4BXIqJX0jzgm8CSiPhS\nXYNLJC0nm/VsakT8uKzuFLKZ446Ef2Rm1uScs80G5zPd1pAi83JE9JZVqRafJ2n8SFYvL4iI3hS/\nk7eZNT3nbLPBudNtDal0fKCkJWRnTADWpvJeSVtL2kvStZK2S/qFpIOSvpnOtpRud2pad5mk96b2\nh4E7Uv0MSaslPS3pQNrWE5KWlm3ni2RnTACeL4npY6m+b3zgJWXrTZJ0p6Tdkl6S9GNJn5Y0saxd\n3/pvk/Tnkv5TUo+kH0q6aqT718xsNDlnO2fb4MbWOwCzQQTwGLACuBn4LPCtVLevpN2Xgd8Hvgqs\nAdqAPwQ2S7o8Ir5ett0u4BxgdVoOpPJ5wBzgQeB5YBzwHuDzkiZHxK2p3WeACcBi4EPA/6TyJyt9\nkTSO8NvAtBTjd4G3AtcCvy3pwog4VLbaivRd7gJeTm2/KOlHEfHtSp9lZlYnztnO2VZJRHjx0nAL\nWSLtBa7Ke1/WdnGqW1pWPgb4d+C5krKpqe1LwJtytvXqnDIBW4GfAWNLypenbZ2Ts86SVHdJSdnf\np7I/KWt7XSr/m5z1t5d95hSgB7in3v9HXrx48dK3OGc7Z3sZfPHwEmsGVwIHgQfTpcBJkiYBrwW+\nAUyVdF7ZOg9FxA/LNxQRh/teS2qTdCZwJrAZOB140wjiXAy8CHyurPyzwP5UX251RLxSEt9/ATuB\nN44gDjOzenLOtpbk4SXWDM4nu2y4r0J9AL8K7Cop25nXUNnNOcvJLk925DR57bCjhHOBx6PsRqOI\nOCppF/CWnHWeyyn7X+DsEcRhZlZPztnWktzptmYgsrMOVwzQ5umy94dzW8E9wKVkZzIeA34KHE1l\nH+bk33x8tEJ5TZ4IYGZ2EjhnW0typ9uKYqDHOO0CFgHfiYhfDPcDJJ0B/A5wd0RcV1b3jipjyvMc\nMEPSmIg4lpgljQWmk3+GxMysiJyzzcp4TLcVRd8d4mfm1N1Ndix/PG9FSe1D/IyjZEn5uN+FpLOA\nqzkxYQ8UU54NwOS0rVLXAJNSvZlZM3DONivjM91WFM+Q3XhzXXpG68+BfRGxNSLul7QW+KCktwIP\nkT0OqgP4LbLHPU0b7AMi4qCkTcCVkn4JPAG8AfgjsjMaF5St8m/p309IuofsDvWnIqL8smiffwDe\nDdyZ4vweMAv4APBsqh8qX6o0s0bmnH0852zzmW5raMfOUkTEL8me6XoAWEk2ju+WkvqlwFVkj2y6\nCVgFvC+1v6mKz7wS+AJwGdnkC+8ke9bsnZSdNYnsmasfIfvj8DlgHfCuvPhT+wPAxWRjDxcBtwML\nyJ7nOjvnMmulS6ExQJ2ZWb04Z+dzzjYAFOHjwMzMzMyslnym28zMzMysxtzpNjMzMzOrMXe6zczM\nzMxqzJ1uMzMzM7Mac6fbzMzMzKzG3Ok2MzMzM6sxd7rNzMzMzGrMnW4zMzMzsxpzp9vMzMzMrMbc\n6TYzMzMzq7H/BxrB1AId2x2SAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xe688d90>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The figure on the left shows the estimates for both $\\mu_a$ and $\\mu_b$ for each iteration and the figure on the right shows the corresponding incomplete likelihood function. The horizontal lines on the left-figure show the true values we are trying to estimate. Notice the EM algorithm converges very quickly, but because each group is equally likely to be chosen, the algorithm cannot distinguish one from the other. The code below constructs a error surface to see this effect. The incomplete likelihood function is monotone which tells us that we have not made a coding error. We're omitting the proof of this monotonicity."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "out, lout = run()\n",
      "mua_step=linspace(0,10,30)\n",
      "mub_step=linspace(0,10,20)\n",
      "z=Lf(xs,mua_step[:,None],mub_step[:,None,None]).sum(axis=2) # numpy broadcasting\n",
      "fig=figure(figsize=(8,5))\n",
      "ax=fig.add_subplot(111)\n",
      "p=ax.contourf(mua_step,mub_step,z,30,cmap=cm.gray)\n",
      "xa,xb=zip(*out) # unpack the container from the previous block\n",
      "ax.plot(xa,xb,'ro')                                    # points per iteration in red\n",
      "ax.plot(mua_true,mub_true,'bs')                        # true values in blue\n",
      "ax.plot(xa[0],xb[0],'gx',ms=15.,mew=2.)                # starting point in green\n",
      "ax.text(xa[0],xb[0],'   start',color='g',fontsize=11.)\n",
      "ax.set_xlabel('$\\mu_a$',fontsize=24)\n",
      "ax.set_ylabel('$\\mu_b$',fontsize=24)\n",
      "ax.set_title('Incomplete Likelihood',fontsize=18)\n",
      "fig.colorbar(p);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAFjCAYAAAD7OehQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXnYHUWZv39/SEISZBMIgkCASACVRXaRdUaDiIrEQYKD\niAzjOALKjN9RFkcBV0DHhR/qjBrAoKJoWMJFWDI4CCI7EYIQloCQRSAGAoSskOf3R/dJTk7Ofnqp\n7n7u6zrXe97q7uqn63T33VVdXS0zw3Ecx3GcYrBO3gE4juM4jtM9Lm7HcRzHKRAubsdxHMcpEC5u\nx3EcxykQLm7HcRzHKRAubsdxHMcpEC5ux2mDpFskPZV3HIMi6VJJKzulpbSecyStlDS6XVpIpFE2\njpMULm5nLSQdGp9U/1/esQRC34MdSDpK0tlJBtOQ/6Xxb7VJh1mN5tuRxkAOjXm2WnfoFDFmpwK4\nuJ12+IlrcI4CUhN3TDe/0yeBkU3SlXAszfL8GjDSzJ5JYV1pkkbZOM7ADM07AMepAGlfAHUUjJm9\nlnIM7db9OvB6Xut3nLLhNW6nKyRtFzfJni3pA5LukbRE0jxJF0ga0mSZHSRdImmOpGWS5kq6WtKe\nDfMdJel2SYskvSLpD5KObJLfXyT9n6TdJE2T9LKk5yV9V9JQSSMl/Ve8niWSfi9p54Y8PhFvx7vj\n+6xPS1oq6QFJE3ooj7GSLpP013jbnorLYb26eW4BPh591cq6z8fr5tlS0o8kPVNXRv8jaVS3sXQZ\nb1f3bOMyvCaO5aN16XtLukrS/Li8Zko6q9nv3iTPdvezR0j6RryPLJX0J0nva5LHUEmnS3o4/m3/\nJulKSbsMOO8ISd+K9+PFku6SdFinbXKcPPEat9MrRwAnAz8CfkrUFPwfwIvAN2szSdobuBkYAkwE\nHgI2BQ4G9gfuj+c7GbgIeAQ4l6j2+AngakmfMrOf1K3bgK2Bm4BfA78B3gucBqwE3ka0T38DGBXH\ndbWkt9rag/KfD6wXr1vAicDlkkaY2c/aFYCkvYDfAS/E5TAXeAfwWeAASYfENdyvAV8CDgI+VpfF\nHXE+o+PvQ+MymgWMBT4N/J2kvc3s5Xax9Ejbmr+kTYFricrxfWb2uzj9/cCVwGPAt4m2+13AV4i2\n+5gBYvoZsBy4ABgO/BvRb7ajmT1dN98vgI8Q/fY/ALYETgHukHSQmf2pz3kvBz4ETAFuBHYAJgNP\n4beKnFAxM//4Z40PcCiRCD9Xl7ZdnPYKMLph/hnAvLr/RSTqxcAuTfJX/PeNwCIiIaxfN30D4Ang\nZWCjuvS/xDH8Q0N+98bpVzWkfyZOP6wu7RNx2lPABnXpG8b5LwBG1KXfAjzZkO8DwMPAGxrSj4rz\nPqEu7VJgZYtyvgZ4FnhzQ/pewArg7C5+q0vjdW7SzXyt0uLfdybRRcjudfOMiGO8BVinYfl/i9d9\nSIf1nBPPN7pJ2pSGefeO079RlzYuTru8Yd7d4nK6tc95D4vnvbhh3g/F6a/nfSz6xz/NPt5U7vTK\n1bZ2J6NbgC3qmonfQVRru8TMHmrMwMxqNZlxRLXeC81sUd30V4ALgfWB9zQsPsfMJjek3R7//f8a\n0v8Q/92hyXb8KF5PbZ0vA/9NdDFxaJP5AZC0K7ArUU1tpKTNap84jsVEQmiLpI2ADxDV9JY35PM0\nUe07kyZbSe8A/khUw3yXmT1QN3kcsDmRkDdpiPP6eJ5B4vx+/T9mdi/RxVz9bzY+/vv1hnkfJGoh\nODBuLeh13qPiv99qmPcaootJxwkSF7fTK082SVsQ/62dEMfGf6d3yGv7+O+fm0x7uGGeGs2eqX6x\nxbRa+qaszSNt0hrXWc9b47/nAs83fJ4juhDZvM3yNXYiapn45yb5PA/s2GU+SfB7os5jB9iazdOw\nensvbhLjI0SyHyTOZvvTC6z5m20fx9fsN2vcT3qZd0w8bzNJN1vecYLA73E7vdKud3AWj8+0W3+r\naUnGVcvr28ANLeZ5sUV6s3wuI7rP24wlPcQ1CL8A/pWor0Djo2u1OP8D+BPNmTfAurP4zRynVLi4\nnTR4NP67R4f5ZsV/dwH+r2Ha2+K/zWpkSfA2oqbTXtdZq52ttLjzVgdaDT7yRJw+vMt80sLM7GRJ\nrwFfkrSumZ1ZN722vYtzjPNJok6ObyPqT1HP24jK8ak+5z2MqPXj4YZ534rjBIo3lTuJE98j/TPw\nT5Le1mbWacCrwGckrV9LlLQBUceyV+J50uDTkjasW+dGRLXOF4majptiZtOJOt79q6S1mtTjR5He\nWJe0KEpeIw0zWwBMBT4sab8m+Si+j9wt3dRQW/aSNrPPAt8FTpf07bpJNxI1i5/RuA1xnCPrf7tO\n6+mTq+K/9RcUxI93HQn8IS7PXue9Ov77+YZ5jyK6VeE4QeI1bictTiR6HOxuSROJRL4xcAhwvZld\nZGYvSfoC0SM7d0m6lNWPg40BPlXfgSxh5sfrvITVj4NtDfyzmS1tmLdRiscTPQ72oKSLiWpr6xF1\nqBoPnAFMiue9g+hRpB9KmkrUs/lOM/sL0WNffwBulTSJqCl6HaJtP5KoCf0rXW7P5yQ1a1q/2czu\naLEda6SZ2f+TtJxI3sPM7DQzWxw/d3418Gi8vbOIfsud4+09Cri1WZ5JYGb/K+kK4Nj44uE6YAui\ncl1M9BheP/PeJOla4ARFQ8beCLwF+Beii7O1nvt2nBBwcTtJsFZzsJndK2kfoueYjwE2Af4G3MXq\n3t6Y2Y8k/ZWo1lO7v/onYLyZTWmynq7W3wWnEz1TfgrwJqLm/ePM7Fed8jazByTtQVSrO5Kopv4K\nURPsJUQXLDUuJ7plcCzRs8W1i4S/mNmc+Jnw04keQfoYsBR4hqi3+RVdbEcttjNbTFtOdPHQrIya\nbduZklYA/ylpqJmdEgtuH6ILko8RPSP/IlFz/3+xZpN0V+tpkda4TfUcR/Ts/yeI+hcsIrq98iUz\na+zc2Mu8E4ietz+OqAf9g0QXI8cBb28Rn+PkSu152mCIr+jfDzxvZrvGaZsQDbixLdGztseY2cLc\ngnQKi6RPEPWQPtTMbu0wu+M4FUXSR4jGG9gZ2MfMaoNGbUf01MHMeNY7zOzkeNpeRI9OjgCmmtlp\ncfpwola4PYmewpnQ5AmOrgnxHvclwOENaWcA08xsR6LazBmZR+U4juNUiRlErS/NLvCfMLM94s/J\ndek/Ak4ys7HAWEk1l50ELIjTv0s0cmPfBCduM7uNtR+nqd3vI/57FI7jOI6TEmY208y6HohH0pZE\nozHeHSdNYrWr6h02GXj3ILEFJ+4WvMnMnou/P0d0T9Jx+iWs+0OO4xSN7SVNl3SLpAPjtK2AOXXz\nzI3TatNmw6o39b0U3wLui8J1TjMzk+QnXqcvzOxSontQjuNUHEnTiJ46aOQsM2sc56HGPGAbM3tR\n0ZsOr5aUaUfGooj7OUlbmNmzcXPE881mcqE7juOEhZmlMgpev+f7+njMbFwfyy8neloDM7tfUu2t\nfnOJHimtsTWra+BzgdHAPElDiV6e9EI/8UNxxD0FOIHohv4JrB44YS2+973vNU0fMWJE2xWMHDmy\n7fROy3ea3u08neLoJa/6eX/4wx9y8sknd565x7zr6Tb2tNY/CBdccAFf+MIXMl9vkixd2vj4ee8s\nWbL2o+AXXnghn/3sZ5vM3f06u5mv0zzNYutm+Xb5tprWal3N5r/22msZN25cx/ma5dlsvsa0bpbr\n9H+zfHrNo135t1p20aJFzWZPjBUrVvQ0/7Bhw/pd1SrZxwMjvWhmr0saQyTtJ81soaSX4wGV7iYa\n7+HCeLGaw+4EjmbNR0Z7Jrh73JIuJ3pT0U6SZks6ETgPGCfpMeDv4/+7Jm1pJxFDWnn1ut4qSttx\nnLWPvcb/2x3jzZYt+rEsabyk2cA7gesk1d6GdwjwgKTpwG+IBoqqPZ58MvBT4HGinue19xlMBDaV\n9DjR63AHejIquBq3mX20xaTG1zsGQ1I7aFLycxyn2owYMSKRFpjGfEaOHNmy5p3UOkPBzK5i9RC6\n9emTiXqGN1vmPqLX/jamLyMaiCoRgqtxJ43XtiP22WefRPOupwy17QMOOCC3dYfOfvutNZR6oUh7\nvxozZkyq+WdJt7fzWh3zZahpF4FSizsLaReltt2NuPMk74Pdxd2aoos7bd7ylrekmn/eLXGtjs1O\nAnfSo9TiTpssxd5rXkW7t+04Tra0OuZ7PRe4vLOntOIOoYm8W6ouvxAP8BBj6kQRY3aKQTfn06qf\nx7IkuM5pSTCotJNYR7fzpPX4Vy94bXtNauXRrlzK1AnHyYZ2HbtCW3+zjmbddD6rnRPy3M4qUEpx\nD4rXXLIjpLJO4uLIhe7kQWg9urO6SAlpm7OkdE3loXRI89p2sUjqAiKkC5GyUbZ9LnSa7cu99Br3\n3ys9SifudoTURO6EUU7++MpgdFt2XsZh0u/vUjtu/HfNh0qJuxMhdkjz2nZ6pPV7F/lkVpbf1kmG\nbisqLvJsqcw9bq9th0We5ZTFuvO65xjavU6nevg5MH0qUePuRtpe2+6OrC6A0iLLdfsJzMmTZvtf\ns+O3U03Z9+PwqIS4O5FlhzQnH7wZz8mStPe1QfoW9DNcqTeFh0XpxR1SL/Ju4uk1v17n7Wf+GkWt\nbVelhl90qlhWWYys2Iu8a/N3qoG3+zjpU2pxhyaaNDr+ZHWgFLXTkp9InKLRSxN3N2m95Nm4TD/H\nj0s8fUrbOS3L+9pJ76Bp7vBVqn2GcuLwDmNOlrTa35qldzNQSv1xFNp+HFo8WVFKcSclbW8ijwit\n5SLE9XUiS3n7hYLTq7xr9CLxRnyfy47SiTtLaTth4r+dUxXaXaT1Iu8avUi8Wb5ONpT6HnfaeG07\n3XWHvq5eCTk2p39C7mfSyz3vZuvzt36FSelq3J0ItYncGZwiiNGbsZ2sSWKfa3Ye8zeA5UelxB1y\nE7nXtouxHscJjW7E3GyeQTudeaUkPyrTVJ6ktPNsIi8SZd2uQfFycfKgl2eznbCphLizGOSgkbSu\nRotU286KIp5oQo65SL99WWm1f7Qb9WyQfJvN5zIPl0o1lbejjE3keeJN5Pni99KdVvS7b/ixFg6l\nF3fITeRpU/badtFPJC5XB/LZD8qy75VhG/qh1OLOo4k8rTxDfuSkkSxiLbq0HSdpepVx4zFUVQkW\nkdKKO+ladLeE0iGtCE3/Ia/D6Z6y1N6yoJshRgfJY5Ce4kkfV75PpEcpO6el0VGjG7xDWvqUTdpl\n256i0WsnsCTy7nbedp3CuhkYJe+OZb5vp0fpxF0EaYfYRJ4ERYrVcfKi14uCdsdVLy18eYu8aEj6\nqqQHJP1J0s2StonTt5O0RNL0+PPDumX2kjRD0uOSvl+XPlzSr+P0OyVtO0hspRN3N+TVjJ5Wnknk\n77Xt8uHllTxJlWm3r+msnzZI7btVfp0+FecCM9vdzN4BXA2cXTftCTPbI/6cXJf+I+AkMxsLjJV0\neJx+ErAgTv8ucP4ggZX2Hncr0pB2KE3k/VKEDmllPon4PeLik9Rv2CmfdtNrx3GSQ5GW+bjrhJm9\nUvfv+sDf2s0vaUtgAzO7O06aBBwF3AAcyWrxTwYuGiS2StW485a2y60/yrpdVaOKv2MSzdzN8uyU\nb+PH6Q9JX5f0DHACcF7dpO3jZvJbJB0Yp20FzKmbZ26cVps2G8DMXgNekrRJv3FVRtxFk3aROqRV\n8YTsOGmR1rnKZb42kqbF96QbPx8EMLMvmtlo4FKiJm6AecA2ZrYH8Dngl5I2yDLuSjSVF00sRYrX\nWxGcUCjiLYdWj3Z1uy2146OqLwnpdFvgrrvu4q677mo53czGdbmqXwJT42WWA8vj7/dLmgWMJaph\nb123zNasroHPBUYD8yQNBTYysxe6XPdalF7cafXgLkMTeegHbpWkXUTpOOHg+09z9ttvP/bbb79V\n/190Ufe3liWNNbPH438/BEyP0zcDXjSz1yWNIZL2k2a2UNLLkvYD7gaOBy6Ml59C1Nx+J3A0cPMg\n21VqcRdR2t5Enn7ejpMlnaTartZdo5fadyMu9L75pqSdgNeBWcCn4/SDga9IWgGsBD5lZgvjaScT\nNauPBKaa2Q1x+kTgMkmPAwuAYwcJrLTiDkHaveKycpKm35pYEiN8VY00a72DNIl3e15xwa+JmR3d\nIv1K4MoW0+4Ddm2Svgw4JqnYSinuUAY4KXMTeZHLLVS8uTM9Qi/bbi+UBhnStJe8nbApXa/yNKUd\nUhN5nri0HSd5BhlExakWpaxxd0PRpV3GDml+AnKKQD+197Rr/M2OnZBbGJzBqKS405R2FpS1idxx\nQm/STpsk+xb4sVpeStdU3om0pV3mJvI08XKI8HJoT5YX0f2sq9/xxhvXG1plwQmLStW4yyDtMta2\nXVZOniTd9N3NeOPQuSm71THrvf1XU9XWmcqIu8rSTgKXtlNF+m267uZioN9HvFzoTunF3Y8YQpT2\nIHizm+MkTyfxZjFsaT1+nFeHUt/jDlHa/eJN5I6TD52OoU73tbvFH+9yuqW04g5V2t5Enm6+ZSDp\nssmirIvwe+Y1/kC/x7xL3GlFKZvKXdoR3nTmON3R73ji3eaR9pClVe2kVVUKJW5JZwIfIxrYfQZw\nYjwG7CrKJO1B8CZyp4rk/Rx4tz3KIVnZ+jFVLQrTVC5pO+CTwJ5mtiswhAHfsBL685J+MDpO8ej2\nuK1vDvdj3emFwogbeBlYAawXv4h8PaKXk/dFv8KuShO517bzpazlVOTtSnPY49oyLnGnGwrTVG5m\nL0j6L+AZYAlwo5n9bz95hS7tQXBpO0lRtVd7Jt3MPkh+fl+7O6paBoURt6S3AP8GbAe8BPxG0nFm\n9ote8kl6GMOQlnMcJ116vZhJ6hntdnk71aMw4gb2Bv5oZgsAJF0JvAtYQ9yXXXbZqu+77bYbu+++\nO5BtLTuP5cBr247Tibw6r+XdaS4rHn/8cZ544om8wyg9RRL3TOBLkkYCS4H3AHc3znT88cevtaBL\nO/3155Gv4yRNN4IdZBjUesoo8rFjxzJ27NhV/9944405RlNeCiNuM3tA0iTgXqLHwe4HftxumUFk\n5s3cTt5UpZbWKyHce08ihk7nCv/tnVYURtwAZnYBcEE38xZN2l7bdpyIQS9YunmWOum3d6VxIVGk\nY8cvMrKlUOLuhkEl5tJOjiKdeJxq0e3bu7oVUtoSDx0/1rOldOLul7zEm7e008IP5LCoWrN7Utvb\nT6/wxuOyiiJ30sXFTTGlnRQhxOA4vZLUPeZuhTzIhUCzC2yXuTMIlRd3UaXtTeSO056km8OTfCa7\n3+PXhb8mVWpFqqey4h5ULkWXdlq4tJ2QSOMWQZ63HUI+9p3sKNJY5Ynh0nbBFoUi/U69vFwjiXny\nPBb8BSFOnlSqxp3EAebSzj5fxxmEbpvMof+m12b7flWbcZ30qUSNO4mr4kHzKLvUyr59TjVIsgbd\nWCv32nmxkHSOpDmSpsefw+P07SQtqUv/Yd0ye0maIelxSd+vSx8u6ddx+p2Sth0ktlLXuJM8APOO\nI/TatuOETK/3peuPkzTukWeB1/gHxoDvmNl3mkx7wsz2aJL+I+AkM7tb0lRJh5vZDcBJwAIzGytp\nAnA+cGy/gZVW3Hk3iyeZR+jS9osBJy96eSys305lab7hK038uEwEdT2jtCWwgZnV3qExCTgKuAE4\nEjg7Tp8MXDRIUKVrKk+qOSqUPLwXqZMnvv+txpu6K8lnJD0gaaKkjevSt4+byW+RdGCcthUwp26e\nuXFabdpsADN7DXhJ0ib9BlXaGvcglFHaXtt2nGQe5arCW76qgqRpwBZNJn2RqNn7K/H/XwX+i6jJ\nex6wjZm9KGlP4GpJb88i3hou7jpCuSeeNC7tYlO14UprpLXdSefby3FQxd8xTTrdJnnggQd48MEH\nW043s3HdrEfST4Fr42WWA8vj7/dLmgWMJaphb1232NasroHPBUYD8yQNBTYysxe6WXczXNwxoUk7\n9PvajjMIScqzn+FP87oYKsrxWJYLjN13353dd9991f8///nPu15W0pZm9tf43/HAjDh9M+BFM3td\n0hgiaT9pZgslvSxpP+Bu4Hjgwnj5KcAJwJ3A0cDNg2xX5cWd5IEUmrTTpCgnIKfYdCvYfuUN5ZFU\nkvjxDcD5kt5B1Lv8KeBTcfrBwFckrQBWAp8ys4XxtJOBS4GRwNS4RznAROAySY8DCxigRzlUXNyh\n1bLB72s76ZFFLTPv4UD7Gcvb71k7zTCzj7dIvxK4ssW0+4Bdm6QvA45JKrZKijvEWjYUQ9pFxE/M\nYdKt5Pt5L/YgL+Nodez4fuOEQqXEnbTMqijtol0QNIvXm0iLR681+SRe+9kshlb4vuRkSSXEXRVh\ng0u7Ri8vsvCTbn6k+U7s+mMr7ddhJn1s+D7ptKPU4k5DNCFL24no9TdygReHfu+hJ9GEniVFukh2\nsqeU4i5Cp6w0pF312naSY8q7xFfTa7NzP+OC9zOO+CACh+JI3HEaKZ24Q69lg0s7DfqJ8Z/+6Ss8\n/fTaQxFvu61x8cVfdnkHzqA92F3ixaeqx2jpxJ0kaQjLpZ0sg8T39NPi1lvPXSv94IPPXpV3VU8M\nRSGp36jdcelSd0LDxd2EtGRVNGmHTpW2vewXEYNsX9pl08tx65J3ssDF3YBLO7v8ByGr2MouzLIQ\nSgfDEDuc+sVE+XBxx6QpApd2soQcmxPRzwVPUm/uylveoRHixYQzGJUXd9GEDdUWV5W33ekOl7dT\ndior7rQFUOSr3FDlmHRc225rqzqiNaY3rtdFkA1JlXWzfcV/Q6csVE7cRRd2lZvIk+bii7+cdwiZ\nM8gAJlm9OjPN93C3wqXuFInKiDsLIaUp7SoJtRlV3/6iMoi8ITuh9rN/ueydvCi9uLM64ZdB2qHK\nMdS4qsYgr80M9VGvQchzvwy1TLKmquVQWnGXQdjg0nacUB71Cgk/XqtN6cSd5Q7t0k6fkGOrInnU\nuuvzqMdF7lSV0ok7C7LoMe7C8jIoG0k3e/tLYZyqsk7eARSNsknb5dgdZSqnPLclzTf31X8cp8x4\njbtLsnou26UdEXJsVaff5vIaWXQ462b/8Vq6U1Rc3B3IciAVl3ZEyLE5yRBCh7NB9zMXv5MXLu4W\nZD3ymcvK6URIj0YNWuuuEYLA+yWEY7aI5eYMjou7gTyGKs36BBDCCacVIcfmpEeRBZ4nfrxUExd3\nTBWEndc6uyXk2Jy1SarWXY/3FHd6oar7SKXFneeLQFzaTp6E1OzeCq+FO05zKinuvN/c5dJem9Dj\nc/LD3/TlOGtSGXHnLWvIT06hSzH0+JzWpNFc3g2d9hkXu1NmSi/uEIQNLicnPJJqLq8dY3kIvBWD\nHG8ufQdA0q+AneJ/NwYWmtkekrYDHgFmxtPuMLOT42X2Ai4FRgBTzey0OH04MAnYE1gATDCzp/uN\nrZTiDkXWNco4UpXjNBKiwPshq2PGLxDCxsyOrX2X9G1gYd3kJ8xsjyaL/Qg4yczuljRV0uFmdgNw\nErDAzMZKmgCcDxzbZPmuKJ24XdphrNspBml0UiuLwNPGLxCKgSQBxwB/12G+LYENzOzuOGkScBRw\nA3AkcHacPhm4aJCYSifuUMhbmnmvv1uKEmeZSauHuQs8DPI4xkp2sXAQ8JyZzapL217SdOAl4D/N\n7A/AVsCcunnmxmnEf2cDmNlrkl6StImZvdBPQIUSt6SNgZ8CbwcM+CczuzPfqNbEReSkSREe42rE\nBV49inIelDQN2KLJpLPM7Nr4+0eBX9ZNmwdsY2YvStoTuFrS21MOdQ0KJW7g+0Q3/I+WNBR4Q94B\n1RPKzhpKHGWiaLLslSwuCFzgTtJ02pceffRRHnvssZbTzWxcu+Vjz4wn6lRWW2Y5sDz+fr+kWcBY\nohr21nWLb83qGvhcYDQwL85zo35r21AgcUvaCDjIzE6AqLmBqJkid0ISZUixOMUiq9p8s34oLnMn\nDXbaaSd22mmnVf9fd911vWbxHuARM5tXS5C0GfCimb0uaQyRtJ80s4WSXpa0H3A3cDxwYbzYFOAE\n4E7gaODmPjcJKJC4ge2B+ZIuAXYH7gNOM7PFeQUUmiRDi8dxuqVVp1IXupMzE4DLG9IOBr4iaQWw\nEviUmdV6nJ9M9DjYSKLW4Rvi9InAZZIeJ3ocrO8e5VAscQ8laq441czukfQ94Azgy/Uzfec731n1\nff/992f//fdPJZjQJBlaPN1QxJjLTmj30Lt9SsQFHwZ33XUXd911V95hJIaZndgk7Urgyhbz3wfs\n2iR9GVHP9ESQmSWVV6pI2oLoQfft4/8PBM4wsw/UzWOzZ89ONY5QZRNqXO0oUswhySyLWELa3jzw\nC4Fk2HHHHTEzpZG3JPvv//7vnpb513/919TiyZLC1LjN7FlJsyXtaGaPEd17+HNW6w9ZMiHH5hST\n0GreWZPFeBB+ceD0S2HEHfMZ4BeS1gVmAWs1YyRN6FIMPb4yEJrAspKqv50rXdK8OPCLgnJTKHGb\n2QPAPmmvpygyLEqczShy7FWj6rXvIuIXBeWmUOJOmyLJpEixOsXHa99OjdCGla4iLm6KJ8GixeuU\nBxe4ExJV3Q/XyTuAPBkxYkThJFi0eFtRlAMu1Djz3g9qx07ecThOFalkjdtPNk4ZCOXec7PjKYS4\nHKesVErcRRd20eNvZOnSpUFvUxHkE4q8G2n1u4YYq+MUjdKLO2Qx9EJZtqMoFEkwocq7Gd3ux0XZ\nHsfJg9KK20VXDEKsdbs08iepfcJ/S6eMlE7coUkgCcq4TaFS1BN9kWrdWZL0seNl7IRA6cRdNqog\n7VBq3UU/Kbu80yeJ/dR/I2dQXNxOENROZnkIvEwnUpd3+PSzj/tv6tTj4g6YEGqhWZN17buMJ0SX\nd/lw2TenCtvYDBd3oFRR2jWyqn2X+aB3eTu9HD++rxQLF7cTLGnWvv1E5Tir6fU48+MnX1zcAVLl\n2nYjadS+q3LS8Vq3kxZ+jsqXRMcql7STpGmSFkt6XtJvJB2c5DrKjh8QzUlKQFUTmY8n7jjlI+ka\n93nAL4GGSBUaAAAgAElEQVQfAzsBBwI3SPoZcLKZWcLrcypEL03nVRN0J+rLzcvGcYpN3zVuSSMl\nqSH5NjO7xMx+Y2ZfM7PDga2BEcBZgwTqOBBJp5uP05r6N3uFVhufOGMi8xfP72mZ+YvnM3HGxK7m\nvfyRy5m1cFY/oQFwwd0XsOL1FX0v7zhJMEhT+R3Ai5Kuk/QFSfsByySNqZ/JzF4wsxOBzQcJ1HGc\ndGgUeV5CnzhjImfcegbjrx7ftbznL57P+KvHc8atZ3Ql71/N/FVf4n595esAfPueb7N85fKel3ec\nJBmkqfxC4LPAbcC7gf+M8/uypAuB24F7zOxVSesS1bodxykIWd+WOPItR3LJjEt49MVHGX/1eK46\n6ipGrTeq5fw1aT/64qPs9MadOPItR66adv2T13Pe3eexjtbh9ZWv882Dv8nTLz/NA/Mf4Iu3fZFv\n3vVNzn3XuYxabxRf+P0XWPzaYpa9tozj3348n9r9UwCcevOpDNVQZi2cxaIVi9h3y30BeP/k9yOJ\na466hg2Hb5jItjtOLwwi7p8Di8zsCuA8SUOBvYFDgEOBLwDrS3oOWAmcNmCslSCU4T8dp1uS2l+3\nGbEN1x97PUdccQQzF8xsK+9GaTfOd/7d5/OdQ7/DXlvshZnx6opXOWCrA7ji0Ss45R2nMG67cQAs\nWr6IyR+azLpD1mXR8kUc/tvDeffod7PDG3cA4OEFD3PN+GsYOXQkAJfMuISp/zCV9Yatl8g2O4NR\n1dtifYvbzJYDV9T9/xpwZ/w5X9IQYCywPvCImb06YKyO45ScUeuNYuoxU1fJ+8NTPszUY6auIeX5\ni+fz4Skf5tEXH2XnTXdeazrAodsdypf++CU+NPZDHPLmQ9h5051XTTNW95Fd/NpiPv/7z/PwgocR\n4tnFz/LQ3x5ihzfugBAffMsHV0nbcUJh4F7lkoYT9SKfBlxVE7SZvQ7MHDT/KuK1bqfKNMr70IsP\n5sjfb8/Gi4yF64sphzzFM8vmtJQ2wHmHnscjf3uEW565hX++6Z85da9T+cRun2CdddZh+LrDVx1f\n5/3+PLbacCsmfmAi62gdjvrtUSx7fdmqfLxm7YRI153TJB0jabykTevTzWwZ8G/AYmCipE8mHKPj\nOBWjJu/Rw7fmmWVz+N2Y25jw2B/43ZjbeGbZHEYP37qltAEee+Ex3rrZW/n0np9mwtsmMP256QBs\nsO4GvLTspVXzvbzsZbbaYCvW0To8/LeH+ePcPzJs2DBGjBjBkCFDGDZ02Bod9jZYdwOWaZlfWDu5\n0kuN+1vANoBJmgncCvweuNXM5gFXAldKmgr8JPFIHcepFKPWG8WRv9+e342Zw8Obw64nR+lvex7e\n/dQYRp3auuPaObedw6yFsxiqoWw8YmN+8N4fAHDibidy1i1n8f17vs/XD/k6X3jnF/jk9Z9k0kOT\n2OGNO3Dg1geukU/jE6+f2fszvP+K97PesPW47pjr2Gj4Ri1jqOr917IgaXfgv4E3AH8BjjOzVyRt\nBzzC6hblO8zs5HiZvYBLiTpjTzWz0+L04cAkYE9gATDBzJ7uO7Zux0SRtDlwFNGgKgcB29ZNfhK4\nD1gIvMPM9us3oEGQZPPn9/YMaMj4Vb1Tdc5/3/uY8NgfVkkbYMYP4dc7Hsjp11+fX2ADUhWpjxo1\nCjNrHO8jESTZOeec09My55xzTtfxSLoH+JyZ3SbpRGB7M/tyLO5rzWzXJsvcDZxqZnfHldgLzewG\nSScDu5jZyZImAOPN7Niegq+j66ZyM3vezH5sZh83s+2Jat//CPwQeBU4Gvg7vPe44zgJsXB9MeHo\nNdMmHA0vbZDoaM2Z0+rZeb9YD4qxZnZb/P1/gX9oN7OkLYENzOzuOGkSUWUX4EjgZ/H3yUSPUPdN\n33u/mc01s1+Z2almtjswDngI+NMgATmrqcpVueM0Y/7i+Uw55Cke3jxqHp/xw+jvw5vDNQc/2fMI\na0XBhR4Mf5b0ofj7R4gqqzW2lzRd0i2SavdXtgLm1M0zN06rTZsNq57AeknSJv0GlthY5Wb2O0mv\nAOcAZySVr+M41WP+4vkcccURqzqivfupMfx6x5W8+6l1WLTNkzyzbA5HXHFE2w5qZaOZvP3ifjAk\nTQO2aDLpLOCfgAslfQmYAtSGzJsHbGNmL0raE7ha0tszCTgm0ZeMmNk9kj6fZJ6O41SLmrRnLpi5\n+pGvuo5o/143vWrybqTqMl+yZEnb6c888wzPPPNMy+lmNq7DKt4LIGlH4P3xMsuJJW5m90uaRTRm\nyVyid3PU2JrVNfC5wGhgXjxY2UZm9kKHdbekl8fBbpN0s6QvSTpI0rAWsw7vNxhnbap0EDpOU2k3\nSLn2qNjOm+68St5lbTbvB29mX83o0aM58MADV316QdKo+O86REN6/yj+f7N4gDHid3OMBZ40s78C\nL0vaL34B1/HANXF2U4AT4u9HAzcPsl293OO+H9gfOJfoMbCFsci/LOlDkg6T9F3AhxlyHKdnupF2\nDZd3b7jI++Kjkh4levRrjpldGqcfDDwgaTrwG+BTZrYwnnYy8FPgceAJM7shTp8IbCrpcaJxTwa6\nndx1U7mZnRY3g+9LNB75wcC7iHqS1/MbSe8G/mhm7dsxnK7wkdScKnDVo1d1Je0ajSOsXfXoVfzL\nHv+SUbTFZsSIEd6a1wEzu5DoZVqN6VcSjVvSbJn7gLUeE4sHKjsmqdh6uscdt+3/If58PW6r34vV\nIj+QqPfdR4AVku4lqp1fbmYzkgq6itQfZC5xp4zUpDt+p/Fd37Ouydul3Tu184gLvHh0PQBLV5lF\n9wLewWqRHwRsAjxrZm9ObEWt11+qAVi6wSXuOM6gpCXvtAdgOf3003ta5vzzz08tnixJulf5SqJ7\n4fcD341v0O8CtB4X0BkIr4k7jjMo3nReLBIVdyMWVee9iTwjXOKO4/SLy7s4pCpuJz9c4o7j9Irf\n9y4GLu4K4BJ3HKcXilL7LkKMaVDskfqdnlm6dOmqj+M4Tiv8Ij9cvMZdYRrl7Qeq4zj1FKXmXTW8\nxu2swmvjjuM04hf04eE1bqcpXht3HKeGd1oLCxe30xUucsdxvOk8DFzcTl+4yB2nmri888fF7SSC\ni7x85P0buhzCxeWdLy5uJxVc5MUixN+nMSYXheNEuLidTGh10g1RGFWiSOXvtbxwCOV3CCWOrCmc\nuCUNAe4lerH5B/OOxxmMdgdekaRSRLx8nX6oqixDonDiBk4DHgY2yDsQJ128lp4OoZbfbVOnMu2i\nixi2bBkrhg9n3KmnctARR6wxj9e688XLPgwKJW5JWwNHAF8HPtdsnm52rFBPXE53+Njr/RFyWd02\ndSrTPvc5vjFr1qq0s558EsDlHQhe5uFQtJHTvgt8Hlg5SCb1I4R183HCxX+nzowYMSJoaQNMu+ii\nNaQN8I1Zs/jfH/yg6fyhb0/Z8OMrLApT45b0AeB5M5su6dBW833nO99Z9X3//fdn//33H3jdXosv\nBt6TfU2KtP3Dli1rmj7UhZE7vUj79ttv5/bbb08xGgcKJG7gXcCRko4ARgAbSppkZh+vn+lzn2va\ngp46LvfwqHKTetG2d8Xw4U3TX2uzHd5kni79lO0BBxzAAQccsOr/b33rW0mG5MQURtxmdhZwFoCk\nQ4D/aJR26Ljc82Pp0qWVKdtQt7Nd57Nxp57KWU8+uUZz+ZljxnDYKae0zdPlnQ5epmFTGHE3wfIO\nIA388aj0qIK8Q92+Tp3PagL/4g9+wNClS3ltxAgOO+WUtTqmNcPlnSxeluEjs/L4T5LNnj077zBy\nI9STdmiUtZxC2q7G2vWz8+dz8Z/+tNZ8XzzsMM697rqB1+eySYaky3HUqFGYmRLNNEaSffzjvTW6\nTpo0KbV4sqTINW6nAa+td0cZa955bU9N0K/Om8ezf/0rG225JUOGDWPZs89y8bPPrprv0yNGcCtw\ncMPySXU+81r34Hj5FQcXd0XwwUzWpEzyzno79n7DG9hy+XLWBxYBs4GtgC8CN77wAu8FboQ1RP2j\npUv5EmuLu13ns15xefePl1uxKNpz3E7CVPnZ9TJsZx7S3mf5cq4HfgNcD+wPzAW+BrwXmEY0QtK0\nhmWfboj1zDFjeE+Hzme9UpaLsSwpw3FQNbzG7TSlKjX0Ite884h7y+XL+UlD2k+A9xGdTKYBQ+L0\nIQ3zrf/Wt/LFUaN67nzmpIdLuzWSPgKcA+wM7GNm98fpmwCTgb2BS83sM3XL3AJsASyJkw4zs/mS\nhgOTgD2BBcAEM3s6XuYEogYrgK+Z2aROsbm4nZ4oo9Br21TkbciK9dukD4k/r8dpr9dNP3PMGD5y\nzjku6oBwaXdkBjAe+J+G9KXAfwK7xJ96DPjHmuTrOAlYYGZjJU0AzgeOjS8CvgzsFc93n6QpZraw\nXWAubicRmp0EXITpkkdrwaI26UOBR4BTgH/fYgte3nJLzt5gg8xq1y6i7vBy6g4zmwkgqTF9MXC7\npLEtFm3Wa/1I4Oz4+2Tgovj7e4GbaqKWNA04HPhVu9hc3E5quMzLx1/XXZdPNjSX/zNRB7VNgBGj\nR3PjzjtzZIbN4C6i7vBySpxWz1L/TNIKYLKZfS1O24roMMHMXpP0kqRNgTcDc+qWnRPP2xYXt5Mp\nIcu8yPe7s+LeV19l7ze8gfc19Cp/85vexJk//nGmTeEuos54GbUnruFu0WTSWWZ2bR9ZHmdm8ySt\nD0yWdLyZXTZYlGtTOnEvWbJk1feRI0fmGInTLX6PuX/yuNi499VXM11fIy6jzlSljDpt5/PPP8/z\nzz/fcrqZjUsyHjObF/9dJOmXwL7AZUQPXowG5kkaCmxkZgskzQUOrctiG+B3ndZTOnHXUy/xNPEL\nhGRwgTvtqIqMBsHLaE0233xzNt9881X/P/zww/1m1ey+9RppkoYAbzSzv0kaBnwQuCmePAU4AbgT\nOBq4OU6/CfiGpI3j/MYBp3cKptTizookLhBc/qvJs8nam8vDw2XUGS+j5JE0HrgQ2Ay4TtJ0M3tf\nPO0vwAbAupKOIhLuM8ANsbSHED0dWesOMhG4TNLjRI+DHQtgZi9I+ipwTzzfuZ16lIOLOxh6lX/Z\nRe+17+4p68WGy6gzXkbpYWZXAVe1mLZdi8X2bjH/MuCYFtMuAS7pJTYXd0HpVvRFF7wLvJq4kNrj\n5VNtXNwlp5PgiyL2LGuVZa3BFgEXUnu8fBxwcVeedmIPTepe+25N0S82XEjt8fJx6nFxOy1pJvUQ\nZF50STmrcSG1x8vHaYaL2+mJUGSetryLeHFQpJhdSO3x8nHa4eJ2BqZR5lmJvEiiciJcSJ0pSxll\nsR1lKatecXE7iZPl6HUu72JQ1RNsL5SljMqyHSFTOnF3s9P4iT47ahJPU+BpybuIFwWhxewn8c6U\noYzKsA1FonTi7oYQdrKQTq5ZsGTJkiA6tlWBPHvfh3BsFYEylFMZtqGoVFLcIZDETl80+bu8s6Vx\nH0tyf/GTdu8UvcyKHn+ZcHEXmG4PpJAEn1bTeWhNxCHSbn9pVnZ+ok6GIpdjkWMvMy7uCtDrCTsL\nilL7rsoFgZ+gk6XI5Vnk2KuCi7viNDtIsxJVFh3XHCdLiiq9osZdVVzczlpkLfOkat9VqR07YVFU\n6RU1bsfF7XRJlccJ9wsCpxlFFV9R425GmbalF0on7kF/SD9BtyctgRflnrdTbYoqiqLG7TSndOIe\nlKx38KJeKKQh8CTk7bVjJw2KKL4ixux0h4s7Z/o5uEISU1Wa0P2CoJoUTX5Fi9fpDxd3AQlxWNf6\nmAZZtzeZO3lTRPkVMWanf1zcJSXPZ7fzrp3mvX6nmBRNfkWL10kOF3cFSap23Gkd/eYdaq3bLwjK\nSdEEGHq8ja/5dZLHxV1x0h7POi95u2SdToQuwEZCjtdlnS0ubmcNkq6Nu0CdkAhZfs0INV4Xdb6U\nTtwh7ehFF1ZVeoz3gl+IFJOQzgudCDXWEGUdYkxZUDpxh0RZXt05qKz6XT7Ue91OMQhVgM0INdaq\nijF0XNyBE8qjX17TdIpAqAJsRoixuqiLgYu7BOT5hq+QSetiwy9iwiNECTYjxDhd1sVjnbwDcNJh\n6dKla3ySyjOPZR2nGUnv32kSWpxLlixZ9XGaI+kjkv4s6XVJe9al7ytpevx5UNKEumm3SJpZN31U\nnD5c0q8lPS7pTknb1i1zgqTH4s/Hu4nNa9wVIane4lnXNv0+t9NISAJsR2hxuqR7ZgYwHvifJul7\nmdlKSVsAD0n6rZm9Dhjwj2Z2f8MyJwELzGxsLPrzgWMlbQJ8Gdgrnu8+SVPMbGG7wEon7rx2ziLJ\nZdDe4t5U7GWQNaFJsBWhxemy7h8zmwkgqTG9vlBHAi/F0q6x5gIRRwJnx98nAxfF398L3FQTtaRp\nwOHAr9rFVjpx50USB0jW8s9aPnnIzgVbbEITYStCijMUWYdUJkkjaV/gEmB74KMNk38maQUw2cy+\nFqdtBcwGMLPXJL0kaVPgzcCcumXnxPO2xcUdEJ0OuDTE3q/YXIheBmlShJN+SDG6rPsjruFu0WTS\nWWZ2bavlzOxu4O2SdgZukHSLmb0EHGdm8yStD0yWdLyZXZZ03C7uAtHq4EziHdbQe9N5P+LqZxm/\nz10dinDiDynGEISdZ3l0WvfLL7/MK6+80nK6mY0bZP1mNlPSLGAH4D4zmxenL5L0S2Bf4DJgLjAa\nmCdpKLCRmS2QNBc4tC7LbYDfdVqvi7sE1B+8WY/v7bVOJwlCkmErQokxb1mHUg7dsOGGG7Lhhhuu\n+n/evHn9ZrXqvrWk7YA5cZP3tsBY4HFJQ4A3mtnfJA0DPgjcFC82BTgBuBM4Grg5Tr8J+IakjeN1\njANO7xRMYcQtaRtgErA5Uc+9H5vZhflGFR6DSryf2neV5V3lbU+C0CUQSnwu6+yRNB64ENgMuE7S\ndDN7H3AQcHp8H3sF8C9m9rKkNxA1mw8DhgDTgJ/E2U0ELpP0OLAAOBbAzF6Q9FXgnni+czv1KAeQ\nmSW2oWkSd7vfwsz+FN8/uA84ysweqZvHrrnmmtRjKdqJut9aeC/b2WuZ9Dr/oE3laf5mRdsfQiB0\nEYQSX57CTqIMdtttN8ysWS/rgZFk++yzT0/L3HPPPanFkyWFqXGb2bPAs/H3RZIeIeqR90jbBVMg\nqYM6qxN+v7XwXmrfadc8/T53OQhFiK0IIT6vXTudKIy464nvMewB3JVvJIPR7gBJS4IuwOTx5vLO\nhCyDEGIres3ayZbCiTtuJv8tcJqZLco7nrRIc/zxtOTdi8BcduUndCGEEF9ewg5h253+KZS445v+\nk4Gfm9nVzea5/PLLV33fZZdd2HXXXTOKLn0aD7asxOeS7YyX0WpCl0Le8ZVZ1vfccw/33HNP5xmd\ngShS5zQBPyMa7/XfW8xjv/71r/teR1FPvP3G3Uutu9t1pNmhLeQOalnkHzp5C7EdecdWZlm3W+e+\n++7rndNSoEg17gOAjwEPSpoep51pZjcktYK0dvK0T+he2wuDqv4OeUuxHXnGVlVZV2G9eVMYcZvZ\nHyjoa0iz6IRWhRHJihBvleQd8kmzasLOentD/u2rQGHEXVaS7ISWpjS6zbtK4mpF2csg5JN2XrG5\nrJ0scXEHyCDvzu5VGkWoxSZF2YWaNiGfuKsibJe1AyUUd787Wqgn9DIPQVqUOPuhbNsW6gnchV3s\ndTn9UTpx90uSO2saJ+xeBd6LOLqtdZdNRmlThvIK9SSeR1wuaycUXNwpkGZntF6a0csgjqJT1N8g\n1BO5C7t463GSx8WdMUl3Ruu0fLfiSPJed8iyyiO2kMujkVBP5mUXdhllHeq+VAZc3AEw6IhoncRQ\nJHF0okqd6bIi9BNs1vG5sMNfT9Upnbj7OehCE0EaHdK6kXc3UizTRUCWhFhuoZ9kyyrsMkk0730o\n7/XnRenE3Q9JHLBpvbQD8h2kJW1CjCktQtnW0E92Lmxfh9MeF3dCtDv4s6rRhyIGpzV5/kZFOOFm\nGWNZhO2yrh4u7gyonSD6FXiWJ/sq3EPO+wIn6/UX4aRbNmGXQaZF2G+qios7QwYVeBLkLS0nO4py\n4i2TtIsu06LsM1WndOJOc8dLSnj1J49uJV60scJDiSNUBhnWttt8Q8eFnX/eWeTvJE/pxJ0mrXbw\nQU68IdTCnXwZ5HHAop50s4o7TWEXWahp5p3X60yrhIs7AZIYVKUbgYdSi807jrLfh++0PxVV1jWy\niL+owi5q3i7rbHFxp0S/TaGdpJS3NMtC0cqx6LIGr2XnkW+aebus86N04h5kJ837XdZFXZ/jdMJr\n2eXINzRZhxZPVpRO3IPQ7c7ejxSTfFvXoGJ2sTtZUmRppxG7y9oZFBd3H/TbDN7LSGh53sct+z1k\nJxtc2OnnmVa+LuuwcXEPSD89gpN4Y1eZXixSpFid7iiqtIsibK9dV5vSibuXHToNWXQr8iTGIXfh\nDYaXXzqkLW0XdrIUpTyzRtJHgHOAnYF9zOz+OH074BFgZjzrHWZ2cjztFmALoFaoh5nZfEnDgUnA\nnsACYIKZPR0vcwLwxXj+r5nZpE6xlU7cvZDmPe36dQxSM/Zma6dIuLSLIewilGMAzADGA//TZNoT\nZrZHk3QD/rEm+TpOAhaY2VhJE4DzgWMlbQJ8Gdgrnu8+SVPMbGG7wCot7m4Z9DntbuTdLk+XtxM6\nRWwaD13YRahdl1DWqzCzmQCSel202QJHAmfH3ycDF8Xf3wvcVBO1pGnA4cCv2q1gnV4jciKWLl3a\n006bR8eUfqclhQ/4UA2yqGWHLJxezwVZ51crv6TKsBZfmaXdBdtLmi7pFkkHNkz7WTztP+vStgJm\nA5jZa8BLkjYF3gzMqZtvTjxvW0pX4261c6ZVY+3lPmmZ7qmWaVuc/ila03jSgk2S0JvDyyjquIa7\nRZNJZ5nZtS0WmwdsY2YvStoTuFrS28xsEXCcmc2TtD4wWdLxZnZZ0nGXTtyt6LQTDyL2pCTmMswe\nL/P+cWmHl1fIZZYGneJbvnw5K1asaDndzMb1uk4zWw4sj7/fL2kWsCNwv5nNi9MXSfolsC9wGTAX\nGA3MkzQU2MjMFkiaCxxal/02wO86xVAZcXei2Q7fi8xdAE6VqKq0qyDs0GXdC+uuuy7rrrvuqv8H\nKKdV960lbQa8aGavSxoDjAWelDQEeKOZ/U3SMOCDwE3xYlOAE4A7gaOBm+P0m4BvSNo4Xsc44PRO\nwbi429D4I3cSeTfydsE7RadI0g5RtEnlE2o5pZFfHkgaD1wIbAZcJ2m6mb0POAQ4V9IKYCXwKTNb\nKOkNwA2xtIcA04CfxNlNBC6T9DjR42DHApjZC5K+CtwTz3dupx7lUEJxJ/GmrlZ007s7LTF7z3Kn\nCoQoo9DySbKTWZKUQdb1mNlVwFVN0icT9QxvTH8V2LtFXsuAY1pMuwS4pJfYSifuZrTbofp5/WY3\n8u4nb8cJnaI8KRCSbMss7LLJuihUQtzt6Gfc8W5rv61q3+1q5aE0pXsN32mkStIOJY/QyiXJfJz+\nqby46+lF4i628hDKxVLIFEHaIYlp0DzKWCZOcpRO3End4+6mudvl7VSBqkg7lDySKBMXdrkpnbib\nMcg97k618H7fnZ1kLa8MNcYybEMZcWlnF8Og5RGarF366VEJcbejl45kreTiNe988HJPF5d2NsuX\nSdhZy7qqFweVF3eNbgXej7x7rU36BYKTN1WQdt7LhyDskKTvdE/pxJ3ECGjQXuBJidWbh50QcWmn\nu3xZhO2yzo/SibsZ/bx4pJPAe5V3XpL2i4Pu8HKKcGmnt2wZhO2yDoNKiLsV3vTsOKtxaae37CDb\nH4JsXdhhUWlxd4PXxBxnMFza/W2/C9tpRenE3atoB6l1u9SdshD6CbqI0s6rlp238JN+s5uzNqUT\nN/Q+Vni/z2IXnbJul9MboTeRV0naLmynG0op7hr+sg/HyY8QTuT9yqhITeN5yj6E37iKlFrcNZIY\nvrSX2mkvze9e63XypMxN5FmLsCrCDknWoe+/aVEJcdfo5SUijlN2ytxE7tJOdjkIS9hVZ528A0ia\nbnfMZvP5jlltqnr1njQu7e7W1c/6sl4Oou3zc2NYlLLGPcjwpb0OXepN3e3x8gmTsl6kZCntLGvZ\nRathl3X/CoVCiVvS4cD3gCHAT83s/HbzuzQcJ1vyrG2HLu0iCNtlXQwK01QuaQhwEXA48Dbgo5Le\n2mm5TjtTVZrMZ8yYkXcIqZDUb3XHHXckkk8R6PUEm2XZFOHtUo373AMPPJDKerJsEh+kOXyQZnin\nPwojbmBf4Akz+4uZrQB+BXyocaZmO1CnHSurnS7Pnfuhhx7Kbd1FoCri7mcf7LZs8hyLOysxNtvG\nBx98MNH19CPCrIVdW58LOx+KJO6tgNl1/8+J09ai1c7Uy05Wxlq346RFnsdLntLutI5+BNwreQjb\nyZciidt6mbnVDtaL1Jvt2N3mGcKbgBynnpD3qSzv5WYl7V7nz2IZF3Y5kFlPPswNSe8EzjGzw+P/\nzwRW1ndQk1SMjXEcx6kIZqY08u33fJ9WPFlSJHEPBR4F3g3MA+4GPmpmj+QamOM4juNkSGEeBzOz\n1ySdCtxI9DjYRJe24ziOUzUKU+N2HMdxHKdYndPaIulwSTMlPS7p9LzjCQVJ20j6P0l/lvSQpM/m\nHVNoSBoiabqka/OOJTQkbSzpt5IekfRw3NfEIepnEx9XMyT9UtLwvGPKE0kXS3pO0oy6tE0kTZP0\nmKSbJG2cZ4xloRTi7ndwloqwAvh3M3s78E7gFC+btTgNeJgen1yoCN8HpprZW4HdAL89BUjaDvgk\nsKeZ7Up0++7YPGMKgEuIzsH1nAFMM7MdgZvj/50BKYW46XJwlipiZs+a2Z/i74uITrxvzjeqcJC0\nNXAE8FOg8L1Nk0TSRsBBZnYxRP1MzOylnMMKhZeJLorXizvOrgfMzTekfDGz24AXG5KPBH4Wf/8Z\ncCBl4PYAAARqSURBVFSmQZWUsoi768FZqkxcS9gDuCvfSILiu8DngZV5BxIg2wPzJV0i6X5JP5G0\nXt5BhYCZvQD8F/AM0VMuC83sf/ONKkjeZGbPxd+fA96UZzBloSzi9ibODkhaH/gtcFpc8648kj4A\nPG9m0/HadjOGAnsCPzSzPYFX8aZOACS9Bfg3YDuiFqz1JR2Xa1CBY1FPaD9XJ0BZxD0X2Kbu/22I\nat0OIGkYMBn4uZldnXc8AfEu4EhJTwGXA38vaVLOMYXEHGCOmd0T//9bIpE7sDfwRzNbYGavAVcS\n7U/OmjwnaQsASVsCz+ccTykoi7jvBcZK2k7SusAEYErOMQWBJAETgYfN7Ht5xxMSZnaWmW1jZtsT\ndSz6nZl9PO+4QsHMngVmS9oxTnoP8OccQwqJmcA7JY2Mj7H3EHVwdNZkCnBC/P0EwCsOCVCYAVja\n4YOztOUA4GPAg5Kmx2lnmtkNOcYUKt6MtzafAX4RXxDPAk7MOZ4gMLMH4taZe4n6R9wP/DjfqPJF\n0uXAIcBmkmYDXwbOA66QdBLwF+CY/CIsDz4Ai+M4juMUiLI0lTuO4zhOJXBxO47jOE6BcHE7juM4\nToFwcTuO4zhOgXBxO47jOE6BcHE7juM4ToFwcTuO4zhOgXBxO47jOE6BcHE7juM4ToFwcTuO4zhO\ngXBxO04OSPq+pNuapG8u6WlJPqaz4zhNcXE7TsbEL+z4JLC8yeSjiF5LuyLToBzHKQwubsfJnn2B\nEcAtTaYdSPSWsluzDMhxnOLg4nac7Dkk/ntLk2kHATPNbEF24TiOUyRc3I6TPQcDS4A76xMlbQVs\ni9e2Hcdpg4vbcTJE0lDgXcCdZtZ4H/ug+O/vs43KcZwi4eJ2nGzZC3gDzeV8YPz3VgBJw7MKynGc\n4uDidpxsqd3ffqDJtIOAJ81snqQ3Aj/OLizHcYqCi9txsuXg+O+S+kRJfw/sAtwXJ/0dMC3DuBzH\nKQhD8w7AcaqCpHVY3Rz+PuCmOH0f4FRgFrAwnn4McFLdsm8GzgX+AqxH1Nz+ZzP7SRaxO44TDi5u\nx8mOdwAbAr8B3i7pGmAZ8CTwj/H070uaClxlZq8CSNqM6NGxU8xsmqStgaeBw7PfBMdx8sbF7TjZ\nUbu//WMzu7nJ9DuB/ZqkX0BUu641nS8FVtLwOJnjONXA73E7TnYcTDSU6R+7XUDSpsBxwOS65AOB\nR8zslWTDcxynCLi4HScDJImo1/i9Zrak0/x1vBMYxpqPjx0E3J5geI7jFAgXt+Nkwy7AJvQ+uMoI\nYImZza5LOxj4o6QdJB2RVICO4xQDF7fjZMNGwLPA5T0udxewsjYYi6STgLcCjxN1TvP73I5TMWRm\necfgOE4bJJ1I1GntOSKR70A0Att0M/tenrE5jpM9Lm7HcRzHKRDeVO44juM4BcLF7TiO4zgFwsXt\nOI7jOAXCxe04juM4BcLF7TiO4zgFwsXtOI7jOAXCxe04juM4BcLF7TiO4zgFwsXtOI7jOAXi/wee\nsfj+7dKI/gAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xf370d50>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The figure shows the incomplete likelihood function that the algorithm is exploring. Note that the algorithm can get to the maximizer but since the surface has symmetric maxima, it has no way to pick between them and ultimately just picks the one that is closest to the starting point. This is because each group is equally likely to be chosen. I urge you to download this notebook and try different initial points and see where the maximizer winds up."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Expectation maximization is a powerful algorithm that is especially useful when it is difficult to de-couple the variables involved in a standard maximum likelihood estimation. Note that convergence to the \"correct\" maxima is not guaranteed, as we observed here. This is even more pronounced when there are more parameters to estimate. There is a nice [applet](http://www.cs.cmu.edu/~alad/em/) you can use to investigate this effect and a much more detailed mathematical derivation [here](http://crow.ee.washington.edu/people/bulyko/papers/em.pdf).\n",
      "\n",
      "As usual, the IPython notebook corresponding to this post can be found [here](https://github.com/unpingco/Python-for-Signal-Processing/blob/master/Expectation_Maximization.ipynb). I urge you to try these calculations on your own. Try changing the sample size and making the choice between the two groups no longer equal to 1/2 (equally likely).  \n",
      "\n",
      "Note you will need at least `sympy` version 0.7.2 to run this notebook.\n",
      "\n",
      "Comments appreciated!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    }
   ],
   "metadata": {}
  }
 ]
}