Population Level Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "import pickle\n",
    "import numpy as np\n",
    "import math\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "from EM import EM\n",
    "from plot import plot\n",
    "from preprocess import preprocess\n",
    "#from EM_individual_params import EM "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling choices (which generates the lowest MSE and the plots in write up)\n",
    "* Cutoff (only consider the data point with at least certain number of observations): 5 --> 2214 data points, 215 bins\n",
    "* Bin size (length of time interval used to bin the time series): 18 hrs\n",
    "* Parameters initialization\n",
    "    * Both treatment and static coefficients: Gaussian with mean 0.1 plus small white noise\n",
    "* Missingness\n",
    "    * After binning the data, only consider the data points with less than 40% missing observations\n",
    "* Chronic conditions and age\n",
    "    * Set to zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coefficients learned with best population level results\n",
    "* Coefficient for nsaid: [0.31291354 0.11150224 0.02927106]\n",
    "* Coefficient for transfusion_plasma: [0.21646024 0.12142273 0.07316288]\n",
    "* Coefficient for transfusion_platelet: [-0.24024057 -0.41378184 -0.28599656]\n",
    "* Coefficient for anticoagulant: [-0.04605939  0.15568198  0.06412978]\n",
    "* Coefficient for aspirin: [0.25202486 0.1102651  0.03499285]\n",
    "\n",
    "Coefficients (jointly) learned with best population level results\n",
    "* Coefficient for nsaid: [ 0.09144572  0.01808053 -0.00736706]\n",
    "* Coefficient for transfusion_plasma: [-0.29370552 -0.15658492 -0.12820925]\n",
    "* Coefficient for transfusion_platelet: [-0.32318352 -0.44511952 -0.43787362]\n",
    "* Coefficient for anticoagulant: [-0.46557922  0.00142221 -0.02525039]\n",
    "* Coefficient for aspirin: [-0.04017403 -0.04024049 -0.00733096]\n",
    "* Coefficient for chronic kidney failure: 0.27480724324859296\n",
    "* Coefficient for sickle cell: 0.238516267304465\n",
    "* Coefficient for age: 0.01752724784756275"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preprocess Setting\n",
    "bin_size = 60 * 18\n",
    "cutoff = 5\n",
    "missing_pct = .4\n",
    "c_zero = False\n",
    "\n",
    "# EM Setting\n",
    "num_past_effects = 3\n",
    "training_pct = .8\n",
    "single_effect = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pickle.load(open('../Data/unimputed_inr_patient_data.pkl', 'rb'), encoding='latin1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "y, X, c = preprocess(data, cutoff, bin_size, missing_pct=missing_pct, c_zero=c_zero)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2214, 215)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reasons for Errors\n",
    "* Assume the a single set of parameters generate all the data. This assumption could be violated since:\n",
    "    1. Each patient respond to treatment differently\n",
    "    2. Different treatment and chronic conditions may have various interaction effects on different individuals\n",
    "* After changing treatment coefficients to be on individual level, the average coefficients are similar to the case of population level coefficients. The fit on some plots improve slightly but not much.  \n",
    "* Assume observed value is well explained by the treatment and chronic variable. whereas in reality there might be other variable that contributes to the observed\n",
    "* Assume linear transition, emission and treatment/chronic effects\n",
    "* Seeing from the data, treatments don't seem to show clear effects on INR, from plots and from calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "em = EM(y, X, c, num_past_effects, 0, train_pct=training_pct, single_effect=single_effect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max iterations: 500 reached\n",
      "CPU times: user 25min 47s, sys: 87.7 ms, total: 25min 47s\n",
      "Wall time: 25min 47s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "em.run_EM(500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.01714340998109"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "em.get_MSE()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient for nsaid: [-0.07408245  0.01390443 -0.02136621]\n",
      "Coefficient for transfusion_plasma: [-0.54602964 -0.20416589 -0.23492377]\n",
      "Coefficient for transfusion_platelet: [-0.78690423 -1.07778153 -0.73084188]\n",
      "Coefficient for anticoagulant: [-0.79181835 -0.1506899  -0.06266534]\n",
      "Coefficient for aspirin: [-0.13027704 -0.1326261  -0.05292395]\n",
      "Coefficient for chronic kidney failure: -1.3121359884570776\n",
      "Coefficient for sickle cell: -2.414278413676523\n",
      "Coefficient for age: -0.1025557457370168\n"
     ]
    }
   ],
   "source": [
    "treatment_types = ['nsaid', 'transfusion_plasma', 'transfusion_platelet', 'anticoagulant', 'aspirin']\n",
    "for i, treatment in enumerate(treatment_types):\n",
    "    print('Coefficient for {}: {}'.format(treatment, em.A[:, i]))\n",
    "static_types = ['chronic kidney failure', 'sickle cell', 'age']\n",
    "for j, static in enumerate(static_types):\n",
    "    print('Coefficient for {}: {}'.format(static, em.b[j]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "random_patients = np.random.randint(0, y.shape[0], 0)\n",
    "for i in random_patients:\n",
    "    plot(em, i, bin_size=bin_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEWCAYAAAC9qEq5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYHVWd//H3Jx0SkCVsQQIJJkoAA4NRWsRtZBMCgwQVMAxK1PwG5QE3nFEYdHBjBtwyosCYMZGASsiATKIGQ1iU0ZElLAJhS8syaYkQSAgIJCSd7++POre70rn31k33vbfT3Z/X89znVp06VXWq06lvn1OnzlFEYGZm1gxD+roAZmY2eDjomJlZ0zjomJlZ0zjomJlZ0zjomJlZ0zjomJlZ0zjomPWQpEMltTfguGMlhaShaf16SVPT8kcl/a6Hxw1Je6fl/5D05bTckOuoUIYnJB3ZjHPZlslBxwaFdLN7RdJfc58fpG0fTTfk73bb54SUfnmfFDqJiGMiYnadj/nJiPh6PY9pVgsHHRtM3hcR2+U+Z+W2/Qn4UKl2kZwGPNrcIpoNbA46Zpm/APcDRwNI2hl4BzC/1gNIeqOk30h6XtISScfntu0i6ReSXpB0p6Rv1NpMlo75/yps+5ak30kakdY/LukhSaskLZT0ugr7XS7pG93SPi/pGUnLJX0slz5C0hWSVkh6UtKXJA1J24ak9SfTvleUypK2fyRte07SebVcrw1sDjpmXa4gq90ATAHmAWtr2VHSVsAvgBuA3YBPAT+VtG/KcgnwErA7MDV9eizd7P8TOBA4KiJWSzoB+GfgA8BI4H+Aq2o85O7ACGBPYBpwiaSd0rbvp22vB95D9jMqBaWPps9haft2QKnZcgJwGfARYA9gF2B0jy7YBgwHHRtM/jvVQkqff+i2/Trg0PSX+mlkQahWh5DdcC+MiFcj4mbgl8ApklqADwLnR8TLEfEg0JtnNFuRBZOdyZoMX07pnwD+LSIeioj1wL8CEyvVdrpZB3wtItZFxALgr8C+qewfAs6NiBcj4gngO2SBBOBU4LsR8VhE/BU4F5iSmilPBH4ZEbdGxFrgy8CGXly3DQBDi7OYDRgnRMSNlTZGxCuSfgV8Cdg1In4v6Zgaj70HsCwi8jfVJ8lqDiPJ/q8ty23LL2+uvYE3AQdHxKu59NcB35P0nVyaUhmeLDjmcylQlbxMFkR3BYZ12790XZBdd/dtQ4HXpm2d1xkRL0l6rqAcNsC5pmO2sSuAzwNXbuZ+TwFjSs86kr2APwMrgPVs3LQ0phdlfIiseev6XPMdZDf4T0TEjrnPNhHxv70417NktaB8bal0XZBdd/dt64GngeXkrlPSa8ia2GwQc9Ax29hvgfeSPcfYHLeTPbP5gqStJB0KvA+YExEdwM+Br0h6jaT96Hp21CMRcRXZ85sbJb0hJf8HcK6k/aGzA8BJvTxPBzAXuEDS9qmp7mzgJynLVcDnJI2TtB1Zk97VqdZ0DXCcpHdJGgZ8Dd9zBj3/Athg8otu7+lc1z1DZG6KiJWbc+DUzHU8cAxZ7eBS4LSIeDhlOYvsYfxfyGpRV1FjJ4Uq55xNdiO/WdLYiLgOuAiYI+kF4IFUnt76FFlAfQz4HfAzYFbaNovsem4FHgfWpPxExBLgzJR/ObAKaMpLqLblkidxM2s+SRcBu0dEr3qxmfU3rumYNYGk/SQdqMzBZN2SN6lpmQ10fRJ0JJ2UXp7bIKm127ZzJbVJekTS0bn0SSmtTdI5ufRxkm6XtFTS1antGEnD03pb2j62WddnVsb2ZM91XiJ7RvIdsveAzAaVvqrpPED2Atut+cT0MtkUYH9gEnCppJb0rsAlZO3TE8jefZiQdrsImB4R48najKel9GnAqojYG5ie8pn1iYi4MyL2jojXRMTYiPi3cNu2DUJ9EnTSy2uPlNk0may3z9qIeBxoAw5On7b0AtqrwBxgsiQBh5P1koHshbsTcscqvYB3DXBEym9mZn1kS3s5dE/gttx6O10voS3rlv42sj7/z+deasvn37O0T0Ssl7Q65X+2+0klnQ6cDrDtttsetN9++9XlYszMBou77rrr2YgYWZSvYUFH0o1k4zl1d15EVGrLLlcTCcrXyKJK/mrH2jQxYgYwA6C1tTUWL15coXhmZlaOpKJRL4AGBp2I6MlETe1s/Kb2aLI3nqmQ/iywo6ShqbaTz186VnsaB2oEsFnvXpiZWX1taV2m55MNFjhc0jhgPHAHcCcwPvVUG0bW2WB+ehB7C9nAgpCN3Dsvd6zSOxAnAjf7wa2ZWd/qqy7T71c2Pe7bgV9JWgidbzDPBR4Efg2cGREdqRZzFrCQbNypuSkvwBeBsyW1kT2zmZnSZwK7pPSzgc5u1mZm1jc8IkE3fqZjZrb5JN0VEa1F+ba05jUzMxvAHHTMzKxpHHTMzKxpHHTq6KaHnuap51/p62KYmW2xHHTqaNrsxRw1/dbijGZmg5SDTp39de364kxmZoOUg46ZmTXNljbgp5kNQhFBx4ZgQ8CGCCJ9Z59s+4aAjg3RuVza3pU3fW+ovj07V/XtGzZAR3Q/V6kc2fZ8+Tba3u06OvLn2dDtPGm5o2D7hs5t3crR7edV+hlW2t7959n9eF89fn8Oet3ODf23dtBpgPUdGxja4krklip/A+v6j9p1oynduDpy/9k7NkTnTahjQ9e++f/onctpvaPzBpQ7R3Tfl9y5ovNcG7rdRErLpZtuR9q3dLMu3bQ2LufGN9dNrrNUttw1dr+m0rlKN9n8zSt/s940GGx64+t+M87fBAcjCYZIDBFIoiUtD5GybUO00fYhghYpWx5S2lcbHaeUVtqu/DEpHROGDBlSdt9hLS0Nv24HnQZYvnoNY3Z+TV2Pub5jA692bODV9V3f6zoifW9gff6GtGHjG0VH7j/8Jjfa/F99G8rnzefpunFVuimXv/F13Sir3/gq3aA7b77db9D56+y+X+exN/5ZxAC4ybUMSTepdHPJltONqfNmpWx5COmGluVpSTebliHZ9uyGR+fy0CFDGD40yyvyx6PzhiZtnCblz8FGN7v8zS07B932zefd+GZbbrtK16Xy5ei69urb8zfu/I09X+b8tZe9piFdwWKT7UPKBINcOQbr9F4OOg1w69IVnPq21xXmiwiefmEtf1rxV556/hWWr17DU8+/wtMvrGH1K+t4Yc16XnhlHS+sWceadRuaUPLNU/qPX7qZdL+Zdb/xdb85du5X4cZXukF13WAr3FSVSy/lyd1gSzexliF026/rplGuDOWuqfQX6SbH6nbe/Pk2LW/+ZpWCQue5cudNwaLz55Arp1l/5aDTAHPuWMZJB41h2NCuJraXX13Po0//lYeXv8DDf3mRh9L36lfWbbTvrtsNZ/cRw9npNcMYNWIbdthmKDtsvRXbDh/K8KFDGDZ0CFu1DGFYS9fyVi1iaEv3m3Dphk7uRlnuJt7tL+YhBTdH3/jMrBccdOrs9SO35f4/r+awb/+G/ffYgVfWdfD4sy/RvqrrpdFth7Ww7+7bc+zfjOKNo7Zn7922Y88dt2H3EVszfGjj21TNzPqKg06dHXfgHkwYtQPX3LWMx599iW2GtfDmvXbiQ61j2Gf37Xnj7jsweqdtGDLENQUzG3wcdOqkNEWEgEkH7M6kA8rN1G1mNri5X2+dlHpE+VGHmVllfTVz6LckPSzpPknXSdoxt+1cSW2SHpF0dC59Ukprk3ROLn2cpNslLZV0dZrOmjTl9dUp/+2Sxjbymkq9cLNOpmZmVk5f1XQWAQdExIHAo8C5AJImAFOA/YFJwKWSWiS1AJcAxwATgFNSXoCLgOkRMR5YBUxL6dOAVRGxNzA95Ws413TMzCrrk6ATETdERGlkzNuA0Wl5MjAnItZGxONAG3Bw+rRFxGMR8SowB5isrN/u4cA1af/ZwAm5Y81Oy9cAR6iB/Xw97beZWbEt4ZnOx4Hr0/KewLLctvaUVil9F+D5XAArpW90rLR9dcq/CUmnS1osafGKFSt6dBFdzWtmZlZJw3qvSboRKNeF67yImJfynAesB35a2q1M/qB8cIwq+asda9PEiBnADIDW1tYeVVnckcDMrFjDgk5EHFltu6SpwHHAEdHVNtUOjMllGw08lZbLpT8L7ChpaKrN5POXjtUuaSgwAljZ8yuqLlI885v6ZmaV9VXvtUnAF4HjI+Ll3Kb5wJTU82wcMB64A7gTGJ96qg0j62wwPwWrW4AT0/5TgXm5Y01NyycCN4cfvJiZ9am+ejn0B8BwYFGqGdwWEZ+MiCWS5gIPkjW7nRkRHQCSzgIWAi3ArIhYko71RWCOpG8A9wAzU/pM4EpJbWQ1nCmNvCCHMzOzYn0SdFI35krbLgAuKJO+AFhQJv0xst5t3dPXACf1rqSbz61rZmaVbQm91waEzo4E7r9mZlaRg06ddHUk6OOCmJltwRx06swxx8ysMgedOnFHAjOzYg46ddL5RqqrOmZmFTno1EnXfDqOOmZmlTjo1IlrOmZmxRx0zMysaSq+HCrpLdV2jIi761+c/ssdCczMilUbkeA76XtroBX4I1mP4AOB24F3NbZo/UznKNNuXzMzq6Ri81pEHBYRhwFPAm+JiNaIOAh4M9nkapbT+XJoH5fDzGxLVssznf0i4v7SSkQ8AExsXJH6J8+nY2ZWrJYBPx+S9CPgJ2SNSB8GHmpoqfoxxxwzs8pqCTofA84APpPWbwUua1iJ+in3IzAzK1YYdCJijaRLgBvJ7q2PRMS6hpesn+l8OdTta2ZmFRU+05F0KLCUbOK1S4FHJf1tb04q6euS7pN0r6QbJO2R0iXpYkltaftbcvtMlbQ0fabm0g+SdH/a52Klu76knSUtSvkXSdqpN2Uu4pdDzcyK1dKR4DvAURHxnoj4W+BoYHovz/utiDgwIiYCvwT+JaUfQzZF9XjgdFIznqSdgfOBt5FN2HZ+LohclvKW9puU0s8BboqI8cBNab1huubTMTOzSmoJOltFxCOllYh4FNiqNyeNiBdyq9vSVVGYDFwRmduAHSWNIgt0iyJiZUSsAhYBk9K2HSLiD5G1b10BnJA71uy0PDuX3liu6piZVVRLR4LFkmYCV6b1U4G7entiSRcApwGrgcNS8p7Asly29pRWLb29TDrAayNiOUBELJe0W5WynE5WW2Kvvfbq0fWEuxKYmRWqpaZzBrAE+DRZD7YHgU8W7STpRkkPlPlMBoiI8yJiDPBT4KzSbmUOFT1I3ywRMSO9/No6cuTIzd19o7O6nmNmVlktvdfWSvoBWZNWzb3XIuLIGsvwM+BXZM9s2oExuW2jgadS+qHd0n+T0keXyQ/wtKRRqZYzCnimxvL0iDsSmJkV66vea+Nzq8cDD6fl+cBpqRfbIcDq1ES2EDhK0k6pA8FRwMK07UVJh6Rea6cB83LHKvVym5pLb4iujgSOOmZmldTyTKfUe+0RAEn7AFcBB/XivBdK2hfYQDa2W6m5bgFwLNnYbi+TvZhKRKyU9HXgzpTvaxGxMi2fAVwObANcnz4AFwJzJU0D/g84qRflrZlrOmZmldUSdDbpvSapt73XPlghPYAzK2ybBcwqk74YOKBM+nPAEb0p5+ZwRwIzs2J91nttoPF7OmZmxWoJOmeQ1T4+TXZPvZXs2Y7luCOBmVmxmnqvAd9NH6ugc+w113XMzCoqDDqS3gl8BXhdPn9EvL5xxerHHHPMzCqqpXltJvA5suc4HY0tTv8V7kdgZlaolqCzOiKuL85m4IqOmVk1FYNOblqBWyR9C/g5sLa0PSLubnDZ+pWu6aoddszMKqlW0/lOt/XW3HIAh9e/OP1X6T0dhxwzs8oqBp2IOKzSNqvMFR0zs8qqNa99OCJ+Iunsctsjwl2oc9yRwMysWLXmtW3T9/bNKEh/55dDzcyKVWte+2H6/mrzitN/+eVQM7Ni1ZrXLq62Y0R8uv7F6f9c0zEzq6xa85oH9dwMfqRjZlasWvPa7Py6pG0j4qXGF6l/ckcCM7Nitcwc+nZJDwIPpfU3SfIo05tIz3TcvmZmVlFh0AH+HTgaeA4gIv4I9Gq66hJJ/ygpJO2a1iXpYkltku7LjYqApKmSlqbP1Fz6QZLuT/tcnKatRtLOkhal/IvSNNcN4/l0zMyK1RJ0iIhl3ZJ6PfCnpDHAe8mmki45BhifPqcDl6W8OwPnA28DDgbOzwWRy1Le0n6TUvo5wE0RMR64Ka03nCs6ZmaV1RJ0lkl6BxCShkn6R1JTWy9NB77Axs/gJwNXROY2YEdJo8hqWosiYmVErAIWAZPSth0i4g9pqusrgBNyxyo9l5qdS28IP9IxMytWS9D5JNnMoXsC7cDEtN5jko4H/pya6vL2BPK1qvaUVi29vUw6wGsjYjlA+t6tSnlOl7RY0uIVK1b04IryzWuu6piZVVLL1AYbIuLUfIKkcaRnPJVIuhHYvcym84B/Bo4qt1uZtOhB+maJiBnADIDW1tYeVVo6B/x0zDEzq6iWms4vJO1QWpH0RuAXRTtFxJERcUD3D/AYMA74o6QngNHA3ZJ2J6upjMkdZjTwVEH66DLpAE+n5jfS9zM1XGuPuSOBmVmxWoLOv5IFnu0kHQRcA3y4pyeMiPsjYreIGBsRY8kCx1si4i/AfOC01IvtELIJ5JYDC4GjJO2UOhAcBSxM216UdEjqtXYaMC+daj5Q6uU2NZfeUK7pmJlVVti8FhG/krQVcAPZ4J8nRMTSBpVnAXAs0Aa8DHwslWGlpK8Dd6Z8X4uIlWn5DOByYBvg+vQBuBCYK2kaWQ+5kxpUZrIyNvLoZmYDQ7Wx177Pxs9HdiBrGvuUpLqNvZZqO6XloEInhYiYBcwqk74YOKBM+nPAEfUoYy2ia5zpZp3SzKzfqVbTWdxt3WOxVdE1XXXflsPMbEtW89hrVhvHHDOzyqo1r82NiJMl3U+ZbsgRcWBDS9ZPeew1M7PKqjWvfSZ9H9eMgvR37khgZlasWvNa6W3+J5tXnP6r8+XQPi6HmdmWrFrz2ouUf7tfZB3NdiizbdByRwIzs2LVajrbN7Mg/V1nh2kHHTOzimqa2sBq5wE/zcwqc9Cpk3BPAjOzQg46ddIZclzRMTOryEGnTjzKtJlZscIBPyv0YltNNkzO5yPisUYUrP8pzafjsGNmVkktk7h9l2yOmp+R/SE/hWxytkfIBuA8tFGF648ccszMKquleW1SRPwwIl6MiBfSLJvHRsTVwE4NLl+/4X4EZmbFagk6GySdLGlI+pyc2+ZbbeL3dMzMitUSdE4FPkI23fMzafnDkrYBzurJSSV9RdKfJd2bPsfmtp0rqU3SI5KOzqVPSmltks7JpY+TdLukpZKuljQspQ9P621p+9ielLVWXR0JHHXMzCopDDoR8VhEvC8idk2f90VEW0S8EhG/68W5p0fExPRZACBpAtkzo/2BScClkloktQCXAMcAE4BTUl6Ai9KxxgOrgGkpfRqwKiL2BqanfA1Tek/HNR0zs8oKg46k0ZKuk/SMpKclXStpdIPKMxmYExFrI+JxsmmrD06fthQAXwXmAJOVdRU7HLgm7T8bOCF3rNKcQNcAR6gJXcscc8zMKqulee3HwHxgD2BP4BcprbfOknSfpFmSSh0S9gSW5fK0p7RK6bsAz0fE+m7pGx0rbV+d8jeEH26ZmRWrJeiMjIgfR8T69LkcGFm0k6QbJT1Q5jMZuAx4AzARWA58p7RbmUNFD9KrHatcWU+XtFjS4hUrVlS5qsqi2lnNzAyo7T2dZyV9GLgqrZ8CPFe0U0QcWUsBJP0n8Mu02g6MyW0eTfaOEBXSnwV2lDQ01Wby+UvHapc0FBgBrKxQ1hnADIDW1tYeVVq65tNx1DEzq6SWms7HgZOBv5DVSk5MaT0maVRu9f3AA2l5PjAl9TwbB4wH7gDuBMannmrDyDobzI/s6f0tqUwAU4F5uWNNTcsnAjdHI0fl9Hw6ZmaFCms6EfF/wPF1Pu83JU0ku1U/AXwinWuJpLnAg8B64MyI6ACQdBawEGgBZkXEknSsLwJzJH0DuAeYmdJnAldKaiOr4Uyp8zWU5ZhjZlZZtZlDv0+V5+MR8emenjQiPlJl2wXABWXSFwALyqQ/Rta7rXv6GuCknpZxc7kjgZlZsWo1ncVNK8UA0DVdtes6ZmaVVJuuenalbbapzo4EjjlmZhV5Pp068Xw6ZmbFHHTqzDUdM7PKHHTqxB0JzMyK1TJz6NZkg2fuD2xdSo+IXr2rM9CEhyQwMytUS03nSrKZQo8Gfkv21v+LjSxUf+T5dMzMitUSdPaOiC8DL6UebX8H/E1ji9UPuSOBmVmhWoLOuvT9vKQDyMYwG9uwEvVzfk/HzKyyWgb8nJGmHvgS2Xhm2wFfbmip+qFwVwIzs0K1BJ2bImIVcCvwesimiG5oqfohv6djZlaslua1a8ukXVMmbVALjzJtZlao2oCf+5F1kx4h6QO5TTuQ6zptma4O0446ZmaVVGte2xc4DtgReF8u/UXgHxpZqP7MNR0zs8qqDfg5D5gn6e0R8YcmlqlfauT8cGZmA0UtHQnukXQmHpGgKoccM7NifTYigaRPSXpE0hJJ38ylnyupLW07Opc+KaW1STonlz5O0u2Slkq6Ok1nTZry+uqU/3ZJY3tb5mrckcDMrFifjEgg6TBgMnBgROwPfDulTyCbVnp/YBJwqaQWSS3AJcAxwATglJQX4CJgekSMB1aRjRNH+l4VEXsD01O+hnNHAjOzyvpqRIIzgAsjYi1ARDyT0icDcyJibUQ8DrSRTUV9MNAWEY9FxKvAHGCystf/D6erC/ds4ITcsUoT0V0DHKGGDhfgBjYzsyK1BJ3SiARfJhuR4EHgm9V3KbQP8O7U7PVbSW9N6XsCy3L52lNapfRdgOcjYn239I2OlbavTvk3Iel0SYslLV6xYkWPLsjNa2ZmxQo7EkTEj9Lib0kjEtRC0o1kz4K6Oy+ddyfgEOCtwFxJr6f8C/1B+eAYVfJTsG3jxIgZwAyA1tbWHlVZPMq0mVmxai+Hnl1tx4j4bsH2I6sc+wzg55H1M75D0gZgV7Kayphc1tHAU2m5XPqzwI6ShqbaTD5/6VjtkoaSNQuurFbm3ugaBsdRx8yskmrNa9unTyvZM5hSM9cnyR7m98Z/kz2LQdI+wDCyADIfmJJ6no0DxgN3AHcC41NPtWFknQ3mp6B1C3BiOu5UYF5anp/WSdtvjia8TOOajplZZdVeDv0qgKQbgLdExItp/SvAf/XyvLOAWZIeAF4FpqaAsETSXLLnRuuBMyOiI533LGAh0ALMiogl6VhfBOZI+gZwDzAzpc8ErpTURlbDmdLLMlflUabNzIrV8nLoXmSBoeRVetl7LfVA+3CFbRcAF5RJXwAsKJP+GFnvtu7pa4CTelPOzeFRps3MitUSdK4ke+5yHdnz8vfT1RXZEnckMDMrVkvvtQskXQ+8OyV9LCLuaWyx+p+ux0WOOmZmldRS0yEi7gbubnBZBgTXdMzMKqvl5VAzM7O6cNCpE3ckMDMr5qBTJ6Uu0w0d3s3MrJ+rNiLBi1QZxTIidmhIifop13TMzIpVezl0ewBJXwP+QtZ1WsCpZCMVWBmu6JiZVVZL89rREXFpRLwYES9ExGXABxtdsP7Gs1WbmRWrJeh0SDo1TaY2RNKpQEejC9bfdL2l46qOmVkltQSdvwdOBp5On5NSmuWUXg5185qZWWW1jEjwBNksnFaFW9fMzIoV1nQk7SPppjQiNJIOlPSlxhetf3JNx8ysslqa1/4TOBdYBxAR99HgaQL6JVd1zMwK1RJ0XhMRd3RLW9+IwvRnfjnUzKxYLUHnWUlvIP0tL+lEYHlDS9UP+eVQM7NitQSdM4EfAvtJ+jPwWbIpq3tM0tWS7k2fJyTdm9t2rqQ2SY9IOjqXPimltUk6J5c+TtLtkpam4w5L6cPTelvaPrY3ZS7i+XTMzIrVEnSejIgjgZHAfhHxroh4sjcnjYgPRcTEiJgIXAv8HEDSBLLnRfsDk4BL0/tBLcAlwDHABOCUlBfgImB6RIwHVgHTUvo0YFVE7A1MT/kazu/pmJlVVkvQeVzSDOAQ4K/1PLmyByAnA1elpMnAnIhYGxGPA21kU1EfDLRFxGNpqus5wOS0/+HANWn/2cAJuWOVZji9BjhCDXzg4hEJzMyK1RJ09gVuJGtme1zSDyS9q07nfzfwdEQsTet7Asty29tTWqX0XYDnI2J9t/SNjpW2r075NyHpdEmLJS1esWJFjy6kqyNBj3Y3MxsUCoNORLwSEXMj4gPAm4EdgN8W7SfpRkkPlPnkXzQ9ha5aDpR/Dh89SK92rE0TI2ZERGtEtI4cObJclkLuSGBmVqym6aolvQf4ENkzlTvJmsSqSs+Bqh1zKPAB4KBccjswJrc+GngqLZdLfxbYUdLQVJvJ5y8dqz2dawSwsqjcPdUZzRx1zMwqqmVEgsfJeqz9D3BARJwcEdfW4dxHAg9HRHsubT4wJfU8GweMB+4gC3TjU0+1YWSdDeZHNuDZLcCJaf+pwLzcsaam5ROBmyMa/+TFHQnMzCqrWtNJvcZ+HBFfa8C5p7Bx0xoRsUTSXOBBshdQz4yIjlSWs4CFQAswKyKWpN2+CMyR9A3gHmBmSp8JXCmpjayG09hRFNyTwMysUNWgExEdkg4D6h50IuKjFdIvAC4ok74AWFAm/TGy3m3d09eQjYjdFH5Px8ysWC3PdP5X0g+Aq4GXSokRcXfDStUPuSOBmVmxWoLOO9J3vrYTZO/HWNI1n47DjplZJbXMp3NYMwoyUDjkmJlVVkvvtddKminp+rQ+QdK0ov0GG3cjMDMrVsuIBJeT9RrbI60/StaF2nI6n+m4qmNmVlEtQWfXiJgLbIDOIWU6GlqqfqhrGARHHTOzSmoJOi9J2oWu+XQOIRvHzHLC3dfMzArV0nvtbLK3+98g6fdkUxycWH2XwcvNa2ZmldXSe+3uNPbavmR/xz8SEesaXjIzMxtwaum9dhKwTRp25gTgaklvaXjJ+hm3rpmZFavlmc6XI+LFNIfO0WQTo13W2GL1P13z6TjsmJlVUkvQKfVU+zvgsoiYBwxrXJE37L8zAAAMvUlEQVT6J9d0zMyK1RJ0/izph2Rz6CyQNLzG/QYlV3TMzCqrJXicTPZy6KSIeB7YGfinhpaqH/KIBGZmxWqZrvpl4AngGEmfAkZFxA2NLlh/09W85qqOmVkltfRe+xeyzgO7ALsCP5b0pd6cVNJESbdJulfSYkkHp3RJulhSm6T78r3kJE2VtDR9pubSD5J0f9rnYqUn+ZJ2lrQo5V8kaafelLlIV0eCRp7FzKx/q6V57RTgrRFxfkScDxwCnNrL834T+GpETAT+Ja0DHEM2RfV44HRSLzlJOwPnA28jm7Dt/FwQuSzlLe03KaWfA9wUEeOBm9J6w3jiUDOzYrUEnSeArXPrw4E/9fK8AeyQlkcAT6XlycAVkbkN2FHSKLKu2osiYmVErAIWAZPSth0i4g+RjUNzBdm7RKVjzU7Ls3PpDeWajplZZRVHJJD0fbLgsBZYImlRWn8v8LtenvezwEJJ3yYLfKWJ4vYEluXytae0auntZdIBXhsRywEiYrmk3XpZZjMz66Vqw+AsTt93Adfl0n9Ty4El3QjsXmbTecARwOci4lpJJwMzgSMp/5pL9CB9s0g6nayJjr322mtzd89OWpo51B0JzMwqqhh0ImI2gKStgb3JbuZ/iog1tRw4Io6stE3SFcBn0up/AT9Ky+3AmFzW0WRNb+3Aod3Sf5PSR5fJD/C0pFGpljMKeKZKWWcAMwBaW1t79HTG8+mYmRWr+ExH0lBJ3yS7sc8GfgIsk/RNSVv18rxPAe9Jy4cDS9PyfOC01IvtEGB1aiJbCBwlaafUgeAoYGHa9qKkQ1KvtdOAebljlXq5Tc2lN5RjjplZZdWa174FbA+Mi4gXASTtAHw7fT5TZd8i/wB8T9JQYA2paQtYABwLtAEvAx8DiIiVkr4O3JnyfS0iVqblM8hmN90GuD59AC4E5qaptf8POKkX5S3kzmtmZsWqBZ3jgH0iujoDR8QLks4AHqYXQScifgccVCY9gDMr7DMLmFUmfTFwQJn058ieHTVFV/Oa6zpmZpVU6zId+YCTS+zAf9hvovPl0D4uh5nZlqxa0HlQ0mndEyV9mKymYznuSGBmVqxa89qZwM8lfZys23QAbyV7dvL+JpStX3LzmplZZdW6TP8ZeJukw4H9yVqOro+Im5pVuP7E7Y1mZsWq1XQAiIibgZubUJb+zYOvmZkV8mRsdRL4eY6ZWREHnTqJcM81M7MiDjp15E4EZmbVOejUSbgrgZlZIQedOnHzmplZMQedOnFHAjOzYg46dZLVdBx1zMyqcdCpJ8ccM7OqHHTqxB0JzMyKOejUizsSmJkVctCpE3ckMDMr1idBR9KbJP1B0v2SfpFmJC1tO1dSm6RHJB2dS5+U0toknZNLHyfpdklLJV0taVhKH57W29L2sY28pohwRwIzswJ9VdP5EXBORPwNcB3wTwCSJgBTyEa1ngRcKqlFUgtwCXAMMAE4JeUFuAiYHhHjgVXAtJQ+DVgVEXsD01O+hnJNx8ysur4KOvsCt6blRcAH0/JkYE5ErI2Ix4E24OD0aYuIxyLiVWAOMFnZuDOHA9ek/WcDJ+SONTstXwMcoQaOU+NBps3MivVV0HkAOD4tnwSMSct7Asty+dpTWqX0XYDnI2J9t/SNjpW2r075NyHpdEmLJS1esWJFjy4ocEcCM7MiDQs6km6U9ECZz2Tg48CZku4CtgdeLe1W5lCV7ufV0qsda9PEiBkR0RoRrSNHjqx2WRVFeMBPM7MihZO49VREHFmQ5SgASfsAf5fS2umq9QCMBp5Ky+XSnwV2lDQ01Wby+UvHapc0FBgBrOzZ1RQLwjUdM7MCfdV7bbf0PQT4EvAfadN8YErqeTYOGA/cAdwJjE891YaRdTaYHxEB3AKcmPafCszLHWtqWj4RuDnlb+CFNfToZmb9Xl890zlF0qPAw2Q1kx8DRMQSYC7wIPBr4MyI6Ei1mLOAhcBDwNyUF+CLwNmS2sie2cxM6TOBXVL62UBnN+tGcEcCM7NiDWteqyYivgd8r8K2C4ALyqQvABaUSX+MrHdb9/Q1ZJ0UmsYVHTOz6jwiQZ1EhDsSmJkVcNCpEw+DY2ZWzEGnjhxzzMyqc9CpE3ckMDMr5qBTJ4Gf6ZiZFemT3msD0QF7jGDdeld3zMyqcdCpkykH78WUg/fq62KYmW3R3LxmZmZN46BjZmZN46BjZmZN46BjZmZN46BjZmZN46BjZmZN46BjZmZN46BjZmZNo0ZPptnfSFoBPNnD3Xclm0J7MPE1Dw6+5sGhN9f8uogYWZTJQaeOJC2OiNa+Lkcz+ZoHB1/z4NCMa3bzmpmZNY2DjpmZNY2DTn3N6OsC9AFf8+Dgax4cGn7NfqZjZmZN45qOmZk1jYOOmZk1jYNOHUiaJOkRSW2Szunr8tSLpFmSnpH0QC5tZ0mLJC1N3zuldEm6OP0M7pP0lr4rec9JGiPpFkkPSVoi6TMpfcBet6StJd0h6Y/pmr+a0sdJuj1d89WShqX04Wm9LW0f25fl7w1JLZLukfTLtD4YrvkJSfdLulfS4pTWtN9vB51ektQCXAIcA0wATpE0oW9LVTeXA5O6pZ0D3BQR44Gb0jpk1z8+fU4HLmtSGettPfD5iHgjcAhwZvr3HMjXvRY4PCLeBEwEJkk6BLgImJ6ueRUwLeWfBqyKiL2B6Slff/UZ4KHc+mC4ZoDDImJi7p2c5v1+R4Q/vfgAbwcW5tbPBc7t63LV8frGAg/k1h8BRqXlUcAjafmHwCnl8vXnDzAPeO9guW7gNcDdwNvI3kwfmtI7f8+BhcDb0/LQlE99XfYeXOvodIM9HPgloIF+zan8TwC7dktr2u+3azq9tyewLLfentIGqtdGxHKA9L1bSh9wP4fUhPJm4HYG+HWnZqZ7gWeARcCfgOcjYn3Kkr+uzmtO21cDuzS3xHXx78AXgA1pfRcG/jUDBHCDpLsknZ7Smvb7PbQ3OxuQ/XXU3WDshz6gfg6StgOuBT4bES9I5S4vy1omrd9dd0R0ABMl7QhcB7yxXLb03e+vWdJxwDMRcZekQ0vJZbIOmGvOeWdEPCVpN2CRpIer5K37dbum03vtwJjc+mjgqT4qSzM8LWkUQPp+JqUPmJ+DpK3IAs5PI+LnKXnAXzdARDwP/IbsedaOkkp/mOavq/Oa0/YRwMrmlrTX3gkcL+kJYA5ZE9u/M7CvGYCIeCp9P0P2B8bBNPH320Gn9+4ExqdeL8OAKcD8Pi5TI80HpqblqWTPPErpp6XeLocAq0vV9f5EWZVmJvBQRHw3t2nAXrekkamGg6RtgCPJHq7fApyYsnW/5tLP4kTg5kgN/v1FRJwbEaMjYizZ/9mbI+JUBvA1A0jaVtL2pWXgKOABmvn73dcPtQbCBzgWeJSsHfy8vi5PHa/rKmA5sI7sL55pZO3YNwFL0/fOKa/IevH9CbgfaO3r8vfwmt9F1nxwH3Bv+hw7kK8bOBC4J13zA8C/pPTXA3cAbcB/AcNT+tZpvS1tf31fX0Mvr/9Q4JeD4ZrT9f0xfZaU7lfN/P32MDhmZtY0bl4zM7OmcdAxM7OmcdAxM7OmcdAxM7OmcdAxM7OmcdAxqzNJf03fYyX9fZ2P/c/d1v+3nsc3azQHHbPGGQtsVtBJo5ZXs1HQiYh3bGaZzPqUg45Z41wIvDvNW/K5NKjmtyTdmeYm+QSApEOVzeHzM7IX8JD032lAxiWlQRklXQhsk47305RWqlUpHfuBNFfKh3LH/o2kayQ9LOmnadQFJF0o6cFUlm83/adjg5IH/DRrnHOAf4yI4wBS8FgdEW+VNBz4vaQbUt6DgQMi4vG0/vGIWJmGpblT0rURcY6ksyJiYplzfYBsLpw3AbumfW5N294M7E82ZtbvgXdKehB4P7BfRERpGByzRnNNx6x5jiIbx+pesukSdiGbHAvgjlzAAfi0pD8Ct5ENuDie6t4FXBURHRHxNPBb4K25Y7dHxAayYX3GAi8Aa4AfSfoA8HKvr86sBg46Zs0j4FORzdg4MSLGRUSppvNSZ6ZsqP0jySYNexPZuGhb13DsStbmljvIJilbT1a7uhY4Afj1Zl2JWQ856Jg1zovA9rn1hcAZaeoEJO2TRvrtbgTZ1MgvS9qPbJqBknWl/bu5FfhQem40EvhbsoEpy0rzBY2IiAXAZ8ma5swazs90zBrnPmB9aia7HPgeWdPW3elh/gqyWkZ3vwY+Kek+sumBb8ttmwHcJ+nuyIbiL7mObHrlP5KNkv2FiPhLClrlbA/Mk7Q1WS3pcz27RLPN41Gmzcysady8ZmZmTeOgY2ZmTeOgY2ZmTeOgY2ZmTeOgY2ZmTeOgY2ZmTeOgY2ZmTfP/Acqpj/9KFCt4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot([i for i in range(1, len(em.log_lik))], em.log_lik[1:])\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('Observed data log likelihood')\n",
    "plt.title('EM log likelihood')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}