Appendix C The E-M Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add lib input sys.path\n",
    "import os\n",
    "import sys\n",
    "import time\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sklearn\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize\n",
    "import math\n",
    "from sklearn.preprocessing import normalize\n",
    "from functools import partial\n",
    "import h5py\n",
    "from scipy.spatial import distance\n",
    "\n",
    "nb_dir = os.path.split(os.getcwd())[0]\n",
    "if nb_dir not in sys.path:\n",
    "    sys.path.append(nb_dir)\n",
    "\n",
    "from matplotlib.colors import ListedColormap\n",
    "import libs.linear_models as lm\n",
    "import libs.data_util as data\n",
    "import libs.nn as nn\n",
    "import libs.plot as myplot\n",
    "from sklearn import decomposition\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise C.1\n",
    "\n",
    "* (a) Apply Lagrange multiplier $\\lambda$, we have $L= \\sum^K_{k=1}N_k\\ln w_k + \\lambda (\\sum_k w_k -1)$, take derivative w.r.t. $w_k$, we have\n",
    "\n",
    "$\\frac{\\partial{L}}{\\partial{w_k}} = N_k\\frac{1}{w_k} + \\lambda = 0$, we thus have\n",
    "\n",
    "$w_k = -\\frac{N_k}{\\lambda}$\n",
    "\n",
    "As we have $\\sum_k w_k = -\\frac{1}{\\lambda}\\sum_k N_k = -\\frac{N}{\\lambda} = 1$, so $\\lambda = -\\frac{1}{N}$ thus $w_k = \\frac{N_k}{N}$\n",
    "\n",
    "* (b) $L_k(X_k; \\theta_k) = \\sum_{x_n \\in X_k} \\ln P_k(x_n; \\theta_k) = \\sum_{x_n \\in X_k} -\\ln (2\\pi)^{\\frac{d}{2}}|\\Sigma_k|^{\\frac{1}{2}} -\\frac{1}{2}(x_n-\\mu_k)^T\\Sigma^{-1}_k (x_n-\\mu_k)$\n",
    "\n",
    "Take derivative w.r.t. $\\mu_k$, and set to 0, we have \n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial{L_k}}{\\partial{\\mu_k}} &= \\sum_{x_n \\in X_k} \\frac{1}{2} 2\\Sigma^{-1}_k (x_n-\\mu_k)\\\\\n",
    "&= \\Sigma^{-1}_k \\sum_{x_n \\in X_k}  (x_n-\\mu_k)\\\\ \n",
    "&= 0\\\\\n",
    "\\end{align*}\n",
    "\n",
    "So we have $\\mu_k^* = \\frac{1}{N_k}\\sum_{x_n \\in X_k}x_n$\n",
    "\n",
    "Let $S_k = \\Sigma^{-1}_k$, then we have $L_k(X_k; \\theta_k) =  \\sum_{x_n \\in X_k} -\\ln (2\\pi)^{\\frac{d}{2}}|\\Sigma_k|^{\\frac{1}{2}} -\\frac{1}{2}(x_n-\\mu_k)^T\\Sigma^{-1}_k (x_n-\\mu_k) = \\sum_{x_n \\in X_k} -\\ln (2\\pi)^{\\frac{d}{2}}+ \\frac{1}{2} \\ln |S| -\\frac{1}{2}(x_n-\\mu_k)^TS_k (x_n-\\mu_k)$\n",
    "\n",
    "\n",
    "Take derivative w.r.t. $S_k$, and set to 0, we have\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial{L_k}}{\\partial{S_k}} &= \\sum_{x_n \\in X_k} \\frac{1}{2}S^{-1} - \\frac{1}{2}(x_n-\\mu_k)(x_n-\\mu_k)^T\\\\\n",
    "&= 0\\\\ \n",
    "\\end{align*}\n",
    "\n",
    "We thus have $(S^{-1})^* =\\frac{1}{N_k} \\sum_{x_n \\in X_k}(x_n-\\mu_k)(x_n-\\mu_k)^T $\n",
    "\n",
    "and $\\Sigma^*_k = S^{-1}_k = \\frac{1}{N_k} \\sum_{x_n \\in X_k}(x_n-\\mu_k)(x_n-\\mu_k)^T $\n",
    "\n",
    "\n",
    "#### Exercise C.2\n",
    "\n",
    "Solve the equations we have $\\hat{\\mu}_1 = \\frac{1}{2}, \\hat{\\mu}_2 = 0$\n",
    "\n",
    "#### Problem C.1\n",
    "\n",
    "The probability to get a green ball is $P(g) = \\frac{1}{2}(1-\\mu_1)$, to get a blue ball $P(b) = \\frac{1}{2}(1-\\mu_2)$ and to get a red ball $P(r) = \\frac{1}{2}\\mu_1 + \\frac{1}{2}\\mu_2$, thus the log-likelihood of the incomplete data is\n",
    "\n",
    "\\begin{align*}\n",
    "\\ln \\prod P &= \\sum \\ln P \\\\\n",
    "&= N_g\\ln P(g) + N_b \\ln P(b) + N_r \\ln P(r) \\\\\n",
    "&= N_g \\ln \\frac{1}{2}(1-\\mu_1) + N_b \\ln \\frac{1}{2}(1-\\mu_2) + N_r \\ln( \\frac{1}{2}\\mu_1 + \\frac{1}{2}\\mu_2) \\\\\n",
    "&= N_g \\ln (1-\\mu_1) + N_b \\ln(1-\\mu_2) + N_r \\ln( \\mu_1 + \\mu_2) - N\\ln 2\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Take derivative w.r.t. $\\mu_1,\\mu_2$ and set them to zero, we have\n",
    "\n",
    "$\\frac{\\partial{L}}{\\partial{\\mu_1}} = -\\frac{N_g}{1-\\mu_1} + \\frac{N_r}{\\mu_1 + \\mu_2} = 0$\n",
    "\n",
    "$\\frac{\\partial{L}}{\\partial{\\mu_2}} = -\\frac{N_b}{1-\\mu_2} + \\frac{N_r}{\\mu_1 + \\mu_2}= 0$\n",
    "\n",
    "Solve for this we have $\\mu_2 = 1-\\frac{N_b}{N_g} + \\frac{N_b}{N_g}\\mu_1$, take this into above equations we have\n",
    "\n",
    "$\\mu_1 = 1 - 2 \\frac{N_g}{N}$ and $\\mu_2 = 1-2\\frac{N_b}{N}$\n",
    "\n",
    "The likelihood under these values is $L_1 = N_g\\ln \\frac{N_g}{N} + N_b \\ln \\frac{N_b}{N} + N_r \\ln \\frac{N_r}{N}$\n",
    "\n",
    "Note, the $\\mu_1,\\mu_2$ are less than 0 when $N_g > \\frac{N}{2}$ or $N_b > \\frac{N}{2}$, which doesn't make sense. So we need solve constrained maximum likelihood with $\\mu_1 \\ge 0, \\mu_2 \\ge 0$. Note, $\\mu_1,\\mu_2$ are implicitly less than 1 from the log-likelihood formula.\n",
    "\n",
    "Set $\\mu_1=0$ and solve for optimal $\\mu_2 = \\frac{N_r}{N_r+N_b}$, which corresponds to a likelihood of $L_2 = N_b \\ln \\frac{N_b}{N_r+N_b} + N_r \\ln \\frac{N_r}{N_r+N_b} - N\\ln 2$\n",
    "\n",
    "similarly, if we let $\\mu_2 = 0$ the optimal $\\mu_1 = \\frac{N_r}{N_r+N_g}$ and \n",
    "\n",
    "$L_3 = N_g \\ln \\frac{N_g}{N_r+N_g} + N_r \\ln \\frac{N_r}{N_r+N_g} - N\\ln 2$\n",
    "\n",
    "We need select the maximum likelihood among $L_1,L_2,L_3$, we compute\n",
    "\n",
    "$L_1 - L_2 = N_g\\ln \\frac{N_g}{N} + N_b \\ln \\frac{N_r+N_b}{N} + N_r \\ln \\frac{N_r + N_b}{N} + N\\ln 2 = N\\left[\\frac{N_g}{N}\\ln \\frac{N_g}{N} + (1-\\frac{N_g}{N})\\ln(1-\\frac{N_g}{N}) + \\ln 2\\right]$\n",
    "\n",
    "Let $x=\\frac{N_g}{N}$ we have $L_1 - L_2 = N[x\\ln x + (1-x)\\ln (1-x) + \\ln 2]$\n",
    "\n",
    "Take derivative w.r.t. $x$ and set it to 0, we see that the derivative is positive when $x>\\frac{1}{2}$ and negative when $x<\\frac{1}{2}$. Note the difference $L_1 - L_2 = 0$ when $x=\\frac{1}{2}$, so we see that $L_1 \\ge L_2$ for all values of $\\frac{N_g}{N}$. \n",
    "Similarly we see that $L_1 \\ge L_3$ for all values of $\\frac{N_g}{N}$. \n",
    "\n",
    "So the maximum likelihhod estimates for $\\mu_1, \\mu_2$ are\n",
    "\n",
    "* $\\mu_1 = 1 - 2 \\frac{N_g}{N}$ and $\\mu_2 = 1-2\\frac{N_b}{N}$ when $N_g \\le \\frac{N}{2}$ and $N_b \\le \\frac{N}{2}$.\n",
    "\n",
    "* $\\mu_1=0, \\mu_2 = \\frac{N_r}{N_r+N_b}$ when $N_g >\\frac{N}{2}$\n",
    "* $\\mu_1=\\frac{N_r}{N_r+N_g}, \\mu_2 = 0$ when $N_b > \\frac{N}{2}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem C.2\n",
    "\n",
    "* (a) Given a red ball, we have \n",
    "\n",
    "\\begin{align*}\n",
    "p_1 &= P[Bag 1|red; \\hat{\\mu_1},\\hat{\\mu_2}] \\\\\n",
    "&= \\frac{P[red \\Cap Bag 1]}{P[red]}\\\\\n",
    "&= \\frac{P[red | Bag 1]P[Bag 1]}{P[red]}\\\\\n",
    "&= \\frac{\\hat{\\mu_1}\\frac{1}{2}}{\\frac{1}{2}\\hat{\\mu_1} + \\frac{1}{2}\\hat{\\mu_2}}\\\\\n",
    "&= \\frac{\\hat{\\mu_1}}{\\hat{\\mu_1} + \\hat{\\mu_2}}\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Similarly we have $p_2 = P[Bag 2|red; \\hat{\\mu_1},\\hat{\\mu_2}]  = \\frac{\\hat{\\mu_2}}{\\hat{\\mu_1} + \\hat{\\mu_2}}$ given a red ball.\n",
    "\n",
    "* (b) The probability to get one green ball is $p_g = \\frac{1}{2}(1-\\mu_1)$, to get $N_g$ green balls, we have $p_g = (\\frac{1}{2}(1-\\mu_1))^{N_g}$, similarly, we have $p_b = (\\frac{1}{2}(1-\\mu_2))^{N_b}$ and the probability to get $N_{r_1}$ red balls from bag 1 is $p_{r_1} = (\\frac{1}{2}\\mu_1)^N_{r_1}$, the probability to get $N_{r_2}$ red balls from bag 2 is $p_{r_2} = (\\frac{1}{2}\\mu_1)^N_{r_2}$\n",
    "\n",
    "The log-likelihood of the complete data is thus\n",
    "\n",
    "$L = N_g\\ln (1-\\mu_1) + N_b\\ln(1-\\mu_2) + N_{r_1}\\ln \\mu_1 + N_{r_2}\\ln \\mu_2 - N\\ln 2$\n",
    "\n",
    "* (c) Take the expectation w.r.t. $p_1, p_2$ of the log-likelihood, note the green balls and blue balls are known to come from bag 1 and bag 2 respectively. So only $N_{r_1}, N_{r_2}$ are random variables w.r.t. $p_1,p_2$.\n",
    "\n",
    "$E[N_{r_1}] = p_1N_r, E[N_{r_2}] = p_2N_r$\n",
    "\n",
    "So the expectation is \n",
    "\n",
    "\\begin{align*}\n",
    "Q_t(\\mu_1, \\mu_2) &=  E[N_g\\ln (1-\\mu_1) + N_b\\ln(1-\\mu_2) + N_{r_1}\\ln \\mu_1 + N_{r_2}\\ln \\mu_2 - N\\ln 2]\\\\\n",
    "&= N_g\\ln (1-\\mu_1) + N_b\\ln(1-\\mu_2) + E[N_{r_1}] \\ln \\mu_1 + E[N_{r_2}]\\ln \\mu_2 - N\\ln 2\\\\\n",
    "&= N_g\\ln (1-\\mu_1) + N_b\\ln(1-\\mu_2) + p_1N_r\\ln \\mu_1 + p_2N_r\\ln \\mu_2 - N\\ln 2 \\\\\n",
    "\\end{align*}\n",
    "\n",
    "* (d) Take derivatives w.r.t. $\\mu_1,\\mu_2$, and set them to 0, we have\n",
    "\n",
    "$\\frac{\\partial{Q_t}}{\\partial{\\mu_1}} = -\\frac{N_g}{1-\\mu_1} + \\frac{p_1N_r}{\\mu_1}$\n",
    "\n",
    "$\\frac{\\partial{Q_t}}{\\partial{\\mu_2}} = -\\frac{N_g}{1-\\mu_2} + \\frac{p_2N_r}{\\mu_2}$\n",
    "\n",
    "Thus $\\mu_1 = \\frac{p_1N_r}{N_g + p_1N_r}, \\mu_2 = \\frac{p_bN_r}{N_b + p_2N_r}$\n",
    "\n",
    "* (e) \n",
    "\n",
    "\\begin{align*}\n",
    "\\mu_1 &= \\frac{p_1N_r}{N_g + p_1N_r} \\\\\n",
    "&= \\frac{N_r}{\\frac{N_g}{p_1} + N_r} \\\\\n",
    "&= \\frac{N_r}{\\frac{N_g}{\\frac{\\mu_1}{\\mu_1 + \\mu_2}} + N_r} \\\\\n",
    "\\end{align*}\n",
    "\n",
    "Re-arrange the above equation we have $\\mu_1(N_r+N_g) + \\mu_2N_g =N_r$\n",
    "\n",
    "Similarly  we have $\\mu_2(N_r+N_b) + \\mu_1N_b =N_r$\n",
    "\n",
    "Solve the two equations we have\n",
    "\n",
    "$\\mu_1 = \\frac{N-2N_g}{N}, \\mu_2 = \\frac{N-2N_b}{N}$\n",
    "\n",
    "\n",
    "#### Problem C.3\n",
    "\n",
    "* (a) We have \n",
    "\n",
    "$P[x=1] = \\frac{1}{2}\\times 1 + \\frac{1}{2}\\times \\pi = \\frac{1+\\pi}{2}$\n",
    "\n",
    "$P[x=0] = \\frac{1}{2}\\times 0 + \\frac{1}{2}\\times (1-\\pi) = \\frac{1-\\pi}{2}$\n",
    "\n",
    "So we have $P[x=x_n] = (\\frac{1+\\pi}{2})^{x_n}(\\frac{1-\\pi}{2})^{1-x_n}$\n",
    "\n",
    "The likelihood is thus\n",
    "\n",
    "\\begin{align*}\n",
    "P[X|\\pi, N] &= \\prod^N_{n=1}P[x=x_n]\\\\\n",
    "&= \\prod^N_{n=1}(\\frac{1+\\pi}{2})^{x_n}(\\frac{1-\\pi}{2})^{1-x_n}\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Take log, we have the log-likelihood $L = \\sum^N_{n=1} x_n \\ln (\\frac{1+\\pi}{2}) + (1-x_n) \\ln (\\frac{1-\\pi}{2})$, take derivative w.r.t. $\\pi$ and let it equal to 0, we have\n",
    "\n",
    "$\\sum^N_{n=1} (\\frac{x_n}{1+\\pi} - \\frac{1-x_n}{1-\\pi}) = \\frac{N_r}{1+\\pi} - \\frac{N-N_r}{1-\\pi}  = 0$, we have $\\pi = \\frac{2N_r-N}{N}$.\n",
    "\n",
    "If $N_r < \\frac{N}{2}$, $\\pi = 0$\n",
    "\n",
    "For $N_r = 600, N = 1000$, we have $\\pi = 0.2$\n",
    "\n",
    "* (b) The appropriate hidden variable $J$ is the bag membership of each $x_n$.\n",
    "\n",
    "* (c) Given $\\pi, N$, we have $P[x_n, j_n] = P[j_n]P[x_n|j_n] = \\frac{1}{2}$ if $j_n \\in bag 1$, and $P[x_n, j_n] = \\frac{1}{2}\\pi $ if $j_n \\in bag 2$\n",
    "\n",
    "Let $N_{r_1} + N_{r_2} = N_r$ where $N_{r_1}$ is the red balls from bag 1 and $N_{r_2}$ is the number of red balls from bag 2. \n",
    "\n",
    "The log-likelihood of the complete data is \n",
    "\n",
    "$N_{r_1}\\ln \\frac{1}{2} + N_{r_2} \\ln (\\frac{\\pi}{2}) + N_b \\ln (\\frac{1-\\pi}{2})  = N_{r_2}\\ln \\pi + N_b \\ln (1-\\pi) - N\\ln 2$\n",
    "\n",
    "  * (iii) We compute the probability of red ball comes from bag 1, as $p_1 = P[bag 1| red] = \\frac{P[bag1, red]}{P[red]} = \\frac{P[red|bag 1]P[bag 1]}{P[red]} = \\frac{\\frac{1}{2}\\times 1}{\\frac{1}{2}\\times 1 + \\frac{1}{2}\\times \\pi_t} = \\frac{1}{1+\\pi_t}$ \n",
    "  \n",
    "  Similarly, we have $p_2 = \\frac{\\pi_t}{1+\\pi_t}$ \n",
    "  \n",
    "  Take expectation w.r.t. $p_1,p_2$, we have\n",
    "  \n",
    "\\begin{align*}\n",
    "Q_t(\\pi)&= E\\left[N_{r_2}\\ln \\pi + N_b \\ln (1-\\pi) - N\\ln 2\\right]\\\\\n",
    "&= E[N_{r_2}]\\ln \\pi + N_b \\ln (1-\\pi) - N\\ln 2\\\\\n",
    "&= p_2N_r\\ln \\pi + (N-N_r) \\ln (1-\\pi) - N\\ln 2\\\\\n",
    "&= \\frac{\\pi_t}{1+\\pi_t}N_r\\ln \\pi  + (N-N_r) \\ln (1-\\pi) - N\\ln 2\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Take derivative of $Q_t(\\pi)$ w.r.t $\\pi$ and set to 0, we solve for the $\\pi_{t+1} = \\frac{\\pi_tN_r}{\\pi_tN_r + (1+\\pi_t)(N-N_r)}$\n",
    "\n",
    "When $N_r \\ge \\frac{N}{2}$, set $\\pi_{t+1} = \\pi_t$, we have $\\pi_t \\to 2\\frac{N_r}{N} - 1$\n",
    "\n",
    "When $N_r < \\frac{N}{2}$, we should have $\\pi = 0$\n",
    "\n",
    "  * (v) After we obtain the convergence of E-M, we can compute the probability of $P[j_n = 1]$ using the $\\pi$ as the $p_1,p_2$ above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xuc3HV97/HXe3MfuUmyIM1lN0A0iRcublFqe8R6aUBPQOVY6HqOcDhuUanl2J4KXaWadtXWU4uPx8F61ha1zQoiLZrS2NRS8FRbKQuGS0giEbPJGpQNl1QJ5Po5f/x+O/nt7MzsJfvbmcm8n4/HPHZ+3/nO7/f5zc7M+3eb308RgZmZGUBLrQswM7P64VAwM7Mih4KZmRU5FMzMrMihYGZmRQ4FMzMrcigcAyR9TNLaKo9vl/Sm6awpL5I+Kena9P4Fkgar9P2MpKvHMc63SPp6Zvjnkk4fx/M+KOlT4629zPM3Sbpgss/PS/Y1Hkffz0v6aN41TSdJp0raLGlOrWupBYdCjaRf1M+nX0A/lfRFScfVuq6jJemlkr4mabekPZIekvQhSTPK9H2DpIclPSvpKUl3SFpYZdytwH8D/u84y/k00C1p9hj9PgEUv9wj4riIeHwc4+8F3i3plEodJIWk59L/84/ToJqRTuflEXHPOKYzJcZaeEj7TOg1joirI+IPxzHtN0raImmvpLsltY2v6vGT1J6Oe286rYoLQpIWSvqGpKclDWYXHiLip8DdQNdU19gIHAq19Z8j4jjgXOAXgY+UdlCiIf5Pks4A7gV2Aq+MiBOB/wJ0AMeXecqjwK9FxEnALwCPAX9eZRJXAOsj4vnx1BMRTwBbgNVVav5F4MSI+N54xlky/heAb5J8iVZzVvp/fiPwG8B7JzqtaXQFE3iNx0PSAuBvgY8CJwP9wFenavwZtwDfB+YD3cDtaciVsxb4EXAq8FbgE5LekHm8D/jNHGqsew3xZXOsi4gfk3y5vAJA0j2SeiR9F9gLnC7pFyStS5dstkkq/WKZK+mrkn4m6QFJZ5WblqQWSddJ+mG6dH6bpJPTx9rTJdsrJe2U9IykqyX9YrrE/6yk/1NlVj4O/GtEfCj9QiYitkbEb0TEs2Xm+6cRsSvTdAg4s8r4LwS+XWaefj9dM9kuqbPk4XtIPvTjHmf6Gpwp6bWSfpJdy5H0dkkPTWD8RRGxBfgXjvyfy27WG2u6ks6T1C/pP9K1zM+MNW1Jq4DfB349XWt5sELXEa+H0k10lV5jSV+S9EdjTP4dwKaI+FoapB8DzpK0fKy6x0vSS0kWrv4gIp6PiL8BHgbeWabvccAFQE9EHIiIB4Hbgf+e6XYvyeduytdo6p1DoQ5IWgxcRLKUM+y/kqy+Hg8MkCwFDZIsUV9KsmTzxkz/i4GvkSyJfQX4uqRZZSb3QeAS4PXpuJ4Bbirp8xpgGfDrwI0kS11vAl4OvEvS6yvMyptIPlzjJmmJpGeB54HfBf6kSvdXAltL2l4CLAAWAu8BeiW9LPP4ZqBsQFYZJwDp2sNzwK9mmn+D5PUd7/iLJK0EfoWR/+fJTPezwGcj4gTgDOC2saYdEf9Aspnsq+nmsUo1T+Y1HsvLgWIIRcRzwA/T9lEk3ZkugJS73VllGo9HxM8ybQ9WmIZK/g7ff0WmxoPANsb5vz2WOBRq6+vpF+J3SJbOPpF57EsRsSl9c74E+GXgwxHxQkRsBP6CJDiG3R8Rt0fEAeAzwFzgtWWm+ZtAd0QMRsQ+kqW2SyXNzPT5w3Q6/0jy5XRLRDyZrtH8C3BOhfmZDzwxkRcgInakm48WkGw+21Kl+0nAz8q0fzQi9kXEt4G/B96Veexn6fMmOs5htwCXA0g6niS8bykZ/4lVng/wgKRngL8j+b99cYz+Y033AHCmpAUR8fPJbPqqYjKv8ViOA/aUtO2h/CZFIuJtEXFShdvbjnYaaXB8F/iopLmSziVZoyiUdB3rvXNMcijU1iXpG70tIt5fsh13Z+b+LwBPlywFDZAsuY3qHxGHObJWUaoNuGN4yYtkSfcQybbVYT/N3H++zHClHeJPAadVeKyqiHga+DLwjZKAynqG0R/yZ9Ilz2EDjJzv44FnobiZ6efp7fNVxpn1FeAdSo5EeQfwQEQMlIy/9Muo1LkR8eKIOCMiPpL+f8ZSbbpXAS8Ftki6T1KlL8rJmMxrPJafAyeUtJ1A9TCeqIlOoxNYSvK5+XOSfQilR7IV3zvNxKFQv7Knr90FnJwuMQ5bAvw4M7x4+E66Y3pR+rxSO4ELS5a+5qZrAUfrnyizDXcCZgKnMPrDPewhki/DrBdLelFmeAkj53sF6aaLiPhEuunkuIgYPtqk3DiLIuJRki/BCxm96WjE+KdStelGxGMRcTnJa/XHJDtUX1R2RCWjHUefybzGY9lEZjNMOq4z0vZRJH0zE96lt29WmcbpJZ+RsypNIyIG0jWS1oh4Dcla7r9naphJsn9ryv+39c6h0AAiYifwr8An09XdV5EsLfZlur1a0jvSN/O1wD6g3GaFzwM9wzvQJLVKuniKSv0D4JckfVrSS9LxnylpraRRq+FpvS9Ld363kmz2+n661lDOepJ9IaU+Lmm2pF8B3kayb2XY60l24ldSaZxZXyHZF/OfSsY9nvEfjbLTlfRuSa3pGsfwkuyhcYzvp0C7qh/NNpnXeCx3AK+Q9E5Jc4EbgIfSHe+jRMSFmfAuvV1Y4Tk/ADYCf5B+Rt4OvAr4m3L9Ja2QdHw6T+8G3kLy/ht2HrC9ZK2wKTgUGsflQDvJEtodJEdZfCvz+DdIdgw/Q7Kv4R3p/oVSnwXWAf8o6WckwfGaqSgwIn4InJ/WuUnSHpIPZT/pany6tPcr6VMWAv+QPvYwcBh4e5VJ/BVwkaR5mbafkMzzLpKQvHr4y0bSacBK4OulI8rU/ACwR1K11+AWkqNV/jkidg83pl9wF5Fs9spD2ekCq0he35+T/D8vS4/qKX19Sw1/kT8l6YEKfSb0Go9HRAyRrEH2pON5DXDZeJ8/AZeRHP78DMnvTi5Np42kTknZtYZfAx5P+14NrBrum+okWYBqOgpfZMcaiKRPAE9GxI3j6PunwA8j4nNj9HsL8P6IuGSCtfwWsDgifm8iz6t32ddYyS+u10bEohqXNW2U/Bjx28A5w2HbTBwKZlZRM4ZCs/PmIzM7aiVHdo1nx7DVKa8pmJlZkdcUzMysqNKPhOrWggULor29vdZlmJk1lPvvv393RFQ6QWBRw4VCe3s7/f39tS7DzKyhSBrXby68+cjMzIocCmZmVuRQMDOzoobbp2BmVg8OHDjA4OAgL7xQXz96njt3LosWLWLWrHKXUxmbQ8HMbBIGBwc5/vjjaW9vR9LYT5gGEcFTTz3F4OAgS5cundQ4ct18JGmVpK1KLh95XZnH2yTdpeRSj/dIyuWn9H19fbS3t9PS0kJ7ezt9fX1jP8nMrIoXXniB+fPn100gAEhi/vz5R7X2klsoKLm+7E0k54NfCVyeXo4w638DfxURrwLWAJ+c6jr6+vro6upiYGCAiGBgYICuri4Hg5kdtXoKhGFHW1OeawrnAdsi4vGI2A/cSnId4ayVwF3p/bvLPH7Uuru72bt374i2vXv30t3dPdWTMjNreHmGwkJGXlJykJGXj4TkqkbDV+p6O3C8pPmlI5LUJalfUv/Q0FDpw1Xt2LFjQu1mZs0sz1Aotw5Teva93wVeL+n7JFd7+jFwcNSTInojoiMiOlpbx/yV9ghLliyZULuZWR4aZd9mnqEwSOa6wZS5ZnBE7IqId0TEOUB32jbWRdAnpKenh0KhMKKtUCjQ09MzlZMxM6uokfZt5nlI6n3AMklLSdYALiO5AHmRpAXA0+m1Zq8Hbp7qIjo7OwG49tpr2b17N6eddhqf/vSni+1mZkfr2muvZePGjRUf/973vse+fftGtO3du5errrqKL3zhC2Wfc/bZZ3PjjWNeYHDK5RYKEXFQ0jXABmAGcHNEbJK0BuiPiHUk15/9pKQA/h/wgTxq6ezs5LjjjuOSSy7hzjvv5Nxzz81jMmZmZZUGwljttZTrj9ciYj2wvqTthsz924Hb86xh2Lx5yXXIn3/++emYnJk1kbGW6Nvb2xkYGH2S0ra2Nu65556cqpqcpjn3kUPBzGqlkfZtNl0olP5mwcwsb52dnfT29tLW1oYk2tra6O3trct9m01z7qPhlPaagpnVQmdnZ12GQKmmW1NwKJiZVeZQMDOzIoeCmdkkRZSepKH2jrampgsF72g2s6kwd+5cnnrqqboKhuHrKcydO3fS42iaHc2zZs1ixowZXlMwsymxaNEiBgcHmehJOvM2fOW1yWqaUJDEvHnzHApmNiVmzZo16aub1bOm2XwEOBTMzMbgUDAzsyKHgpmZFTVdKPjoIzOzypoqFAqFgtcUzMyqaKpQ8OYjM7PqHApmZlbkUDAzs6KmCwXvaDYzq6ypQsE7ms3Mqss1FCStkrRV0jZJ15V5fImkuyV9X9JDki7Ksx5vPjIzqy63UJA0A7gJuBBYCVwuaWVJt48At0XEOcBlwOfyqgccCmZmY8lzTeE8YFtEPB4R+4FbgYtL+gRwQnr/RGBXjvUwb948Dh06xIEDB/KcjJlZw8ozFBYCOzPDg2lb1seAd0saBNYDv1VuRJK6JPVL6j+a09T6QjtmZtXlGQoq01Z6NYrLgS9FxCLgIuCvJY2qKSJ6I6IjIjpaW1snXZAvtGNmVl2eoTAILM4ML2L05qGrgNsAIuLfgLnAgrwKKhQKgNcUzMwqyTMU7gOWSVoqaTbJjuR1JX12AG8EkLSCJBRyu4yRNx+ZmVWXWyhExEHgGmADsJnkKKNNktZIWp12+x3gvZIeBG4BrogcL3jqUDAzqy7Xy3FGxHqSHcjZthsy9x8FXpdnDVkOBTOz6prqF83e0WxmVl1ThYJ3NJuZVddUoeDNR2Zm1TkUzMysyKFgZmZFDgUzMytqylDw0UdmZuU1VSjMmjWLmTNnek3BzKyCpgoF8DUVzMyqcSiYmVmRQ8HMzIqaMhS8o9nMrLymC4VCoeA1BTOzCpouFLz5yMysMoeCmZkVORTMzKzIoWBmZkVNFwqFQsFHH5mZVdB0oeA1BTOzyhwKZmZWlGsoSFolaaukbZKuK/P4n0namN5+IOnZPOuBI6EQEXlPysys4czMa8SSZgA3AW8GBoH7JK2LiEeH+0TE/8z0/y3gnLzqGTZv3jwOHz7MgQMHmD17dt6TMzNrKHmuKZwHbIuIxyNiP3ArcHGV/pcDt+RYD5DsaAZfU8HMrJw8Q2EhsDMzPJi2jSKpDVgK/HOFx7sk9UvqHxoaOqqifPU1M7PK8gwFlWmrtCH/MuD2iDhU7sGI6I2IjojoaG1tPaqiHApmZpXlGQqDwOLM8CJgV4W+lzENm47AoWBmVk2eoXAfsEzSUkmzSb7415V2kvQy4MXAv+VYS5FDwcysstxCISIOAtcAG4DNwG0RsUnSGkmrM10vB26NaTpG1KFgZlZZboekAkTEemB9SdsNJcMfy7OGUj76yMyssqb8RTN4TcHMrByHgpmZFTkUzMysyKFgZmZFTRcK3tFsZlZZ04WC1xTMzCprulCYMWMGs2bNciiYmZXRdKEAvtCOmVklDgUzMytyKJiZWVFThkKhUPDRR2ZmZTRlKHhNwcysPIeCmZkVORTMzKzIoWBmZkVNGQre0WxmVl5ThoLXFMzMynMomJlZkUPBzMyKcg0FSaskbZW0TdJ1Ffq8S9KjkjZJ+kqe9QwbDoWImI7JmZk1jJl5jVjSDOAm4M3AIHCfpHUR8WimzzLgeuB1EfGMpFPyqierUCgQEezfv585c+ZMxyTNzBpCnmsK5wHbIuLxiNgP3ApcXNLnvcBNEfEMQEQ8mWM9RcPXVPARSGZmI+UZCguBnZnhwbQt66XASyV9V9L3JK0qNyJJXZL6JfUPDQ0ddWG+0I6ZWXl5hoLKtJVuxJ8JLAMuAC4H/kLSSaOeFNEbER0R0dHa2nrUhTkUzMzKyzMUBoHFmeFFwK4yfb4REQci4kfAVpKQyJVDwcysvDxD4T5gmaSlkmYDlwHrSvp8HXgDgKQFJJuTHs+xJsChYGZWSW6hEBEHgWuADcBm4LaI2CRpjaTVabcNwFOSHgXuBv5XRDyVV03DCoUC4B3NZmalcjskFSAi1gPrS9puyNwP4EPpbdp4TcHMrLym/UUzOBTMzEo5FMzMrMihYGZmRU0ZCsM7mh0KZmYjNWUo+DQXZmblNXUoeE3BzGykpgyFlpYWZs+e7VAwMyvRlKEAvtCOmVk5VX+8JumPI+LDkr7G6JPZERHvyq2ynBUKBYeCmVmJsX7R/B1JJwP3AvuBPfmXND3mzZvnHc1mZiXGCoUCyZlMfwbMAd4ZEXflXtU08OYjM7PRxtqn0A38UkScCrwd+Gj+JU0Ph4KZ2WhjhcLhiNgIEBF3AyfmX9L0cCiYmY021uaj2ZJWcOQqanOywxHxaJ7F5WnevHns2XPM7CIxM5sS49mnsL6kbXg4gNOnvKJpUigU+MlPflLrMszM6krVUIiI9mmqY9r56CMzs9H84zUzMytyKJiZWZFDwczMinINBUmrJG2VtE3SdWUev0LSkKSN6e1/5FlP1vBpLpLLRJuZGYx99NGkSZoB3AS8meRX0fdJWlfmMNavRsQ1edVRyfDps1944YXifTOzZpfnmsJ5wLaIeDwi9gO3AhfnOL0J8TUVzMxGyzMUFgI7M8ODaVupd0p6SNLtkhaXG5GkLkn9kvqHhoampDiHgpnZaHmGgsq0lW7A/zugPSJeBfwT8OVyI4qI3ojoiIiO1tbWKSnOoWBmNlqeoTAIZJf8FwG7sh0i4qmI2JcOfgF4dY71jOBQMDMbLc9QuA9YJmmppNnAZcC6bAdJp2UGVwObc6xnhEKhADgUzMyycjv6KCIOSroG2ADMAG6OiE2S1gD9EbEO+KCk1cBB4GngirzqKTW8puBTXZiZHZFbKABExHpKTqgXETdk7l8PXJ9nDZV485GZ2WhN/YtmcCiYmWU5FBwKZmZFTRsK3tFsZjZa04aCdzSbmY3W9KHgNQUzsyOaNhTmzp0LOBTMzLKaNhRaWlqYM2eOQ8HMLKNpQwGOXFPBzMwSTR0KvvqamdlITR8KPvrIzOyIpg8FrymYmR3hUHAomJkVORQcCmZmRU0dCj76yMxspKYOBe9oNjMbqWlDoa+vjw0bNvDII4/Q3t5OX19frUsyM6u5pgyFvr4+urq6imsJAwMDdHV1ORjMrOk1ZSh0d3eP2my0d+9euru7a1SRmVl9aMpQ2LFjx4TazcyaRa6hIGmVpK2Stkm6rkq/SyWFpI486xm2ZMmSCbWbmTWL3EJB0gzgJuBCYCVwuaSVZfodD3wQuDevWkr19PQUr7w2rFAo0NPTM10lmJnVpTzXFM4DtkXE4xGxH7gVuLhMvz8E/gR4IcdaRujs7KS3t5fFixcDcOKJJ9Lb20tnZ+d0lWBmVpfyDIWFwM7M8GDaViTpHGBxRNxZbUSSuiT1S+ofGhqakuI6OzvZsWMHp59+OqtWrXIgmJmRbyioTFsUH5RagD8DfmesEUVEb0R0RERHa2vrFJYIK1asYPPmzVM6TjOzRpVnKAwCizPDi4BdmeHjgVcA90jaDrwWWDddO5uHLV++nK1bt3Lo0KHpnKyZWV3KMxTuA5ZJWippNnAZsG74wYjYExELIqI9ItqB7wGrI6I/x5pGWbFiBfv27WNgYGA6J2tmVpdyC4WIOAhcA2wANgO3RcQmSWskrc5ruhO1fPlyAG9CMjMDZuY58ohYD6wvabuhQt8L8qylkhUrVgCwZcsW3vrWt9aiBDOzutGUv2jOOvnkkznllFO8pmBmhkMBSDYhbdmypdZlmJnVnEOBI4elRsTYnc3MjmEOBZI1haeffprdu3fXuhQzs5pyKHBkZ7P3K5hZs3Mo4MNSzcyGORSAxYsXUygUvLPZzJqeQwFoaWlh+fLlXlMws6bnUEj5sFQzM4dC0YoVKxgYGOC5556rdSlmZjXjUEgN72z+wQ9+UONKzMxqx6GQ8mGpZmYOhaIzzzyTlpYWh4KZNTWHQmrOnDmcccYZ3tlsZk3NoZDhS3OaWbNzKGQsX76cxx57jIMHD9a6FDOzmnAoZKxYsYL9+/fzox/9qNalmJnVhEMhY/iwVO9XMLNm5VDIePjhhwFYvXo17e3t9PX11bgiM7Pp5VBI9fX1ce211xaHBwYG6OrqcjCYWVPJNRQkrZK0VdI2SdeVefxqSQ9L2ijpO5JW5llPNd3d3ezdu3dE2969e+nu7q5RRWZm0y+3UJA0A7gJuBBYCVxe5kv/KxHxyog4G/gT4DN51TOWHTt2TKjdzOxYlOeawnnAtoh4PCL2A7cCF2c7RMR/ZAZfBNTsIslLliyZULuZ2bEoz1BYCOzMDA+mbSNI+oCkH5KsKXyw3IgkdUnql9Q/NDSUS7E9PT0UCoURbYVCgZ6enlymZ2ZWj/IMBZVpG7UmEBE3RcQZwIeBj5QbUUT0RkRHRHS0trZOcZmJzs5Oent7aWtrK7Z9+MMfprOzM5fpmZnVozxDYRBYnBleBOyq0v9W4JIc6xlTZ2cn27dvZ8+ePcybN48nnniiluWYmU27PEPhPmCZpKWSZgOXAeuyHSQtywy+FXgsx3rG7YQTTuDSSy/llltu4fnnn691OWZm0ya3UIiIg8A1wAZgM3BbRGyStEbS6rTbNZI2SdoIfAh4T171TNSVV17Jnj17uOOOO2pdipnZtFFEzQ74mZSOjo7o7+/PfTqHDx/mjDPO4Mwzz+Rb3/pW7tMzM8uTpPsjomOsfv5FcwUtLS1cccUV3HXXXQwMDNS6HDOzaeFQqOKKK64A4Mtf/nJtCzEzmyYOhSra2tpYsWIFa9asoaWlxSfJM7Nj3sxaF1DP+vr62LZtG4cOHQKOnCQP8O8XzOyY5DWFKrq7u9m/f/+INp8kz8yOZQ6FKnySPDNrNg6FKnySPDNrNg6FKsqdJK+lpYU1a9bUqCIzs3w5FKrIniRPEgsWLODw4cNs3ry51qWZmeXCoTCG4ZPkHT58mKGhId773vfyqU99ilNPPdWHqZrZMcehMEHnn38+knjyySeJCF/L2cyOKQ6FCfr4xz9O6fmifJiqmR0rHAoTVOlw1IGBAdrb271JycwamkNhgqodjjowMOBNSmbW0BwKE1TuMNVyvEnJzBqRQ2GCSg9TzV7TuZQ3KZlZo3EoTEL2MNXt27ePGQzZTUrvf//7HRRmVrccClNgIpuUPv/5zzsozKxuORSmwEQ2KZU7nLU0KK688koWLFgwIiT6+vpGBUe5NjOzoxIRud2AVcBWYBtwXZnHPwQ8CjwE3AW0jTXOV7/61dEI2traApiS26xZs2L27NljthUKhXjf+94XbW1tISna2tpi7dq1sXbt2rpoq+faPF+er0ZrmyigP8bzvT2eTpO5ATOAHwKnA7OBB4GVJX3eABTS++8DvjrWeBslFNauXRuFQmHEl7akKQuKiQTKrFmzJhUyU9022bCr9zbPV2O1HQvzVSgUJhwM1EEonA9syAxfD1xfpf85wHfHGm+jhEJElF1CqYeg8M033xr/1tbWNqHvI8YZCnnuU1gI7MwMD6ZtlVwFfLPcA5K6JPVL6h8aGprCEvNVepTS5z73uVH7Hq6++upx7aQ2M8vK62JfeYaCyrRF2Y7Su4EO4NPlHo+I3ojoiIiO1tbWKSxx+o0VFPPnz2f27NkjnjNr1qxxtUnlXnIzOxbldbGvPENhEFicGV4E7CrtJOlNQDewOiL25VhP3coGxe7du7n55ptHrE188YtfHFdbubWO8QbKdLSVqqfaPF+er0ZqKxQK9PT0VJ2vSRvPNqbJ3ICZwOPAUo7saH55SZ9zSHZGLxvveBtpn0It1NPREcfqUR+er8ZuO1bma6IY5z4FJX3zIeki4EaSI5FujogeSWvS4tZJ+ifglcAT6VN2RMTqauPs6OiI/v7+3Go2MzsWSbo/IjrG6jczzyIiYj2wvqTthsz9N+U5fTMzmxj/otnMzIocCmZmVuRQMDOzIoeCmZkV5Xr0UR4kDQEDk3z6AmD3FJZTC56H+nEszIfnoT5Mxzy0RcSYv/5tuFA4GpL6x3NIVj3zPNSPY2E+PA/1oZ7mwZuPzMysyKFgZmZFzRYKvbUuYAp4HurHsTAfnof6UDfz0FT7FMzMrLpmW1MwM7MqHApmZlbUNKEgaZWkrZK2Sbqu1vWMh6SbJT0p6ZFM28mSviXpsfTvi2tZ41gkLZZ0t6TNkjZJ+u20vWHmQ9JcSf8u6cF0Hj6eti+VdG86D1+VVP3E/XVA0gxJ35d0ZzrcUPMgabukhyVtlNSftjXMewlA0kmSbpe0Jf1cnF9P89AUoSBpBnATcCGwErhc0sraVjUuXwJWlbRdB9wVEcuAu9LhenYQ+J2IWAG8FvhA+to30nzsA341Is4CzgZWSXot8MfAn6Xz8AzJJWXr3W8DmzPDjTgPb4iIszPH9TfSewngs8A/RMRy4CyS/0f9zMN4LrrQ6DfgfGBDZvh64Ppa1zXO2tuBRzLDW4HT0vunAVtrXeME5+cbwJsbdT6AAvAA8BqSX6DOTNtHvMfq8UZy9cO7gF8F7iS5ZG6jzcN2YEFJW8O8l4ATgB+RHuRTj/PQFGsKwEJgZ2Z4MG1rRKdGxBMA6d9TalzPuElqJ7na3r002Hykm102Ak8C3yK5YuCzEXEw7dII76kbgd8DDqfD82m8eQjgHyXdL6krbWuk99LpwBDwxXQz3l9IehF1NA/NEgrlrmjvY3GnkaTjgL8Bro2I/6h1PRMVEYci4mySpe3zgBXluk1vVeMn6W3AkxFxf7a5TNe6nYfU6yLiXJJNwR+Q9J9qXdAEzQTOBf48Is4BnqPONnc1SygMAoszw4uAXTWq5Wj9VNJpAOnfJ2tcz5gkzSIJhL6I+Nu0ueHmAyAingXuIdk/cpKk4asX1vt76nXAaknbgVtJNiHdSGPNAxGxK/37JHAHSUA30ntpEBiMiHvT4dtJQqJu5qFZQuE+YFl6pMVs4DJgXY1rmqx1wHsFVaFGAAACIElEQVTS++8h2UZftyQJ+Etgc0R8JvNQw8yHpFZJJ6X35wFvItk5eDdwadqtruchIq6PiEUR0U7y/v/niOikgeZB0oskHT98H3gL8AgN9F6KiJ8AOyW9LG16I/Ao9TQPtd7xMo07eC4CfkCyLbi71vWMs+ZbgCeAAyRLGFeRbAe+C3gs/XtyrescYx5+mWSTxEPAxvR2USPNB/Aq4PvpPDwC3JC2nw78O7AN+Bowp9a1jnN+LgDubLR5SGt9ML1tGv4cN9J7Ka33bKA/fT99HXhxPc2DT3NhZmZFzbL5yMzMxsGhYGZmRQ4FMzMrciiYmVmRQ8HMzIpmjt3FrPlIuheYA8wGXkpyKCokJ43bFckx/mbHHB+SalZFer6m/ohYUONSzKaFNx+ZTYCkCzLn8W+XtFvSJ9OTm22R9GpJX5D0UHqdgpdknvt76XUZHpD0d9nHzOqFQ8Hs6MwHvhPJyc3+kuTXqDdFxKuA+4FrACS9GzgTeG0kJ3RbD/xpbUo2q8z7FMyOzs8j4u/T+w+QnOxsYzp8P8m1IwBWAx3AA8npoJgJ7JnOQs3Gw6FgdnT2Ze4fAl4oGR7+jAn4o4i4eboKM5sMbz4ymx7rgPcPX3tX0hxJZ9W4JrNRvKZgNg0i4q8lLQC+nW4+agE+R3LGT7O64UNSzcysyJuPzMysyKFgZmZFDgUzMytyKJiZWZFDwczMihwKZmZW5FAwM7Oi/w/IoR4lll4D6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#### Problem C.3\n",
    "\n",
    "def calc_pis(pi0, Nr, N, maxit = 1000, tol = 1.0e-12):\n",
    "    pi = math.inf\n",
    "    itx = 1\n",
    "    pis = [pi0]\n",
    "    while itx < maxit:\n",
    "        r1 = pi0 * Nr\n",
    "        r2 = r1 + (1+pi0)*(N-Nr)\n",
    "        pi = r1/r2\n",
    "        pis.append(pi)\n",
    "        itx += 1\n",
    "        if np.abs(pi-pi0) <= tol:\n",
    "            break\n",
    "        pi0 = pi\n",
    "    return pis, itx\n",
    "\n",
    "Nr = 600\n",
    "N = 1000\n",
    "pis_1, itx1 = calc_pis(0.9, Nr, N)\n",
    "myplot.plt_plot([np.arange(itx1)], [pis_1], 'plot', ['k', 'g', 'r', 'b', 'y'], ['o', 'x', '.', ',', 's'], \n",
    "                [' ', '', '', '', ''], \n",
    "        title = 'Problem C.3 (b)-(iv) Pi vs. t (pi_0 = 0.9)', yscale = None, ylb = 0, yub = 1.2,\n",
    "        xlb = None, xub = None, xlabel = 'Time', ylabel = 'Pi',\n",
    "        legends = ['']*6, legendx = None, legendy = None, marker_sizes=[10]*5) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+YXVV97/H3hySQpiCEJFaawEyUWA0aAxwCtlcF+ZWgJv5ASQwVemNT9aE+Xm6tVGrRaLRAW9CntBAqSkv4IbTgqCAgAqWV0ExCiCSQmzHNjxEqIQlBiAQj3/vHXhN2Ts7MOTOZNT8yn9fznCdnr73W2mudkznfvdbe5yxFBGZmZr3tgP5ugJmZ7Z8cYMzMLAsHGDMzy8IBxszMsnCAMTOzLBxgzMwsCwcY24OkL0q6oYv96yWd1pdtykXS1yR9Jj0/WVJ7F3n/TtInGqjzDEl3lLZfkPT6Bsp9WtJfN9r2GuVXSTq5p+VzKb/GDeS9WtIXcrepL0n6HUlPSDqov9vSHxxg9gPpQ/9X6cPsF5K+Jeng/m7XvpL0Rkm3SnpW0nZJKyVdKGlYjbynSPqppOckbZF0u6TxXdQ9DvgYcE2DzbkcuFjSgXXyfRXYHSgi4uCIWNdA/YuAcyW9trMMkkLSi+l9/nkKesPScY6JiAcaOE6vqHcikvJ06zWOiE9ExJcbOPapkp6UtEPS/ZKaGmt14yQ1p7p3pGN1elIl6W8krZX0y5T3Yx37IuIXwP3A/N5u42DgALP/eF9EHAwcB5wA/GV1BhUGxXsu6Q3AI8Am4K0RcSjwYaACHFKjyGrgzIg4DPhdYC3wj10c4nzgzoj4VSPtiYingSeBmV20+QTg0IhY0kidVfW/BNxF8YHclbel9/lU4KPAH3f3WH3ofLrxGjdC0ljg34AvAIcDrcAtvVV/yU3Ao8AY4GLgthQwa3kReB9wKHAe8HVJv1/avxj4kwxtHPAGxYeNNS4ifk7xQfUWAEkPSFoo6T+BHcDrJf2upBZJWyW1Sar+kBop6ZZ0RrZc0ttqHUvSAZIukvSzNGr4jqTD077mdMb9R5I2Sdom6ROSTkgjkeck/X0XXfkS8JOIuDB9uBMRayLioxHxXI1+/yIiniol/QY4uov6ZwAP1ujT59OIab2kuVW7HwDe050602twtKSTJP1PefQl6QOSVnaj/t0i4kngIV59n2tOXdY7rqRpklolPZ9Gv39X79iSpgOfB85Jo6nHOsm6x+uhNA3Z2Wss6duSvlLn8B8EVkXErSkofxF4m6Q31Wt3oyS9keJE7ZKI+FVE/CvwU+BDtfJHxCUR8WREvBIRj1C8L28vZXmE4u+u10daA50DzH5G0pHAWRRnXx3+kGKIfgiwgeLsrJ3iTP9s4KuSTi3lnwXcSnGGeCNwh6QRNQ73aeD9wLtSXduAq6rynAhMAs4BrqQ4GzwNOAb4iKR3ddKV04Db6vf4VZKOkvQc8Cvgz4DLusj+VmBNVdrrgLHAeIoz0UWSfq+0/wmgZrDtok4A0qjmReDdpeSPUry+jda/m6TJwDvY833uyXG/Dnw9Il4DvAH4Tr1jR8QPKaYCb0lTgJ21uSevcT3HALsDWkS8CPwspe9F0vfTyUytx/e7OMa6iPhlKe2xzo5RdbzfophBWFVq4y6gjQbf2/2JA8z+44704fofFGeNXy3t+3ZErEr/0V8H/C/gcxHxUkSsAP6JIgh1WBYRt0XEr4G/A0YCJ9U45p8AF0dEe0TspDibPFvS8FKeL6fj3EPxQXdTRDyTRloPAcd20p8xwNPdeQEiYmOaIhtLMUX4ZBfZDwN+WSP9CxGxMyIeBH4AfKS075epXHfr7HATMAdA0iEUJwI3VdV/aBflAZZL2gZ8j+J9+1ad/PWO+2vgaEljI+KFnkzvdaEnr3E9BwPbq9K2U3valIh4b0Qc1snjvb1xjCpXUwSju6vS6/3f2S85wOw/3p/+aJoi4lNV896bSs9/F9hadXa2geKMcq/8EfEKr452qjUBt3ecEVKcgf8G+J1Snl+Unv+qxnZnNyNsAY7oZF+XImIrcD3w3apgV7aNvT8wtqUz4g4b2LPfhwDPwe6ptBfS4+ou6iy7EfigijuKPggsj4gNVfVXf7BVOy4iRkfEGyLiL9P7U09Xx50HvBF4UtJSSZ196PZET17jel4AXlOV9hq6Duzd1aNjSLqcYsryI7H3rwjv/r8zlDjADA3l/+xPAYenM9kORwE/L20f2fEk3RQwIZWrtgmYUXVWODKNTvbVj+hkzrtBw4HXsvcHRYeVFB+sZaMl/XZp+yj27PebSdMzEfHVND10cER03L5cq87dImI1xQfqDPaeHtuj/t7U1XEjYm1EzKF4rS6luJj92zUrqqq2gTw9eY3rWUVpqinV9QZKU1Jlku4qnQhUP+7q4hivr/obeVtnx0jH+RLF63tGRDxftW84xfXAXn9vBzoHmCEmIjYBPwG+JmmkpCkUZ7GLS9mOl/TB9IfxGWAnUGvq5GpgYcfFS0njJM3qpaZeAvy+pMslvS7Vf7SkGyTtNdWQ2vt76caDcRRTe4+m0Uwtd1JcO6r2JUkHSnoH8F6Ka1Ed3kVxA0VnOquz7EaKa1fvrKq7kfr3Rc3jSjpX0rg0Euo4w/5NA/X9AmhW13cl9uQ1rud24C2SPiRpJPBXwMp008NeImJG6USg+jGjkzL/D1gBXJL+Rj4ATAH+tVZ+SX9BEbhPj4gtNbJMA9ZXjVaHBAeYoWkO0Exx5ng7xd0y95b2f5fiovw2imszH0zXY6p9HWgB7pH0S4ogdGJvNDAifkZxJ04zsErSdoo/8FbSVEU6C31HKjIe+GHa91PgFeADXRzin4Gz0kXZDv9D0eenKALuJzo+uCQdAUwG7qiuqNTm5cB2SV29BjcBJwM/johnOxLTh+VZFFN7OdQ8LjCd4vV9geL9nJ3uzqp+fat1BIUtkpZ3kqdbr3EjImIzxch2YarnRGB2o+W7YTbFLfHbKL7XdHY6NpLmSiqPZr5KMRJbWxodfb60fy7FydiQo72nCs2GBklfBZ6JiCsbyPu3wM8i4h/q5DsD+FREvL+bbflT4MiI+PPulBvoyq+xil8auCEiJvRzs/qMii/OPggc2xG4hxIHGDPrE0MxwAx1niIzswGl6g69Ri7K2wDlEYyZmWXhEYyZmWXR2ZfQhoSxY8dGc3NzfzfDzGxQWbZs2bMR0dmPf+42pANMc3Mzra2t/d0MM7NBRVJD3+nxFJmZmWXhAGNmZlk4wJiZWRZD+hqMmdlA8Otf/5r29nZeemlgfdl/5MiRTJgwgREjai0HVZ8DjJlZP2tvb+eQQw6hubkZSf3dHAAigi1bttDe3s7EiRN7VEfWKTJJ0yWtUbEs70U19l8oabWKJXTvKy8pKumHtVadkzRR0iOS1qpY1vfAlH5Q2m5L+5tz9s3MrLe89NJLjBkzZsAEFwBJjBkzZp9GVdkCjIo1wK+iWCNhMjBHxTKvZY8ClYiYQrE8bnmJ28vZc5XFDpcCV0TEJIpfOp2X0udRLGZ0NHBFymdmNigMpODSYV/blHMEMw1oi4h1EfEycDPFWu+7RcT9EbEjbS6hWNiqY999VK0gp6K37+bVtdqvp1gTnlR3x0+d3wacqoH4jpmZDRE5A8x49lyqt509l+WtNo/6iy2NAZ5La8tX17n7eGn/9pR/D5LmS2qV1Lp58+a6nTAzs57JGWBqjR5q/rKmpHMpFve5fB/qbOh4EbEoIioRURk3ru4vHZiZDTiLFy+mubmZAw44gObmZhYvXly/UD/IeRdZO6W13elkXXdJpwEXA++KiJ116nwWOEzS8DRKKdfZcbz2tNTvoUBny+WamQ1KixcvZv78+ezYUVxd2LBhA/Pnzwdg7ty5/dm0veQMMEuBSZImAj+nWIL0o+UMko4FrgGmR8Qz9SqMiJB0P3A2xTWd8yiW94Vi6d7zgIfT/h+H1yIws0HmM5/5DCtWrOh0/5IlS9i5c89z8R07djBv3jyuvfbammWmTp3KlVfWXbi112WbIksjjAuAu4EngO9ExCpJCyTNTNkuBw4GbpW0QlJLR3lJD1Gs+32qpHZJZ6ZdnwMulNRGcY3lmyn9m8CYlH4hsNdt0WZmg111cKmX3p+G9IJjlUol/GvKZtbfnnjiCd785jc3lLe5uZkNG/b+MeOmpibWr1/fyy2r3TZJyyKiUq+sf4vMzGwQWbhwIaNGjdojbdSoUSxcuLCfWtQ5Bxgzs0Fk7ty5LFq0iKamJiTR1NTEokWLBtwFfvBvkZmZDTpz584dkAGlmkcwZmaWhQOMmZll4QBjZjYADMQ7eve1TQ4wZmb9bOTIkWzZsmVABZmO9WBGjhzZ4zp8kd/MrJ9NmDCB9vZ2BtoP8HasaNlTDjBmZv1sxIgRPV41ciDzFJmZmWXhAGNmZlk4wJiZWRYOMGZmloUDjJmZZeEAY2ZmWTjAmJlZFlkDjKTpktZIapO01wqTki6UtFrSSkn3SWoq7TtP0tr0OC+lHZJWvux4PCvpyrTvfEmbS/s+nrNvZmbWtWxftJQ0DLgKOB1oB5ZKaomI1aVsjwKViNgh6ZPAZcA5kg4HLgEqQADLUtltwNTSMZYB/1aq75aIuCBXn8zMrHE5RzDTgLaIWBcRLwM3A7PKGSLi/ojYkTaXAB2/SXAmcG9EbE1B5V5germspEnAa4GHMvbBzMx6KGeAGQ9sKm23p7TOzAPu6kbZORQjlvKvw30oTbfdJunIWgeRNF9Sq6TWgfa7P2Zm+5OcAUY10mr+VKikcymmwy7vRtnZwE2l7e8BzRExBfgRcH2tY0XEooioRERl3LhxXTTfzMz2Rc4A0w6URxETgKeqM0k6DbgYmBkROxspK+ltwPCIWNaRFhFbSuWvBY7vjU6YmVnP5AwwS4FJkiZKOpBixNFSziDpWOAaiuDyTGnX3cAZkkZLGg2ckdI6zGHP0QuSjihtzgSe6LWemJlZt2W7iywidkm6gCIwDAOui4hVkhYArRHRQjEldjBwqySAjRExMyK2SvoyRZACWBARW0vVfwQ4q+qQn5Y0E9gFbAXOz9U3MzOrTwNpBbW+VqlUorW1tb+bYWY2qEhaFhGVevn8TX4zM8vCAcbMzLJwgDEzsywcYMzMLAsHGDMzy8IBxszMsnCAMTOzLBxgzMwsCwcYMzPLwgHGzMyycIAxM7MsHGDMzCwLBxgzM8vCAcbMzLJwgDEzsywcYMzMLIusAUbSdElrJLVJuqjG/gslrZa0UtJ9kppK+86TtDY9ziulP5DqXJEer03pB0m6JR3rEUnNOftmZmZdyxZgJA0DrgJmAJOBOZImV2V7FKhExBTgNuCyVPZw4BLgRGAacImk0aVycyNiano8k9LmAdsi4mjgCuDSTF0zM7MG5BzBTAPaImJdRLwM3AzMKmeIiPsjYkfaXAJMSM/PBO6NiK0RsQ24F5he53izgOvT89uAUyWpF/phZmY9kDPAjAc2lbbbU1pn5gF3NVj2W2l67AulILK7TETsArYDY6oPImm+pFZJrZs3b+5Of8zMrBtyBphao4eomVE6F6gAlzdQdm5EvBV4R3r8YXeOFxGLIqISEZVx48Z10XwzM9sXOQNMO3BkaXsC8FR1JkmnARcDMyNiZ72yEfHz9O8vgRsppuL2KCNpOHAosLWX+mJmZt2UM8AsBSZJmijpQGA20FLOIOlY4BqK4PJMadfdwBmSRqeL+2cAd0saLmlsKjsCeC/weCrTAnTcbXY28OOIqDliMjOz/Ibnqjgidkm6gCJYDAOui4hVkhYArRHRQjEldjBwa7qUsjEiZkbEVklfpghSAAtS2m9TBJoRqc4fAdemPN8E/kVSG8XIZXauvpmZWX0ayif5lUolWltb+7sZZmaDiqRlEVGpl8/f5DczsywcYMzMLAsHGDMzy8IBxszMsnCAMTOzLBxgzMwsCwcYMzPLwgHGzMyycIAxM7MsHGDMzCwLBxgzM8vCAcbMzLJwgDEzsywcYMzMLAsHGDMzyyJrgJE0XdIaSW2SLqqx/0JJqyWtlHSfpKbSvvMkrU2P81LaKEk/kPSkpFWS/rqU/3xJmyWtSI+P5+ybmZl1LVuAkTQMuAqYAUwG5kiaXJXtUaASEVOA24DLUtnDgUuAE4FpwCVp6WSAv4mINwHHAn8gaUapvlsiYmp6/FOuvpmZWX05RzDTgLaIWBcRLwM3A7PKGSLi/ojYkTaXABPS8zOBeyNia0RsA+4FpkfEjoi4P5V9GVheKmNmZgNIzgAzHthU2m5PaZ2ZB9zVaFlJhwHvA+4rJX8oTbfdJunIWgeRNF9Sq6TWzZs3N9YTMzPrtpwBRjXSomZG6VygAlzeSFlJw4GbgG9ExLqU/D2gOU23/Qi4vtaxImJRRFQiojJu3LiGOmJmZt2XM8C0A+VRxATgqepMkk4DLgZmRsTOBssuAtZGxJUdCRGxpVT+WuD4fe6BmZn1WM4AsxSYJGmipAOB2UBLOYOkY4FrKILLM6VddwNnSBqdLu6fkdKQ9BXgUOAzVXUdUdqcCTzRy/0xM7NuGJ6r4ojYJekCisAwDLguIlZJWgC0RkQLxZTYwcCtkgA2RsTMiNgq6csUQQpgQUqbQDHaeRJYnsr8fbpj7NOSZgK7gK3A+bn6ZmZm9Smi5mWRIaFSqURra2t/N8PMbFCRtCwiKvXy+Zv8ZmaWhQOMmZll4QBjZmZZOMCYmVkWDjBmZpaFA4yZmWXhAGNmZlk4wJiZWRYOMGZmloUDjJmZZeEAY2ZmWTjAmJlZFg4wZmaWhQOMmZll4QBjZmZZOMCYmVkWWQOMpOmS1khqk3RRjf0XSlotaaWk+yQ1lfadJ2ltepxXSj9e0k9Tnd9QWtZS0uGS7k35701LLfe6xYsX09zczAEHHEBzczOLFy/OcRgzsyz69DMsIrI8KJZJ/hnweuBA4DFgclWeU4BR6fkngVvS88OBdenf0en56LTvv4C3AwLuAmak9MuAi9Lzi4BL67Xx+OOPj+644YYbYtSoUQHsfowaNSpuuOGGbtVjZtYfeuszjGLZ+7pxINuSyZLeDnwxIs5M23+RAtrXOsl/LPD3EfEHkuYAJ0fEn6R91wAPpMf9EfGmlL47n6Q16fnTko4AHoiI3+uqjd1dMrm5uZkNGzbslX7QQQdx0kknNVyPmVl/WLJkCTt37twrvampifXr1zdcT6NLJg+vU8mlEfE5SbdSRLs9RMRHuig+HthU2m4HTuwi/zyKEUlnZcenR3uNdIDfiYinU7uelvTaWgeRNB+YD3DUUUd10Zy9bdy4sWZ6rTfMzGyg6eyzqrPPtn3VZYAB/kPS4cAjwMvA9m7UrRppNYdLks4FKsC76pRtuM7ORMQiYBEUI5julD3qqKNqjmCampp44IEHulOVmVmf62wWprsn242qd5F/FMUo4bPAAqA9Iq7veNQp2w4cWdqeADxVnUnSacDFwMyI2FmnbHt6XqvOX6SpMdK/z9RpX7ctXLiQUaNG7ZE2atQoFi5c2NuHMjPrdX3+GdbVBRpgJTA1Xr0g/0AjF3ZS/uEUF+cn8upF/mOq8hxLcSPApKr0w4H/prjAPzo9PzztWwqcxKsX+c9K6Zez50X+y+q1sbsX+SOKi2RNTU0hKZqamnyB38wGld74DKM3LvJLWhERU0vbj0bEsY0GL0lnAVdS3FF2XUQslLQgNa5F0o+AtwJPpyIbI2JmKvu/gc+n9IUR8a2UXgG+DfxWCjB/GhEhaQzwHeAoYCPw4YjY2lX7unuR38zMGr/IXy/ArAY+xKvXPm4rb0fE6n1vav9xgDEz675euYuM4hrMnVVpHdtB8R0XMzOzvXQZYCKiuY/aYWZm+xn/FpmZmWXhAGNmZlk4wJiZWRYOMGZmloUDjJmZZeEAY2ZmWTjAmJlZFg4wZmaWhQOMmZll4QBjZmZZOMCYmVkWDjBmZpaFA4yZmWXhAGNmZllkDTCSpktaI6lN0kU19r9T0nJJuySdXbXvUkmPp8c5pfSHJK1Ij6ck3ZHST5a0vbTvr3L2zczMulZvwbEekzQMuAo4HWgHlkpqqVoFcyNwPvBnVWXfAxwHTAUOAh6UdFdEPB8R7yjl+1fgu6WiD0XEe3P0x8zMuifnCGYa0BYR6yLiZeBmYFY5Q0Ssj4iVwCtVZScDD0bEroh4EXgMmF7OIOkQ4N3AHbk6YGZmPZczwIwHNpW221NaIx4DZkgaJWkscApwZFWeDwD3RcTzpbS3S3pM0l2SjqlVsaT5kloltW7evLnB5piZWXdlmyIDVCMtGikYEfdIOgH4CbAZeBjYVZVtDvBPpe3lQFNEvCDpLIqRzaQadS8CFgFUKpWG2mNmZt2XcwTTzp6jjgnAU40WjoiFETE1Ik6nCFZrO/ZJGkMxBfeDUv7nI+KF9PxOYEQa/ZiZWT/IGWCWApMkTZR0IDAbaGmkoKRhKYggaQowBbinlOXDwPcj4qVSmddJUno+jaJvW3qlJ2Zm1m3ZpsgiYpekC4C7gWHAdRGxStICoDUiWtI02O3AaOB9kr4UEccAI4CHUrx4Hjg3IspTZLOBv6465NnAJyXtAn4FzI4IT4GZmfUTDeXP4EqlEq2trf3dDDOzQUXSsoio1Mvnb/KbmVkWDjBmZpaFA4yZmWXhAGNmZlk4wJiZWRYOMGZmloUDjJmZZeEAY2ZmWTjAmJlZFg4wZmaWhQOMmZll4QBjZmZZOMCYmVkWDjBmZpaFA4yZmWWRNcBImi5pjaQ2SRfV2P9OScsl7ZJ0dtW+SyU9nh7nlNK/Lem/Ja1Ij6kpXZK+kY61UtJxOftmZmZdy7aipaRhwFXA6UA7sFRSS0SsLmXbCJwP/FlV2fcAxwFTgYOAByXdFRHPpyyfjYjbqg45A5iUHicC/5j+NTOzfpBzBDMNaIuIdRHxMnAzMKucISLWR8RK4JWqspOBByNiV0S8CDwGTK9zvFnAP0dhCXCYpCN6pSdmZtZtOQPMeGBTabs9pTXiMWCGpFGSxgKnAEeW9i9M02BXSDqoO8eTNF9Sq6TWzZs3N9oXMzPrppwBRjXSopGCEXEPcCfwE+Am4GFgV9r9F8CbgBOAw4HPded4EbEoIioRURk3blwjzTEzsx7IGWDa2XPUMQF4qtHCEbEwIqZGxOkUwWNtSn86TYPtBL5FMRW3z8czM7PelTPALAUmSZoo6UBgNtDSSEFJwySNSc+nAFOAe9L2EelfAe8HHk/FWoCPpbvJTgK2R8TTvdkhMzNrXLa7yCJil6QLgLuBYcB1EbFK0gKgNSJaJJ0A3A6MBt4n6UsRcQwwAnioiCE8D5wbER1TZIsljaMY1awAPpHS7wTOAtqAHcAf5eqbmZnVp4iGLovslyqVSrS2tvZ3M8zMBhVJyyKiUi+fv8lvZmZZOMCYmVkWDjBmZpaFA4yZmWXhAGNmZlk4wJiZWRYOMGZmloUDjJmZZeEAY2ZmWTjAmJlZFg4wZmaWhQOMmZll4QBjZmZZOMCYmVkWDjBmZpaFA4yZmWWRNcBImi5pjaQ2SRfV2P9OScsl7ZJ0dtW+SyU9nh7nlNIXpzofl3SdpBEp/WRJ2yWtSI+/ytk3MzPrWrYAI2kYcBUwA5gMzJE0uSrbRuB84Maqsu8BjgOmAicCn5X0mrR7MfAm4K3AbwEfLxV9KCKmpseC3u2RmZl1R84RzDSgLSLWRcTLwM3ArHKGiFgfESuBV6rKTgYejIhdEfEi8BgwPZW5MxLgv4AJGftgZmY9lDPAjAc2lbbbU1ojHgNmSBolaSxwCnBkOUOaGvtD4Iel5LdLekzSXZKOqVWxpPmSWiW1bt68udG+mJlZNw3PWLdqpEUjBSPiHkknAD8BNgMPA7uqsv0D8O8R8VDaXg40RcQLks4C7gAm1ah7EbAIoFKpNNQeMzPrvpwjmHb2HHVMAJ5qtHBELEzXUk6nCFZrO/ZJugQYB1xYyv98RLyQnt8JjEijHzMz6wc5A8xSYJKkiZIOBGYDLY0UlDRM0pj0fAowBbgnbX8cOBOYExGvlMq8TpLS82kUfdvSi/0xM7NuyDZFFhG7JF0A3A0MA66LiFWSFgCtEdGSpsFuB0YD75P0pYg4BhgBPJTixfPAuRHRMUV2NbABeDjt/7d0x9jZwCcl7QJ+BcxONwKYmVk/0FD+DK5UKtHa2trfzTAzG1QkLYuISr18/ia/mZll4QBjZmZZOMCYmVkWDjBmZpaFA4yZmWXhAGNmZlk4wJiZWRYOMGZmloUDjJmZZeEAY2ZmWTjAmJlZFg4wZmaWhQOMmZll4QBjZmZZOMCYmVkWDjBmZpZF1gAjabqkNZLaJF1UY/87JS2XtEvS2VX7LpX0eHqcU0qfKOkRSWsl3ZKWY0bSQWm7Le1vztk3MzPrWrYAI2kYcBUwA5gMzJE0uSrbRuB84Maqsu8BjgOmAicCn5X0mrT7UuCKiJgEbAPmpfR5wLaIOBq4IuUzM7N+knMEMw1oi4h1EfEycDMwq5whItZHxErglaqyk4EHI2JXRLwIPAZMlyTg3cBtKd/1wPvT81lpm7T/1JTfzMz6Qc4AMx7YVNpuT2mNeAyYIWmUpLHAKcCRwBjguYjYVaPO3cdL+7en/HuQNF9Sq6TWzZs3d7NLZmbWqJwBptboIRopGBH3AHcCPwFuAh4GdtWps6HjRcSiiKhERGXcuHGNNMfMzHogZ4Bppxh1dJgAPNVo4YhYGBFTI+J0iuCxFngWOEzS8Bp17j5e2n8osHWfemBmZj2WM8AsBSalu74OBGYDLY0UlDRM0pj0fAowBbgnIgK4H+i44+w84LvpeUvaJu3/ccpvZmb9IFuASddBLgDuBp4AvhMRqyQtkDQTQNIJktqBDwPXSFqVio8AHpK0GlgEnFu67vI54EJJbRTXWL6Z0r8JjEnpFwJ73RZtZmZ9R0P5JL9SqURra2t/N8PMbFCRtCwiKvXy+Zv8ZmaWhQOMmZll4QBjZmZZOMCYmVkWQ/oiv6TNwIYeFh9L8b2cocR9Hhrc56FhX/rcFBF1v6k+pAPMvpDU2shdFPsT93locJ+Hhr7os6fIzMwsCwcYMzPLwgGm5xb1dwP6gfs8NLgsRodmAAAEvElEQVTPQ0P2PvsajJmZZeERjJmZZeEAY2ZmWTjA1CFpuqQ1ktok7fULzZIOknRL2v+IpOa+b2XvaqDPF0paLWmlpPskNfVHO3tTvT6X8p0tKSQN+ltaG+mzpI+k93qVpBv7uo29rYH/20dJul/So+n/91n90c7eIuk6Sc9IeryT/ZL0jfR6rJR0XK82ICL86OQBDAN+BrweOJBiKefJVXk+BVydns8GbunvdvdBn08BRqXnnxwKfU75DgH+HVgCVPq73X3wPk8CHgVGp+3X9ne7+6DPi4BPpueTgfX93e597PM7geOAxzvZfxZwF8WijicBj/Tm8T2C6do0oC0i1kXEy8DNwKyqPLOA69Pz24BTJdVavnmwqNvniLg/InakzSUUK4sOZo28zwBfBi4DXurLxmXSSJ//GLgqIrYBRMQzfdzG3tZInwN4TXp+KN1YhXcgioh/p+uVfWcB/xyFJRQrBh/RW8d3gOnaeGBTabs9pdXME8WiaNspFkIbrBrpc9k8ijOgwaxunyUdCxwZEd/vy4Zl1Mj7/EbgjZL+U9ISSdP7rHV5NNLnLwLnpoUQ7wT+tG+a1m+6+/feLcPrZxnSao1Equ/rbiTPYNJwfySdC1SAd2VtUX5d9lnSAcAVwPl91aA+0Mj7PJximuxkilHqQ5LeEhHPZW5bLo30eQ7w7Yj4W0lvB/4l9fmV/M3rF1k/vzyC6Vo7cGRpewJ7D5l355E0nGJY3dWQdKBrpM9IOg24GJgZETv7qG251OvzIcBbgAckraeYq24Z5Bf6G/2//d2I+HVE/DewhiLgDFaN9Hke8B2AiHgYGEnxo5D7q4b+3nvKAaZrS4FJkiZKOpDiIn5LVZ4W4Lz0/Gzgx5Gung1SdfucpouuoQgug31eHur0OSK2R8TYiGiOiGaK604zI2Iwr7fdyP/tOyhu6EDSWIops3V92sre1UifNwKnAkh6M0WA2dynrexbLcDH0t1kJwHbI+Lp3qrcU2RdiIhdki4A7qa4A+W6iFglaQHQGhEtwDcphtFtFCOX2f3X4n3XYJ8vBw4Gbk33M2yMiJn91uh91GCf9ysN9vlu4AxJq4HfAJ+NiC391+p902Cf/y9wraT/QzFVdP5gPmGUdBPFFOfYdF3pEmAEQERcTXGd6SygDdgB/FGvHn8Qv3ZmZjaAeYrMzMyycIAxM7MsHGDMzCwLBxgzM8vCAcbMzLLwbcpmmUl6BDiI4gcW3wh0/LLtNuCpiJjbX20zy8m3KZv1kbSUQ2tE7M/fDDfbzVNkZv1E0smSWtPzZknPSvpaWovkSUnHS7o2rdPxiKTXlcr+uaT/krRc0vfK+8wGCgcYs4FjDPAfEXEsxS9E3Efxc/lTgGXABbD7R0aPBk6KiOMovo39t/3TZLPO+RqM2cDxQkT8ID1fDrRHxIq0vQw4PT2fSfEr1svTT/UMp1gmwmxAcYAxGzjKv0r9G/Zc2Ow3vPr3KuArEXFdXzXMrCc8RWY2+LQAn5I0GkDSQZLe1s9tMtuLRzBmg0xE/Ev6+fwH0xTZAcA/UKwxbzZg+DZlMzPLwlNkZmaWhQOMmZll4QBjZmZZOMCYmVkWDjBmZpaFA4yZmWXhAGNmZln8f78uDJyuFThcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pis_1, itx1 = calc_pis(0.2, Nr, N)\n",
    "myplot.plt_plot([np.arange(itx1)], [pis_1], 'plot', ['k', 'g', 'r', 'b', 'y'], ['o', 'x', '.', ',', 's'], \n",
    "                [' ', '', '', '', ''], \n",
    "        title = 'Problem C.3 (b)-(iv) Pi vs. t (pi_0 = 0.2)', yscale = None, ylb = 0, yub = 1.2,\n",
    "        xlb = None, xub = None, xlabel = 'Time', ylabel = 'Pi',\n",
    "        legends = ['']*6, legendx = None, legendy = None, marker_sizes=[15]*5) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem C.4\n",
    "\n",
    "* (a) In our assumed model for the temperature, we have $w$ as our parameters. And $y=h(x) + \\epsilon = w_0 + w_1x+\\epsilon$, $y_n = w_0 +w_1x_n + epsilon ~ N(w_0 + w_1x_n, 1)$ so \n",
    "\n",
    "$P[y=y_n|w] = \\int P[y=y_n|x=x_n,w]dP[x=x_n] =\\int^1_0 dx \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}$\n",
    "\n",
    "Thus the log-likelihood is \n",
    "\n",
    "$\\ln P[y|w] = \\sum^N_{n=1}\\ln\\left(\\int^1_0 dx \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}\\right)$\n",
    "\n",
    "* (b) The natural hidden variable $j_n$ will be the time that $x_n$ is measured.\n",
    "\n",
    "* (c) $P[y, j_n|w] = \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x_n)^2}$, thus the log-likelihood for the complete data is\n",
    "\n",
    "$L = \\ln P[y, J|w] =  N\\ln \\frac{1}{\\sqrt{2\\pi}}  + \\sum^N_{n=1} -\\frac{1}{2}(y_n -w_0 - w_1x_n)^2 $\n",
    "\n",
    "* (d) \n",
    "\n",
    "\\begin{align*}\n",
    "\\gamma_n(x|w) &= P(x_n = x|y_n, w) \\\\\n",
    "&= \\frac{P(y_n, x_n }{P(y_n}\\\\\n",
    "&= \\frac{P(y_n| x_n)P(x_n) }{P(y_n}\\\\\n",
    "&= \\frac{ \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x_n)^2}}{\\int^1_0 dx \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}}\\\\\n",
    "&= \\frac{ e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}}{\\int^1_0 dx e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}}\\\\\n",
    "\\end{align*}\n",
    "\n",
    "We then compute $Q_t(w)$, \n",
    "\n",
    "\\begin{align*}\n",
    "Q_t(w) &= E[L] \\\\\n",
    "&=E\\left[N\\ln \\frac{1}{\\sqrt{2\\pi}}  + \\sum^N_{n=1} -\\frac{1}{2}(y_n -w_0 - w_1x_n)^2\\right]\\\\\n",
    "&=N\\ln \\frac{1}{\\sqrt{2\\pi}} -\\frac{1}{2}\\sum^N_{n=1} E\\left[(y_n -w_0 - w_1x_n)^2\\right]\\\\\n",
    "&=N\\ln \\frac{1}{\\sqrt{2\\pi}}  -\\frac{1}{2} \\sum^N_{n=1} E\\left[(y_n -w_0)^2 - 2(y_n-w_0)w_1x_n+ w^2_1x_n^2\\right]\\\\\n",
    "&=N\\ln \\frac{1}{\\sqrt{2\\pi}}  -\\frac{1}{2}\\sum^N_{n=1}(y_n -w_0)^2 + \\sum^N_{n=1}(y_n-w_0)w_1E[x_n] -\\frac{1}{2} \\sum^N_{n=1}w^2_1E[x_n^2]\\\\\n",
    "&=N\\ln \\frac{1}{\\sqrt{2\\pi}}  -\\frac{1}{2}\\sum^N_{n=1}(y_n -w_0)^2 + w_1\\sum^N_{n=1}(y_n-w_0)\\alpha_n -\\frac{1}{2} w^2_1\\sum^N_{n=1}\\beta_n\\\\\n",
    "&=N\\ln \\frac{1}{\\sqrt{2\\pi}}  -\\frac{1}{2}\\sum^N_{n=1}(y_n -w_0)^2 + w_1\\sum^N_{n=1}(y_n-w_0)\\alpha_n -\\frac{1}{2} w^2_1N\\bar{\\beta}\\\\\n",
    "\\end{align*}\n",
    "\n",
    "* (e) \n",
    "Take derivative w.r.t. $w_0,w_1$ and set to 0, we have\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial{Q_t}}{\\partial{w_0}} &= \\sum^N_{n=1}(y_n-w_0) - w_1\\sum^N_{n=1}\\alpha_n \\\\\n",
    "\\frac{\\partial{Q_t}}{\\partial{w_1}} &= \\sum^N_{n=1}(y_n-w_0)\\alpha_n -w_1N\\bar{\\beta}\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Solve the equations we have \n",
    "\n",
    "$w_1(t+1) = \\frac{\\frac{1}{N}\\sum^N_{i=1}(y_i-\\bar{y})(\\alpha_i - \\bar{\\alpha})}{\\bar{\\beta}-\\bar{\\alpha}^2}$\n",
    "\n",
    "$w_0(t+1) = \\bar{y} - w_1(t+1)\\bar{\\alpha}$\n",
    "\n",
    "* (f) If the temperature measurement is distrubted according to $P(x)$, then clearly we have\n",
    "\n",
    "\\begin{align*}\n",
    "\\gamma_n(x|w) &= P(x_n = x|y_n, w) \\\\\n",
    "&= \\frac{P(y_n, x_n }{P(y_n}\\\\\n",
    "&= \\frac{P(y_n| x_n)P(x_n) }{P(y_n}\\\\\n",
    "&= \\frac{ P_t(x)\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x_n)^2}}{\\int^1_0 dx P_t(x)\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}}\\\\\n",
    "&= \\frac{ P_t(x)e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}}{\\int^1_0 dx P_t(x)e^{-\\frac{1}{2}(y_n -w_0 - w_1x)^2}}\\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "The updates in (e) are unchanged except that they use this new $\\gamma_n(x|w_t)$\n",
    "\n",
    "With $P(x)$, the $P[y, j_n|w] = P(x)\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}(y_n -w_0 - w_1x_n)^2}$, \n",
    "Let $\\hat{N}_k$ be the number of occurance of bump $k$ in the random variable $J$. Let $z_{nk}=1$ if $x_n$ is from bump $k$ and 0 otherwise. \n",
    "\n",
    "thus the log-likelihood for the complete data is\n",
    "\n",
    "$L = \\ln P[y, J|w] =  N\\ln \\frac{1}{\\sqrt{2\\pi}}  + \\sum^K_{k=1}\\hat{N}_k\\ln P(x) + \\sum^N_{n=1} -\\frac{1}{2}(y_n -w_0 - w_1x_n)^2 $\n",
    "\n",
    "The expectation will have an extra term $E[\\sum^K_{k=1}\\hat{N}_k\\ln P(x)] = \\ln P(x) \\sum^N_{n=1}E[\\hat{N}_k] =\\sum^N_{n=1} E[\\sum^N_{n=1}z_{nk}] \\ln P(x) = \\sum^N_{n=1} \\sum^N_{n=1}\\gamma_{nk}\\ln P(x) $, \n",
    "\n",
    "This achieves maximum when $P(x) = \\frac{N_k}{N} = \\frac{1}{N}\\sum^N_{n=1}\\gamma_{n}$\n",
    "\n",
    "If we maximize the log-likelihood for the incomplete data, we won't be able to separate $P(x)$ from $w$ during our maximization, thus making the problem very hard to solve."
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}