Lecture13_01 - EM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The EM algorithm\n",
    "\n",
    "##### Keywords: maximum likelihood, mixture model, ELBO, full-data likelihood, x-likelihood, latent variables, log-likelihood, training set, normal distribution, z-posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib as mpl\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "pd.set_option('display.width', 500)\n",
    "pd.set_option('display.max_columns', 100)\n",
    "pd.set_option('display.notebook_repr_html', True)\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_context(\"poster\")\n",
    "import pymc3 as pm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\newcommand{\\isum}{\\sum_{i}}$$\n",
    "$$\\newcommand{\\zsum}{\\sum_{k=1}^{K}}$$\n",
    "$$\\newcommand{\\zsumi}{\\sum_{\\{z_i\\}}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A toy problem\n",
    "\n",
    "This example is taken from Efron and Hastie, Computer Age Statistical Inference.\n",
    "\n",
    "Assume we have data drawn from a bi-variate Normal\n",
    "\n",
    "![](images/gform1.png)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1.85,  1.  ]), array([[ 1.    ,  0.615 ],\n",
       "        [ 0.615 ,  0.5625]]))"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sig1=1\n",
    "sig2=0.75\n",
    "mu1=1.85\n",
    "mu2=1\n",
    "rho=0.82\n",
    "means=np.array([mu1, mu2])\n",
    "cov = np.array([\n",
    "    [sig1**2, sig1*sig2*rho],\n",
    "    [sig2*sig1*rho, sig2**2]\n",
    "])\n",
    "means, cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the samples below as blue circles. Now say we lose the y-values of the last-20 pieces of data. We are left with a missing data or hidden data or latent-variables problem. We plot both datasets below, with the y-values of the lost points set to 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "samples=np.random.multivariate_normal(means, cov, size=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1147ef1d0>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEHCAYAAACgHI2PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wHPd93/H37t4dHkiAJEDSBAhLICnqR1AUJIquY4uW5cie2pZHHqel6j/cOJmO1c50HGkk20k9HTt2MnHSJJ44qe16Mm3Gjd2ZJqOmiR8647p2XVpPtoSQAkWRS1IASBEgQREgQBCP97D943Dg4XAPe3eLe8LnNcMhube3+9s74Lu//e1vv1/L8zxERKRx2dVugIiIrC8FehGRBqdALyLS4BToRUQanAK9iEiDU6AXEWlwoaA2ZIyJAF8Efh3YDvwC+Kzruv8Y1D5ERKR4Qfbo/wx4Evgj4GPAHPB/jTF3BrgPEREpUiCB3hizBXgC+JLruv/Jdd0fA48DYZI9fBERqZKghm5mgV8BRtKWRQEPaApoHyIiUoJAAr3rujHgBIAxxgZ6gS+RDPTfDWIfIiJSmsBuxqb5AskgD/BF13XdYjcwMDCgBDwiIkU6cuSIlW35egT6/wn8DPhV4IvGmIjrul8odiNHjhwJul0bxpkzZwDo6+urckskH31Pta+evqOBgYGcrwUe6F3XHVz+5/8zxrQBnzPG/J7rutGg9yUiIoUFEuiNMbuADwPPuq47k/bSCZI3YzuBq0HsS0REihPUPPqtwF8BxzKW/1Pg2vIfERFJE40lmJlbIhpLrOt+gpp1c9YY8z+Ary4/ITsE/DOSc+j/leu663sUIiJ1ZHhsmuMnRjl14TrxRALHtunfv533Hu6ht6s98P0FOUb/SeB3gc8DXcDrwOOu6z4b4D5EROraC4NjPPvT86RX94snEpxwr3Hy3Fsce2Q/D/Z3B7rPwAK967pzwO8s/xERkQzDY9Nrgnw6z/N49qfn6d6xOdCevbJXiohUyPETozmDfIrneRw/cTnQ/SrQi4hUQDSW4NSF677WHTx/PdAbtAr0IiIVsLAUI57wF7zjiQSLS7HA9q1ALyJSAc2REI7tL+Q6tk1TJLi5Mgr0IiIVEA7Z3HvXdl/r9u/fTjgUXHhWoBcRqZCHH+jBsrLmHVthWRbvPdwT6H4V6EVEKqS3q51jj+zPGewty+Lx9+8P/KGp9cheKSIiOTzY3033js0cP3GZwfP192SsiIj40NvVTm/XQaKx5Oyapkgo0DH5TAr0IiJVEg7ZhEORdd+PxuhFRBqcAr2ISINToBcRaXAK9CIiDU6BXkQCV6nKSeKPZt2ISGAqXTmpXr0yOrhmWTzusRSLEwk5OM7qB6resbu/rP0p0ItIIKpROakRjE/OcnpokpErN0kkEti2zZ7udu7Z28nOba2B7EOBXkTKVq3KSfXu7Mgkzw+OrfrcEokEb1yeYmh0mqP93Rzo7Sh7PxqjF5GyVatyUj0bn5xdE+TTeZ7H84NjXLsxV/a+1KMXkbKUUjkp83H/So9Z14LTQ5O+To6nhyZ49FB5+1KgF5GylFI5Kd9j/5UYs662eNxj5MpNX+sOj93MenIshgK9iJQlVTnJT7AvVDmpUmPW1bYUi5PweXJM+Dg5FqIxehEpS1CVkyo5Zl1tkZCD7bOsoB1AWUEFehEpWxCVk4oZs653jmP5nn20p7u97BTGCvQiUrZyKyeVMmZd7w7t6/R1crxnb2fZ+9IYvYgEopzKSZUes64FO7e1crS/O+dwlWVZHL2vO5Ab0Ar0IhKYUisnpcas/QT7IMasa8WB3g46tjRzemiC4TE9GSsidaTYykmpMeuh0amC6/bs3FxO02rCqucAdsOjh1jXsoIK9CJSEw7t62R4bDrrMEY0lmB+KU40Gica9/j8N55ruGRp61lWUDdjRaQmpMasM29QLizFmb61xNJSnE0tYcKOtZIs7S/+5gQvDI5VqcX1Qz16EakZmWPWi0sxZuejNEUcmptChDNSIShZmj8K9CJSdbnGrL/9g9OcHprAsoDY8p8MqWRpvV0HK9XcuqOhGxGpWe7FGxSYag7cTpYm2SnQi0hNKiVZmmSnQC8iNSmVLM2PQsnSNjoFehGpSUElSxMFehGpYUEkS5MAZ90YYxzgKeAJ4A7gIvBN4Buu6+ZPSScikkUqWVquerSFkqVJUpCDWl8A/h3w+8BLwEPA14BW4I8D3I+IbCDlJEuTpEAC/XJv/hngT1zX/YPlxT8xxuwAPosCvYiUodRkaZIU1CfVDvw18HcZy11ghzFmU0D7EZENLByy2dwaUZAvUiA9etd1bwCfzvLSY8Bl13Vng9iPiNSuV0YH1yyLxz2WYnEiIQcnI33BqqdhZV2t28RTY8yngA8AT5by/jNnzgTboA1kfn4e0GdY6xrte7p8882Vf0/ejDI8Ps+ViUUSJIcOujqb2LOrhY62MACbboar09AiNMp3tC6B3hjzCeBbwLPA19djHyJSm0bG5zk1fIv0OTIJYHRikbGJRe7ds5net7VUq3kbUuCB3hjzDPCnwPeAT5Q6tbKvry/Qdm0kqd6HPsPa1mjf0+xolPHJWdzRYcLh3L11d3SJu/f21MVx19N3NDAwkPO1QO9oGGO+AnwV+A5wzHXdpSC3LyK17fTQZNb57uk8z+P00ESFWiQQYKA3xjwFfB74c+A3XddVhiGRDSQe9xi5ctPXusNjN5VtsoKCmkffBfwH4BTw34FfMcakr/KKAr9IY1uKxX0V9wZILGebXK/SebJaUGP0HwSagHuBF7O8vgO4HtC+RKQGRUIOtm37Cva2sk1WVFDz6L8NfDuIbYlIfYjGEiwsxWhefkrVcSx6u9oZGp0q+N493e166KmCdEoVkaIMj01z/MQopy6szjuz/Y45Du3rZHhsOu8NWcuyuGdvZwVbLAr0IuLbC4NjazJJxhMJTrjXWLg4xtH+bo72d/P84FjObJNH7+tm57bWSjZ7w9O1k0gNiMYSzMwt1fRMlOGx6ZzpgiE5bfL5wTE6tjTz2EN72dezFXu5QpRt2+zr2cpjD+3lwJ0dlWy2oB69SFVdmVjklz98fc0wSC2m3z1+YjTvkEzL4m4Arl9q55OPHuTRQyjbZI1QoBepklPDM/xscJKmpuaVZalhkJPn3uLYI/t5sL+7ii28LRpLcOqCv4lzg+evE40lCIfs5T+aQlltOsWKVMHw2DQ/G5wkkaOD7Hkez/70vO8HkNbbwlKMuM858vHlOfJSOxToRarg+InRnEE+xfM8jp+4XJkGFdAcCeHY/sKFoznyNUeBXqTCShkGqbZwyObeu7b7Wrd//3aNx9cYfRsiFVavwyAPP9CDZVl517Esi/ce7qlQi8QvBXqRCqvXYZDernaOPbI/Z7C3LIvH37+/5mYLiWbdiFRcahjkxcFLBdddr2GQUsv+PdjfTfeOzRw/cZnB87U/JVSSFOhFquDhB3r4xalLeW/IVmoYZHxyltNDk4xcuUkikcC2bfZ0t9O+Y44tm5pWrZs6GZhDDubQFqKxOOGQg+PEFORrmAK9SBX0drXzvv4OfjY4mfX1Sg2DnB2ZXJOuIJFI8MblKaLXJunr7aBnZ1vOk8E9ezuVzqAOKNCLVMm9e9rYviXC5ekIr557i2g8TthxuO/uHRUZBhmfnM2ZkwaS0zvPjEwyuxDj6pu3sp4MhkanOdrfzYHe1WkNSh0akvWhQC9STR54HsuFtC1KKrBcIj9l/2KxBGcvTrLF2pn19fT8NuzOvg1dDVSfAr1IlWSmQHBsC8/zKpICIR73GBmbJuF52HmmTC7GEsTiCchd63ulBuyjh9a+lm9oKNfVgARPgV5kWTnDDZlFOArxmwKhe8fmwIdwhsem+cnLl7g2NZ+8lLCgKezQ3OQQdla3PRZPgAcJD+w8U+hTNWDTj93P0FChqwEJhgK9SBZ+hxtyFeHINsaefiL56StvstDyVjKIhsNZg2jL4m6On7hMb9fBwI4rlU9+LpKWQ8eDxaU4i0txNreGaY44qcUrJ4J8QR6y14D1MzSU72pAgqNAL5LB73BDviIc+YZe4nGP82/eYG4hTjSewCKe1qsOEU67ckjPBFmuzHzyTWGHxaX4qnVuzUVxHIuwY2MBWKzp5WeTWQM2Hvd8J2TLdjUgwVKgF0njd7jBA068fCPvermGXl4bus6Nm4t4pKVByNGrjmfpKZcqM598c1NoTaAHWFiME25NBt2wYxMJO7B2tVUya8AuxeK+ioRD9qsBCZZOoSJp/A43vHTqiq/1MrNPjk/O8vLrVyHPUMituSjReHLbQaVAyJZILexYbG5de5d1MZqM6pZlcaB3GyGncH6bzBqwkZCzUl2qkMyrAQmeAr3IsmKGG67dmKNAnAfWZp88PZR8QCoSdvK+b2ExmcgsqBQIuRKpNUcctrQ10RRxVp187tjVzjvv2YW5o4O+3o68+W2y1YB1HMv3TeTMqwEJnk6jIsv8DjckPA/PS/7tFMjmmD70kn4iaWlyWFyK5nzfYjQOTnApEFKJ1FLBPlX2D6AFaHfAs5PHFHJsPvfRh3h1/DUA7uqEg7vmOD00wfCY/7nwh/Z1Mjw2nffKJ9vVgARPgV5kWWq4oVCwty0LyyLv/POU9KGX9BNJ2LFpiTjMZxkjB8CDjz60J7CplalEaifPXcu5jmWBY1ncf/cOwiF79dOquym6BuzOba0c7e/Oec8j19WABE+BXmRZarhhaHSq4Lo7t7VizRfeZvrQS+aJJBK2cWyIJexkDz5tTvum5jDvuT/YhGYPP9DDq+ffKtjDzncVUWwN2AO9HXRsaS76akCCpUAvksbvcMO77+3iH1+2igqa2U4kjmPT0hymjfCqh5L29WwNfNw6lU8+c0poenuDSqRW7tWABEuftkia1HBDoZuPB3o7SirCcWhfZ873pIL8eo5bP9jfzZMfP8xhs3Ol+Ilj2xw2O3ny44d5973rk3IBklcDm1sjCvJVoB69SAa/ww2lFOFIH7fOphLj1r1d7fR2HVQPewNRoBdZVspwQylBM3UiefHkMFcmFgGqMm5d7Hi71C8FepE8/AZDP+tlnkj2OFuIxT327L0r6wlCOd0lKAr0IlUUciw2txY+kSinu5RDgV6kximnu5RLd2BEapjfJGvXbsxVuGVSTxToRWpYMTndRXJRoBepUaXkdC9FNJZgZm6p5PdL7dMYvUiNWu+c7sVUx5L6pkAvDakRpib6TbIGxed0L7U6ltQnBXppePU6NbGYJGvF5HTPLCmYaT0Lk0t1rMsYvTHmo8aYmfXYtkgxzo5M8oPnhhkanVrpGaemJn7/50OcHZmscgvzy5cbJ6XY3DiZJQWzyVYdS+pX4IHeGPMg8F3yFksTWX/1NDUx1w1Rv0nW/F6ZZCspmEtmdSypX4EN3RhjmoCngN8HZgEl0ZCqKmZq4qOHKtSoDFcmFvnlD1/PekM0Jcic7rlKCmYTZGFyqa4gx+g/DHwe+BzQCXwmwG2LFKWUqYmVzuB4aniGnw1O0tTUvLIs5w3RgHK6Z5YUzCeowuRSfUF+iy8De1zXnTLGfKncjZ05c6b8Fm1Q8/PJ0kcb+TMcuX6JhcVF3+sPvvY6rU35C3aXwr05tGbZUtTj6o1FTlyYhiZYsG8H7MjcrpV/f/d/DbIwc42ujqZA29S9Dc6NLhRc7+7dm7hw3g103/WmUX6XAgv0ruuOBrUtkXKFQjY24GeQwgYiFejND1+d5/TFW9y4FSWxXDnQcaAlkiwrmCnhwck3ZgIP9Pfva+fC2ByJPKNatgX372sLdL9SPTV7XdbX11ftJtStVO9jI3+Gs6NR7r6Gr6mJ+3q2cu+hg+vWDoD/88uLvHJ2OhndsfCS/yAeh1vzcVo9m/ZNYVqam1e9f+yGx137TaDDSn1Ac9vaefQpqepY61ltql7U0+/SwMBAztdqNtCLlMtv/df1KtuX8vrwBK+cvbYc5LObW4gSCdu0ZCxfrxuipVTHkvqlQC8NK71sX66e63qX7QN46bWreYN8yux8lG2rO/TrekNUJQU3DgV6aWhBTk3MFI0lWFiK0ZwnQC5G47yVMU/fwsKyPDLPPdFYgkQC0u7N0r9/+7oHX5UUbHwK9NKQ3rG7fyUQ3/e2EI8esgPruRaTDGxuPromoAPYtkU8vvaF+PKJCJJXHOnz6UVKpUAvDWc9szIWmwystSWMZbEm2Ofq1TtpQf7x9+/XWLkEYl0Cveu6XwK+tB7blsZSbpbJzPefHZnkuVfHSHgedvj28hdG3uTFiyf59aPvKzkrYynJwJrCDju2tXJtcm2aBdtOngBSf5JDKLohKsFTj15qRrlZJl8fnuBHL11kMbo8Sd1KBtrmphBhxyo7K2MxycB6u25P13zw3i7+/vgbWW/IWhbYlo2Hx4fe3csn3/Me3RCVwOknSmpCuVkmz45MJoP8Uvx2QPVgcSnO9MwiC0vx5KISszKWkwzsQG8H7ziwM3eaPwveeXAX9+3foSAv60I9eqk6v1kmO7Y0w+7s73/u1dFkTz6HW3NRHCc5Tz0ViHMF1WyzacpNBvaBd95Jz842Xjx1hWs35lbG5rdtDvO+I70c6O3wtW2RUijQS9WVm2Xy9NAkCc8rOFd9YTFGu5P7IaR8N3F379hcUjKw9HsK79gN//IozC3EuDm7yJXLIzRH7Lp46lLqmwK9VFW5WSZT77ctKzk0kifYL0bjeDaEnLUPIfmZTXPvXds5ee5awXYWmvve2hyitTnEjWsappHK0E+aVFUpBbBzvT8SLpB90oOE560JxH5n0+x/+1Zf1Z40911qjQK9VFWqALYf2Qpgp7+/pVCaYSs5rJIZiP3OprlweYpjj+zPW+1Jc9+lFinQS1WlCmD7ka0Advr7w47N5tZwtrcC0BRx+BcfWB2Ii51N808O7uLJjx/msNm58nCTY9scNjt58uOHlfFRapLG6KXqys0ymf7+5kjyIauFxfia+fQfflfvmkBcymwaJQOTeqOfTqm6cgtgZ74/7Ni0tYbZvqWZji3NbN/SzIcfzD6FMVVaz4/MTJLhkM3m1oiCvNQ89eilJpSSZTIaSzC/GCMScnK+f1+BJ2vDITuw2TQitUqBXqpqVe4anwWwc893P8ijh9qLHlJ5+IEeXj3/VsGhI82mkXql7onUnHxDIi8MjvEf//YkJ89dWxlbT813/4u/OcELg2NFD6n0drVrNo00NPXopW6Ukj3SL5XWk0amQC91o9TskX5pNo00KgV6qQulZI8sNUirtJ40GnVXpC6UMt9dRJIU6KUulDPfXWSjU6CXupCa7+6H5ruLrKZuT5WUWyt1I9J8d5HSKNDXgHJrpW4UqfnuuaZYar67SHYK9FV2dmRyTRm9VK3UodFpjvZ3q8xcmqDmu+uKSjYSBfoqKrdW6kYV9Hx3XVFJo9MdqyoqplaqrBVE9sizI5P84LlhhkanVipVpa6ovv/zIc6OTAbVXJGqUaCvklJqpUqw/F5RXbsxV+GWiQRLgb5Kyq2VKuXTFZVsFAr0VVJurVRZKxpLMDO35OvqR1dUspEoelRJqtbp0OhUwXWz1UqV23Lnp889C6eUKyrlv5F6pUBfReXWSpVkfvrMefWp/PQnz73FsUf282D/2oLdqSsqP8FeV1RS79RNrKJya6VudH7z02cbokldUfmhKyqpd/rprbIDvR089tBe9vVsXRmzt22bfT1beeyhvRy4Uw9L5VJMfvpsDu3rzHmSTdEVlTQCXY9WSSm1UuW2IPLTp66ock2x1BWVNAoF+hqighf+lZKfPttne6C3g44tzZwemmB4TE/GSmNSoK+iaCzBwlKM5gbtwa9nPplUfno/wT5bfnpdUclGokBfBaVMB6x3QeeTSeWnP3nuWsF1/ean1xWVNCoF+gordTpgPVuvDJ3KTy/ijwJ9hbwyOsj45Cw/eH4YL5I7MH3n+VGm7L08euhdWbfhx0tXXwFgtj3KhYkRABIJj1giQci2se3VQyZ3dfau/LuUdLy52nVhYoSpWwu8fGYcL3T7mCPR27NY/GbozNxH6ri6D81wZmSSqLN6CmUo1oZlWfTt6eClt6a5nuhd9XrqOP1+pqn93dXZu/LvfFKf6YWJkVWfL8A/nPkRAPPz8wC0jLasvNbRsnXNtp4++oSvNuaTOk4/bYfb7ffzOWXbZuYxp/P7M+b3uyl2u9VSzeMJNNAbY54AfhvoAU4Cz7iu+2KQ+6hnxeRWefRQMPucurXApau3GL8xh5dIYNk2uzpauWNXG1s2NQWzkzwuXb21rsfcs7ONtk0RXh+LMjWziOclsCybru2bK3aMIrUusEBvjPkN4FvA7wEvA78F/MgYc5/rusNB7adelZJbpdwbgmdHJpO96bRA6yUSXLl+i6sTs/T1dtCzs62sfeSTSHiM+8z8WM4xb9nURO+uNhJva1sZ/9/d7q++rB8XJka4MjO+8v8rM+MkPFb2lXGBRFfb2wLbt0gQAgn0xhgL+DLwl67rfnl52Y8BF3gaeDKI/dSzSudWuTUX4/TZsVVDJuk8z+PMyCRtmyKwTs8DxRIJvAoes22B7azvbJmpWwuMXJlh6tbtq4etbU3s2NrCpmaNhEptCuq34i7gTuB7qQWu60aBHwIfCmgfda3S2SrHp5Z8DZlcujpT1n7yCdk2VgNl6Lw+Pc/Lr49zY2Yez0uewDwvwY2b85y/dIPr0/NVbqFIdkH9Zt29/PeFjOVDwD5jjOO6bryYDZ45cyaQhtWKKzcv07kZRieiBdfd3dnEhfPumuWXb77pa19LS1Emby7hRKPE7Pwf++XxaS5eepNUR3jTzbCvffhp1+TcBO3NMDGT0Ybo2s8g1zHn2sf1udVPxd5cXH3CCi85q/7fPL/6/6nj9POZXp+7zvjNKS5fX8CObSZuZb9KuXhlmnhskfCSs7K/63PX1+w7dRM2Hk+s+j/AzdjqdSGY34XUcWZ+brmk2uznc8q2zcxjTuf3Z8zvz3ux2y1G6rsJ8jvwK8jjCapHn5r8ndk9nFnex6aA9lPX9uxqIX9mFbCW1ytHPO6RyN+ZX+EBsfj65Vp/27ZIRY55vU3Nxij0kXrA1C0ViJHaE1SPPvW7nOt3oehI0tfXV3pratDsaJQeINKydk55Siq3yoE7O+jbvfb4Z0cLXw0AXLo1jm3PEwmHwcndswKwbJveO+5YeUo1234LydWuhYk424FQJDkFMnXMkfDtnkqhY861j4WJ1VcJ0Yyrhu1tq2/G9nS+fdX/U/vy85nOvRVj/uI4jm0TchzIMxw1vwQd2zrp2fH2lXZm7rtl5vXkuqnplS23T3LtLWtvjgfxu5A6zszPLZdUm/18Ttm2mXnM6fz+jPn9eS92u8VI9eSD/A78KvZ4BgYGcr4WVKCfXv67DRhPW94GxF3XvRXQfupe0LlVsqUUsG3oaAszP13gzcCujtY1qQiClpoCeenqDFcn55bbWD/5ZGKJxMqYfCGelyDm8wa0SKUEFejPL/+9l9Xj9HuBcwHto2Hs3NbKziOtvPd+j2gsTjhL3pdC8qUUgOSQycWbhVPw3rFr/aZXrnpophOO3JE8Md2zva+u8smEbBvLsn0Fe8uyCfm8AS1SKUH9RJ4H3gQ+llpgjAkDHwF+EtA+Go7jWDQ3hYoO8mdHJvnBc8MMjU6tTNlMpRT4/s+HuDa1yOaWUMGiJn17Oir+QJHjWGxujdRNkAewbYutm/19TlvbmtY8eSxSbYH06F3X9YwxfwR83RhzA3ge+DSwHfizIPYhSeOTsznH+CE5ZfLi+AKtTQ7vuruDWXvXypBJNZ6MbRQ7trUwNbOQdx0Lix1ba/umsmxMVqG51sUwxnwGeIpkgD8JfKaUFAgDAwPekSNHAmtXI/mvP3y9YMbG+YUF7t69iWc++dDKMqXgLV0qR0m25Gwp6TeV/eYoCfJGn6yPevqOBgYGOHLkSNbLyUCfUHFd96vAV4PcptxWTFWlN8bmVqUUUAre8qlIidSr2n4UUVYpqqqS55WdUkCSVKRE6p0CfR0pqqqSZdV8SoF6piskqSfqitSRVFUlP/Z1t6qnKSKAAn3defiBnpxTJlNsC+7ft37z40Wkvujavs70drVz7JH9a8oRpliWxa/e10FXR21NnVzPQuEikp8CfR16sL+b7h2bOX7iMoPn1xYYn58arXYT8wq6ULiI5KdAX6d6u9rp7TqYdfbHmRoO9OtVKFxEclOgr3P1NPvDz1O9fgqFi0hxdDNWKqaY4ugiEhwFeqmIUoqji0gwFOilIkopji4iwVCgl4qodHF0EblNgV4qwnEservaC68I7Olu11O9IgHSb5NUzKF9nQWf6rUsa6VKlogEQ4FeKmbnttaCVa+O3teth6ZEAqaBUFkX0ViChaUYzRlpfJXTXaTyFOglUMNj0xw/McqpC5mpGXpvj9Erp7tIRSnQb0C5etvlemFwbE2ytXgiwQn3GifPvcWxR/bzYH/3ymv19FSvSD1ToN9Acve2e3zPiMm37VwZNSH5xOuzPz1P947NZe9LRIqjQL9BFNvbLtbxE6O+0hscP3GZ3q6DJe9HRIqngdENwG9v22+KgkzFFC0fPH9d6Q1EKkyBfgMoprddiqKKliu9gUjFKdA3uEr0tlNFy/1wlN5ApOIU6BtcJXrbxRQt79+/XVMpRSpMv3ENrlK9bT9Fyy3L4r2He0ravoiUToG+wVWqt50qWp4vvcHj79+vqZUiVaDB0g3g4Qd6ePX8W3lvyAbR2y5UtFxBXqQ6FOg3gFRvO9cUyyB72/mKlotIdSjQbxCV7m0rvYFI7VCg30DU2xbZmBToNyD1tkU2FnXnREQanAK9iEiDU6AXEWlwCvQiIg1OgV5EpMEp0IuINLjAA70xps0Yc9EYcyzobYuISPECDfTGmDbgH4A7gtyuiIiULrBAb4x5GPglcH9Q2xQRkfIF2aP/e+AU8KEAtykiImUKMgXCQ67rvmaM6Q1iY2fOnAliMxvS/Pw8oM+w1ul7qn2N8h0VDPTGmDCwL88q467r3nBd97XgmiUiIkHx06PfDeQ7nT0NfC2Y5tzW19eX9/VXRgfXLIvHPZY20cvTAAAEeElEQVRicSIhB8dZXenoHbv7A21fLUv1Pgp9hlJd+p5qXz19RwMDAzlfKxjoXdcdAfIXA62y8clZTg9NMnLlJolEAtu22dPdzj17O9m5rbXazRMRqaq6T1N8dmSS5wfHVlVOSiQSvHF5iqHRaY72d3Ogt6OKLRQRqa66fjJ2fHJ2TZBP53kezw+Oce3GXIVbJiJSO+o60J8emsxb8BqSwf700ESFWiQiUnvqNtDH4x4jV276Wnd47CbRWGKdWyQiUpsCH6Ov1M3bpVicRMJf8E4kkjVSVT5PRDaiuu3RR0IOtu2v+bZt0xSp+/vOIiIlqdtA7zgWvV3tvtbd091OOFS3hyoiUpa6jn6H9nViWflHiSzL4p69nRVqkYhI7anrQL9zWytH+7tzBnvLsjh6X7cemhKRDa3uB64P9HbQsaWZ00MTDI/pyVgRkUx1G+hX5a7ZDY8egmgsObumKRLSmLyIyLK6DfTZhEO2plCKiGRQt1dEpMEp0IuINDirUK6YahgYGKi9RomI1LgjR45knYJYk4FeRESCo6EbEZEGp0AvItLgFOhFRBqcAr2ISINToBcRaXAK9CIiDU6BXkSkwSnQi4g0OAV6EZEGp0AvItLgGipNsSQZY54AfhvoAU4Cz7iu+2J1WyUpxhgHeAp4ArgDuAh8E/iG67rKSVJDjDFNJH+HfuG67m9WuTklU4++wRhjfgP4FvBd4J8DU8CPjDF7qtowSfcF4Cskv6OPAn8LfA34XDUbJVn9LnCg2o0ol3r0DcQYYwFfBv7Sdd0vLy/7MeACTwNPVrF5wkpv/hngT1zX/YPlxT8xxuwAPgv8cdUaJ6sYYw6T/J25Xu22lEs9+sZyF3An8L3UAtd1o8APgQ9Vq1GySjvw18DfZSx3gR3GmE2Vb5JkMsaEgL8C/gQYrXJzyqYefWO5e/nvCxnLh4B9xhjHdd14hdskaVzXvQF8OstLjwGXXdedrXCTJLvfASLAHwK/VuW2lE09+sbSvvz3TMbyGZLftXqLNcgY8yngA2jYpiYYY/qAfw98ynXdpWq3JwgK9I0lVV0m18yNRKUaIv4YYz5B8ub5s8DXq9ycDc8YYwP/GfgvjTRTTUM3jWV6+e82YDxteRsQd133VuWbJLkYY54B/pTkPZVPaGplTfgtklNeP7I8Tp9iGWNCruvGqtSusqhH31jOL/+9N2P5XuBchdsieRhjvgJ8FfgOcKxRhggawK+RfP7kBhBd/nMf8EkgaozprV7TSqcefWM5D7wJfAz43wDGmDDwEZIzb6QGGGOeAj4P/DnwtHryNeXfkLwCTvffSHaUvgyMVbxFAVBx8AZjjPm3JMd6/xB4nuQMj/cA97uuO1TNtgkYY7qAYZKB419nWeWVeh0eaFTGmJPAyXp+MlY9+gbjuu43jTEtJB+xf5rk49sfVJCvGR8EmoB7gWw3+3bQAA/oSG1Rj15EpMHpZqyISINToBcRaXAK9CIiDU6BXkSkwSnQi4g0OAV6EZEGp0AvItLgFOhFRBrc/wepDTbgGJy04wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1147ef2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples_censored=np.copy(samples)\n",
    "samples_censored[20:,1]=0\n",
    "plt.plot(samples[:,0], samples[:,1], 'o', alpha=0.8)\n",
    "plt.plot(samples_censored[:,0], samples_censored[:,1], 's', alpha=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We would use MLE if we had all the data. \n",
    "\n",
    "![](images/gform2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu1 = lambda s: np.mean(s[:,0])\n",
    "mu2 = lambda s: np.mean(s[:,1])\n",
    "s1 = lambda s: np.std(s[:,0])\n",
    "s2 = lambda s: np.std(s[:,1])\n",
    "rho = lambda s: np.mean((s[:,0] - mu1(s))*(s[:,1] - mu2(s)))/(s1(s)*s2(s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we dont. So we shall follow an iterative process to find them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu1s=[]\n",
    "mu2s=[]\n",
    "s1s=[]\n",
    "s2s=[]\n",
    "rhos=[]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bur remember our data are missing in the y-direction. Assume 0 and go. Since we are using the MLE of the **full-data** likelihood, with this assumption we can use the MLE formulae. This is called M-step or maximization step since we used the MLE formulae."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1.7547657491036741],\n",
       " [0.44394554402533365],\n",
       " [1.2718969990037354],\n",
       " [0.77180100755764103],\n",
       " [0.52348014627786521])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu1s.append(mu1(samples_censored))\n",
    "mu2s.append(mu2(samples_censored))\n",
    "s1s.append(s1(samples_censored))\n",
    "s2s.append(s2(samples_censored))\n",
    "rhos.append(rho(samples_censored))\n",
    "mu1s,mu2s,s1s,s2s,rhos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now having new estimates of our parameters due to our (fake) y-data, lets go in the other direction. Using these parameters let us calculate new y-values.\n",
    "\n",
    "One way we might do this is to replace the old-missing-y values with the means of these fixing the parameters of the multi-variate normal and the non-missing data. In other words:\n",
    "\n",
    "$$E_{p(z \\vert \\theta, x)}[z]$$\n",
    "\n",
    "where we have used the notation $z$ to refer to the missing values.\n",
    "\n",
    "This **posterior** distribution (in the sense of bayes theorem, not bayesian analysis) for the multi-variate gaussian is a gaussian..see [wikipedia](https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Conditional_distributions) for the formulae\n",
    "\n",
    "$$\\bar{y}(t+1) - \\hat{\\mu_2}(t) = \\hat{\\rho}(t)\\frac{\\hat{\\sigma_2}(t)}{\\hat{\\sigma_1}(t)} \\left( \\bar{x} - \\hat{\\mu_1}(t) \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ynew(x, mu1, mu2, s1, s2, rho):\n",
    "    return mu2 + rho*(s2/s1)*(x - mu1)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.25149771,  0.71551635,  0.81675223,  0.85921124,  0.08168773,\n",
       "       -0.14712658,  0.7300238 ,  0.42993898,  0.51504132,  0.28024916,\n",
       "       -0.06713977,  0.44793596,  0.54224183,  1.30632573,  0.23563191,\n",
       "        0.5396655 ,  0.25878933,  0.36945243,  0.99094696,  0.53980901])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newys=ynew(samples_censored[20:,0], mu1s[0], mu2s[0], s1s[0], s2s[0], rhos[0])\n",
    "newys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is called the E-step as it computes an expectation for us. Lets run this iteratively and see if we converge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mu1</th>\n",
       "      <th>mu2</th>\n",
       "      <th>rho</th>\n",
       "      <th>s1</th>\n",
       "      <th>s2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.443946</td>\n",
       "      <td>0.523480</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.771801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.673782</td>\n",
       "      <td>0.822228</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.718523</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.792335</td>\n",
       "      <td>0.904408</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.749225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.853302</td>\n",
       "      <td>0.925737</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.775802</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.884575</td>\n",
       "      <td>0.932120</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.790958</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.900583</td>\n",
       "      <td>0.934331</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.798740</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.908762</td>\n",
       "      <td>0.935193</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.802589</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.912935</td>\n",
       "      <td>0.935558</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.804467</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.915062</td>\n",
       "      <td>0.935723</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.805378</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.916144</td>\n",
       "      <td>0.935799</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.805821</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.916695</td>\n",
       "      <td>0.935835</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806036</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.916974</td>\n",
       "      <td>0.935853</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917116</td>\n",
       "      <td>0.935861</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917189</td>\n",
       "      <td>0.935866</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806218</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917225</td>\n",
       "      <td>0.935868</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917244</td>\n",
       "      <td>0.935869</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917253</td>\n",
       "      <td>0.935869</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917258</td>\n",
       "      <td>0.935869</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917260</td>\n",
       "      <td>0.935870</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806242</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>1.754766</td>\n",
       "      <td>0.917261</td>\n",
       "      <td>0.935870</td>\n",
       "      <td>1.271897</td>\n",
       "      <td>0.806242</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mu1       mu2       rho        s1        s2\n",
       "0   1.754766  0.443946  0.523480  1.271897  0.771801\n",
       "1   1.754766  0.673782  0.822228  1.271897  0.718523\n",
       "2   1.754766  0.792335  0.904408  1.271897  0.749225\n",
       "3   1.754766  0.853302  0.925737  1.271897  0.775802\n",
       "4   1.754766  0.884575  0.932120  1.271897  0.790958\n",
       "5   1.754766  0.900583  0.934331  1.271897  0.798740\n",
       "6   1.754766  0.908762  0.935193  1.271897  0.802589\n",
       "7   1.754766  0.912935  0.935558  1.271897  0.804467\n",
       "8   1.754766  0.915062  0.935723  1.271897  0.805378\n",
       "9   1.754766  0.916144  0.935799  1.271897  0.805821\n",
       "10  1.754766  0.916695  0.935835  1.271897  0.806036\n",
       "11  1.754766  0.916974  0.935853  1.271897  0.806141\n",
       "12  1.754766  0.917116  0.935861  1.271897  0.806192\n",
       "13  1.754766  0.917189  0.935866  1.271897  0.806218\n",
       "14  1.754766  0.917225  0.935868  1.271897  0.806230\n",
       "15  1.754766  0.917244  0.935869  1.271897  0.806236\n",
       "16  1.754766  0.917253  0.935869  1.271897  0.806239\n",
       "17  1.754766  0.917258  0.935869  1.271897  0.806241\n",
       "18  1.754766  0.917260  0.935870  1.271897  0.806242\n",
       "19  1.754766  0.917261  0.935870  1.271897  0.806242"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for step in range(1,20):\n",
    "    samples_censored[20:,1] = newys\n",
    "    #M-step\n",
    "    mu1s.append(mu1(samples_censored))\n",
    "    mu2s.append(mu2(samples_censored))\n",
    "    s1s.append(s1(samples_censored))\n",
    "    s2s.append(s2(samples_censored))\n",
    "    rhos.append(rho(samples_censored))\n",
    "    #E-step\n",
    "    newys=ynew(samples_censored[20:,0], mu1s[step], mu2s[step], s1s[step], s2s[step], rhos[step])\n",
    "df=pd.DataFrame.from_dict(dict(mu1=mu1s, mu2=mu2s, s1=s1s, s2=s2s, rho=rhos))\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voila. We converge to stable values of our parameters.\n",
    "\n",
    "But they may not be the ones we seeded the samples with. The Em algorithm is only good upto finding local minima, and a finite sample size also means that the minimum found can be slightly different."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The EM algorithm\n",
    "\n",
    "** Expectation-maximization (EM)** method is an iterative method for maximizing difficult \n",
    "likelihood (or posterior) problems. It was first introduced by Dempster, Laird, and Rubin (1977). \n",
    "\n",
    "EM recognizes that if the data were fully observed, then  ML/ MAP estimates would be easy to compute. It thus alternates between inferring the missing values given the parameters (E step), and then optimizing the parameters given the “filled in” data (M step). \n",
    "The idea is to find a lower-bound \n",
    "on the log-likelihood $\\ell$ (E-step) and the optimize the lower-bound (M-step). \n",
    "\n",
    "**We want to use EM wherever we dont know how to optimize simply the x-data likelihood.**\n",
    "\n",
    "The EM algorithm is naturally useful for missing data, truncated data, censored data and grouped data problems. But it is also useful at places where you have clusters like mixture models, and you are working in an unsupervised or semi-supervised mode. Basically, any place you can set up a latent variable model or a data augmentation procedure, EM is useful.\n",
    "\n",
    "There are applications in astronomical image analysis, genetics, engineering, etc\n",
    "\n",
    "The basic idea is:\n",
    "\n",
    "**calculate MLE estimates for the incomplete data problem by using the complete-data likelihood. To create complete data, augment the observed data with manufactured data**\n",
    "\n",
    "Sorta like, just assign points to clusters to start with and iterate.\n",
    "\n",
    "Then, at each iteration, replace the augmented data by its conditional expectation given current observed data and parameter estimates. (E-step) Maximize the full-data likelihood to do the M-step.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why does EM work?\n",
    "\n",
    "To understand why EM works we will start with our old friend the KL-divergence. All images in this section are stolen from Bishop\n",
    "\n",
    "We cast the problem as a MLE problem, but in general any problem in which there is a $\\theta, x, z$ triad will work.\n",
    "\n",
    "We start with:\n",
    "\n",
    "$$p(x) = \\sum_z p(x, z)$$\n",
    "\n",
    "Specifically, lets cast these as likelihoods\n",
    "\n",
    "$$p(x \\vert \\theta) = \\sum_z p(x, z \\vert \\theta)$$\n",
    "\n",
    "where the $x$ and $z$ range over the multiple points in your data set.\n",
    "\n",
    "Then $$\\ell( x \\vert \\theta) = log\\, p(x \\vert \\theta) = log\\, \\sum_z p(x, z \\vert \\theta)$$\n",
    "\n",
    "Assume $z$ has some normalized distribution:\n",
    "\n",
    "$$z \\sim q(z)$$.\n",
    "\n",
    "We wish to compute conditional expectations of the type:\n",
    "\n",
    "$$E_{p(z \\vert x, \\theta)}[z]$$\n",
    "\n",
    "but we dont know this \"posterior\". Lets however say that we somehow know $q$ and thus can define our loss function as the KL-divergence between $q$ and the posterior $p(z \\vert x, \\theta)$, henceforth abbreviated $p$.\n",
    "\n",
    "Then (notice new notation to be consistent with Bishop):\n",
    "\n",
    "$$ \\mathrm{KL}(q \\vert\\vert p) = D_{KL}(q, p) = E_q[log \\frac{q}{p}] = -E_q[log \\frac{p}{q}]$$\n",
    "\n",
    "$$D_{KL}(q, p) = -E_q[log \\frac{p(x, z \\vert \\theta)}{q\\,p(x \\vert \\theta)}]$$\n",
    "\n",
    "where we have used Bayes Theorem.\n",
    "\n",
    "Now at first blush this does not seem to have bought us anything as we dont know the full data likelihood and we dont know how to optimize the x-data likelihood. But lets persevere:\n",
    "\n",
    "$$D_{KL}(q, p) = - \\left( E_q[log \\frac{p(x, z \\vert \\theta)}{q}] - E_q[log\\,p(x \\vert \\theta)] \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second term does not depend on $q$, and $q$ is normalized so the expectation just gives 1 and we can rewrite the equation in terns of the x-data log likelihoood thus:\n",
    "\n",
    "$$log\\,p(x \\vert \\theta) = E_q[log \\frac{p(x, z \\vert \\theta)}{q}] + D_{KL}(q, p)$$\n",
    "\n",
    "If we define the ELBO or Evidence Lower bound as:\n",
    "\n",
    "$$\\mathcal{L}(q, \\theta) = E_q[log \\frac{p(x, z \\vert \\theta)}{q}]$$\n",
    "\n",
    ", then\n",
    "\n",
    "$log\\,p(x \\vert \\theta)$ = ELBO + KL-divergence\n",
    "\n",
    "a situation made clear in the diagram below where the ELBO goues upto the blue line and the divergence from the blue to the red. Be careful with the p's, the one on the right is for the x-data and the one on the left for the latent-variable posterior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](images/klsplitup.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now recall that the Kullback Liebler divergence is 0 only if the distributions as a function of $z$ are the same at every poiny; it is otherwise **ALWAYS** greater than 0. This tells us that the quantity $\\mathcal{L}(q, \\theta)$, is **ALWAYS** smaller than or equal to the log-likelihood of $p(x  \\vert  \\theta)$, as illustrated above. In other words, $\\mathcal{L}(q, \\theta)$ is a lower bound on the log-likelihood. This is why its called the ELBO, with the \"evidence\" aspect of it coming from Variational calculus (next lecture).\n",
    "\n",
    "The ELBO itself is the expected value of the full-data log-likelihood minus the entropy of $q$:\n",
    "\n",
    "$$\\mathcal{L}(q, \\theta) = E_q[log \\frac{p(x, z \\vert \\theta)}{q}] = E_q[log p(x, z \\vert \\theta)] - E_q[log\\, q]$$\n",
    "\n",
    "### E-step\n",
    "\n",
    "These observations set up the EM algorithm for us. If we choose, in the **E-step**, at some (possibly initial) value of the parameters $\\theta_{old}$,\n",
    "\n",
    "$$q(z) = p(z  \\vert  x, \\theta_{old}),$$ \n",
    "\n",
    "we then set the Kullback Liebler divergence to 0, and thus $\\mathcal{L}(q, \\theta)$ to the log-likelihood at $\\theta_{old}$,  and maximizing the lower bound. \n",
    "\n",
    "Using this missing data posterior, conditioned on observed data, and $\\theta_{old}$, we compute the expectation of the missing data with respect to the posterior and use it later.\n",
    "\n",
    "![](images/klsplitestep.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now the **M-step**. \n",
    "\n",
    "Since after the E-step, the lower bound touches the log-likelihood, any maximization of this ELBO from its current value with respect to $\\theta$ will also “push up” on the likelihood itself. Thus M step guaranteedly modifies the parameters $\\theta$ to increase (or keep same) the likelihood of the observed data.\n",
    "\n",
    "Thus we hold now the distribution $q(z)$ fixed at the hidden variable posterior calculated at $\\theta_{old}$, and maximize $\\mathcal{L}(q, \\theta)$ with respect to $\\theta$ to obtain new parameter values $\\theta_{new}$. This is a regular maximization.\n",
    "\n",
    "The distribution $q$, calculated as it is at $\\theta_{old}$ will not in general equal the new posterior distribution $p(z \\vert x,\\theta_{new})$, and hence there will be a nonzero KL divergence. Thus the increase in the log-likelihood will be greater than the increase in the lower bound $\\mathcal{L}$, as illustrated below.\n",
    "\n",
    "The M in “M-step” and “EM” stands for “maximization”."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](images/klsplitmstep.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The process\n",
    "\n",
    "Note that since $\\mathcal{L}$ is maximized with respect to $\\theta$, one can equivalently maximize the expectation of the full-data log likelihood $\\mathrm{E_q[\\ell( x,z  \\vert  \\theta)]}$ in the M-step since the difference is purely a function of $q$. Furthermore, if the joint distribution $p(x, z \\vert  \\theta)$ is a member of the exponential family, the log-likelihood will have a particularly simple form and will lead to a much simpler maximization than that of the incomple-data log-likelihood $p(x \\vert \\theta)$.\n",
    "\n",
    "We now set $\\theta_{old} = \\theta_{new}$ and repeat the process. This **EM algorithm** is presented and  illustrated below:\n",
    "\n",
    "![](images/emupdate.png)\n",
    "\n",
    "1. We start with the log-likelihood $p(x  \\vert  \\theta)$(red curve) and the initial guess $\\theta_{old}$ of the parameter values\n",
    "2. Until convergence (the $\\theta$ values dont change too much):\n",
    "    1. E-step: Evaluate the hidden variable posterior $q(z, \\theta_{old}) = p(z  \\vert  x, \\theta_{old})$ which gives rise to a lower bound function of $\\theta$: $\\mathcal{L}(q(z, \\theta_{old}), \\theta)$(blue curve) whose value equals the value of $p(x  \\vert  \\theta)$ at $\\theta_{old}$.\n",
    "    2. M-step: maximize the lower bound function with respect to $\\theta$ to get $\\theta_{new}$.\n",
    "    3. Set $\\theta_{old} = \\theta_{new}$\n",
    "    \n",
    "One iteration more is illustrated above, where the subsequent E-step constructs a new lower-bound function that is tangential to the log-likelihood at $\\theta_{new}$, and whose value at $\\theta_{new}$ is higher than the lower bound at $\\theta_{old}$ from the previous step.\n",
    "\n",
    "Thus\n",
    "\n",
    "$$\\ell(\\theta_{t+1}) \\ge \\mathcal{L}(q(z,\\theta_t), \\theta_{t+1}) \\ge \\mathcal{L}(q(z,\\theta_t), \\theta_{t}) = \\ell(\\theta_t)$$\n",
    "\n",
    "The first equality follows since $\\mathcal{L}$ is a lower bound on $\\ell$, the second from the M-step's maximization of $\\mathcal{L}$, and the last from the vanishing of the KL-divergence after the E-step. As a consequence, you **must** observe monotonic increase of the observed-data log likelihood $\\ell$ across iterations. **This is a  powerful debugging tool for your code**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that as shown above, since each EM iteration can only improve the likelihood, you are guaranteeing convergence to a local maximum. Because it **IS** local , you must try some different initial values of $\\theta_{old}$ and take the one that gives you the largest $\\ell$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The EM algorithm with indices laid out\n",
    "\n",
    "I often find it confusing as to what the indices are actually doing if I dont write them out explicitly. So lets visit the EM derivation once more, focussing on mixtures, and explicitly writing out indices. The derivation does not need mixtures, but I find it helpful to imagine that we are fitting such a model.\n",
    "\n",
    "Suppose we have an estimation problem in which we have data consising of $m$ independent examples $\\{x_1,\\ldots,x_m\\}$ . \n",
    "The goal is to fit the parameters of the model, where the log-likelihood is given by \n",
    "\\begin{eqnarray}\n",
    "\\ell(x  \\vert  \\theta)&=& \\log \\prod_{i=1}^{m} p(x_i \\vert  \\theta) =   \\sum_{i=1}^{m} \\log \\,p(x_i \\vert  \\theta)  \\\\ \n",
    "   &=& \\sum_{i=1}^{m} \\log \\zsumi \\,p(x_i,z \\vert  \\theta)  \\\\ \n",
    "\\end{eqnarray}\n",
    "\n",
    "where the $z$ are the latent random variables. If $z$ were observed then the maximum likelihood estimation would be easy. \n",
    "\n",
    "Indeed then, let us start with the full data log-likelihood, \n",
    "\n",
    "$$\\ell(x, z  \\vert  \\theta) = \\sum_{i=1}^{m}  \\log \\,p(x_i, z_i  \\vert  \\theta),$$\n",
    "\n",
    " which is the log-likelihood we'd calculate if we knew all the $z_i$. But we do not know thse, so lets assume the $\\{z_i\\}$ have some normalized distribution $q(z)$, and calculate the expected value of the full data log likelihood with respect to this distribution:\n",
    "\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\mathrm{E_q[\\ell( x,z  \\vert  \\theta)]}  &=& \\sum_i \\zsumi q_{i}(z_i) \\log \\,p(x_i, z_i  \\vert  \\theta)\\\\\n",
    "    &=& \\sum_i \\zsumi q_{i}(z_i) \\log \\,\\frac{p(x_i, z_i  \\vert  \\theta)}{q_{i}(z_i)} +  \\sum_i \\zsumi q_{i}(z_i) \\log \\,q_{i}(z_i)\n",
    "\\end{eqnarray}\n",
    "\n",
    "The second term only involves $q$ and is independent of $\\theta$. Looking only at the first term inside the i-summation:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\mathcal{L}(i, q, \\theta) &=&  \\zsumi q_{i}(z_i) \\log \\,\\frac{p(x_i, z_i  \\vert  \\theta)}{q_{i}(z_i)} \\\\\n",
    "&=& \\zsumi  q_i(z_i) \\left( \\log \\frac{p(z_i \\vert  x_i,  \\theta)}{ q_i(z_i)} + \\log p(x_i  \\vert  \\theta)\\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "we can see that, since $\\zsumi q_i(z_i) = 1$:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\mathcal{L}(i, q, \\theta) &=& \\zsumi  q_i(z_i) \\left( \\log \\frac{p(z_i \\vert  x_i,  \\theta)}{ q_i(z_i)} + \\log p(x_i  \\vert  \\theta)\\right)\\\\\n",
    "    &=& -\\mathrm{KL}\\left(q_i  \\vert  \\vert  p_i \\right) + \\log p(x_i  \\vert  \\theta)\\\\\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $\\mathrm{KL}$ is the Kullback-Leibler divergence between $q(x)$ and the hidden variable posterior distribution $p(z \\vert x,\\theta)$ at the poin $i$.\n",
    "\n",
    "Since the sum over the data-points of the second term is just the log-likelihood we desire, it can then can be written as:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\ell(x  \\vert  \\theta) &=& \\sum_i \\left(\\mathcal{L}(i, q, \\theta) +\\mathrm{KL}\\left(q_i   \\vert  \\vert  p_i \\right)\\right)\\\\\n",
    "&=& \\mathcal{L}(q, \\theta) + \\mathrm{KL}\\left(q  \\vert  \\vert  p \\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "where we are defining:\n",
    "\n",
    "$$\\mathrm{KL}(q  \\vert  \\vert  p) = \\sum_i \\mathrm{KL}\\left(q_i   \\vert  \\vert  p_i \\right)$$\n",
    "\n",
    "as the sum of the KL-divergence at each data point, and $\\mathcal{L}(q, \\theta)$ as the sum of $\\mathcal{L}$ at each data point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Gaussian Mixture model using EM\n",
    "\n",
    "We dont know how to solve for the MLE of the unsupervised problem. The EM algorithm comes to the rescue. As described above here is the algorithm:\n",
    "\n",
    "\n",
    "* Repeat until convergence \n",
    "*  E-step: For each $i,j$ calculate \n",
    "\n",
    "$$ w_{i,j} = q_i(z_i=j)=p(z_i=j \\vert  x_i, \\lambda, \\mu, \\Sigma) $$\n",
    "     \n",
    "* M-step: We need to maximize, with respect to our parameters the\n",
    "  \n",
    "$$\n",
    "\\begin{eqnarray}\n",
    " \\mathcal{L} &=& \\sum_i \\sum_{z_i} q_i(z_i) \\log \\frac{p(x_i,z_i  \\vert \\lambda, \\mu, \\Sigma)}{q_i(z_i)} \\nonumber \\\\\n",
    " \\mathcal{L} &=& \\sum_i \\sum_{j=i}^{k}  q_i(z_i=j) \\log \\frac{p(x_i \\vert z_i=j , \\mu, \\Sigma) p(z_i=j \\vert \\lambda)}{q_i(z_i=j)} \\\\\n",
    " \\mathcal{L} & =&  \\sum_{i=1}^{m} \\sum_{j=i}^{k} w_{i,j}  \\log \\left[   \\frac{ \\frac{1}{ (2\\pi)^{n/2} \\vert \\Sigma_j \\vert ^{1/2}} \\exp \\left(    -\\frac{1}{2}(x_i-\\mu_j)^T \\Sigma_j^{-1} (x_i-\\mu_j) \\right)  \\, \\lambda_j   }{w_{i,j}}\\right]\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Taking the derivatives yields the following updating formulas:\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    " \\lambda_j &=& \\frac{1}{m} \\sum_{i=1}^m w_{i,j} \\nonumber \\\\ \n",
    " \\mu_j&=& \\frac{ \\sum_{i=1}^m  w_{i,j} \\, x_i}{ \\sum_{i=1}^m  w_{i,j}} \\nonumber \\\\ \n",
    " \\Sigma_j &=& \\frac{ \\sum_{i=1}^m  w_{i,j} \\, (x_i-\\mu_j)(x_i-\\mu_j)^T}{ \\sum_{i=1}^m  w_{i,j}}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "To calculate the E-step we basically calculating the posterior of the  $z$'s given the $x$'s and the\n",
    "current estimate of our parameters. We can use Bayes rule \n",
    "\n",
    "$$ w_{i,j}= p(z_i=j \\vert  x_i, \\lambda, \\mu, \\Sigma) = \\frac{p( x_i \\vert  z_i=j,  \\mu, \\Sigma)\\, p(z_i=j \\vert \\lambda)}{\\sum_{l=1}^k p(x_i  \\vert  z_i=l,  \\mu, \\Sigma) \\, p(z_i=l \\vert \\lambda)} $$\n",
    "\n",
    "Where $p(x_i  \\vert  z_i =j,  \\mu, \\Sigma)$ is the density of the Gaussian with mean $\\mu_j$ and covariance \n",
    "$\\Sigma_j$ at $x_i$ and $p(z_i=j \\vert  \\lambda)$ is simply $\\lambda_j$. \n",
    "If we to compare these formulas in the M-step with the ones we found in GDA we can see\n",
    "that are very similar except that instead of using $\\delta$ functions we use the $w$'s. Thus the EM algorithm corresponds here to a weighted maximum likelihood and the weights are interpreted as the 'probability' of coming from that Gaussian instead of the deterministic \n",
    "$\\delta$ functions. Thus we have achived a **soft clustering** (as opposed to k-means in the unsupervised case and classification in the supervised case)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAECCAYAAAAciLtvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAEeRJREFUeJzt3X2U3FV9x/H3ykNMSfBwJB6fagnWfF1agppja6WttkABaQiKnhaxooeIaMVUEiE8SSQFEYlVQIoaPZDiqQ8UhViLQLQWRHpwPYjW5RsURBRqpQYJccWQbP/4/bYOw27Ymbk7k519v87ZM+y9v7vzzewwn7339zQwOjqKJGlme0qvC5Ak9Z5hIEkyDCRJhoEkCcNAkoRhIEkCdu11Ae0YGhryeFhJasOiRYsGxmuflmEAsGjRopbHDA8PAzA4OFi6HDXxte4uX+/umq6v99DQ0IR9LhNJkgwDSZJhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJIlpfAayumPx8ms6Gr9+zfQ6Q1OaqZwZSJIMA0mSYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSSJNs5AjogjgU9l5tyGttnAmcBfAc8E7gLOz8zPNGxzNHDVOD/ypMy8pNU6JEnltBQGEfFy4EpgoKnrH4GjqALhTuBI4NMRMZqZn623OQD4PvA3TWPvabVoSVJZkwqDiJgFLANWA1uA3Rv6ngEcByzNzE/UzTdGxPOBFcBYGCwEhjLz1kK1S5IKmew+g8OB04B3Axc39c0BLgOub2pPYH7D9wuBO9qoUZI0xSa7THQbMD8zH4qIVY0dmXk38LbGtojYhSpA7qy/nwvsA7w4IjZShcQwsDIzv9TJP0CS1LlJhUFm/qTFn/te4IVU+w4A9qfazzAfOBl4DHg7sD4iDs7Mr7b48xkeHm51CCMjI22PVXt8rbvD93Z39ePrXfx+BhFxKnAGsCYz19fN3wOOAG7OzIfr7W4Avk2107nlMND0cMrajW2PvWDpgoKVSNqRYmEQEQPAGuBdwKVU+xcAyMyHgMctB2XmtjoQmo8umpTBwdZvmjKW4u2Mnbna/zDvlL+nyfO93V3T9fUeGhqasK9IGETEU4ArgDcA52XmGU39LwYWZebapqGzgQdL1CBJal+pM5DXUAXB8uYgqL0I+HgdCsD/n6j2KuBrhWqQJLWp45lBRLyE6hyEG4BbIuJlDd3bMvM24HNUh6Z+LiLOAEaolpHmAH/faQ2SpM6UWCY6kupIoUPqr0ZbgDmZ+UhEHARcAFxEFQI3A3+amfcVqEGS1IGWwyAzVwGrJvp+B+PuA45p9fkkSVPPq5ZKkgwDSZJhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCdi11QERcSTwqcyc29A2AJwOvBXYG/g6cFJm3tmwzSzgfOAYYA/gy8A7M/P+jv4FkqSOtTQziIiXA1cCA01d7wHOBC4E/hp4GrAhIp7WsM1lwBuBlcCbgQOAL0XELu2VLkkqZVIzg/qv+mXAamALsHtD31xgBbAqMy+q224C7gWOBz4YEc+nCoLXZ+Zn6m2+DSSwBLi61D9IktS6yc4MDgdOA94NXNzU9zJgDnDtWENmbgK+BhxWN/15/fjFhm3uAv6rYRtJUo9MNgxuA+bXf/mPNvUtqB9/0NR+d0PfAuC/M3PLDraRJPXIpJaJMvMnO+jeE3g0M3/d1L657hvbZvM4YzcDvz2ZGpoNDw+3PGZkZKTtseo+f0+T53u7u/rx9S5xaOkAT5wtjNnewjaSpB5p+dDScfwCmBURu2Xm1ob2uXXf2DZznzDy8du0ZHBwsOUxYyneztiZa2PPntnf0+T53u6u6fp6Dw0NTdhXYmZwF9Vf/vOb2velOlpobJtnRsTsHWwjSeqREmFwC/Ar4KixhojYC3gFsKFu2gDsAixu2OYFwO81bCNJ6pGOl4ky85GIuBhYHRHbqdYVzgAeBtbW2/wgIj4HfLw+EW0T8D7gDuALndYgSepMiX0GUF2KYjvVyWdzqGYLx2Vm4/6ANwP/ALyfakZyI9XlKLYVqkGS1KaWwyAzVwGrmtoeo7rMxModjNsCnFB/SZJ2Il61VJJkGEiSDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJlLschVTc4uXXtD12/ZolBSuR+p8zA0mSYSBJcplIUoc6Wc4Dl/R2Fs4MJEmGgSTJZSJJPeZRYzsHZwaSJMNAkuQy0YzQ6dEekvqfMwNJkmEgSTIMJEkYBpIkDANJEoaBJIlCh5ZGxCuBr+5gk32AvYFvjtO3JjNXlKhDktSeUucZfAv4o6a2pwJXAUPAfcBBwBbg4Kbt7i9UgySpTUXCIDMfBm5tbIuIDwGjwBsyc3tELAS+m5m3jvczJEm9MyVnIEfEfsA7gL/NzJ/VzQuBO6bi+SRJnZmqy1GcC2wEPt7Qtj/waETcDuwH/AhYnZlXTFENkqRJKh4GEbEvcCRwQmZur9ueTbUD+QXAacAm4Bjg8ogYzcx1rT7P8PBwy7WNjIy0PVbTy0z7Hc/U93av/r39+HpPxcxgKdWH/ZUNbZuAQ4HvZOYDdduNdUicDbQcBpKkcqYiDI4CvpCZj441ZOYIcP04214HHBYRczLzkVaeZHBwsOXCxlK8nbHT28ZeF9B1M+133Nv3du/eX736PU/Xz5KhoaEJ+4qGQUQ8DxgEVjS1L6A6tPSTjSEBzAZGqA45lST1SOkzkP+gfmw+fPQ5wKXAq8YaImIAeA1wU2aOFq5DktSC0stEvw88mJk/b2r/D+Bm4LKI2At4ADiB6nDTAwvXIElqUemZwTOAh5obM3MbsAT4PHAOcHW97SGZOfEiliSpK4rODDLz7Tvo+zlwYsnnkySV4VVLJUmGgSTJMJAkYRhIkjAMJEkYBpIkDANJElN3PwNJmnKLl1/T9tj1a5YUrGT6c2YgSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEgVvexkRTwceHKfrXzLztRExAJwOvBXYG/g6cFJm3lmqBklSe0reA/mA+vEvgM0N7f9bP74HWAmcCvwQOBPYEBH7ZeYvCtYhqUWd3EtY/aFkGCwEfpqZNzR3RMRcYAWwKjMvqttuAu4Fjgc+WLAOSVKLSu4zWAjcMUHfy4A5wLVjDZm5CfgacFjBGiRJbSg9M/hVRNwCvIRq/8GHgQuBBfU2P2gaczewpJ0nGx4ebnnMyMhI22M1vXSy7HHB0gVPvtFOxvd26zp5rfrx9S4SBhGxC7AfsIVqOehe4AjgfGA2sBV4NDN/3TR0M7BniRokSe0rOTP4S+BHmfn9+vt/j4g5VDuMzwVGJxi3vZ0nGxwcbHnMWIq3M3Z629jrAqaV6fj+6Py9PfPeI538nqfrZ8nQ0NCEfUXCIDO3AV8Zp+s64ESqGcOsiNgtM7c29M8FPJJIknqs1DLRs6lmBp/PzJ81dM2uHzcBA8B8Hv8nyL5AlqhBktS+UkcTzQI+Cryhqf1oqg//q4FfAUeNdUTEXsArgA2FapAktanUMtE9EfHPwOqI2A4MA6+jCoOjMvORiLi4oX8jcAbwMLC2RA2SpPaV3IF8PHAW8HfAs6gC4ejMHDu34HSqncUrqM45uAU4zrOPJan3ioVBZo5QfeCfPkH/Y1SXo1hZ6jlnEi8XIGkqedVSSZJhIEkyDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkibL3M5DUQ6es3chMvLG9ynBmIEkyDCRJhoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEgUvRxERuwDLgLcAzwPuBS4FPpKZoxGxCPjmOEPXZOaKUnVIklpX8tpEZwErgdXArcCfAB8Cfgu4ADgA2AIc3DTu/oI1SNKkLF5+TdtjL1i6oGAlO4ciYVDPCk4GPpCZ59bNGyJiHrCCKgwWAt/NzFtLPKckqZxS+wz2BNYBVze1JzAvIvagCoM7Cj2fJKmgIjODzNwEvGOcrsXAjzNzS0TsDzwaEbcD+wE/AlZn5hUlapBK6WT5YP2aJQUrkbpnyu5nEBFLqfYPvDMing3sDbwAOA3YBBwDXB4Ro5m5rtWfPzw83HJNIyMjbY+VJqNX762x97a6ox8/S6YkDCLiWOAy4CrgEuCpwKHAdzLzgXqzG+uQOJtqiUmS1CPFwyAiTgYuBK4Fjs3MUWAEuH6cza8DDouIOZn5SCvPMzg42HJtYyneztje8w5W00Gv3lv99BfqdHD2p+5re2wvlxKHhoYm7CsaBhFxHtUy0Drg+Mx8rG5fABwEfDIzH20YMpsqKLaUrEOS1JpiZyBHxDKqIPgw8KaxIKg9h+oEtFc1bD8AvAa4qZ49SJJ6pNR5Bs8C3g98B/g08IcR0bjJLcDNwGURsRfwAHAC1eGmB5aoQZLUvlLLRIcCs4D9gW+M0z8PWAKcB5wDPB34FnBIZk68iCVJ6opS5xlcDlw+iU1PLPF8kqSyvGqpJMkwkCRN4RnIeqJOLnMgSVPJMJB2Iv7BoF5xmUiSZBhIkgwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiS8HIVUlJeT0HTlzECSZBhIkgwDSRKGgSQJw0CShEcTSVJXdXrE2fo1SwpV8njODCRJhoEkyTCQJGEYSJLowQ7kiHgLcArwXOB24OTM/Ea365Ak/UZXwyAijgMuA84BbgNOAr4cEQdk5j3dqOGUtRuBjd14KkmaNrq2TBQRA8B7gY9l5nsz80vAkcCDwLu6VYck6Ym6uc/gd4HfAa4da8jMrcC/Aod1sQ5JUpNuhsGC+vH7Te13A8+PiF26WIskqUE39xnsWT9ubmrfTBVKewAPT/aHDQ8Pt1zAyMhIy2MkaWfSzmffZHQzDAbqx9EJ+re38sN++ctftlXEqtc/t61xkrQzaPez78l0Mwx+UT/OBX7a0D4X2JaZj0z2By1atGjgybeSJE1WN/cZ3FU/7tvUvi8e6ylJPdXtMLgPOGqsISJ2A44ANnSxDklSk4HR0YmW8MuLiLcDlwDvA74OvAP4Y+BFmXl31wqRJD1OV8MAICKWA8uAvakuR7Hcy1FIUm91PQwkSTsfr1oqSTIMJEmGgSQJw0CSRA9ubtMr3lSne+qLDi4D3gI8D7gXuBT4SGZ6xMIUiYhZVO/t/8zMN/W4nL4VEQcB5wELgf8BLgfOycxtvayrUzNiZtBwU50rgaOBh6huqjO/p4X1r7Oo/me5kuqeFZ8FPgS8u5dFzQBnAy/sdRH9LCIOBP4NGKY6YfYS4FTgzF7WVULfzwyab6pTt90AJNVNdd7Zw/L6Tj0rOBn4QGaeWzdviIh5wArggp4V18ci4sVU7+UHe11LnzsfuL5h5vWViHg68GdUnzPT1kyYGXhTne7aE1gHXN3UnsC8iNij+yX1t4jYFfgk8AHgJz0up2/Vf9AcCHyssT0zV2bmK3tSVEF9PzNgEjfVme5rfTuTzNxEdZmRZouBH2fmli6XNBOcCuxOdZmXV/e4ln62P9Wl+LdExHrgEKp7sFxKtc+gpcvw72xmwsxgMjfV0RSKiKXAwbhEVFxEDAJnAEsz89e9rqfPzasf1wF3AodTBcGZ9MH+sJkwMyh6Ux21JiKOpdp5fxXVzjYVEhFPAdYCn/DIuK7YrX78cmaOffh/NSL2Bs6MiAun8yrDTJgZNN5Up1HLN9VRayLiZOCfgC8Cx3pYaXEnUR26e1ZE7FrvOwAYaPhvlTP2WXFdU/sNwBxgn65WU9hMCANvqtMDEXEesIYqDF7rEsaUeDXVeTObgK311wHAG4GtEbFP70rrS2P7HXdvah+bMUzrP3ZmShh4U50uiohlwGnAh4E3ZeZjPS6pX70VeGnT10aqmdhLgft7V1pf+h7V0Vqva2o/guq1/mG3CyppRlzC2pvqdE9EPAu4h+pD6YRxNvmm4TB1IuJ24HbPQJ4aEfFG4Ap+sx/sYKqjud6WmR/tZW2dmhHripl5aUTMprpEwruoTtk/1CCYEocCs6gOwxtvp+Y8PDFK01RmrouIrcDpwJupVh1OzMyP7Xjkzm9GzAwkSTs2E/YZSJKehGEgSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkScD/AVzF0dyO+i+rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x117c46e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#In 1-D\n",
    "# True parameter values\n",
    "mu_true = [2, 5]\n",
    "sigma_true = [0.6, 0.6]\n",
    "lambda_true = .4\n",
    "n = 1000\n",
    "\n",
    "# Simulate from each distribution according to mixing proportion psi\n",
    "z = np.random.binomial(1, lambda_true, n)\n",
    "x = np.array([np.random.normal(mu_true[i], sigma_true[i]) for i in z])\n",
    "\n",
    "plt.hist(x, bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from Bios366 lecture notes\n",
    "from scipy.stats.distributions import norm\n",
    "\n",
    "def Estep(x, mu, sigma, lam):\n",
    "    a = lam * norm.pdf(x, mu[0], sigma[0])\n",
    "    b = (1. - lam) * norm.pdf(x, mu[1], sigma[1])\n",
    "    return b / (a + b)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Mstep(x, w):\n",
    "    lam = np.mean(1.-w) \n",
    "    \n",
    "    mu = [np.sum((1-w) * x)/np.sum(1-w), np.sum(w * x)/np.sum(w)]\n",
    "    \n",
    "    sigma = [np.sqrt(np.sum((1-w) * (x - mu[0])**2)/np.sum(1-w)), \n",
    "             np.sqrt(np.sum(w * (x - mu[1])**2)/np.sum(w))]\n",
    "    \n",
    "    return mu, sigma, lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4 [2, 5] [0.6, 0.6]\n",
      "Initials, mu: [  0.72500527 -20.77035111]\n",
      "Initials, sigma: [ 4.59386658  3.6262629 ]\n",
      "Initials, lam: 0.6261111131564271\n",
      "Iterations 95\n",
      "A: N(5.0083, 0.6288)\n",
      "B: N(2.0261, 0.5936)\n",
      "lam: 0.4116\n"
     ]
    }
   ],
   "source": [
    "print(lambda_true, mu_true, sigma_true)\n",
    "# Initialize values\n",
    "mu = np.random.normal(4, 10, size=2)\n",
    "sigma = np.random.uniform(0, 5, size=2)\n",
    "lam = np.random.random()\n",
    "print(\"Initials, mu:\", mu)\n",
    "print(\"Initials, sigma:\", sigma)\n",
    "print(\"Initials, lam:\", lam)\n",
    "\n",
    "# Stopping criterion\n",
    "crit = 1e-15\n",
    "\n",
    "# Convergence flag\n",
    "converged = False\n",
    "\n",
    "# Loop until converged\n",
    "iterations=1\n",
    "\n",
    "\n",
    "while not converged:\n",
    "    # E-step\n",
    "    if np.isnan(mu[0]) or np.isnan(mu[1]) or np.isnan(sigma[0]) or np.isnan(sigma[1]):\n",
    "        print(\"Singularity!\")\n",
    "        break\n",
    "        \n",
    "    w = Estep(x, mu, sigma, lam)\n",
    "\n",
    "    # M-step\n",
    "    mu_new, sigma_new, lam_new = Mstep(x, w)\n",
    "    \n",
    "    # Check convergence\n",
    "    converged = ((np.abs(lam_new - lam) < crit) \n",
    "                 & np.all(np.abs((np.array(mu_new) - np.array(mu)) < crit))\n",
    "                 & np.all(np.abs((np.array(sigma_new) - np.array(sigma)) < crit)))\n",
    "    mu, sigma, lam = mu_new, sigma_new, lam_new\n",
    "    iterations +=1           \n",
    "\n",
    "print(\"Iterations\", iterations)\n",
    "print('A: N({0:.4f}, {1:.4f})\\nB: N({2:.4f}, {3:.4f})\\nlam: {4:.4f}'.format(\n",
    "                        mu_new[0], sigma_new[0], mu_new[1], sigma_new[1], lam_new))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why is EM important?\n",
    "\n",
    "We have motivated the EM algorithm using mixture models and missing data, but that is not its only place of use. \n",
    "\n",
    "Since MLE's can overfit, we often prefer to use MAP estimation. EM is a perfectly reasonable method for MAP estimation in mixture models; you just need to multiply in the prior.\n",
    "\n",
    "Basically the EM algorithm has a similar setup to the data augmentation problem and can be used in any problem which has a similar structure. Suppose for example you have two parameters $\\phi$ and $\\gamma$ in a posterior estimation, with daya $y$. Say that we'd like to estimate the posterior $p(\\phi  \\vert  y)$. It may be relatively hard to estimate this, but suppose we can  work with $p(\\phi  \\vert  \\gamma, y)$ and $p(\\gamma  \\vert  \\phi, y)$. Then you can use the structure of the EM algorithm to estimate the marginal posterior of any one parameter. Start with:\n",
    "\n",
    "$$log p(\\phi  \\vert  y) = log p(\\gamma, \\phi  \\vert  y) - log p(\\gamma  \\vert  \\phi, y)$$\n",
    "\n",
    "Notice the similarity of this to the above expressions with $\\phi$ as $x$, $y$ as $\\theta$, and $\\gamma$ as $z$. Thus the same derivations apply toany problem with this structure.\n",
    "\n",
    "This structure can also be used in type-2 likelihood or emprical bayes estimation."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}