ExpectationMaximization(EM).ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "ExpectationMaximization(EM).ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "The EM algorithm is a popular tool for simplifying difficult maximum likelihood problems. We first describe it in the context of a simple mixture model\n",
        "\n",
        "<font color='blue'>**Two-Component Mixture Model**</font>\n",
        "\n",
        "We are given a small toy dataset, as plotted below "
      ],
      "metadata": {
        "id": "Rb908BqgeUjp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import random\n",
        "import math\n",
        "\n",
        "x = np.array([-0.39,0.12,0.94,1.67,1.76,2.44,3.72,4.28,4.92,5.53,0.06,0.48,1.01\n",
        "              ,1.68,1.80,3.25,4.12,4.60,5.28,6.22])\n",
        "#x = np.sort(x)\n",
        "sns.histplot(x,stat=\"probability\",bins=15)\n",
        "# Plotting the density curve"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        },
        "id": "qDJB-xgEVHM0",
        "outputId": "228f8ab5-bae7-4611-faac-ea0dd6fe7f85"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7f609ab08ad0>"
            ]
          },
          "metadata": {},
          "execution_count": 67
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAD4CAYAAAD7CAEUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYpklEQVR4nO3df7RdZX3n8feHhATG8jvprKybYOIQdQLWIJfg6EhbKBjWIKFtkGQhP1zU4Ghau5x2jNMxaKSrpbNaOs6iDqkJPxRMmFjKXdPQyBSonTJgbjASAkYvEckNdHElCCgFjHzmj/OkHg735p4dzr73XPJ5rbXX3fvZ+3nud5MsPtnP3mcf2SYiIqJdh4x3ARERMbEkOCIiopIER0REVJLgiIiIShIcERFRyeTxLmAsTJs2zbNnzx7vMiIiJpQtW7b80Pb01vaDIjhmz55Nf3//eJcRETGhSPrBcO2ZqoqIiEoSHBERUUmCIyIiKklwREREJQmOiIioJMERERGV1BockhZK2iFpQNKKYfZ/UtLDkh6U9HeS3ty071JJ3yvLpU3tp0jaVsb8giTVeQ4REfFqtQWHpEnAtcA5wDxgqaR5LYd9C+i1/UvABuBPSt9jgSuB04AFwJWSjil9vgh8BJhbloV1nUNERLxWnVccC4AB2zttvwysAxY1H2D7btsvlM37gJll/f3Anbb32H4GuBNYKGkGcKTt+9z4IpGbgPNrPIeIiGhRZ3D0ALuatgdL20guB+4YpW9PWR91TEnLJPVL6h8aGqpYenRSz6zjkdSxpWfW8eN9ShEHta545YikDwG9wC93akzbq4HVAL29vfmaw3H0xOAuLrzu3o6Nt/6K93RsrIiors4rjt3ArKbtmaXtVST9GvAHwHm2Xxql725+Pp014pgREVGfOoNjMzBX0hxJU4AlQF/zAZJOBq6jERpPNe3aBJwt6ZhyU/xsYJPtJ4HnJL27PE11CXB7jecQEREtapuqsr1X0nIaITAJWGt7u6RVQL/tPuC/Ab8A/K/yVO3jts+zvUfS52mED8Aq23vK+seAG4DDadwTuYOIiBgztd7jsL0R2NjStrJp/df203ctsHaY9n7gpA6WGRERFeST4xERUUmCIyIiKklwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEoSHBERUUmCIyIiKklwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEpqDQ5JCyXtkDQgacUw+0+X9ICkvZIWN7X/qqStTcuLks4v+26Q9P2mffPrPIeIiHi12r4BUNIk4FrgLGAQ2Cypz/bDTYc9DlwG/F5zX9t3A/PLOMcCA8DXmw75fdsb6qo9IiJGVudXxy4ABmzvBJC0DlgE/Etw2H6s7HtlP+MsBu6w/UJ9pUZERLvqnKrqAXY1bQ+WtqqWAF9taftDSQ9KukbS1AMtMCIiquvqm+OSZgDvADY1NX8aeDtwKnAs8KkR+i6T1C+pf2hoqPZaIyIOFnUGx25gVtP2zNJWxQeB22z/dF+D7Sfd8BJwPY0psdewvdp2r+3e6dOnV/y1ERExkjqDYzMwV9IcSVNoTDn1VRxjKS3TVOUqBEkCzgce6kCtERHRptqCw/ZeYDmNaaZHgFttb5e0StJ5AJJOlTQIXABcJ2n7vv6SZtO4Yvn7lqFvlrQN2AZMA66q6xwiIuK16nyqCtsbgY0tbSub1jfTmMIaru9jDHMz3fYZna0yIiKq6Oqb4xER0X0SHBERUUmCIyIiKklwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEoSHBERUUmCIyIiKklwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEpqDQ5JCyXtkDQgacUw+0+X9ICkvZIWt+z7maStZelrap8j6f4y5vryfeYRETFGagsOSZOAa4FzgHnAUknzWg57HLgMuGWYIf7Z9vyynNfUfjVwje0TgGeAyztefEREjKjOK44FwIDtnbZfBtYBi5oPsP2Y7QeBV9oZUJKAM4ANpelG4PzOlRwREaOpMzh6gF1N24OlrV2HSeqXdJ+kfeFwHPAj23tHG1PSstK/f2hoqGrttemZdTySOrr0zDp+vE8rIg4ik8e7gP14s+3dkt4C3CVpG/Bsu51trwZWA/T29rqmGit7YnAXF153b0fHXH/Fezo6XkTE/tR5xbEbmNW0PbO0tcX27vJzJ3APcDLwNHC0pH2BV2nMiIh4/eoMjs3A3PIU1BRgCdA3Sh8AJB0jaWpZnwa8F3jYtoG7gX1PYF0K3N7xyiMiYkS1BUe5D7Ec2AQ8Atxqe7ukVZLOA5B0qqRB4ALgOknbS/d/C/RL+jaNoPhj2w+XfZ8CPilpgMY9jzV1nUNERLxWrfc4bG8ENra0rWxa30xjuqm1373AO0YYcyeNJ7YiImIc5JPjERFRSYIjIiIqSXBEREQlCY6IiKgkwREREZUkOCIiopIER0REVJLgiIiIShIcERFRSYIjIiIqSXBEREQlCY6IiKgkwREREZUkOCIiopIER0REVJLgiIiIShIcERFRSVvBIekDkiqHjKSFknZIGpC0Ypj9p0t6QNJeSYub2udL+n+Stkt6UNKFTftukPR9SVvLMr9qXRERceDaDYMLge9J+hNJb2+ng6RJwLXAOcA8YKmkeS2HPQ5cBtzS0v4CcIntE4GFwJ9LOrpp/+/bnl+WrW2eQ0REdEBbwWH7Q8DJwKPADeVqYJmkI/bTbQEwYHun7ZeBdcCilnEfs/0g8EpL+3dtf6+sPwE8BUxv96QiIqI+bU8/2X4O2EAjAGYAvw48IOm3R+jSA+xq2h4sbZVIWgBMoRFa+/xhmcK6RtLUEfotk9QvqX9oaKjqr42IiBG0e49jkaTbgHuAQ4EFts8B3gn8p7qKkzQD+DLwYdv7rko+DbwdOBU4FvjUcH1tr7bda7t3+vRcrEREdMrkNo/7DeAa299obrT9gqTLR+izG5jVtD2ztLVF0pHA3wB/YPu+pt/5ZFl9SdL1wO+1O2ZERLx+7U5V/VNraEi6GsD2343QZzMwV9IcSVOAJUBfO7+sHH8bcJPtDS37ZpSfAs4HHmrzHCIiogPaDY6zhmk7Z38dbO8FlgObgEeAW21vl7RK0nkAkk6VNAhcAFwnaXvp/kHgdOCyYR67vVnSNmAbMA24qs1ziIiIDtjvVJWk/wh8DPg3kh5s2nUE8I+jDW57I7CxpW1l0/pmGlNYrf2+AnxlhDHPGO33RkREfUa7x3ELcAfwR0DzB/iet72ntqoiIqJrjRYctv2YpI+37pB0bMIjIuLg084Vx7nAFsCAmvYZeEtNdUVERJfab3DYPrf8nDM25URERLcb7eb4u/a33/YDnS0nIiK63WhTVX+6n30G8oRTRMRBZrSpql8dq0IiImJiGG2q6gzbd0n6jeH22/6resqKiIhuNdpU1S8DdwEfGGafgQRHRMRBZrSpqivLzw+PTTkREdHt2n2t+nGSvlC+5nWLpP8u6bi6i4uIiO7T7ksO1wFDwG8Ci8v6+rqKioiI7tXu93HMsP35pu2rJF1YR0EREdHd2r3i+LqkJZIOKcsHabwuPSIiDjKjPY77PD9/R9Xv8vNXnR8C/Jh8+15ExEFntKeqjhirQiIiYmJo9x4Hko4B5gKH7Wtr/TrZiIh442v3cdzfAr5B477G58rPz7bRb6GkHZIGJK0YZv/p5RHfvZIWt+y7VNL3ynJpU/spkraVMb9Qvns8IiLGSLs3xz8BnAr8oLy/6mTgR/vrIGkScC2N7yafByyVNK/lsMeBy2h870dz32OBK4HTgAXAleWKB+CLwEdoXP3MBRa2eQ4REdEB7QbHi7ZfBJA01fZ3gLeN0mcBMGB7p+2XaXwWZFHzAbYfs/0g8EpL3/cDd9reY/sZ4E5goaQZwJG277Nt4Cbg/DbPISIiOqDd4BiUdDTw18Cdkm4HfjBKnx5gV/MYpa0dI/XtKeujjilpmaR+Sf1DQ0Nt/toA6Jl1PJI6tsTBodN/b3pmHX9Q1jgRtHVz3Pavl9XPSrobOAr429qq6gDbq4HVAL29vR7nciaUJwZ3ceF193ZsvPVXvKdjY0X3mgh/byZCjRNBu1ccSHqXpN8BfgkYLNNP+7MbmNW0PbO0tWOkvrvL+oGMGRERHdDuU1UrgRuB44BpwPWS/uso3TYDcyXNkTQFWAL0tVnXJuBsSceUm+JnA5tsPwk8J+nd5WmqS4Db2xwzIiI6oN3PcVwEvLPpBvkfA1uBq0bqYHuvpOU0QmASsNb2dkmrgH7bfZJOBW4DjgE+IOlztk+0vUfS52mED8Aq23vK+seAG4DDgTvKEhERY6Td4HiCxgf/XizbU2ljisj2RmBjS9vKpvXNvHrqqfm4tcDaYdr7gZParDsiIjpstHdV/Q8a76p6Ftgu6c6yfRbwzfrLi4iIbjPaFUd/+bmFxpTSPvfUUk1ERHS90V5yeOO+9XKD+61lc4ftn9ZZWEREdKe27nFI+hUaT1U9RuMV67MkXZqXHEZEHHzavTn+p8DZtncASHor8FXglLoKi4iI7tTuBwAP3RcaALa/CxxaT0kREdHN2r3i2CLpS/z8GwAv4uc3ziMi4iDSbnB8FPg48Dtl+x+Av6ilooiI6GqjBkf5Xo1v23478Gf1lxQREd1s1Hsctn8G7JB0cL4/OCIiXqXdqapjaHxy/JvAT/Y12j6vlqoiIqJrtRscn6m1ioiImDBGe1fVYTRujJ8AbAPW2N47FoVFRER3Gu0ex41AL43QOIfGBwEjIuIgNtpU1Tzb7wCQtIa8ETci4qA32hXHv7zIMFNUEREBo19xvFPSc2VdwOFlW4BtH1lrdRER0XX2e8Vhe5LtI8tyhO3JTeujhoakhZJ2SBqQtGKY/VMlrS/775c0u7RfJGlr0/KKpPll3z1lzH37fvHATj0iIg5Euy85rKx84vxaGjfV5wFLJc1rOexy4BnbJwDXAFcD2L7Z9nzb84GLge/b3trU76J9+20/Vdc5RETEa9UWHMACYMD2TtsvA+uARS3HLKLx5BbABuBMSWo5ZmnpGxERXaDO4OgBdjVtD5a2YY8pN9+fBY5rOeZCGt/90ez6Mk31mWGCBgBJyyT1S+ofGho60HOIiIgWdQbH6ybpNOAF2w81NV9UHhF+X1kuHq6v7dW2e233Tp8+fQyqjYg4ONQZHLuBWU3bM0vbsMdImgwcBTzdtH8JLVcbtneXn88Dt9CYEouIiDFSZ3BsBuZKmiNpCo0Q6Gs5pg+4tKwvBu6ybQBJhwAfpOn+hqTJkqaV9UOBc4GHiIiIMdPuSw4rs71X0nJgEzAJWGt7u6RVQL/tPmAN8GVJA8AeGuGyz+nALts7m9qmAptKaEwC/g/wl3WdQ0REvFZtwQFgeyOwsaVtZdP6i8AFI/S9B3h3S9tPgFM6XmhERLStq2+OR0RE90lwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEoSHBERUUmCIyIiKklwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEoSHBERUUmCIyIiKqk1OCQtlLRD0oCkFcPsnyppfdl/v6TZpX22pH+WtLUs/7OpzymStpU+X5CkOs8hIiJerbbgkDQJuBY4B5gHLJU0r+Wwy4FnbJ8AXANc3bTvUdvzy/LRpvYvAh8B5pZlYV3nEBERr1XnFccCYMD2TtsvA+uARS3HLAJuLOsbgDP3dwUhaQZwpO37bBu4CTi/86VHRMRI6gyOHmBX0/ZgaRv2GNt7gWeB48q+OZK+JenvJb2v6fjBUcYEQNIySf2S+oeGhg78JGYdj6SOLdGdOv3n3DPr+K6uTxKTpxzW3X+3D5nc8XPu9ho7/WfS6b+H+0yuZdTX70ngeNtPSzoF+GtJJ1YZwPZqYDVAb2+vD7SQJwZ3ceF19x5o99dYf8V7OjZWdE63/zl3uj5o1NjN58wre2s5547qcI1d/2dS1HnFsRuY1bQ9s7QNe4ykycBRwNO2X7L9NIDtLcCjwFvL8TNHGTMiImpUZ3BsBuZKmiNpCrAE6Gs5pg+4tKwvBu6ybUnT1bi5jqS30LgJvtP2k8Bzkt6txnXnJcDtNZ5DRES0qG2qyvZeScuBTcAkYK3t7ZJWAf22+4A1wJclDQB7aIQLwOnAKkk/BV4BPmp7T9n3MeAG4HDgjrJERMQYqfUeh+2NwMaWtpVN6y8CFwzT72vA10YYsx84qbOVRkREu/LJ8YiIqCTBERERlSQ4IiKikgRHRERUkuCIiIhKEhwREVFJgiMiIipJcERERCUJjoiIqCTBERERlSQ4IiKikgRHRERUkuCIiIhKEhwREVFJgiMiIipJcERERCUJjoiIqKTW4JC0UNIOSQOSVgyzf6qk9WX//ZJml/azJG2RtK38PKOpzz1lzK1l+cU6zyEiIl6ttq+OlTQJuBY4CxgENkvqs/1w02GXA8/YPkHSEuBq4ELgh8AHbD8h6SQa31ve09TvovIVshERMcbqvOJYAAzY3mn7ZWAdsKjlmEXAjWV9A3CmJNn+lu0nSvt24HBJU2usNSIi2lRncPQAu5q2B3n1VcOrjrG9F3gWOK7lmN8EHrD9UlPb9WWa6jOSNNwvl7RMUr+k/qGhoddzHhER0aSrb45LOpHG9NUVTc0X2X4H8L6yXDxcX9urbffa7p0+fXr9xUZEHCTqDI7dwKym7ZmlbdhjJE0GjgKeLtszgduAS2w/uq+D7d3l5/PALTSmxCIiYozUGRybgbmS5kiaAiwB+lqO6QMuLeuLgbtsW9LRwN8AK2z/476DJU2WNK2sHwqcCzxU4zlERESL2oKj3LNYTuOJqEeAW21vl7RK0nnlsDXAcZIGgE8C+x7ZXQ6cAKxseex2KrBJ0oPAVhpXLH9Z1zlERMRr1fY4LoDtjcDGlraVTesvAhcM0+8q4KoRhj2lkzVGREQ1XX1zPCIiuk+CIyIiKklwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEoSHBERUUmCIyIiKklwREREJQmOiIioJMERERGVJDgiIqKSBEdERFSS4IiIiEoSHBERUUmtwSFpoaQdkgYkrRhm/1RJ68v++yXNbtr36dK+Q9L72x0zIiLqVVtwSJoEXAucA8wDlkqa13LY5cAztk8ArgGuLn3nAUuAE4GFwF9ImtTmmBERUaM6rzgWAAO2d9p+GVgHLGo5ZhFwY1nfAJwpSaV9ne2XbH8fGCjjtTNmRETUSLbrGVhaDCy0/Vtl+2LgNNvLm455qBwzWLYfBU4DPgvcZ/srpX0NcEfptt8xm8ZeBiwrm28DdnT8JKuZBvxwnGs4EKl7bKXusZW69+/Ntqe3Nk4eg188LmyvBlaPdx37SOq33TvedVSVusdW6h5bqfvA1DlVtRuY1bQ9s7QNe4ykycBRwNP76dvOmBERUaM6g2MzMFfSHElTaNzs7ms5pg+4tKwvBu5yY+6sD1hSnrqaA8wFvtnmmBERUaPapqps75W0HNgETALW2t4uaRXQb7sPWAN8WdIAsIdGEFCOuxV4GNgLfNz2zwCGG7Ouc+iwrpk2qyh1j63UPbZS9wGo7eZ4RES8MeWT4xERUUmCIyIiKklwjIGJ+JoUSWslPVU+azNhSJol6W5JD0vaLukT411TOyQdJumbkr5d6v7ceNfUrvJWh29J+t/jXUsVkh6TtE3SVkn9411PuyQdLWmDpO9IekTSvxvzGnKPo17lNSnfBc4CBmk8GbbU9sPjWtgoJJ0O/Bi4yfZJ411PuyTNAGbYfkDSEcAW4PwJ8N9bwJts/1jSocD/BT5h+75xLm1Ukj4J9AJH2j53vOtpl6THgF7bE+oDgJJuBP7B9pfK06X/yvaPxrKGXHHUb0K+JsX2N2g86Tah2H7S9gNl/XngEaBnfKsanRt+XDYPLUvX/6tO0kzgPwBfGu9aDgaSjgJOp/FEKrZfHuvQgATHWOgBdjVtDzIB/kf2RlDetnwycP/4VtKeMuWzFXgKuNP2RKj7z4H/DLwy3oUcAANfl7SlvKJoIpgDDAHXl+nBL0l601gXkeCINyRJvwB8Dfhd28+Ndz3tsP0z2/NpvBFhgaSuniKUdC7wlO0t413LAfr3tt9F423bHy/Ts91uMvAu4Iu2TwZ+Aoz5fdMER/3ympQxVu4RfA242fZfjXc9VZWph7tpfKVAN3svcF65V7AOOEPSV8a3pPbZ3l1+PgXcRmNaudsNAoNNV6MbaATJmEpw1C+vSRlD5SbzGuAR23823vW0S9J0SUeX9cNpPEzxnfGtav9sf9r2TNuzafy9vsv2h8a5rLZIelN5eIIy1XM20PVPENr+J2CXpLeVpjNpvGFjTL1h347bLUZ69co4lzUqSV8FfgWYJmkQuNL2mvGtqi3vBS4GtpX7BQD/xfbGcaypHTOAG8tTeIcAt9qeUI+3TjD/Grit8e8MJgO32P7b8S2pbb8N3Fz+IboT+PBYF5DHcSMiopJMVUVERCUJjoiIqCTBERERlSQ4IiKikgRHRERUkuCIiIhKEhwREVHJ/wd/SNc72fXJogAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "${EM \\text{ }Algorithm\\text{ } for\\text{ } Two-component\\text{ } Gaussian \\text{ }Mixture}$\n",
        "\n",
        "1. Take initial guesses for the parameters $\\hat{\\mu_1},\\hat{σ_1^2},\\hat{\\mu_2},\\hat{σ_2^2},\\hat{π}$\n",
        "\n",
        "2. $Expectation\\text{ } Step$: Compute the responsibilities\n",
        "\\begin{align*}\n",
        "\\hat{\\gamma_i} = \\frac{\\hat{π}𝛷_{\\hat{θ_2}}(y_i)}{(1-\\hat{π})𝛷_{\\hat{θ_1}}(y_i)+ \\hat{π}𝛷_{\\hat{θ_2}}(y_i)}  \\text{   }i = 1,2,\\ldots,N\n",
        "\\end{align*}\n",
        "\n",
        "3.$Maximization\\text{ } Step$: compute the weighted means and variances:\\\n",
        "$\\hat{\\mu}_1 = \\frac{∑_{i=1}^N((1-\\hat{\\gamma}_i)y_i}{∑_{i=1}^N(1-\\hat{\\gamma}_i)}$ $\\text{  }$  $\\hat{σ_1^2} = \\frac{∑_{i=1}^N(1-\\hat{\\gamma}_i)(y_i-\\hat{\\mu_1})^2}{∑_{i=1}^N(1-\\hat{γ}_i)}$ $\\hat{\\mu}_2 = \\frac{∑_{i=1}^N\\hat{\\gamma}_iy_i}{∑_{i=1}^N\\hat{\\gamma}_i}$ $\\text{  }$  $\\hat{σ_2^2} = \\frac{∑_{i=1}^N\\hat{\\gamma}_i(y_i-\\hat{\\mu_2})^2}{∑_{i=1}^N\\hat{γ}_i}$\n",
        "\n",
        "and the mixing probability $\\hat{π} = ∑_{i=1}^N \\hat{γ_i}/N$\n",
        "\n",
        "4. Iterate steps 2 and 3 until convergence."
      ],
      "metadata": {
        "id": "EjujrczFe1aw"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Univariate Gaussiam Mixture Implementaion of EM"
      ],
      "metadata": {
        "id": "NrsmiHvuP_g3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def univariateGaussian(x,mu,var):\n",
        "\n",
        "  sigma = np.sqrt(var)\n",
        "\n",
        "  u = (x - mu)/abs(sigma)\n",
        "  density = (1/(np.sqrt(2*np.pi)* abs(sigma)))* np.exp(-u*u/2)\n",
        "\n",
        "  #density = (1/(np.sqrt(2*np.pi)* abs(sigma)))* np.exp(-(np.square(x-mu)/(2*np.square(sigma))))\n",
        "\n",
        "  return density"
      ],
      "metadata": {
        "id": "XyCI78EbbGnA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "######################################################################\n",
        "# Implementation of EM algorithm,for a simple Univarite data list.\n",
        "# The Data used, is the same used to build the histogram above\n",
        "######################################################################\n",
        "class EM1:\n",
        "\n",
        "  ######################################################################\n",
        "  # Initializing the parameters \n",
        "  # Mean-1,Mean-2, Variance-1, Variance-2 & mix \n",
        "  ######################################################################\n",
        "\n",
        "  def __init__(self,x):\n",
        "    self.x = x\n",
        "    self.mean_1 = 4.12 # Random Initialization of means \n",
        "    self.mean_2 = 5.28\n",
        "    \n",
        "    self.var_1 = self.var_2 =  1/len(self.x) * np.sum(np.square(self.x-np.mean(self.x)))\n",
        "    self.mix = 0.5\n",
        "\n",
        "  ######################################################################\n",
        "  # Below function executes the Expectation Step of the EM algorithm\n",
        "  # i.e calculate the responsibility value \"gamma\"\n",
        "  ######################################################################\n",
        "\n",
        "  def expectation(self):\n",
        "\n",
        "    self.phi_1 = univariateGaussian(self.x,self.mean_1,self.var_1)\n",
        "    self.phi_2 = univariateGaussian(self.x,self.mean_2,self.var_2)\n",
        "   \n",
        "    # Responsibility of Model-2\n",
        "    self.gamma = (self.mix* self.phi_2)/(((1-self.mix)*self.phi_1) + (self.mix*self.phi_2))\n",
        "\n",
        "    return self.gamma\n",
        "\n",
        "  ######################################################################\n",
        "  # Below function executes the maximization Step of the EM algorithm\n",
        "  # Implements the closed form equations for updating the \n",
        "  # algorithm parameters. \n",
        "  ######################################################################\n",
        "\n",
        "  def maximization(self):\n",
        "\n",
        "    self.mean_1 = np.sum((1-self.gamma)*self.x)/np.sum(1-self.gamma)\n",
        "    self.mean_2 = np.sum(self.gamma*self.x)/np.sum(self.gamma)\n",
        "\n",
        "\n",
        "    self.var_1 = np.sum((1-self.gamma)*((self.x-self.mean_1)**2))/np.sum(1-self.gamma)\n",
        "    self.var_2 = np.sum(self.gamma*((self.x-self.mean_2)**2))/np.sum(self.gamma)\n",
        " \n",
        "    self.mix = np.sum(self.gamma)/len(self.x)\n",
        "\n",
        "    return self.mix,self.mean_1,self.mean_2,self.var_1,self.var_2\n",
        "\n",
        "  ######################################################################\n",
        "  # Calculates the log-LikeLihood after each iteration\n",
        "  # Basically to view the grapy of how the liklihood converges.\n",
        "  ######################################################################\n",
        "\n",
        "  def logLikeLihood(self):\n",
        "\n",
        "    # Calculating the log Likelihood\n",
        "    self.loglikehood = np.sum(((1-self.gamma)*np.log(self.phi_1)) + (self.gamma*np.log(self.phi_2)) +\\\n",
        "    ((1-self.gamma)*np.log(1-self.mix))+(self.gamma*np.log(self.mix)))\n",
        "\n",
        "    return self.loglikehood\n",
        "\n",
        "  ######################################################################\n",
        "  # EM is an iterative algorithm\n",
        "  # EMIteration facilitates the execution of the E & M steps,\n",
        "  # until convergence.\n",
        "  ######################################################################\n",
        "\n",
        "  def EMIteration(self):\n",
        "\n",
        "    Newlog = 10\n",
        "    _ = self.expectation()\n",
        "    loglikehoodlst = []\n",
        "    Oldlog = self.logLikeLihood()\n",
        "    loglikehoodlst.append(Oldlog)\n",
        "    while abs(Newlog - Oldlog) > 0.00001:\n",
        "\n",
        "      Oldlog = Newlog\n",
        "      _ = self.maximization()\n",
        "      _ = self.expectation()\n",
        "      Newlog = self.logLikeLihood()\n",
        "      loglikehoodlst.append(Newlog)\n",
        "\n",
        "    finalResponsibilityAssng =  self.expectation()\n",
        "\n",
        "    return   finalResponsibilityAssng , loglikehoodlst"
      ],
      "metadata": {
        "id": "IbFLHvhVbvI4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Execution #\n",
        "EM = EM1(x)\n",
        "gamma,loglst = EM.EMIteration()\n",
        "mix,mean_1,mean_2,var_1,var_2 = EM.maximization()"
      ],
      "metadata": {
        "id": "wmSmlUVsdRMl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "t = np.arange(0,len(loglst))\n",
        "plt.plot(t,loglst,color='green',alpha=0.5, linestyle='dashed',linewidth=2)\n",
        "plt.title(\"LogLikeLihood\")\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 281
        },
        "id": "cTwGyRe1nKf-",
        "outputId": "6c6047be-8caf-47dc-965b-da6e188d8509"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3Qe9X3n8fdX97stWbJkS76BLcDExoDMJRACOYQSoIeQk2VJN5Cy25BwkrQ9m542l26bZkO3TdOQ5vT0QnYpbNM2UG4hKeE0pCwQSExssI2NbZDBxhfZutm63x493/1jRg8PsuTbI3kezfN5nfMczTMzz8x35lgfj37PzO9n7o6IiMRTXtQFiIjI7FHIi4jEmEJeRCTGFPIiIjGmkBcRiTGFvIhIjCnkJaeYWb+ZnRVOP2Bm35hL259mn1eb2f7Z3o/MTQp5iZyZ7TGza2dwe9OGq7tXuPtb2bx9kZmkkBcRiTGFvGQlMys2s++Y2cHw9R0zK05b/vtm1hYu+y0zczNbeRLbnXI9M6s0s2fN7LsWONfMfmpm3Wa2y8xuPcm6J2+/2sz+zcz6zGyDmZ2dtu77zexXZtYT/nx/2rLFZvZkuP9WM/t02rLS8K+JI2b2OrD+ZGqT3KSQl2z1VeAyYB1wAXAJ8IcAZnY98N+Ba4GVwNWZ7MjMFgA/A150998GyoCfAv8MLARuA/7GzFafxuZvA/4EqAZagXvCfdYA/wZ8F1gAfBv4t7AWgB8A+4HFwMeBPzWzD4XL/hg4O3z9GvCp06hLcoRCXrLVfwG+7u7t7t5BEJS3h8tuBf7B3be7+yDwtQz2sxh4DvhXd//DcN5NwB53/wd3T7j7q8CjwH86je0/7u4vu3sC+CeC/7QAbgTedPd/DPfxL8BO4NfNbAlwBfAH7j7s7puB/w3cEX72VuAed+92930E/1GITKkg6gJEprEY2Jv2fm84b2LZxrRl+zLYz41AP/B3afOWAZea2dG0eQXAP57G9g+lTQ8CFeH05OMjfN8YLut2975Jy1rSPrtv0jKRKelKXrLVQYKwnbA0nAfQBjSlLVuSwX6+BzwNPGVm5eG8fcBz7j4/7VXh7ndnsJ/JJh8fBMd4IFxWY2aVUyyD4PiXTFomMiWFvGSLQjMrmXgB/wL8oZnVmVkt8EfA98N1HwbuNLPzzKwM+B9TbC8/fXtmVnScfX8e2AX8yMxKgR8DzWZ2u5kVhq/1ZnbeaW5/Kk+F+/gNMysws/8MrAZ+HDbBvAT8r3Dba4H/Nun4v2xm1WbWBHzhFPctOUQhL9niKWAo7VVC0CSzFXgNeAX4BoC7/4SgHfpZgi8zfxluYyRte1+atL3/mG7HHgyqcBfBF50/BMaA6wi+ND1I0OTy50Bx2sdOevvT7LOLoO3/i0AX8PvATe7eGa7yCWB5uP/HgT9292fCZX9C0ETzNvDvnF4zkuQI06AhMteFV9jbgOLwC04RCelKXuYkM7slvJe+muAq+0cKeJFjKeRlrvoM0A7sBsaBmfxSVCQ21FwjIhJjupIXEYmxrHoYqra21pcvXx51GSIic8qmTZs63b1uqmVZFfLLly9n48aNJ15RRERSzGzap57VXCMiEmMKeRGRGFPIi4jEmEJeRCTGZizkzeyL4ag4teF7C0fZaTWzrWZ20UztS0RETs6MhHw4yMF1wDtpsz8CrApfdwF/OxP7EhGRkzdTV/L3EvSil/747M3A//XAL4H5ZrZohvYnIiInIeP75M3sZuCAu28xs/RFjbx39Jr94by2TPcpEjdHho4wMDbA2PgYiWSCsWT4c3yM8qJyzq09F4BEMsFze54j6UmSnmTcx1PTSU9y8aKLWTIvGE/kja432HJoC47j7qmfAPl5+dx6/rtjk/9o14/oGekhvZsTD6/Zzqs9j/WNwVjhh/sP83Tr08fUP7HuR8/9KPNL5gPwwt4X2H1k95THu7B8ITesuiF1TP+4Zfrekj+w7AOsrAnGRt/evp2XD7w85XoFeQXcfsHtqfeP73ico8NHp1x3dd1qLm26NHVMT7351LT7v+W8W1LH9Pze59ndPfUx1VfUn/Yxvdb+Gre977Zp18/ESYW8mT0DNEyx6KvAVwiaak6Lmd1F0JzD0qUa4EbiYyQxwv7e/XQOdtI52En/aD9DiSGGxoYYSgxx57o7qS6tBuCZt55he8f2KbezdN7SVMi7Oy+888K0+1wxf0Uq5LuHuqfdZmFe4Xve7+3ZS+dg55TrLixfmJoeTgzz9tG3p93/2PhYarpjsIM9R/dMuV7Sk6lpd2dvz/QjGF40+u7XeX2jfdOuO/mYDvYdpGOwY8p1F1W+26gwnBg+7v7Tj6lzsHPadT2tIeNUj2l/7/5p181URh2UmdkaglHuB8NZTQSDHFxCMLDB/wsHKMbMdgFXu/u0V/ItLS2uJ15lrkl6krePvE1bfxs1pTWsrlsNwDs973D/q/dP+7lPX/RpGqsageAKcVfnLgrzCynIK6AwL/yZX8iC0gVcsfSK1L5+/s7PybM88i2fPMt7z2vZ/GXUlNYAQSAd6j9EngWtsoYx8dd2nuXRvKA5Vcveo3sZS46ll4cRrDuvZB61ZbUADI0N0dbfdsw6ExqrGinKDwbJ6hjoYGBsYMpjL84vTgVt0pO80/POlOsB1JbVUlEUDI3bO9JL91D3lOsZxrL5746oeKD3wDHHNKGquCp1noYTwxzqP3TMsUxYVLkodUydg50MjE5zTAXFNFQ0pI5pX8/0Qw8vKFvwnmPqGe5J/ed8Osxsk7u3TLlsJnuhNLM9QIu7d5rZjQTDqt0AXAp8190vOd7nFfIylySSCbYc2sJL+16ia6gLgHNrz0392T04NshD2x6itqyWuvI6KosqKS0spaywjNKCUiqKKsjPy4/yECQmjhfys9l3zVMEAd9KcKV/5yzuS+SMGU4M86sDv2LDgQ30j/YDML9kPucsOIel895tciwrLOPOC/XPXqI1oyHv7svTph343ExuXyQbbD28lZ+9/TMAFlUs4oqlV7C6bnWqWUQkm2RVL5Qi2WhsfIy3jrzFObXnAHBhw4XsObqHixddzFnVZzHprjKRrKJLD5ETeLr1aR7f+TgjiREACvMLufX8Wzm75mwFvGQ9hbzIcbze8Tqb2jaRSCboHemNuhyRU6aQF5lGz3APT+56EoDrzr6OuvIpB94RyWoKeZEpJD3JYzseYzgxTPOCZtYvXh91SSKnRSEvMoUX9r7A3p69VBZVcvM5N6vtXeYshbzIJL0jvTy/93kM45bzbqG8qDzqkkROm26hFJmkqriKOy64g/29+zmr+qyoyxHJiEJeZArL5i97Tz8oInOVmmtEQtvat9Ha3Rp1GSIzSiEvQtC74A93/pDvb/0+B/sORl2OyIxRyEvOG0+O8+jrjzKWHGNt/VoWVWgAM4kPhbzkvK2Ht9LW30Z1STU3rrpRt0tKrCjkJee90vYKAFctu4riguKIqxGZWQp5yWntA+3s691HcX4x5y88P+pyRGacQl5y2sRV/Nr6takh3kTiRPfJS067atlVzC+Zr4eeJLYU8pLTygrLuKzpsqjLEJk1aq6RnJX0ZNQliMw6hbzkpK7BLu79xb08v/f5qEsRmVUKeclJm9o20Tfax5GhI1GXIjKrFPKSc8aT42w5tAWAixdfHHE1IrNLIS85Z2fnTgbGBqgvr6exsjHqckRmlUJecs7EvfEXLbpIXRhI7CnkJaccGTrC7iO7KcgrYG392qjLEZl1CnnJKa93vA7A6rrVlBaWRlyNyOzTw1CSU96/5P0srlyscVslZyjkJaeYGSuqV0RdhsgZo+YayRlDY0NRlyByxinkJSf0jvTyrZe+xUPbHsLdoy5H5IxRyEtOeLXtVcZ9nDzL022TklMU8hJ7SU+m7o3XE66SaxTyEntvHXmLnpEeqkuqWTFfX7pKbpmRkDezL5qZm1lt+P5qM+sxs83h649mYj8ip2Pzoc0ArGtYp6YayTkZ30JpZkuA64B3Ji16wd1vynT7IpkYTgyzs3MnABc0XBBxNSJn3kxcyd8L/D6gWxYk6+zv3U/Skyyfv5z5JfOjLkfkjMvoSt7MbgYOuPuWKf4MvtzMtgAHgd9z9+3TbOMu4C6ApUuXZlKOyDFW1qzki5d/kcGxwahLEYnECUPezJ4BGqZY9FXgKwRNNZO9Aixz934zuwF4Alg11fbd/T7gPoCWlhb9NSAzrryoXN0YSM46Yci7+7VTzTezNcAKYOIqvgl4xcwucfdDaZ9/ysz+xsxq3b1zhuoWOaGjw0epKq4iz3QTmeSu026ucffXgIUT781sD9Di7p1m1gAcdnc3s0sI2v67Mi1W5GS5Ow9sfoCkJ7lz3Z1Ul1ZHXZJIJGarg7KPA3ebWQIYAm5zPUsuZ9Denr0cHT7KvOJ5+sJVctqMhby7L0+b/mvgr2dq2yKnamIM17X1a3VvvOQ0NVZK7IyOj7K9I7iZS/fGS65TyEvs7Ozcyej4KE1VTdSW1UZdjkikFPISO+ndGIjkOoW8xMrY+BjdQ90U5BVwft35UZcjEjkN/yexUphfyO9c+ju0D7RroG4RdCUvMWRm1FfUR12GSFZQyEts9I700j/aH3UZIllFIS+x8cLeF/j2L76dGgVKRBTyEhOJZIJt7dtIepLGysaoyxHJGgp5iYU3ut5gKDFEQ0WD2uNF0ijkJRZ0b7zI1BTyMuf1j/bT2t1KnuWxZuGaqMsRySoKeZnzNh7cSNKTNC9o1uAgIpMo5GXOGxgdwDAub7o86lJEso6eeJU578bmG7li6RXMK54XdSkiWUchL7GggUFEpqbmGpmz3ux6k9buVjTomMj0FPIyJyU9ydOtT/P9rd9nV9euqMsRyVoKeZmTdnbupGuoi/kl82le0Bx1OSJZSyEvc4678/N3fg7A+5e8nzzTP2OR6ei3Q+act468xcG+g5QXlnNhw4VRlyOS1RTyMudMXMVf1nQZhfmFEVcjkt0U8jKnHOg9wNtH36Y4v5j1jeujLkck6+k+eZlTakpruGb5NeTn5VNSUBJ1OSJZTyEvc0ppYSkfXP7BqMsQmTPUXCNzhh56Ejl1CnmZE3qGe/irDX/Fhv0boi5FZE5RyMuc8NK+lzg6fJT9vfujLkVkTlHIS9YbGB1IDc595dIrI65GZG5RyEvWe7r1acaSYzQvaNb4rSKnSCEvWW1Hxw5ea3+NwrxCrl95fdTliMw5CnnJWgOjA/z4jR8D8OGzP0xNaU3EFYnMPRmFvJl9zcwOmNnm8HVD2rIvm1mrme0ys1/LvFTJNSPjI1QVV7Fi/grWL9bTrSKnYyYehrrX3b+VPsPMVgO3AecDi4FnzKzZ3cdnYH+SI2pKa/iti36L4cQwZhZ1OSJz0mw119wM/MDdR9z9baAVuGSW9iUxk0gmUtP5efmUF5VHWI3I3DYTIf95M9tqZvebWXU4rxHYl7bO/nDeMczsLjPbaGYbOzo6ZqAcmcvcnUdef4R/3f6vDI4NRl2OyJx3wpA3s2fMbNsUr5uBvwXOBtYBbcBfnmoB7n6fu7e4e0tdXd0pH4DEy2vtr7Gzcyet3a2MjY9FXY7InHfCNnl3v/ZkNmRm3wN+HL49ACxJW9wUzhOZVt9IH0+9+RQA16+8nnkl8yKuSGTuy/TumkVpb28BtoXTTwK3mVmxma0AVgEvZ7IviTd350dv/IjhxDCralaxrmFd1CWJxEKmd9d808zWAQ7sAT4D4O7bzexh4HUgAXxOd9bI8Ww5vIU3ut6gpKCEXz/n13U3jcgMySjk3f324yy7B7gnk+1Lbugf7ecnb/4EgI+s/AhVxVURVyQSHxo0RCJXXljO2vq19I32sbZ+bdTliMSKQl4iZ2bcsOoGEsmEmmlEZpj6rpFI9Az38NC2h1L3wpsZhfmFEVclEj+6kpczrneklwe3PEj3UDfFBcV89NyPRl2SSGzpSl7OqL6RPh7cHAT84srF6j5YZJYp5OWM6R/t58EtD9I11EVDRQO3r72dkoKSqMsSiTWFvJwRA6MDPLj5QToHO6kvr+eOC+6gtLA06rJEYk8hL2fEK22v0DHYwcLyhdxxwR2UFZZFXZJITtAXr3JGXLn0SsyMdQ3r1HWwyBmkK3mZFcOJYX6484d0D3UDwS2SVy69koqiiogrE8ktupKXGbfn6B4e3/E4PSM99Iz0cMcFd0RdkkjOUsjLjEkkEzz79rO8tO8lHKexspEbV90YdVkiOU0hLzOiY6CDR3c8yqH+QxjGB5d9kKuWXUV+Xn7UpYnkNIW8ZGw8Oc4Dmx9gYGyA6pJqPnbex1gyb8mJPygis04hL6elf7SfkoISCvIKyM/Lp2VxC70jvVy/8nqKC4qjLk9EQgp5OSUDowO8uO9FfnXgV1x39nWsb1wPwNXLr1YPkiJZSCEvJyU93MeSwQDbhwcOp5Yr4EWyk0Jejqutr41f7P8FOzp2pML9nAXn8MHlH2Rx5eKIqxORE1HIyzESyQQFecE/jaPDR9l6eCsAzQuauXr51Qp3kTlEIS8ADI4Nsq19G5sPbaa2rJaPnfcxAFYtWMWHVnyI9y18HzWlNRFXKSKnSiGfo9ydtv42Wrtb2d29m329+0h6Egj6fB9PjpOfl09BXgFXLbsq4mpF5HQp5HOIu6e+IH35wMv8pPUnqWV5lseqmlVc0HAB5yw4Rw8xicSEQj6m3J2uoS4O9B5gf+9+9vXu45wF53DNimsAWFG9gnnF81hZs5KVNStZUb1CA3iIxJBCPmY2HdzEjs4dHOg9wFBi6D3LivKLuIYg5OvK6vjdy35Xtz6KxJxCfo4ZHR+lY6CD9oH21OuGVTewoGwBAAf7DtLa3QpAZVElTVVNNFY10lTVRFNVU2o7CneR3KCQz0LjyXGGE8OpwTUGxwZ5YucTdAx0cGT4yDHrHx44nAr5ixZdxNk1Z9NY2UhVcZXCXCTHKeQj1D7QzqH+Q/QMB/2uHxk6QvdQNz0jPSypWsKdF94JQHF+MW92vYnj5Fs+tWW11JXXsbB8IQvLF7Kk6t3OwBqrGmmkMapDEpEso5CfYYlkgvaBdvpH+1OvvpE++kf76R3p5abmm1hUuQgI2s83HNhwzDYMI5FMpN7n5+XziTWfoLqkmprSGt35IiInTSE/jaQnMSzV3HGw7yCdg50Mjg0yODbIwOhAarqhooGPrPoIEPTxct+m+6bd7pHhI6mQb6pqYmBsgKriKuYVz6O6NAjx+SXzU0+cTmhe0DxLRyoicZYTIe/uHOo/RO9IL/NK5tFQ0QBA12AXL+57kZHECMOJYYYSQwyNDTGUGGIkMcJvX/rbVJdWA/DiOy+yvWP71NvHU9MVRRU0VDRQUVRBZVElFUUVqVdVcRV15XWpddfUr2FN/ZpZPHIRyXU5EfLbO7bzyOuPAHBJ4yXcsOoGAIYSQ7zS9sq0nxtODKeml85biplRVlhGeWE5ZYVlwXRROZVFlan18vPy+WzLZ2fpSERETk1OhPyh/kMALCxfSF3Zu1fSNaU13NR8EyUFJZQUlFBaUBr8LAx+5lleat1Lmy7lUi4947WLiGQio5A3s68BnwY6wllfcfenzGw5sAPYFc7/pbtHdnnbN9IHwOVNl3PhogtT88sKy2hZ3BJVWSIis24mruTvdfdvTTF/t7uvm4HtZ6xvNAj5yuLKE6wpIhIveSdeZe6buJJPbzsXEckFMxHynzezrWZ2v5lVp81fYWavmtlzZvaB6T5sZneZ2UYz29jR0THdahnRlbyI5Cpz9+OvYPYM0DDFoq8CvwQ6AQf+J7DI3f+rmRUDFe7eZWYXA08A57t77/H21dLS4hs3bjyNwzi+0fFR+kf7qS6p1mP+IhI7ZrbJ3af8gvGEbfLufu1J7uR7wI/Dz4wAI+H0JjPbDTQDM5/gJ6Eov0ijGolITsqoucbMFqW9vQXYFs6vM7P8cPosYBXwVib7EhGRU5fp3TXfNLN1BM01e4DPhPOvAr5uZmNAEvisu3dnuK/T8mbXm7y470XOrT2Xy5oui6IEEZHIZBTy7n77NPMfBR7NZNszpXOwkz1H91BfXh91KSIiZ1zsb6HUnTUiksviH/LhPfIVRRURVyIicubFP+RH9SCUiOSu+If8iJprRCR3xT/kdSUvIjks1l0NJz3J2vq1DIwOUFJQEnU5IiJnXKxDPs/yuKn5pqjLEBGJTOyba0REclmsQ75vpI+2vjaGxoaiLkVEJBKxDvnX2l/j7zf9Pc/tfS7qUkREIhHrkNdgISKS6+Id8urSQERyXLxDXlfyIpLj4h3yupIXkRwX25B3d13Ji0jOi23Ij4yPMJYcoyi/iKL8oqjLERGJRGyfeC3KL+ILl3yBocSQBu8WkZwV25DPszwWlC2IugwRkUjFtrlGRERiHPLb2rfx2I7HeKPrjahLERGJTGxDfn/vfrYe3krXYFfUpYiIRCa2Ia8RoURE4hzyGhFKRCTGIR9eyVcUVURciYhIdGIZ8u6uLg1ERIhpyA8nhkkkExTnF+tpVxHJabF8GCqRTNC8oJl8y4+6FBGRSMUy5CuLK/mNNb8RdRkiIpGLZXONiIgEYnklPzA6AEBZYZk6JxORnBbLK/nn9j7HX7z0F2w4sCHqUkREIhXLkNdgISIigYxD3sy+YGY7zWy7mX0zbf6XzazVzHaZ2a9lup9ToXvkRUQCGbXJm9k1wM3ABe4+YmYLw/mrgduA84HFwDNm1uzu45kWfDJ0JS8iEsj0Sv5u4M/cfQTA3dvD+TcDP3D3EXd/G2gFLslwXyfF3ekf7Qd0JS8ikmnINwMfMLMNZvacma0P5zcC+9LW2x/OO4aZ3WVmG81sY0dHR4blwODYIOM+TmlBKQV5sbx5SETkpJ0wBc3sGaBhikVfDT9fA1wGrAceNrOzTqUAd78PuA+gpaXFT+WzU1F7vIjIu04Y8u5+7XTLzOxu4DF3d+BlM0sCtcABYEnaqk3hvFlXU1rDb677TYKSRERyW6bNNU8A1wCYWTNQBHQCTwK3mVmxma0AVgEvZ7ivk1KUX8Ty+ctZUb3iTOxORCSrZdpofT9wv5ltA0aBT4VX9dvN7GHgdSABfO5M3VkjIiLvyijk3X0U+OQ0y+4B7slk+6dj86HNdAx0sKZ+DQ0VU32VICKSO2L3xOuOjh28uO9FjgwdiboUEZHIxS7kdXeNiMi74hfyetpVRCQlViGf9CQDY0E3wxrAW0QkZiE/ODZI0pOUF5aTn6eh/0REYhXyE001uooXEQnEqnOXpCdZVLGI2rLaqEsREckKsQr5xqpGPtPymajLEBHJGrFqrhERkfeKVciPjo+qYzIRkTSxCvlHXn+Ebzz/DVq7W6MuRUQkK8Qq5PtG+lIDhoiISNxCXl0aiIi8R2xCPulJBkYHMIzywvKoyxERyQqxCfn+0X4cp7xIT7uKiEyITcirYzIRkWPFJ+TVHi8icozYPPG6uHIxHzvvY5QVlkVdiohI1ohNyFcVV7G2fm3UZYiIZJXYNNeIiMixYnMl/2rbq4yOj7K6brXa5UVEQrEJ+V/u/yWHBw6zdN5ShbyISCg2zTW6u0ZE5FixCPlEMsHg2CB5lqe7a0RE0sQi5PtH+4Fg2L88i8UhiYjMiFgkYnrIi4jIu2IR8urSQERkarEI+UQyQXlhOVXFVVGXIiKSVWJxC+Wa+jWsqV+jof9ERCaJxZX8BDOLugQRkawSq5AXEZH3UsiLiMRYxiFvZl8ws51mtt3MvhnOW25mQ2a2OXz9XealiojIqcroi1czuwa4GbjA3UfMbGHa4t3uvi6j6kREJCOZXsnfDfyZu48AuHt75iWJiMhMyTTkm4EPmNkGM3vOzNanLVthZq+G8z8w3QbM7C4z22hmGzs6OjIsR0RE0p2wucbMngEaplj01fDzNcBlwHrgYTM7C2gDlrp7l5ldDDxhZue7e+/kjbj7fcB9AC0tLbrRXURkBp0w5N392umWmdndwGMePIX0spklgVp37wAmmnA2mdlugqv+jTNTtoiInIxMn3h9ArgGeNbMmoEioNPM6oBudx8Pr+xXAW+daGObNm3qNLO9GdRTC3Rm8PnZpvoyo/oyo/oyk831LZtuQaYhfz9wv5ltA0aBT7m7m9lVwNfNbAxIAp919+4Tbczd6zIpxsw2untLJtuYTaovM6ovM6ovM9le33QyCnl3HwU+OcX8R4FHM9m2iIhkTk+8iojEWNxC/r6oCzgB1ZcZ1ZcZ1ZeZbK9vSqbueUVE4ituV/IiIpJGIS8iEmOxCHkzu97MdplZq5l9Kep6JjOzPWb2WtgjZ+QPhJnZ/WbWHt76OjGvxsx+amZvhj+rs6y+r5nZgbSeTW+IsL4lZvasmb0e9r76O+H8rDiHx6kvK86hmZWY2ctmtiWs70/C+SvCLlJazewhMyvKsvoeMLO3087f3OiA0d3n9AvIB3YDZxE8jLUFWB11XZNq3EPwJHDktYT1XAVcBGxLm/dN4Evh9JeAP8+y+r4G/F7U5y6sZRFwUThdCbwBrM6Wc3ic+rLiHAIGVITThcAGgq5RHgZuC+f/HXB3ltX3APDxqM/fqb7icCV/CdDq7m95cN/+Dwi6P5ZpuPvzwOSH024GHgynHwQ+ekaLSjNNfVnD3dvc/ZVwug/YATSSJefwOPVlBQ/0h28Lw5cDHwIeCedHef6mq29OikPINwL70t7vJ4v+QYcc+Hcz22Rmd0VdzDTq3b0tnD4E1EdZzDQ+b2Zbw+acyJqT0pnZcuBCgqu9rDuHk+qDLDmHZpZvZpuBduCnBH+NH3X3RLhKpL/Hk+tz94nzd094/u41s+Ko6jsVcQj5ueBKd78I+AjwubDbh6zlwd+p2Xbl8rfA2cA6gl5O/zLacsDMKgie7P5dn9TDajacwynqy5pz6O7jHgwq1ETw1/i5UdUylcn1mdn7gC8T1LmeoPfdP4iwxJMWh5A/ACxJe98Uzssa7n4g/NkOPE7wjzrbHDazRQDhz6waAMbdD4e/eEnge0R8Ds2skCBA/8ndHwtnZ805nKq+bDuHYU1HgWeBy4H5ZjbR1UpW/B6n1Xd92AzmHgyS9A9kwfk7GXEI+V8Bq8Jv5ouA24AnI64pxRBHW2MAAAEhSURBVMzKzaxyYhq4Dth2/E9F4kngU+H0p4AfRljLMSbCM3QLEZ5DMzPg/wA73P3baYuy4hxOV1+2nEMzqzOz+eF0KfBhgu8NngU+Hq4W5fmbqr6daf+BG8H3Bdn4e3yMWDzxGt4K9h2CO23ud/d7Ii4pJexq+fHwbQHwz1HXZ2b/AlxN0HXqYeCPCbqNfhhYCuwFbvWT6Dn0DNZ3NUEzgxPcrfSZtPbvM13flcALwGsEvawCfIWg3Tvyc3ic+j5BFpxDM1tL8MVqPsGF5sPu/vXwd+UHBE0hrwKfDK+as6W+/wDqCO6+2UzQu27/9FvKDrEIeRERmVocmmtERGQaCnkRkRhTyIuIxJhCXkQkxhTyIiIxppAXEYkxhbyISIz9f1sJ1ozuJXxuAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "density1 = []\n",
        "for i in range(0,100):\n",
        "  val = np.random.normal(mean_1,var_1,1)\n",
        "  density1.append(val)\n",
        "\n",
        "density2 = []\n",
        "for i in range(0,100):\n",
        "  val = np.random.normal(mean_2,var_2,1)\n",
        "  density2.append(val)\n",
        "\n",
        "density1 = np.array(density1).ravel()\n",
        "density2 = np.array(density2).ravel()\n"
      ],
      "metadata": {
        "id": "kYJ2xNDE0H4Q"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import matplotlib.mlab as mlab\n",
        "from scipy.stats import norm\n",
        "arr = np.random.randn(100)\n",
        "\n",
        "plt.figure(figsize=(8,4))\n",
        "plt.hist(density1,bins=30, density=True,color='skyblue',alpha=0.5)\n",
        "plt.hist(density2,bins=30, density=True,color='skyblue',alpha=0.5)\n",
        "\n",
        "\n",
        "\n",
        "sigma1 = np.sqrt(var_1)\n",
        "sigma2 = np.sqrt(var_2)\n",
        "x1 = np.linspace(min(density1), max(density1), 50)\n",
        "x2 = np.linspace(min(density2), max(density2), 50)\n",
        "sns.histplot(x,stat=\"probability\",bins=20,edgecolor='white',label='Original Unlabelled Data')\n",
        "plt.plot(x1, norm.pdf(x1, mean_1, sigma1),color='blue',linewidth=2,label='Gaussian One')\n",
        "plt.plot(x2, norm.pdf(x2, mean_2, sigma2),color='blue',linewidth=2,label='Gaussian Two')\n",
        "plt.title(\"To Illustrate the accuracy of EM Model parameters to given data\")\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 281
        },
        "id": "UER_9Ftu7V4S",
        "outputId": "d5d4eabe-9fb9-4934-feac-e4b0df6ec6f3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 576x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAEICAYAAABPr82sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZxN9f/A8dd77p251xaSJHtSsg0hFUopSUohWUNKttAmS2VPi1KJRIulRZZCffm1SdayJUW2JFQKWRuzf35/fM6Ma5qZe4e5c+fOvJ+Pxzzm3rPd9zn33PM+53M+5/MRYwxKKaWUCj8RoQ5AKaWUUmdHk7hSSikVpjSJK6WUUmFKk7hSSikVpjSJK6WUUmFKk7hSSikVpjSJh5CI7BGRm5zXI0Tk3VDHFAwiUlFEjIi4Qx1LXiEil4vIJhE5ISL9Qx3PuRCRbiKyMsBpp4vImGDHpM6eiCwRka4hjiHfHHM0iQdIRE76/CWLyCmf950CXEbAB6uzjDFbdlzfk4tQzK8CMgj42hhTxBjzatqRIrJMRGLT7LefOOOaOPvJx2nmiXaGL0vvA332r+/TDL9AROJFZE+2rV0+k5MnJ873vz9YyzfG3GqMmRGs5We3YG+PYNMkHiBjTOGUP2AvcLvPsPdCHV+g8sOZaTDlou1XAdjiZ5p+vvutMeZ2n3EHgWtEpITPsK7AjgA+u6CI1PB53xH4NaCow0Qu+p4DEm7xquyjSfwciYhHRF4WkT+cv5dFxHOOy/zPmWGaoverRGS9iBwXkb9E5CVnsuXO/6POldc1ztX/KhGZICKHgREiUllElorIYRE5JCLviUgxZ9mzgPLAJ84yBjnDrxaR1SJyVER+EJEmGcSe7vyOTiKy1/nMYT7zRIjIYBH5xYlpjoicn8Hyi4vIpyJyUESOOK/L+ow/X0Tecb6LIyKywGdcK6cI+rjzWc3TblvnfeqtDZ+rzx4ishdY6gyfKyIHROSYiCwXkeo+8xcQkRdF5Ddn/Epn2P9E5KE067NZRO7KYF3vEJEtzjZfJiJXOMOXAjcArznb+LL05vcjHlgAtHeW6QLuAQI5IZ2FTfgp7gVmpon9Cifmo8463OEzroSILHK+h7VA5TTzVhWRL0TkHxHZLiLtAlkhn339NWe7bxORpj7ju4vIz2JvQewWkQd9xjURkf0i8oSIHADeCWBfWyYiY5zfxUkR+cRZt/ecdVsnIhX9rZeI9AQ6AYPkzBKTi0VkvvP5v4rPbRNnH50nIu+KyHGgm2R8XPDdRoWAJcDFcrqE5mLJwnFMRFzO/n3Iiauf+JQAOtvlfmeZR8XnhE9ESootxbzQed9S7G/yqLMda/lMu0dEHnN+I8dE5EMR8WYS03gnpt3AbWnGp/vdZ7I9rhKRNU5cfzr7VFR6nx1yxhj9y+IfsAe4yXk9CvgWuBAoCawGRmcwXzdgZQbLGQG867xuAuzP5DPXAF2c14WBq53XFQEDuNN8ZiLwEOAGCgCXAjcDHifm5cDL6X2W874McBhogT3xu9l5X9Lf9kkT1zTn86OBOOAKZ/wAZxuWdWJ6A/ggg2WXANoABYEiwFxggc/4/wEfAsWBSOB6Z/hVwDEn9ghnnapmEK/vd5ES+0ygEFDAGX6f8/ke4GVgk8/8k4Blzme4gGud6doB3/lMF+1sx6h01vMy4F8n3khs8fmulGmd5d+fyT6a4Xic/cuJ6ztnWAvgM+B+YFkG86Vsi4rAPmfdqgHbgJuAPc50kU6sQ4Eo4EbgBHC5M342MMfZnjWA33F+F86wfUB37P5aBzgEVHPGTwfGZPL7SgQedmK4x/nOz3fG34Y9YRDgeiAGuNJnmyQCzznfVQH872vLnPWsDBQFtmJLMm5yYp8JvHM264XdRzcATzvb8BJgN3CLzz6aANzpTFuADI4LGX3/aYZl5TjWy1nXstjf2Zf4HHfw2feAt4GxPvP2Bf7PeV0H+BtogN2XumJ/ix6f3+Va4GLgfOBnoFcmMW0DyjnTfp0mJn/ffdrtURe42vmuKjqfPTCQ/JDTfyEPIBz/ODOh/gK08Bl3C87BLJ35upE9SXw5MBK4IM00FUk/ie/1sz53At+n91nO+yeAWWnm+Qzo6m/7pImrrM+wtUB75/XPQFOfcaWxByh3ZnE709YGjvjMlwwUT2e6N4AJAcbr+12kxH5JJjEUc6Ypij2gngKi05nOCxwBqjjvxwOTM1jmU8Acn/cR2GTXxHm/DP9JPAY46vM3Ou3+BewELscm1k4ElsTd2AP3LcCzwDDOTOKNgQNAhM+8Hzjb1eV8t1V9xj3D6SR+D7Aine9uuPN6Opkn8T8ASbOfdclg+gXAAJ9tEg94A9nXfLbxMJ/3LwJLfN7fjnNyl9X1wia2vWmmH8Lpk4IRwPI049M9LqSzHqnfv8+wrBzHlgIP+ry/iYyT+E3ALz7TrgLudV6/TpoTBWA7p0+89wCdfcY9D0zJJKZePu+bkeZY6Oe735/edD7TDwQ+zmyaUP1pcfq5uxj4zef9b86wYOqBvVLb5hTZtfQz/T7fNyJSSkRmi8jvTlHcu8AFmcxfAbjbKVo6KiJHgUbYpJkVB3xex2CvFlKW/7HPsn8GkoBSaRcgIgVF5A2xRdXHsQeuYmKLg8sB/xhjjqTz2eWwB6qzlboNnaK7Z8UWyR/HHmzAbsMLsMn6P59ljInFlhJ0FpEIoAO2aDo9Z+xXxphkJ4YyWYi5vzGmmM/fU+lMMwvohy2e/zid8RmZiU2a6a3DxcA+J+YUv2FjL4k9CdiXZlyKCkCDNPtaJ+CiAOP63ThHXZ9lXwwgIreKyLdOcfZRbOmD735/0PmOcKbPbF9L8ZfP61PpvPfdx7OyXhWwRby+0w/lzN/EvjTzZPW44Csrx7GL03x22jh8fY2tQ9HAubVQm9P7WQXg0TTrWC7N52Z0zPAXk++6BPLdk2b6y8TePjngfPfPZDZ9KGkSP3d/YHfGFOWdYefiX2wRHpB6v7JkyntjzE5jTAds0ddzwDzn3o5Ju6CUWdK8f8YZVtMYcx7QGVvMlNH0+7BX4r4JoZAx5tkAP8+ffcCtaZbvNcb8ns60j2KvHBs4sV/nDBdnOeeLc38/nc+onM5wSLO9Sf/A6rtOHYFW2KuMotgr1JQYDgGxmXzWDOzBuykQY4xZk8F0Z+xXIiLYA1x62+RczAL6AIuNMTFZmG8+tohytzFmb5pxfwDlnBOVFOWxsR/EFluXSzMuxT7gmzT7QmFjTO8A4yrjbCvfZf/h3N+djy39KGWMKQYsJvP9PrN9Lav8rVd6v7lf00xfxBjTIqN4MzkupJXe7zMrx7E/sUXpKcplMB3GmCTsrZMOzt+nxpgTzuh92KJ233UsaIz5IKPlZeJPMtinAvju09ser2OL56s43/1Qzu57DzpN4ufuA+BJp8LGBdh7WOf6vPcOwCsit4lIJPAk9j4dACLSWURKOlc6R53BydgDZDL2/llmigAngWMiUgZ4PM34v9Is413gdhG5xbkK9YqtCFSW9KWd358pwFgRqQCplV9aZRL7KWzlvfOB4SkjjDF/YiupTBZbKSlSRFIOvG8B3UWkqdiKdGVEpKozbhPQ3pm+HtDWT7xFsPf0D2OT/zM+MSRj7wO+5FSQcYmtYOhxxq/BfkcvkvFVONgD321OvJHYhBKHvVeZbYwxv2LvEQ7zN22a+f7F3uu+P53R32GvmgY527QJtmh5tnNQ/whbwbKgiFTjzEpynwKXiUgXZ95IEakvTqW+AFwI9Hfmuxu4AnvAjsL+hg4CiSJyK7bINTMZ7mtnwd96pf3NrAVOiK1oV8DZj2qISP2MPiCT40JafwElRKSoz7CsHMfmAAOc31Ax7O22zLyPvZ3QyXmdYhrQy7lKFxEp5BzzivhZXkYx9ReRsiJSHBjsM87fd5/e9igCHAdOOseJQE8ic5wm8XM3BlgPbAZ+BDY6w86aMeYY9uroTezVy7/YikgpmgNbROQk8Ar23vIp50pqLLDKKZ66OoOPGAlcia308z/sQdXXOOwP+qiIPGaM2Ye98hyK/SHswyb+jPafM+YPYJVfARYBn4vICWwFmwYZTPsythLPIWe6/0szvgv2nus2bKWZgQDGmLXYSkUTnPX+htNXHk9hr5yPYLfN+2RuJra47ndsBZ9v04x/DLsvrAP+wV4VRaSZvyaZnOwZY7ZjS0gmOut6O/axxng/sflKqb2e8rchg89aaYzJcumRMWa9MSa92wbxTry3OrFPxt4H3eZM0g9bLHoAey/4HZ95T2APsO2xV4IHOF3ZLBDfAVWczx0LtDXGHHaW2x97sD+CLU1Z5GdZ/va1gAWwXm8B1ZzfzALnZKcltvj5VyeGN7ElPxlJ97iQTizbsEl7t/N5F5O149g04HNn2u+xJ0mJ2Ftg6a37d9hj2MXYk+yU4euBB4DXsN/JLuwtmrMxDVtP5wcn9tRjmr/vPoPt8Zgz3Qln2R+eZVxBJ2fePlJKBZuI3Av0NMY0CnUseYmIdMNWqNLtmoOcK9spxpgKfidW2U6vxJXKQSJSEFvKMjXUsSh1Npzi/RYi4nZuxw0na5UiVTbSJK5UDhGRW7C3I/7Cf5G9UrmVYG87HcEWp/+MvYeuQkCL05VSSqkwpVfiSimlVJgKu0bzL7jgAlOxYsVQh6GUUkrliA0bNhwyxpRMb1zYJfGKFSuyfv36UIehlFJK5QgR+S2jcVqcrpRSSoUpTeJKKaVUmNIkrpRSSoWpsLsnrpRSOSEhIYH9+/cTGxvrf2KlsoHX66Vs2bJERkYGPI8mcaWUSsf+/fspUqQIFStW5MyO0ZTKfsYYDh8+zP79+6lUqVLA82lxulJKpSM2NpYSJUpoAlc5QkQoUaJElkt+NIkrpVQGNIGrnHQ2+5smcaWUUipMaRJXSqlc6q+//qJjx45ccskl1K1bl2uuuYaPPw5+h2Hr16+nf//+2bKs+Ph4Bg4cyKWXXkqVKlVo1aoV+/fvz5ZlK63YFnZW/Pnvf4Y1Ll0oBJEopYLJGMOdd95J165def992+ndb7/9xqJFi4L+2fXq1aNevXrZsqyhQ4dy4sQJtm/fjsvl4p133qF169Z89913ersiG+iVuFJK5UJLly4lKiqKXr16pQ6rUKECDz30EAB79uyhcePGXHnllVx55ZWsXr0agGXLltGyZcvUefr168f06dMBGDx4MNWqVaNWrVo89thjAMydO5caNWoQHR3Ndddd959lrF27lmuuuYY6depw7bXXsn37dgCmT59O69atad68OVWqVGHQoEH/WYeYmBjeeecdJkyYgMvlAqB79+54PB6WLl3Knj17uOKKK3jggQeoXr06zZo149SpUwD88ssvNG/enLp169K4cWO2bduWbds2L9ErcaWU8iNYF4yZ9QS9ZcsWrrzyygzHX3jhhXzxxRd4vV527txJhw4dMu1X4vDhw3z88cds27YNEeHo0aMAjBo1is8++4wyZcqkDvNVtWpVVqxYgdvt5ssvv2To0KHMnz8fgE2bNvH999/j8Xi4/PLLeeihhyhXrlzqvLt27aJ8+fKcd955ZyyzXr16bNmyhcqVK7Nz504++OADpk2bRrt27Zg/fz6dO3emZ8+eTJkyhSpVqvDdd9/Rp08fli5dmvEGy6c0iSulVBjo27cvK1euJCoqinXr1pGQkEC/fv3YtGkTLpeLHTt2ZDp/0aJF8Xq99OjRg5YtW6ZeaTds2JBu3brRrl07Wrdu/Z/5jh07RteuXdm5cyciQkJCQuq4pk2bUrRoUQCqVavGb7/9dkYSD0SlSpWoXbs2AHXr1mXPnj2cPHmS1atXc/fdd6dOFxcXl6Xl5hdanK6UUn4YE5y/zFSvXp2NGzemvp80aRJfffUVBw8eBGDChAmUKlWKH374gfXr1xMfHw+A2+0mOTk5db6U547dbjdr166lbdu2fPrppzRv3hyAKVOmMGbMGPbt20fdunU5fPjwGXE89dRT3HDDDfz000988sknZzzH7PF4Ul+7XC4SExPPmLdy5crs3buXEydOnDF8w4YNVK9ePcNlJCcnU6xYMTZt2pT69/PPP2e+wfIpTeJKKZUL3XjjjcTGxvL666+nDouJiUl9fezYMUqXLk1ERASzZs0iKSkJsPfNt27dSlxcHEePHuWrr74C4OTJkxw7dowWLVowYcIEfvjhB8Dee27QoAGjRo2iZMmS7Nu374w4jh07RpkyZQBS760HqlChQnTt2pVHHnkkNb6ZM2cSExPDjTfemOF85513HpUqVWLu3LmAreSXEq86kyZxpZTKhUSEBQsW8M0331CpUiWuuuoqunbtynPPPQdAnz59mDFjBtHR0Wzbto1ChexTKuXKlaNdu3bUqFGDdu3aUadOHQBOnDhBy5YtqVWrFo0aNeKll14C4PHHH6dmzZrUqFGDa6+9lujo6DPiGDRoEEOGDKFOnTr/udIOxLhx4/B6vVx22WVUqVKFuXPn8vHHH/utmf7ee+/x1ltvER0dTfXq1Vm4cGGWPzs/EOOvTCeXqVevnsms8kZep4+YKZUzfv75Z6644opQh6HymfT2OxHZYIxJ95k/vRJXSimlwpQmcaWUUipMaRJXSimlwpQmcaWUUipMaRJXSimlwpQmcaWUUipMaRJXSqlcKty7Ir3rrruoXbs2l156KUWLFqV27drUrl07tbMWde6C2na6iDQHXgFcwJvGmGfTjC8PzACKOdMMNsYsDmZMSikVDvJCV6QpJxzLli1j/PjxfPrpp+e8THWmoF2Ji4gLmATcClQDOohItTSTPQnMMcbUAdoDk4MVj1JKhZO80BVpemrWrMnRo0cxxlCiRAlmzpwJwL333ssXX3xBbGws3bt3p2bNmtSpU4evv/76bDdhvhDMK/GrgF3GmN0AIjIbaAVs9ZnGACl91BUF/ghiPEopdVa0K9Kz64o0PQ0bNmTVqlVUqFCBSy65hBUrVnDvvfeyZs0aXn/9dSZNmoSI8OOPP7Jt2zaaNWvGjh078Hq9mS43vwrmPfEygG9L+vudYb5GAJ1FZD+wGHgovQWJSE8RWS8i61N68FFKqfykb9++REdHU79+fQASEhJ44IEHqFmzJnfffTdbt27NdH7frkg/+ugjChYsCJzuinTatGmpnZT4OnbsGHfffTc1atTg4YcfZsuWLanjUroi9Xq9qV2R+tO4cWOWL1/O8uXL6d27Nz/++CO///47xYsXp1ChQqxcuZLOnTsD9gSiQoUKfrtZzc9CXbGtAzDdGFMWaAHMEpH/xGSMmWqMqWeMqVeyZMkcD1Iplb9pV6Rn1xVpeq677jpWrFjBihUraNKkCSVLlmTevHk0btzY77zqv4KZxH8HfMtVyjrDfPUA5gAYY9YAXuCCIMaklFJhIS90RZqecuXKcejQIXbu3Mkll1xCo0aNGD9+fOr9+MaNG/Pee+8BsGPHDvbu3cvll19+zp+bVwUzia8DqohIJRGJwlZcS1utci/QFEBErsAmcS0vV0rle3mlK9L0NGjQgMsuuwywSfv333+nUaNGqeuVnJxMzZo1ueeee5g+ffoZV/zqTEHtilREWgAvYx8fe9sYM1ZERgHrjTGLnNrq04DC2Epug4wxn2e2TO2KVLsiVSonaFekKhSy2hVpUJ8Td575Xpxm2NM+r7cCDYMZgzo36Z00QPacOOgJiVJKnZtQV2xTSiml1FnSJK6UUkqFKU3iSimlVJjSJK6UUkqFKU3iSimlVJjSJK6UUgGIT0z2P1E2L2///v20atWKKlWqULlyZQYMGJDaMltaf/zxB23btvW7zBYtWqTbRnogRowYwfjx4/8zvFu3bsybN++MYYULF/a7vCZNmmTa3jtAxYoVOXToUMAxTp8+nX79+mUab2YyitvlclG7dm2qV69OdHQ0L7744hkt46Vnz549qT3QBUtQHzFTSqm8IsodQfupa7JtebN7XpPpeGMMrVu3pnfv3ixcuJCkpCR69uzJsGHDeOGFF86YNjExkYsvvvg/iTQ9ixdrb89no0CBAmzatAmAv//+m44dO3L8+HFGjhyZ4TwpSbxjx45Bi0uvxJVSKhdaunQpXq+X7t27A/ZKcMKECbz99tvExMQwffp07rjjDm688UaaNm3Knj17qFGjBmCbZ23Xrh3VqlXjrrvuokGDBqlXvClXtnv27OGKK67ggQceoHr16jRr1oxTp04BMG3aNOrXr090dDRt2rQ5o7nXrFq2bBlNmjShbdu2VK1alU6dOpFeI2O9e/emXr16VK9eneHDh58x7vnnn6dmzZpcddVV7Nq1C4CDBw/Spk0b6tevT/369Vm1alWmcfzyyy80b96cunXr0rhxY7Zt2wbAr7/+yjXXXEPNmjV58sknA1qnCy+8kKlTp/Laa69hjMmwW9jBgwezYsUKateuzYQJEzKc7lxoEldKqVxoy5Yt1K1b94xh5513HuXLl09NZBs3bmTevHl88803Z0w3efJkihcvztatWxk9ejQbNmxI9zN27txJ37592bJlC8WKFUvtYrR169asW7eOH374gSuuuIK33nrrnNbl+++/5+WXX2br1q3s3r073YQ7duxY1q9fz+bNm/nmm2/YvHlz6riiRYvy448/0q9fPwYOHAjAgAEDePjhh1m3bh3z58/n/vvvzzSGnj17MnHiRDZs2MD48ePp06dP6nJSelMrXbp0wOt0ySWXkJSUxN9//53aLezGjRv58MMP6d+/PwDPPvssjRs3ZtOmTTz88MMZTncutDhdKaXC1M0338z555//n+ErV65kwIABANSoUYNatWqlO3+lSpWoXbs2AHXr1mXPnj0A/PTTTzz55JMcPXqUkydPcsstt2Qah6TT4brvsKuuuoqyZcsCULt2bfbs2ZPaVnqKOXPmMHXqVBITE/nzzz/ZunVratwdOnRI/f/www8D8OWXX57R/erx48c5efJkuvGdPHmS1atXc/fdd6cOi4uLA2DVqlWpJy9dunThiSeeyHRd05OQkEC/fv3YtGkTLpcrw65TA50uKzSJK6VULlStWrX/3OM+fvw4e/fu5dJLL2Xjxo2pnZ6crbRdiaYUp3fr1o0FCxYQHR3N9OnTWbZsWabLKVGiBEeOHEl9/88//3DBBac7pPTXZemvv/7K+PHjWbduHcWLF6dbt25ndHnqe0KQ8jo5OZlvv/0Wr9frdz2Tk5MpVqxY6j3ttNI7CfFn9+7duFwuLrzwQkaOHJnaLWxycnKGMfl2H5vZdFmhxelKKZULNW3alJiYGGbOnAlAUlISjz76KN26daNgwYKZztuwYUPmzJkDwNatW/nxxx+z9NknTpygdOnSJCQkpHYLmpkmTZrw4Ycfptacnz59OjfccEPAn3f8+HEKFSpE0aJF+euvv1iyZMkZ4z/88MPU/9dcYysENmvWjIkTJ6ZOk1GCBnsbolKlSsydOxewlQZTumJt2LAhs2fPBghoXcHej+/Vqxf9+vVDRDLsFrZIkSKcOHEidb6MpjsXeiWulFIBiE9M9lujPKvLi3JnfB0lInz88cf06dOH0aNHk5ycTIsWLXjmmWf8LrtPnz507dqVatWqUbVqVapXr07RokUDjm306NE0aNCAkiVL0qBBgzMSUXpatmzJhg0bqFu3Li6Xi8qVKzNlypSAPy86Opo6depQtWpVypUrR8OGZ/aLdeTIEWrVqoXH4+GDDz4A4NVXX6Vv377UqlWLxMRErrvuukw/87333qN3796MGTOGhIQE2rdvT3R0NK+88godO3bkueeeo1WrVhnOf+rUKWrXrk1CQgJut5suXbrwyCOPAHZ7t2nThpkzZ9K8efPUEpJatWrhcrmIjo6mW7duGU53LoLaFWkwaFekOdvzl/ZipvKrcO6KNCkpiYSEBLxeL7/88gs33XQT27dvJyoqKtShKT9yVVekSimlcl5MTAw33HADCQkJGGOYPHmyJvA8SpO4UkrlMUWKFPHbEprKG7Rim1JKZSDcbjeq8HY2+5smcaWUSofX6+Xw4cOayFWOMMZw+PDhLD92psXpSimVjrJly7J//34OHjwY6lBUPuH1elMbxQmUJnGllEpHZGQklSpVCnUYSmVKi9OVUkqpMKVJXCmllApTWpyuwkIwG51RSqlwpVfiSimlVJjSK3GllApARqVBKbRUSIWCXokrpZRSYUqTuFJKKRWmNIkrpZRSYUqTuFJKKRWmNIkrpZRSYUqTuFJKKRWmNIkrpZRSYUqTuFJKKRWmNIkrpZRSYUqTuFJKKRWmNIkrpZRSYUqTuFJKKRWmNIkrpZRSYSqoSVxEmovIdhHZJSKDM5imnYhsFZEtIvJ+MONRSiml8pKgdUUqIi5gEnAzsB9YJyKLjDFbfaapAgwBGhpjjojIhcGKRymllMprgtmf+FXALmPMbgARmQ20Arb6TPMAMMkYcwTAGPN3EONRSuVR/vr6Bu3vW+VNwSxOLwPs83m/3xnm6zLgMhFZJSLfikjz9BYkIj1FZL2IrD948GCQwlVKKaXCS6grtrmBKkAToAMwTUSKpZ3IGDPVGFPPGFOvZMmSORyiUkoplTsFM4n/DpTzeV/WGeZrP7DIGJNgjPkV2IFN6koppZTyI5hJfB1QRUQqiUgU0B5YlGaaBdircETkAmzx+u4gxqSUUkrlGUGr2GaMSRSRfsBngAt42xizRURGAeuNMYuccc1EZCuQBDxujDkcrJhUxgKpGKSUUip3CWbtdIwxi4HFaYY97fPaAI84f0oppZTKglBXbFNKKaXUWdIkrpRSSoUpTeJKKaVUmAooiYvI7SKiCV8ppZTKRQJNzPcAO0XkeRGpGsyAlFJKKRWYgJK4MaYzUAf4BZguImucplCLBDU6pZRSSmUo4CJyY8xxYB4wGygN3AVsFJGHghSbUkoppTIR6D3xViLyMbAMiASuMsbcCkQDjwYvPKWUUkplJNDGXloDE4wxy30HGmNiRKRH9oellFJKKX8CTeIH0iZwEXnOGPOEMearIMSllFJhxV/TxdqfuQqGQO+J35zOsFuzMxCllFJKZU2mV+Ii0hvoA1QWkc0+o4oAq4IZmFJKKaUy5684/X1gCTAOGOwz/IQx5p+gRaWUUiGMO+4AACAASURBVEopv/wlcWOM2SMifdOOEJHzNZErpZRSoRPIlXhLYANgAPEZZ4BLghSXUkoppfzINIkbY1o6/yvlTDhKKaWUCpS/im1XZjbeGLMxe8NRSimlVKD8Fae/mMk4A9yYjbEopZRSKgv8FaffkFOBKKWUUipr/BWn32iMWSoirdMbb4z5KDhhKaWUUsoff8Xp1wNLgdvTGWcATeJKKaVUiPgrTh/u/O+eM+HkPxm1t6ztLJ89f21Yp6XbWikVrgLtirSEiLwqIhtFZIOIvCIiJYIdnFJKKaUyFmgHKLOBg0AboK3z+sNgBaWUUkop/wLtirS0MWa0z/sxInJPMAJSSimlVGACvRL/XETai0iE89cO+CyYgSmllFIqc/4eMTvB6TbTBwLvOqMigJPAY0GNTimllFIZ8lc7vUhOBaKUUkqprAn0njgiUhyoAnhThhljlgcjKKWUUkr5F1ASF5H7gQFAWWATcDWwBm07XSmllAqZQK/EBwD1gW+NMTeISFXgmeCFpXKjxAT4eVMEv+2I4LOj8NtvcOAAuN0QFWX/SpeG6tWhRg37V7RoqKNWKmccPQzH/hEMYJIhORlMslCiVDLFLwh1dCqvCjSJxxpjYkUEEfEYY7aJyOVBjUzlCv8chK8/cbPuGxeb1riIOSkBzxsRAQ0bwp132r9LLglioErloLhY2PFjBD9/H8HWjS62fh/BgX0ZP+xTsnQy114Fdevav4YN9QRXZY9Ak/h+ESkGLAC+EJEjwG/BC0uFkjGwaU0EC2dGsnyJi8SE04m7/KXJXFEniQZXRFKxor3yNgbi4yE21l6db9kCP/0EP/4IK1bYv0cfhWuvhccegzvuAJcrdOun1Nnas0NYMCOS/5vr/s8JbYGChhIXGSIi7AmsODn9wD7h4J8RLFwICxfaYV4v3HMPPPggXH01SODnxkqdIaAkboy5y3k5QkS+BooC/xe0qFRIGAPffe1i6jOR7Npqs2xEhKFhs0QaN0+ibuMkSpUxADQuHel3ecePw5IlsGAB/O9/sHo1tG4Nl14KjzwCVW8Dt//FKBVSCQmwaBGMneDl+1Wnzz4rXZ5MtSuTqH5lMldcmUTFy0y6J6fJybDvF0H2FWTDBvj2W1izBmbMsH+1atlk3rkznHdeDq6YyhOyUjv9SqAR9rnxVcaY+KBFpXLcth8ieH1MVOpBqkSpZG7vmMhtHRNTE3dWnXeevdq45x44eRLeeQcmTIBdu6BPH6hQpQD9R8VR//rk7FwVpbKFMfbK+eGHYc8eABcFChpubpPIXV0TqFwtsN9FRARUqGJofB106mSH7doF06bB22/D5s3Qty8MHw7jxsF999l5lApEoB2gPA3MAEoAFwDviMiTwQxM5YyTJ+0BpOetBfh+lYsixQx9no5j9ppT3Pd4wlkn8LQKF4aHHoIdO+DDD6FKFfhtZwSPdijA0O4e/tyn5Ykq99i9G1q2hLvusgn8ssug/6g45m+M4bHn4gNO4Bm59FJ47jnYvx8++MDeajp0CB54wBavr12bPeuh8r5Az/c6AfWNMcOd7kmvBroELyyVE775xhblTZ4M7khDh97xzF4dQ/teiXi8/uc/G243tGtn75f3GhZPgUKGlZ+56d60AJ/Nc2Oy55xBqbMSGwujRkG1arB4sa189tprsHUrtL0/kcLZXNzt8UD79rByJcyeDWXKwLp10KAB9OgBf/+dvZ+n8p5Ak/gf+DTyAniA37M/HJUTEhLsPekmTeDXX6F2bZi25BS9n0qgSLGcicHjgY59E3hvxSka35pIzElhbH8PHTvC0aM5E4NSvn74AaKjbbF2XBx06QLbt9uSqmBXxBSxt522bYPBgyEy0ha116gBX34Z3M9W4S3TJC4iE0XkVeAYsEVEpovIO8BPgN9DrYg0F5HtIrJLRAZnMl0bETEiUi+rK6Cy5s8/4cYb7b1ptxtGjLBFd+daPHi2LrjIMObNOJ54MY4CBQ2zZ9sD6caNIQlH5VMzZ9pi7B077FX4smV2WKlSORtH4cL2vvhPP8ENN8DBg9CsGYwZYyvIKZWWvyvx9cAG4GNgKPA1sAwYBizMbEYRcQGTgFuBakAHEamWznRFsI3JfJfF2FUWrV5tn1FduRIuvtg++jV8uD3rDyURuK1DIm9+for69WHvXmjUyN47VyqY4uKgd2/o2tUWpd93H6xfD9dfH9q4LrsMvvjC/j4BnnrK3qM/fDi0cancx18HKDNSXotIFHCZ83a7MSbBz7KvAnYZY3Y7888GWgFb00w3GngOeDwLcassmjkT7r/fFqVfdx3MmZPzVxn+lLvEsGKFPai+8469V7hliy0tUOFtxZ//Zjq+celCORTJaXv3wt1325IojwcGjImjZadE1h8lgHLGrPO3DeDM7eBy2X3/6qvt42dLlsCVV8LcuXDVVcGJIRTfgzo3gdZObwLsxF5ZTwZ2iMh1fmYrA+zzeb/fGea73CuBcsaY//n5/J4isl5E1h88eDCQkJXDGFsLtmtXm8AHDLD32HJbAk/h8cBbb8FLL9nHbEaPho4dIUEfaFTZaP16Wyq1di1UqGBLp1p2Sgx1WOlq3tzeXmrQwJ54XH89fPJJqKNSuUWgFdteBJoZY643xlwH3AJMOJcPFpEI4CXgUX/TGmOmGmPqGWPqlSxZ8lw+Nl9JSrJJe/BgW2T9yivw8suhLz73R8Q+m7t4MRQpYovVh3b3cCom1JGpvODrr+395kOH4KabYMMGqJfLa+OULw/Ll9vStNhY++jbjBn+51N5X6BJPNIYsz3ljTFmB+AvFfwOlPN5X5Yza7QXAWoAy0RkD/axtUVauS17JCTYIriJE23HJLNnQ//+oY4qa265xR5wL7gAvvvazaPtvZzQmuvqHCxYALfeattH6NDBtiRYokSoowpMVBRMnQrDhtkT9G7d4MUXQx2VCrVAW2zbICJvAu867zthK71lZh1QRUQqYZN3e6BjykhjzDFswzEAiMgy4DFjjL/lhkRG95Kycg8pkHtiZyPtchMTYFRfD8s+dVOkiG0yskmToHx0QM5lvevWtRXwrm+azE/rXfRv62XCh7EUC5MDr8o9pk+3z14nJ9sWAydODL+W0URsTfWSJWHgQNsXwcGDtka7tr+ePwW6C/fCVkjr7/xtBXpnNoMxJhHoB3wG/AzMMcZsEZFRInLH2YesMpMQDyN62wRetKi9/x3KBJ4dqlaFSQtiKV85mV+2unjkHi/Hj4Q6KhVO5kxz0727TeBPPWUbcAm3BO5rwACYNcs+Jvrcc7btdX0ELX/yeyXuPCr2gzGmKvYedsCMMYuBxWmGPZ3BtE2ysmz1XwnxMKKXhxX/56ZwUcMXXwj164c6quxRqqzh5XmxDGjjZddWF4+0t1fklA51ZCq3mzPNzWvDPYBtH2HgwBAHlE06d4bzz4e2bW077G43TJqkV+T5jd9zUWNMErBdRMrnQDzqLCUlwTMDbAIvUsww4cPYPJPAU1xQyvDy3FjKVEpmx48uHu3o5dixUEelcrNP3zudwN98M+8k8BQtWtia6h4PvP46DBqENl2czwRaoFQc22LbVyKyKOUvmIGpwBkDE4ZG8dVCNwULG178IJbLa+XNsrWSpQ0vz4mldPlktm1ycfvtcOpUqKNSudGXH7t4YVAUAANGx9GjR4gDCpKmTWH+fPvUyfjxMHJkqCNSOSnQim1PBTUKdU6mPRvJolmRRHkMz86IpWp03kzgKUqVMbwyL5a+rbysWBFB+/b2IOYOuGNdldet+D8XY/t7MEboOSSeNj0SsV0+5E233Qbvv2/bXx85EgoVgse1+ax8wV/b6V4RGQjcDVTF9iP+TcpfjkSoMjV7ipt3J0bhchlGvhFH7WvydgJPcVFZW+JQvLitff/gg1qMqKx130QwopeHpCSh80PxdH7IX+OSeUPbtrYGvogtVp80KdQRqZzgrzh9BlAP+BHbBro+lZiLfLXAxeRR9upiyMtxNGyWFOKIclbFywz/+x8UKGB7fBo6NNQRqVDb+n0Ew+7zkhAvtLkvgQcG548EnqJLF3tvHOChh2DevNDGo4LPXxKvZozpbIx5A2gLNM6BmFQANq2J4JmBNoH3fiqOZm3yVwJPcc019kDlcsGzz8KUKaGOSIXKH3uFIV29xJ4SmrdL4KFR8fmypvaDD9rnxo2xNdhXrw51RCqY/CXx1NNY57lvlQvs2SkMda42WndPoH2v/P3VtGhhH7EB2/fzt0uD3PmzynWOH4FBnbwcOSTUuy6JQS/Eh/Vz4OfqiSegVy/bS9sdd9guVlXe5G83jxaR487fCaBWymsROZ4TAaoz/XPQHqxOHhMaN0/Mt1cbaXXvbpujTE6G4Q962LUlHx/B85m4OBjWw8veXyK45IpkRk+LxZ3L+wcINhHbIt1tt9nuS2+9Ff7+O9RRqWDI9EhnjHEZY85z/ooYY9w+r8/LqSCVFRsLw+7zcmB/BNWuTOKp1+Jw6UVnqtGjbXvYp/4VnrjXw8E/9ewmr0tOtm2I//CtiwsuSub5WbEUKhLqqHIHt9v2mVC3LuzeDbffDrHaiVCeo5crYcIY24PRlg0uLrw4mWfeicVbMNRR5S4ith/yWlclcfDPCIZ29+hBK48bNswmqgKFDM/NjOPCi/URBV+FC8Onn0LFirbb1VF9PSTlz+ozeZYm8TAxbhy89x4UKGh4dkYc52uPrOnyeGDMW7YxmO2bXTz7iEcfPcujZsywlRldLhg9LY4qNfLH45VZddFFsGQJFC8OKz9z8+Zz+fxeQx6jSTwMfPyxveIQgacmxXFpdT1YZaZYCXh2RiwFChmWLnIz8xU9aOU1q1dDz5729aRJcFUTvbzMTNWqMHcuuFyG916L4vOP9D5cXqFJPJf78Uf77CfYq/FGt+jBKhCVLjcMnxyHiOGt56NYvkQPWnnFvn3QujXEx9unER58MNQRhYemTeGhUfEAPP+oh60b9fCfF+i3mIsdPgytWsG//9rnPQcNCnVE4eXam5N4cKh9SnJMPw+//KwV3cJdTAzceSf89RfceKPtlUwF7q5uidzRJYH4OGHofR7+/kN/E+FOk3gulZho20H+9Vdbu3TqVO1i8Gx06JPAza0TiT0lDLtP+yEPZ8bYRwk3boTKlW3xcKTeKckSERg4Jp461ybxz98RDLtPK3+GO03iudSgQfDVV3DhhfaeeIECoY4oPInAoBfiuKxmEn/8FsHIPlo7N1yNHQtz5kCRIra9/PPPD3VE4ckdCSPf0MqfeYUm8Vzo8/kuJkywVxnz50O5cqGOKLx5CsDYt+Ioer5h3Tdupo7Ty7dws2gRPPWUPSn74AOoVi3UEYW3tJU/P5isv4lwpUk8l9n5UwQvPG7bRH/lFWjUKMQB5RGlyhpGTY3F5TJ8MDmKrxZqRbdwsW2brRMC8MwzthUyde4qXW548tU4AKaOi2TtMv1NhCNN4rnIsX9gWA8PcbFC9+627WOVfepcm0zf4bZ27nOPeNi9TSsZ5HbHj8Ndd8GJE7arzSeeCHVEeUvjW5PoOjCe5GRhZB8Pu3eHOiKVVZrEc4mkJNua0oF9EVSNTmLyZK3IFgxteiTSrI1T0a2HlxPHQh2RykhyMtx7r70Sr1HDtsanv4ns1/2xBK69KZETR4W77rJPw6jw4Q51ALnNij9Dswe/9Xwk675xU/R8w+hpcXiD2KZqdqxjsLZTsLe/CDz2XBy7fxZ2bXUx9iEPzT/jPz1eZRRH49KFghqfOm3sWFi4EIoVs5U7CxcOdUR5U0QEDJsYx4O3RbB5cwQ9eth6B3rCFB70SjwXWLHExbsTo4iIMIyYEkupslpVNJi8BWHMW3EUKWZY/aWb0aNDHZFK69NPYfhwm0jefx8uvTTUEeVtRYrCM2/HUrgwfPghjB8f6ohUoDSJh9jeXcLYAbYi24NDE6jbSJtUzQkXVzA8Pcm26DZyJCxeHOqIVIodO2xFNmNgzBjbjaYKvoqXGWbOtK8HD4YvvwxtPCowmsRDKOZfePJ+LzEnhSYtE2nfOyHUIeUrDW5I4r7HEzAGOnVCK/XkAidP2opsx47Z/0OGhDqi/OWuu2w/DcnJ0L49/PZbqCNS/mgSDxFjbPvFe3ZEUKFKMoNfitN7UCHQpX8Ct98OR49Cmza2WU8VGiktsm3dajvsmDFD78uGwsiR0Ly5bfa5dWs4dSrUEanMaMW2EJkz1c3SRW4KFjaMeTOWglppJyQiImDmTKhfHzZtgt69Yfr0UEeVP40fD/Pm2RbZFiyw/1Pk5sqY2clfjDlRsdLlst0e169vm7jt3TtrTwac63bWyqNZo1fiIfD96gimjIkCYMjLcVSoohXZQqlYMfjoIyhY0Cb0KVNCHVH+8+WX9j4s2O/g8stDG09+d/759jdRoIAtEXn99VBHpDKiSTyHHfxTGNHLS1KS0LFvPNe30Ia8c4OaNeHNN+3rAQPgp/X608gpB/YL7dvb+7DDhtleylToRUef+ZtYvTq08aj06ZEqByXEw9M9PRw5JNRtlMT9T2hFttykQwd7sEpIsN/T4b/1hmywxZ2Cp+73cPiwvQ87cmSoI1K+OnaEgQNtr4pt2sAff4Q6IpWWJvEcNHF4FFs2uLjw4mSGvx6LW2sk5DovvACNG8OhAxGM6OUhUc+zgsYYGD84iu2bXVSqZO/DurT57lzn+efh+uvhwAHb9G18fKgjUr40ieeQ6dNhwYxIIqMMo6bFUaxEqCNS6YmMtN1dliiVzA/funjdqbugst/8t9x8NjcSbwHDggXatWhulfKbKFsW1qyB/v1DHZHypUk8B6xff7ozk4efiadaHW3QJTe76CIYNTUOl9swd1okX3ykl4fZbdOaCCaNtCdIgyfEUatWiANSmbrwQtv0rccDb7wB06aFOiKVQpN4kB08aJ+1jIuDO7ok0LJjYqhDUgGoWT+Z/qNsueHzj3nY+ZP+VLLLX/uFp3vayp0d+sRz4x1auTMc1Kt3+smNfv3g229DG4+y9MgURImJttWjffvg6qtJTQoqPNzZNZEW9yQQFys8eb+HY/+EOqLwF3cKnrzfw9HDQr3rkug5RCsdhJNu3WwCj4+3Fd0OHAh1REqTeBANGgRLl0KpUrYRiyhPqCNSWSECD4+Lp2p0En/ujWBkHy9JetF41oyBF4fYimylyycz4vVYrcgWhl56yVb+/OMPm8jj4kIdUf6mSTxIZs6ECRNspZD586FMmVBHpM6Gxwuj34yjWAnD+uUubcv7HMyd5ub/5tiKbM+8Hct5xUMdkTobkZEwd66t6LZ6NfTta0/QVGhoEg+C9euhZ0/7euJEaNgwtPGoc1OqjGHU1FhcbsMLL9hHoVTWrF3mYvKo0xXZKlfTo344K1XKNo3r9cJbb9njnAoNTeLZ7K+/bE9AcXHw4IP2T4W/2tecruh2//32RE0FZvt2GNHLQ3Ky0HWgVmTLK+rWtW2qAzz8MHzxRWjjya+CmsRFpLmIbBeRXSIyOJ3xj4jIVhHZLCJfiUiFYMYTbHFxNoHv32+vvl99NdQRqex0Z9dEHngAYmPt96yVevw7cgTuuANOHhca35pI98e0Ilte0r49DB1qm8xt1w527gx1RPlP0JK4iLiAScCtQDWgg4hUSzPZ90A9Y0wtYB7wfLDiCTZjbBH6mjVQrpxTkU3bCclTROC11+wJ2v799tHB2NhQR5V7pTydsWMHVK6WxLBX44jQsr88Z/RoaNXKdudrT9hCHVH+Esyf1FXALmPMbmNMPDAbaOU7gTHma2NMSg/O3wJlgxhPUL34oq3MVrAgLFxoGwxReU9UlK2oWK6cPWHr2VMr9WTk8cfh88+hZEkYNz2OgtrDZJ4UEQGzZkGNGrBtG4zs4yFRm8PIMcFsvbsMsM/n/X6gQSbT9wCWBDGeoFm82D5OBjaR16kT2nhUcJUqBZ98Yq/IZ82CatVOd6MZDMHuYzoY/WwvnOXm5Zc9uCMNT0+N5aKymZ/phENf3+EgVNuxSBFYtMj2Qf7dUjeTRhgGjNF2MXJCruiCQ0Q6A/WA6zMY3xPoCVC+fPkcjMy/zZttkaExMGqUfW5SBV+dEl4KR6X/kHFWktrJ+LOrZBUdDe++a4vUhwyx/V/fdddZLSrP+e5rFy8PtfeSHn02nugG2sxwflCpkm2atelNhvlvR1KmYjJt79dL8mALZhL/HSjn876sM+wMInITMAy43hiTbrMBxpipwFSAevXq5ZrCy0N/CW1vS+bEiQiatkrk+h5xrPgz1FHlD4WjXNz62iqKRv33jtCx+MCTxpJ+Z//83513wrhx9iq8c2dYvtzW2E2R1auic72izg1+2SoMf9BDUpLQpX88t3XQg3h+0rgxDH4pjtH9vEwcHsVF5QyNbtGnEYIpmPfE1wFVRKSSiEQB7YFFvhOISB3gDeAOY8zfQYwl252KgSFdPfz9RwQ16iUxeEIcot1P5zuDBkHXrhATAy1bwt69oY4odA4dEAZ18RJzUmjaKpEeg7Qmen50c+skejwejzHCqD4etm/W2ozBFLSta4xJBPoBnwE/A3OMMVtEZJSI3OFM9gJQGJgrIptEZFEGi8tVkpJgdF8P2ze7uLhCMs+8E4vHG+qoVCiIwNSpcMMN9pGzFi1sLd38JuZfeOJeDwf/jKBmfXtSqzXR8697ByZwy90JxJ4SnrjXw1/79QonWIJ6T9wYsxhYnGbY0z6vbwrm5weDMfDa8ChWfuamSDHD87NitW/wfC6lxnrDhrBli60XsSQsq2iencREGNnbw86fXJSplMzYt/WkNr8TgUEvxPP37xF8v9rFE129TPzoFEWKhjqyvEfPlbPo/UmRzH87ksgow5g3Yyl/aa65Ra9CqHhx+5RCqVK205v7788fj54ZA+MHRbHmSzdFi+tJrTotMgrnGJnM7p8jGNLNS9ypUEeV92gSz4LP5rl545koRAzDXo2jzrVa61adVrEifPqpbStg1iyYMiYy1CEF3ZSxkSyebTs1GTcjlnKX5IMzFxWwIsVg/PuxlCydzObvXIzorc+QZzdN4gFau8zFs4/Yx2b6jdT2n1X66tWDjz4Ctxs+eD2K2VNyxVOcQfHB624+mByFy20Y/WYcNerpSa36r4vKGsa/H0uRYoZVn7sZ/3hUviilyimaxAPw0/oInuzhISlR6NA7nrv12UeViVtugenT7evJozx8Pj/vdZq95EM3r4/2ADD05Tga3KAntSpjlS43PDczFo/XsPjDSN54Ju+XUuUUTeJ+7NoSwaAuXmJPCbfcncCDw/SxGeVfp07Qb4Rt9mDcwx5WfpZ3EvmK/3Px/GO2VGrA6Dhubq0JXPlXo14yo9+Mw+U2vD8pb5dS5SRN4pnY94vwaAcvJ4/ZHpieeDFeH5tRAWvXM5FO/eJJSrQNoKz7Jvx3nvXLIxjZ2zbm0nVgPG16aKmUCtzVNyYxZII9uZ08ysOCGZrIz1X4H1WC5MB+4eF7vBw5JNS7Lonhk+Nw6/6msqjnkARad08gIV4Y2t3LD9+F709uw4oIBnfzEh8n3Nk1gfse11IplXXN2iQxcKxN5C8N8bDoXT2wnovwPaIE0V/7hf5tvKmtsY19O5YoT6ijUuFIBPqPjqfFPQnExQpPdPGy9fvw+9ltWBnB4K5e4mOF2zslMHBsvLZQqM5a6+6JPDTKJvLxgzx8+p4m8rMVfkeTIPtrv9C/rZcD+yK4ok4Sz78bS4GCoY5KhbOICHh8fDxNWyUSc1J4tL2XLRvC56f3/eoIBt/rJS5WuK1DAo8+p7eV1Lm7+/5E+g63ifyFQVEsnq2J/GzoT9HHvn0w4G4vf+6NoGrtJMa/H0vh80IdlcoLXC4Y9mocN9yeyL8nbF2LH9fl/p/fpjURPNHFJvAW7RN4/AVN4Cr73PNgIr2G2XbWn3s0iv+bq4k8q3SLOU6etO1f//GbTeAvfhCrTQSqbOWOhKcm2TbFv1ro5rGOXp5/N5bi1ezVyApyV5/a65dHMOw++2TGre0SGDReE3ioBbu/8ECWn9297XXsm0ByMkwdF8W4gVHExsRxZ9eMK0zmhd7+spP+JB2FC0PPnlA1WhO4Ch63G4ZNjKNZm0RO/Ss83tHL5hW571z6609cDOri5VSM0LxdAoP0yQwVRJ0fSqDnEHtF/tIQDzNfidQGYQKkP0sfgwbBaws0gavgcrthyMtxtLjH9vL0Yu9CfLs49zR+sXCmmxG9PCQmCHc/kMDgl+Jx5Z3H3FUu1fmhBB57Pg4Rw5vPRTFpZBTJ2gigX5rE09Ba6ConuFww6MV42j2YQFKCMOmRgiF/ZtYYmPFyJC8O9mCM8MDgePqN0CtwlXPu6JzI8NfjcEca5kyN5NlHorStdT/056lUiEREQN+n42n36KnUYsS3ng9NMWJyMkwcHsVbz9sOfh57Po4u/RP0MTKV4268I4lnZ8ThLWD4vzmRPNnDw6mYUEeVe2kSVyqEROCOnnHcNyqGiAjDjJejGPOQh/i4nIsh5iQMu8/DvDdtF7sj34jjjs56+aNC56omSUyYYztNWf2Fm36tvPz1u55RpkeTuFK5wI33xDNuehwFChm++MjNI/d4OXo4+J/766/Q+44CrPrcTZFihuffjaVJS20LXYVe9brJTF54ijKVktm5xcWDLcKrfYWcoltEqVzimpuSeO1jp+/ltS56316A3duCd/XxzTdQvz78ui2C8pcmM+XTU9RtpDWJVO5RoYrhjU9PcWWjJP45GMGAtl7eey/UUeUuue/ZljBWp4SXwlHnXo33ZHzuuhLKzvX6/nBsNkSUd1WpkcyUT2MZ0s3Djh9d9LqtAIMnxGV7//XTpkGfPpCYCA1uTGT45Dht2CgPyUu/2fOKw/j3YnnlqSgWzoykc2fYsgVGj0afmkCTeLYqHOXi1tdWnfNylvRrmC2NOmRXwxDZuV65UbAb0AjU3hNOhyKF89fQIAAACwBJREFU4YlZ8bz9dEFWLYpiRC8va9fGcs+jsbicX2z5Imf3SNrJk9C3L8ycad8/9hjcNiBOD4Z5TF77zboj4dFn46l0eTITn/YwbhysXg3vvQdlyoQ6utDS4nSlciFPAej1fAz3PhmDy21Y/LaXsV0Kc+gcKvd8/z3UrWsTeIECMH06vPCCXs2o8NG6eyKffw4XXWRvB9WuDUuWhDqq0NIkrlQuJQLNusQzZMZJil+YzI6NbobeWYS1n2XtKtwYePVVuPpq2LEDataE9euha9cgBa5UEN14I2zaBM2awaFD0KIFPP44xMeHOrLQ0CSuVC5XtV4Szyw8QZ0bEog5HsGr/Qsx/okoYgN4dvbAAWjVCgYMsAe53r3hu++gWrXgx61UsJQqZa/Ax42zJUnjx0OjRrB1a6gjy3maxJUKA0XONzzy+r/c+2QM7kjDqs9dxJ7KeHpj4J13bLL+5BMoVgzmzYPJk21RulLhLiICBg+G5cuhfHlYtw7q1IFRo/LXVbkmcaXCRErx+sh5Jxg5JY5iJdKf7pdf4Oab4b774MgRuOUWW/zYpk3OxqtUTrj2Wti82XZgFR8Pw4fbuh/r1oU6spyhSVypMFOhajK1Gvz3ee6EBFusWLMmfPUVlCgB775rix0rVAhBoErlkKJF4Y03YOlSqFwZfvrJ1gF55BE4dizU0QWXPmKmVJgzBhYsgCeegJ077bBOnWDCBChZMrSxqbwntzySmZ4bbrBX5cOHw0sv2d/ArFkwcqS9UncHmPHOdR1zss9zvRJXKoxt3RjBdddB69Y2gV92GSxebK/ANYGr/KhgQfvo5Nq1trLboUO2bYRateB//yPP9VOuSVypMLRnhzCil4deLQuwciVccAG89potRrz11lBHp1To1a1rK73Nm2eL2H/+GVq2hJtuss+Y55VkrklcqTDyy2YXE/oW5N4mBVm6yE2UxzB4MOzaZa82Is+uITel8iQRW6FzyxZ48UX7lMbSpdCkCTRubOuLhHsyz/f3xOMTk4lynz6Xycl7GRmJT0zOljhyQ7vH+UXa/QiysC+VPv3yWGwin+49ecZoY2DLt24+ecPDljU2S0d5DC3aJ9KpbwJt6xc8p9j9ya52uOMSk/C4895ysnNZ2RlTbpNd+1F6vzV/PB5bya1bN5g4EV55BVatsg3F1KljH1W7667wPAnO90k8yh1B+6lrUt8fiz/7Xpyyq53hKHdEnmr3OD9Iux9B4PtSbOLp6b4e2Dj19Yl/hBULo/j6wyj+/NUe/LyFDP/f3v3HVlWfcRx/f3pvS0cpKBYI0gbJRroYnOgu4FJSfuow/mB/uAA6XRamsMiC4Y9l2z9uSzRziRuGLDOmuKBzmrViwsyy4kKNWwjIrxL5URAYrCU4cG5aXCy9t8/+OGdb18gocs49Pcfnldz0nvbknuebpn3u/Z7n+3wX39vHyocHuGZieT5CRNmHO4uvE+VrRfk6I81I6Oc+fnxQ9LZuXVDN/uSTQTviZcuCVq4rV8KDDwJVVxxm2Xzqk7hzI8nAABzakaejtYpd7ZUU+4Ne6VdPHGDxfX0svvcCNWONaz7hBijOOaitDTb/WbMm2ENgw4ag29tjj8Hjj8PsBaO4674icxaUGFWddLT/nydx5xJmA/D+n8dyrnMC9T+BM2fGACAZN87rZ+GyPmbOK/5nBzPnXDSqq2H1ali1Kphef/ppaG2Fndvy7NyWp6bWmPvlIgvuKjFrXonKEfgJ3f8tOJeQ3tOj6d4+iXf21XGh97//HSZMKdG0tJ/5X+2j7tqUV904lwJSsBxt7lxYvx4e3dDHa5vzvH0gR3tbJe1tlYwZZzTdVmTO/BKF5tJFOyaWmydx5xLyfvcY/vLGtQBUj/+ICTPPsfmJBrpqe9En33HUOXcF6upg+eoiy1cX6T4uOl7N0/HbHMcP5WhvraS9tRLJaPzCALPml5jVXOLzNw5QHW996UV5EncuIRNveI9pi3u4+oZz1DacR4JCoYEjR5OOzDkH0PBZ44G1/Tywtp9Tb4vtf8ix6/U8+3dW0LU/R9f+HM8/Bbm8MX3GADMKA8wolJh2B9TXlydGT+LOJaSqpkjj0lP/U53unBuZpk43pk4vsuJbRT76J3TuyPHm6zk6t1dwoquCrs4cXZ052loq+QGwdWuwEVHcPIk755xzl6F6NNyysMQtC0sAfNgLh/dVcGB3joN7KjjamadQKE8ssSZxSUuAp4Ac0GJmPx7y81HAc8AXgb8By8zsZJwxOeecc1GqqYVC8wCF5mBWrWlSnooy9UON7TKScsDPgduB64EVkq4fctpK4O9m9jngZ8ATccXjnHPOlUO5EjjE2zt9NnDMzE6Y2QXgJWDpkHOWApvC523AIsnrcp1zzrnhkMXU/V3SPcASM/tmeHw/MMfM1gw650B4Tk94fDw8590hr/UQ8FB42AgciSXoaNUB717yrHTysaWTjy29sjw+H9ulTTWzj91cOBWFbWb2DPBM0nFcDkm7zaxMpQ3l5WNLJx9bemV5fD62KxPndPppoGHQcX34vY89R1IeGEdQ4Oacc865S4gzie8CpkuaJqkKWA5sGXLOFuDr4fN7gG0W1/y+c845lzGxTaebWVHSGqCdYInZs2Z2UNKPgN1mtgXYCDwv6RjwHkGiz4pUTf9fJh9bOvnY0ivL4/OxXYHYCtucc845F68yrmZzzjnnXJQ8iTvnnHMp5Uk8YpKWSDoi6Zik7yYdT5QkPSvpbLi+P1MkNUjqkHRI0kFJa5OOKSqSqiW9KWl/OLYfJh1T1CTlJO2T9GrSsURJ0klJb0nqlLQ76XiiJOkqSW2SuiQdlvSlpGOKgqTG8Pf178cHkh6J7Xp+Tzw6YavZo8CtQA9Bhf4KMzuUaGARkdQMnAeeM7MZSccTJUmTgclmtldSLbAH+EoWfndhF8QaMzsvqRL4E7DWzHYkHFpkJK0DCsBYM7sz6XiiIukkUBjaACsLJG0C/mhmLeEKptFm9o+k44pSmBNOEzQxOxXHNfyTeLSG02o2tczsDYJVBJljZmfMbG/4vBc4DExJNqpoWOB8eFgZPjLz7l1SPXAH0JJ0LG54JI0DmglWKGFmF7KWwEOLgONxJXDwJB61KUD3oOMeMpIIPk0kXQfcBOxMNpLohNPNncBZ4DUzy8zYgPXAd4AsbsxuwFZJe8L201kxDTgH/DK8DdIiqSbpoGKwHHgxzgt4EnduEEljgJeBR8zsg6TjiYqZlcxsJkHnxNmSMnE7RNKdwFkz25N0LDGZa2Y3E+wG+XB4SysL8sDNwC/M7CbgQyBrNURVwN1Aa5zX8SQereG0mnUjVHi/+GXgBTPbnHQ8cQinLDuAJUnHEpEm4O7w3vFLwEJJv0o2pOiY2enw61ngFYJbdlnQA/QMmhFqI0jqWXI7sNfM/hrnRTyJR2s4rWbdCBQWf20EDpvZT5OOJ0qSJki6Knz+GYLCy65ko4qGmX3PzOrN7DqCv7dtZva1hMOKhKSasMiScKr5NiATK0PM7B2gW1Jj+K1FQOqLSIdYQcxT6ZCSXczS4mKtZhMOKzKSXgTmA3WSeoBHzWxjslFFpgm4H3grvHcM8H0z+12CMUVlMrAprJStAH5jZplaipVRk4BXgveX5IFfm9nvkw0pUt8GXgg/8JwAvpFwPJEJ33TdCqyK/Vq+xMw555xLJ59Od84551LKk7hzzjmXUp7EnXPOuZTyJO6cc86llCdx55xzLqU8iTvnnHMp5UncOeecS6l/AarFGnoONJN2AAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Multivariate Gaussian Implementation of EM\n"
      ],
      "metadata": {
        "id": "qaQkPtnkQJzg"
      }
    }
  ]
}