ModelingExamples2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation of a Categorical distribution\n",
    "\n",
    "## Maximum Likelihood Estimation\n",
    "\n",
    "We observe a dataset $\\{x^{(n)}\\}_{n=1\\dots N}$. The model for a single observation is a categorical distribution with parameter $\\pi = (\\pi_1, \\dots, \\pi_I)$ where \n",
    "\n",
    "\\begin{eqnarray}\n",
    "x^{(n)} & \\sim & p(x|\\pi) = \\prod_{i=1}^{I} \\pi_i^{\\ind{i = x^{(n)}}}\n",
    "\\end{eqnarray}\n",
    "where $\\sum_i \\pi_i  = 1$.\n",
    "\n",
    "The loglikelihood of the entire dataset is\n",
    "\n",
    "\\begin{eqnarray}\n",
    "{\\cal L}(\\pi_1,\\dots,\\pi_I) & = & \\sum_{n=1}^N\\sum_{i=1}^I \\ind{i = x^{(n)}} \\log \\pi_i\n",
    "\\end{eqnarray}\n",
    "This is a constrained optimisation problem.\n",
    "Form the Lagrangian\n",
    "\\begin{eqnarray}\n",
    "\\Lambda(\\pi, \\lambda) & = & \\sum_{n=1}^N\\sum_{i'=1}^I \\ind{i' = x^{(n)}} \\log \\pi_{i'}  + \\lambda \\left( 1 - \\sum_{i'} \\pi_{i'} \\right ) \\\\\n",
    "\\frac{\\partial \\Lambda(\\pi, \\lambda)}{\\partial \\pi_i} & = & \\sum_{n=1}^N \\ind{i = x^{(n)}} \\frac{1}{\\pi_i} - \\lambda = 0 \\\\\n",
    "\\pi_i & = & \\frac{\\sum_{n=1}^N \\ind{i = x^{(n)}}}{\\lambda}\n",
    "\\end{eqnarray}\n",
    "\n",
    "We solve for $\\lambda$\n",
    "\\begin{eqnarray}\n",
    "1 & = & \\sum_i \\pi_i = \\frac{\\sum_{i=1}^I \\sum_{n=1}^N \\ind{i = x^{(n)}}}{\\lambda} \\\\\n",
    "\\lambda & = & \\sum_{i=1}^I \\sum_{n=1}^N \\ind{i = x^{(n)}} =  \\sum_{n=1}^N 1 = N\n",
    "\\end{eqnarray}\n",
    "\n",
    "Hence\n",
    "\\begin{eqnarray}\n",
    "\\pi_i & = & \\frac{\\sum_{n=1}^N \\ind{i = x^{(n)}}}{N}\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table width=100%><tr><td align=\"center\">\\begin{eqnarray}\\mathcal{D}(\\pi_{1:I}; a_{1:I} )\\end{eqnarray}</td><td>\\mathcal{D}(\\pi_{1:I}; a_{1:I} )</td></tr><tr><td align=\"center\">\\begin{eqnarray}\\frac{\\Gamma(\\sum_{i} a_{i})}{\\prod_{i} \\Gamma(a_{i})} \\prod_{{i}=1}^{I} {\\pi}_{i}^{a_{i} - 1} \\end{eqnarray}</td><td>\\frac{\\Gamma(\\sum_{i} a_{i})}{\\prod_{i} \\Gamma(a_{i})} \\prod_{{i}=1}^{I} {\\pi}_{i}^{a_{i} - 1} </td></tr><tr><td align=\"center\">\\begin{eqnarray}\\exp\\left(\\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})} + \\sum_{{i}=1}^{I} (a_{i} - 1) \\log{\\pi}_{i} \\right)\\end{eqnarray}</td><td>\\exp\\left(\\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})} + \\sum_{{i}=1}^{I} (a_{i} - 1) \\log{\\pi}_{i} \\right)</td></tr><tr><td align=\"center\">\\begin{eqnarray}\\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})} + \\sum_{{i}=1}^{I} (a_{i} - 1) \\log{\\pi}_{i} \\end{eqnarray}</td><td>\\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})} + \\sum_{{i}=1}^{I} (a_{i} - 1) \\log{\\pi}_{i} </td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# %load template_equations.py\n",
    "from IPython.display import display, Math, Latex, HTML\n",
    "import notes_utilities as nut\n",
    "from importlib import reload\n",
    "reload(nut)\n",
    "Latex('$\\DeclareMathOperator{\\trace}{Tr}$')\n",
    "\n",
    "L = nut.pdf2latex_dirichlet(x=r'\\pi', a=r'a',N=r'I', i='i')\n",
    "\n",
    "display(HTML(nut.eqs2html_table(L)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Maximum A-Posteriori Estimation\n",
    "\n",
    "$$\n",
    "\\pi \\sim \\mathcal{D}(\\pi_{1:I}; a_{1:I} )\n",
    "$$\n",
    "where $\\sum_i \\pi_i  = 1$. For $n = 1\\dots N$\n",
    "\\begin{eqnarray}\n",
    "x^{(n)} & \\sim & p(x|\\pi) = \\prod_{i=1}^{I} \\pi_i^{\\ind{i = x^{(n)}}}\n",
    "\\end{eqnarray}\n",
    "$X = \\{x^{(1)},\\dots,x^{(N)} \\}$\n",
    "\n",
    "The posterior is \n",
    "\\begin{align}\n",
    "\\log p(\\pi_{1:I}| X) & =^+  \\log p(\\pi_{1:I}, X) \\\\\n",
    "& = \\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})} \n",
    "+ \\sum_{{i}=1}^{I} (a_{i} - 1) \\log{\\pi}_{i}  \n",
    "+ \\sum_{i=1}^I\\sum_{n=1}^N \\ind{i = x^{(n)}} \\log \\pi_i \\\\\n",
    " & =^+    \\sum_{i=1}^I \\left(a_i - 1 +  \\sum_{n=1}^N \\ind{i = x^{(n)}}\\right) \\log \\pi_i \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finding the parameter vector $\\pi_{1:I}$ that maximizes the posterior density is a constrained optimisation problem. After omitting constant terms that do not depend on $\\pi$, we form the Lagrangian\n",
    "\\begin{eqnarray}\n",
    "\\Lambda(\\pi, \\lambda) & = & \\sum_{i=1}^I \\left(a_i - 1 +  \\sum_{n=1}^N \\ind{i = x^{(n)}}\\right) \\log \\pi_i   + \\lambda \\left( 1 - \\sum_{i'} \\pi_{i'} \\right ) \\\\\n",
    "\\frac{\\partial \\Lambda(\\pi, \\lambda)}{\\partial \\pi_i} & = & \\left(a_i - 1 + \\sum_{n=1}^N \\ind{i = x^{(n)}}\\right) \\frac{1}{\\pi_i} - \\lambda = 0 \\\\\n",
    "\\pi_i & = & \\frac{a_i - 1 + \\sum_{n=1}^N \\ind{i = x^{(n)}}}{\\lambda}\n",
    "\\end{eqnarray}\n",
    "\n",
    "\n",
    "We solve for $\\lambda$\n",
    "\\begin{eqnarray}\n",
    "1 & = & \\sum_i \\pi_i = \\frac{- I + \\sum_{i=1}^I \\left( a_i + \\sum_{n=1}^N \\ind{i = x^{(n)}  \\right) }}{\\lambda} \\\\\n",
    "\\lambda & = &  N - I + \\sum_{i=1}^I a_i\n",
    "\\end{eqnarray}\n",
    "\n",
    "Setting the count of observations equal to $i$ as $C_i \\equiv \\sum_{n=1}^N \\ind{i = x^{(n)}}$, we obtain \n",
    "\n",
    "Hence\n",
    "\\begin{eqnarray}\n",
    "\\pi_i & = & \\frac{C_i + a_i - 1}{N + \\sum_{i=1}^I a_i - I}\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full Bayesian Inference\n",
    "\n",
    "Setting the count of observations equal to $i$ as $C_i \\equiv \\sum_{n=1}^N \\ind{i = x^{(n)}}$, we obtain \n",
    "\n",
    "\n",
    "The posterior is \n",
    "\\begin{eqnarray}\n",
    "\\log p(\\pi_{1:I}, X) & = & \\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})} +  \\sum_{i=1}^I \\left(\\left(a_i + \\sum_{n=1}^N \\ind{i = x^{(n)}}\\right) - 1  \\right) \\log \\pi_i \\\\\n",
    "& =  & \\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})}\\\\\n",
    "& &  + \\sum_{i=1}^I (a_i + C_i - 1) \\log \\pi_i \\\\\n",
    " & & + \\log{\\Gamma(\\sum_{i} (a_{i} + C_i) )} - {\\sum_{i} \\log \\Gamma(a_{i} + C_i)} \\\\\n",
    " & & - \\log{\\Gamma(\\sum_{i} (a_{i} + C_i) )} + {\\sum_{i} \\log \\Gamma(a_{i} + C_i)} \\\\\n",
    " & = & \\log \\mathcal{D}(\\pi_{1:I}, a_{1:I} ) + \\log p(X) \n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\log p(X) & = & \\log{\\Gamma(\\sum_{i} a_{i})} - {\\sum_{i} \\log \\Gamma(a_{i})} - \\log{\\Gamma(\\sum_{i} (a_{i} + C_i) )} + {\\sum_{i} \\log \\Gamma(a_{i} + C_i)}\n",
    "\\end{eqnarray}\n",
    "\n",
    "$a_+ = \\sum_{i} a_{i}$ and $C_+ = \\sum_{i}  C_i$\n",
    "\\begin{eqnarray}\n",
    "p(X) & = & \\frac{\\Gamma(a_+) \\prod_{i} \\Gamma(a_{i} + C_i)}{\\Gamma(a_+ + C_+ ) \\prod_{i} \\Gamma(a_{i})} = \\frac{ \\left(\\prod_{i} \\Gamma(a_{i} + C_i)\\right)/ \\Gamma(a_+ + C_+ )}{ \\left(\\prod_{i} \\Gamma(a_{i}) \\right)/ \\Gamma(a_+)} = \\frac{B(a + C)}{B(a)}\n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Are $x$ and $y$ independent?\n",
    "Suppose we observe a dataset of $(x, y)$ pairs where\n",
    "$x \\in \\{1,\\dots,I_1\\}$ and\n",
    "$y \\in \\{1,\\dots,I_2\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|  |  |\n",
    "| --- | --- |\n",
    "| $x^{(1)}$ | $y^{(1)}$ | \n",
    "| $\\vdots$ | $\\vdots$ | \n",
    "| $x^{(n)}$ | $y^{(n)}$ | \n",
    "| $\\vdots$ | $\\vdots$ | \n",
    "| $x^{(N)}$ | $y^{(N)}$ | "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\ind}[1]{\\left[#1\\right]}$\n",
    "We are given the counts of observations where $x = i$ while $y = j$. These counts can be stored as an array, that is known as a contingency table\n",
    "$\n",
    "S_{i,j} = \\sum_{n=1}^N \\ind{x^{(n)} = i}\\ind{y^{(n)} = j}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "hide_input": true,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<TABLE cellpadding=\"4\" border=\"1\" style=\"border: 1px solid #000000; border-collapse: collapse;\">\n",
       " <TR>\n",
       "  <TD width=\"65px\">&nbsp;</TD>\n",
       "  <TD width=\"65px\">$y=1$</TD>\n",
       "  <TD width=\"65px\">$y=j$</TD>\n",
       "  <TD width=\"80px\">$y=I_2$</TD>\n",
       " </TR>\n",
       " <TR>\n",
       "  <TD>$x=1$</TD>\n",
       "  <TD>$S_{1,1}$</TD>\n",
       "  <TD>&nbsp;</TD>\n",
       "  <TD>$S_{1,I_2}$</TD>\n",
       " </TR>\n",
       " <TR>\n",
       "  <TD>$x=i$</TD>\n",
       "  <TD>&nbsp;</TD>\n",
       "  <TD>$S_{i,j}$</TD>\n",
       "  <TD>&nbsp;</TD>\n",
       " </TR>\n",
       " <TR>\n",
       "  <TD>$x=I_1$</TD>\n",
       "  <TD>$S_{I_1,1}$</TD>\n",
       "  <TD>&nbsp;</TD>\n",
       "  <TD>$S_{I_1,I_2}$</TD>\n",
       " </TR>\n",
       "</TABLE>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Math, Latex, HTML\n",
    "import html_utils as htm\n",
    "wd = '65px'\n",
    "L = [[htm.TableCell('', width=wd), htm.TableCell('$y=1$', width=wd), htm.TableCell('$y=j$', width=wd), htm.TableCell('$y=I_2$', width='80px')],\n",
    "     [r'$x=1$',r'$S_{1,1}$',r'',r'$S_{1,I_2}$'],\n",
    "     [r'$x=i$',r'',r'$S_{i,j}$',r''],\n",
    "     [r'$x=I_1$',r'$S_{I_1,1}$',r'',r'$S_{I_1,I_2}$']]\n",
    "\n",
    "t = htm.make_htmlTable(L)\n",
    "display(HTML(str(t)))\n",
    "#print(str(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is deciding if the random variables $x$ and $y$ are independent or dependent, given some observations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| |  || |\n",
    "|-|-    |-|-|\n",
    "| |$y$|| |\n",
    "|$x$ |$3$|$5$|$9$|\n",
    "| |$7$|$9$|$17$|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Independent model $M_1$\n",
    "\n",
    "\\begin{equation}\n",
    "p(x, y) = p(x) p(y) \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{align}\n",
    "\\pi_1 & \\sim  \\mathcal{D}(\\pi_1; a_1) &\n",
    "\\pi_2 & \\sim  \\mathcal{D}(\\pi_2; a_2) \\\\\n",
    "x^{(n)} & \\sim  \\mathcal{C}(x; \\pi_1) &\n",
    "y^{(n)} & \\sim  \\mathcal{C}(y; \\pi_2) \n",
    "\\end{align}\n",
    "\n",
    "We let \n",
    "\n",
    "* $X_{i+} = \\sum_j X_{i,j}$ \n",
    "* $X_{+j} = \\sum_i X_{i,j}$\n",
    "\n",
    "The marginal likelihood can be found as \n",
    "\\begin{eqnarray}\n",
    "\\log p(X|M_1) & = & \\log{\\Gamma(\\sum_{i} a_{1}(i))} - {\\sum_{i} \\log \\Gamma(a_{1}(i)} - \\log{\\Gamma(\\sum_{i} (a_{1}(i) + \\sum_{j} C(i,j)) )} + {\\sum_{i} \\log \\Gamma(a_{1}(i) + \\sum_{j} C(i,j))} \\\\\n",
    "& & +\\log{\\Gamma(\\sum_{j} a_{2}(j)} - {\\sum_{j} \\log \\Gamma(a_{2}(j))} - \\log{\\Gamma(\\sum_{j} (a_{2}(j) + \\sum_{i} C(i,j)) )} + {\\sum_{j} \\log \\Gamma(a_{2}(j) + \\sum_{i} C(i,j))} \\\\\n",
    "& = & \\log{\\Gamma(A_1)} - \\sum_{i} \\log \\Gamma(a_{1}(i)) - \\log{\\Gamma(A_1+ N)} + {\\sum_{i} \\log \\Gamma(a_{1}(i) + C_1(i))} \\\\\n",
    "& & + \\log{\\Gamma(A_2)} - {\\sum_{j} \\log \\Gamma(a_{2}(j))} - \\log{\\Gamma(A_2 + N )} + {\\sum_{j} \\log \\Gamma(a_{2}(j) +  C_2(j))} \\\\\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dependent model $M_2$\n",
    "\\begin{equation}\n",
    "p(x_1, x_2)\n",
    "\\end{equation}\n",
    "\n",
    "$\\pi_{1,2}$ is a $S_1 \\times S_2$ matrix where the joint distribution of entries is Dirichlet $\\mathcal{D}(\\pi_{1,2}; a_{1,2})$ with $S_1 \\times S_2$ parameter matrix $a_{1,2}$. Then, the probability that $p(x_1 = i, x_2 = j|\\pi_{1,2}) = \\pi_{1,2}(i,j)$.\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\pi_{1,2} & \\sim & \\mathcal{D}(\\pi_{1,2}; a_{1,2}) \\\\\n",
    "(x_1, x_2)^{(n)} & \\sim & \\mathcal{C}((x_1,x_2); \\pi_{1,2}) \\\\\n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\log p(X|M_2) & = & \\log{\\Gamma(A_{1,2})} - {\\sum_{i,j} \\log \\Gamma(a_{1,2}(i,j))} - \\log{\\Gamma(A_{1,2}+ N)} + {\\sum_{i,j} \\log \\Gamma(a_{1,2}(i,j) + C(i,j))} \n",
    "\\end{eqnarray}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dependent model $M_3$\n",
    "\\begin{equation}\n",
    "p(x_1, x_2) = p(x_1) p(x_2|x_1) \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\pi_1 & \\sim & \\mathcal{D}(\\pi_1; a_1) \\\\\n",
    "\\pi_{2,1} & \\sim & \\mathcal{D}(\\pi_2; a_2) \\\\\n",
    "\\vdots \\\\\n",
    "\\pi_{2,S_1} & \\sim & \\mathcal{D}(\\pi_2; a_2) \\\\\n",
    "x_1^{(n)} & \\sim & \\mathcal{C}(x_1; \\pi_1) \\\\\n",
    "x_2^{(n)} & \\sim & \\mathcal{C}(x_2; \\pi_{2}(x_1^{(n)},:)) \n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\log p(x_1^{(1:N)}|\\pi_1) & = & \\sum_n \\sum_i \\sum_j \\ind{x_1^{(n)} = i} \\ind{x_2^{(n)} = j} \\log \\pi_{1}(i)  =  \\sum_i \\sum_j C(i,j) \\log \\pi_{1}(i)\n",
    "\\end{eqnarray}\n",
    "\\begin{eqnarray}\n",
    "\\log p(x_2^{(1:N)}|\\pi_2, x_1^{(1:N)} ) & = & \\sum_n \\sum_i \\sum_j \\ind{x_1^{(n)} = i} \\ind{x_2^{(n)} = j} \\log \\pi_{2}(i,j) =  \\sum_i \\sum_j C(i,j) \\log \\pi_{2}(i,j)\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{eqnarray}\n",
    "\\log p(\\pi_1) & = & \\log{\\Gamma(\\sum_{i} a_{1}(i))} - {\\sum_{i} \\log \\Gamma(a_{1}(i))} + \\sum_{{i}=1}^{S_1} (a_{1}(i) - 1) \\log{\\pi_1}(i) \n",
    "\\end{eqnarray}\n",
    "\\begin{eqnarray}\n",
    "\\log p(\\pi_2) & = & \\sum_i \\left(\\log{\\Gamma(\\sum_{j} a_{2}(i,j))} - {\\sum_{j} \\log \\Gamma(a_{2}(i,j))} + \\sum_{{j}=1}^{S_2} (a_{2}(i,j) - 1) \\log{\\pi_2}(i,j) \\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "The joint distribution is\n",
    "\\begin{eqnarray}\n",
    "\\log p(X, \\pi| M_2)&= & \\log{\\Gamma(\\sum_{i} a_{1}(i))} - {\\sum_{i} \\log \\Gamma(a_{1}(i))} + \\sum_{{i}=1}^{S_1} (a_{1}(i) + C_1(i) - 1) \\log{\\pi_1}(i)  \\\\\n",
    "& & + \\sum_i \\left(\\log{\\Gamma(\\sum_{j} a_{2}(i,j))} - {\\sum_{j} \\log \\Gamma(a_{2}(i,j))} + \\sum_{{j}=1}^{S_2} (a_{2}(i,j) + C(i,j) - 1) \\log{\\pi_2}(i,j) \\right) \n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will assume $a_2(i,j) = a_2(i',j)$ for all $i$ and $i'$.\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\log p(X| M_2) & = & \\log{\\Gamma(\\sum_{i} a_{1}(i))} \n",
    "- {\\sum_{i} \\log \\Gamma(a_{1}(i))}\n",
    "- \\log{\\Gamma(\\sum_{i} a_{1}(i) + C_1(i))} \n",
    "+ \\sum_{i} \\log \\Gamma(a_{1}(i) + C_1(i)) \\\\ \n",
    "& & + \\sum_i \\left( \\log\\Gamma(\\sum_{j} a_{2}(i,j)) - \\sum_{j} \\log \\Gamma(a_{2}(i,j)) - \\log\\Gamma( \\sum_{j} a_{2}(i,j) + C(i,j)) + \\sum_j \\log\\Gamma( a_{2}(i,j) + C(i,j) ) \\right) \\\\\n",
    "& = & \\log{\\Gamma(A_1)} - {\\sum_{i} \\log \\Gamma(a_{1}(i))} - \\log{\\Gamma(A_1+ N)} + {\\sum_{i} \\log \\Gamma(a_{1}(i) + C_1(i))} \\\\\n",
    "& & + \\sum_i \\left( \\log{\\Gamma(A_2)} - {\\sum_{j} \\log \\Gamma(a_{2}(j))} - \\log{\\Gamma(A_2 + C_1(i) )} + {\\sum_{j} \\log \\Gamma(a_{2}(j) + C(i,j))} \\right) \n",
    "\\end{eqnarray}\n",
    "\n",
    "\n",
    "### Dependent model $M_3b$\n",
    "The derivation is similar and corresponds to the factorization:\n",
    "\n",
    "\\begin{equation}\n",
    "p(x_1, x_2) = p(x_2) p(x_1|x_2) \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M1: -25.8064501902\n",
      "M2: -24.7950874021\n",
      "M3: -24.7950874021\n",
      "M3b: -24.7950874021\n",
      "Log Odds, M1-M2\n",
      "-1.01136278809\n",
      "[ 0.26671324  0.73328676]\n",
      "Log Odds, M1-M3\n",
      "-1.01136278809\n",
      "[ 0.26671324  0.73328676]\n",
      "Log Odds, M1-M3b\n",
      "-1.01136278809\n",
      "[ 0.26671324  0.73328676]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from notes_utilities import randgen, log_sum_exp, normalize_exp, normalize\n",
    "\n",
    "import scipy as sc\n",
    "from scipy.special import gammaln\n",
    "\n",
    "#C = np.array([[3,1,9],[7,9,17]])\n",
    "C = 1*np.array([[4,1,1],[1,1,6]])\n",
    "\n",
    "#C = np.array([[0,1,1],[1,0,2]])\n",
    "\n",
    "C_i = np.sum(C, axis=1)\n",
    "C_j = np.sum(C, axis=0)\n",
    "\n",
    "N = np.sum(C)\n",
    "S_1 = C.shape[0]\n",
    "S_2 = C.shape[1]\n",
    "\n",
    "#M1 Parameter\n",
    "M1 = {'a_1': S_2*np.ones(S_1), 'a_2': S_1*np.ones(S_2), 'A_1': None, 'A_2': None}\n",
    "M1['A_1'] = np.sum(M1['a_1'])\n",
    "M1['A_2'] = np.sum(M1['a_2'])\n",
    "#p(x_1) p(x_2)\n",
    "log_marglik_M1 = gammaln(M1['A_1']) - np.sum(gammaln(M1['a_1'])) - gammaln(M1['A_1'] + N) + np.sum(gammaln(M1['a_1'] + C_i)) \\\n",
    "            + gammaln(M1['A_2']) - np.sum(gammaln(M1['a_2'])) - gammaln(M1['A_2'] + N) + np.sum(gammaln(M1['a_2'] + C_j))\n",
    "\n",
    "# p(x_1, x_2)\n",
    "M2 = {'a_12': np.ones((S_1,S_2)), 'A_12':None}\n",
    "M2['A_12'] = np.sum(M2['a_12'])\n",
    "log_marglik_M2 = gammaln(M2['A_12']) - np.sum(gammaln(M2['a_12'])) - gammaln(M2['A_12'] + N) + np.sum(gammaln(M2['a_12'] + C)) \n",
    "\n",
    "    \n",
    "    \n",
    "M3 = {'a_1': S_2*np.ones(S_1), 'a_2': np.ones(S_2), 'A_1': None, 'A_2': None}\n",
    "M3['A_1'] = np.sum(M3['a_1'])\n",
    "M3['A_2'] = np.sum(M3['a_2'])\n",
    "\n",
    "#p(x_1) p(x_2|x_1)\n",
    "log_marglik_M3 = gammaln(M3['A_1']) - np.sum(gammaln(M3['a_1'])) - gammaln(M3['A_1'] + N) + np.sum(gammaln(M3['a_1'] + C_i)) \n",
    "for i in range(S_1):\n",
    "    log_marglik_M3 +=  gammaln(M3['A_2']) - np.sum(gammaln(M3['a_2'])) - gammaln(M3['A_2'] + C_i[i]) + np.sum(gammaln(M3['a_2'] + C[i,:]))\n",
    "\n",
    "# Beware the prior parameters\n",
    "M3b = {'a_1': np.ones(S_1), 'a_2': S_1*np.ones(S_2), 'A_1': None, 'A_2': None}\n",
    "M3b['A_1'] = np.sum(M3b['a_1'])\n",
    "M3b['A_2'] = np.sum(M3b['a_2'])\n",
    "\n",
    "#p(x_2) p(x_1|x_2)\n",
    "log_marglik_M3b = gammaln(M3b['A_2']) - np.sum(gammaln(M3b['a_2'])) - gammaln(M3b['A_2'] + N) + np.sum(gammaln(M3b['a_2'] + C_j)) \n",
    "for j in range(S_2):\n",
    "    log_marglik_M3b +=  gammaln(M3b['A_1']) - np.sum(gammaln(M3b['a_1'])) - gammaln(M3b['A_1'] + C_j[j]) + np.sum(gammaln(M3b['a_1'] + C[:,j]))\n",
    "\n",
    "\n",
    "print('M1:', log_marglik_M1)\n",
    "print('M2:', log_marglik_M2)\n",
    "print('M3:', log_marglik_M3)\n",
    "print('M3b:', log_marglik_M3b)\n",
    "\n",
    "print('Log Odds, M1-M2')\n",
    "print(log_marglik_M1 - log_marglik_M2)\n",
    "print(normalize_exp([log_marglik_M1, log_marglik_M2]))\n",
    "\n",
    "print('Log Odds, M1-M3')\n",
    "print(log_marglik_M1 - log_marglik_M3)\n",
    "print(normalize_exp([log_marglik_M1, log_marglik_M3]))\n",
    "\n",
    "print('Log Odds, M1-M3b')\n",
    "print(log_marglik_M1 - log_marglik_M3b)\n",
    "print(normalize_exp([log_marglik_M1, log_marglik_M3b]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conceptually $M_2$, $M_3$ and $M_3b$ should have the same marginal likelihood score, as the dependence should not depend on how we parametrize the conditional probability tables. However, this is dependent on the choice of the prior parameters. \n",
    "\n",
    "How should the prior parameters of $M_2$, $M_3$ and $M_3b$ be chosen such that we get the same evidence score?\n",
    "\n",
    "The models \n",
    "$M_2$, $M_3$ and $M_3b$ are all equivlent, if the prior parameters are chosen appropriately. For $M_2$ and $M_3$, we need to take $a_1(i) = \\sum_j a_{1,2}(i,j)$.\n",
    "\n",
    "For example, if in $M_2$, the prior parameters $a_{1,2}$ are chosen as\n",
    "\\begin{eqnarray}\n",
    "a_{1,2} & = & \\left(\\begin{array}{ccc} 1 & 1 & 1\\\\ 1 & 1 & 1 \\end{array} \\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "we need to choose in model $M_3$\n",
    "\\begin{eqnarray}\n",
    "a_{1} & = & \\left(\\begin{array}{c} 3 \\\\ 3  \\end{array} \\right)\n",
    "\\end{eqnarray}\n",
    "\\begin{eqnarray}\n",
    "a_{2} & = & \\left(\\begin{array}{ccc} 1 & 1 & 1\\\\ 1 & 1 & 1 \\end{array} \\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "and in model $M_3b$\n",
    "\\begin{eqnarray}\n",
    "a_{2} & = & \\left(\\begin{array}{ccc} 2 & 2 & 2  \\end{array} \\right)\n",
    "\\end{eqnarray}\n",
    "\\begin{eqnarray}\n",
    "a_{1} & = & \\left(\\begin{array}{ccc} 1 & 1 & 1\\\\ 1 & 1 & 1 \\end{array} \\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "This is due to fact that the marginals of a Dirichlet distribution are also Dirichlet. In particular,\n",
    "if a probability vector $x$ and corresponding parameter vector $a$ are partitioned as $x = (x_\\iota, x_{-\\iota})$ \n",
    "and $a = (a_\\iota, a_{-\\iota})$, the Dirichlet distribution\n",
    "$$\n",
    "\\mathcal{D}(x_\\iota, x_{-\\iota}; a_\\iota, a_{-\\iota})\n",
    "$$\n",
    "has marginals \n",
    "$$\n",
    "\\mathcal{D}(X_{\\iota}, X_{-\\iota}; A_{\\iota}, A_{-\\iota})\n",
    "$$\n",
    "where $X_\\iota = \\sum_{i \\in \\iota}  x_i$ and $X_{-\\iota} = \\sum_{i \\in -\\iota}  x_i$, where $A_\\iota = \\sum_{i \\in \\iota}  a_i$ and $A_{-\\iota} = \\sum_{i \\in -\\iota}  a_i$. The script below verifies that the marginals are indeed distributed according to this formula.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAFkCAYAAAC0KZhSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAGWBJREFUeJzt3X2QZFd53/Hvswhrs7I1qmSwwIEtIQsmTeySmJGFFUXC\nGIKCAzJy+YUWE8CUylFklfE4FSMCKmOrHPNSYhABB2I7SLCmMVQSEGZLwsayIiQQQROBbZo1tpEa\nvSDnIjSKtIze9uSP7o3njGd35/bc0z3T8/1UdUlz5p4+T5/t6f71ubfvjZQSkiRJh+0adwGSJGlr\nMRxIkqSM4UCSJGUMB5IkKWM4kCRJGcOBJEnKGA4kSVLGcCBJkjKGA0mSlDEcSJKkTO1wEBHnRsR1\nEXFPRByKiAvW2aYVEZ+MiAcj4uGIuC0intlMyZIkqaRhVg5OAO4ALgX+3oUZIuIHgZuBrwLnAT8M\nXAmsDF+mJEkaldjMhZci4hDwypTSdavaOsBjKaXXNlCfJEkasUaPOYiIAP4V8PWIuD4i7o+IL0TE\nTzY5jiRJKue4hu/v+4HvBd4IvBn4VeBlwH+PiB9LKd28tkNE/CPgfOBO3PUgSVIdu4FTgBtSSt9u\n6k6bDgeHVyI+kVJ6z+D/vxIR/wy4hP6xCGudD/x+w3VIkrSTvBr4SFN31nQ4qIAngO6a9i5wzhH6\n3Amwb98+Wq1Ww+XoSBYWFlhcXBx3GTuKcz56zvnoOeej1e12mZ+fh8F7aVMaDQcppccj4n8BM2t+\n9VzgriN0WwFotVrMzs42WY6OYmpqyvkeMed89Jzz0XPOx6bR3fK1w0FEnACcBsSg6dSIOB14IKX0\nTeCdwEcj4mbgRvrHHLwceGEzJUuSpJKGWTk4k/6bfhrcrhq0Xwu8PqX0iYi4BPgPwNXAAeCnUkqf\nb6BeSZJUWO1wkFK6iWN8BTKldA1wzXAlSZKkcfLaCjtUu90edwk7jnM+es756Dnnk2FTZ0hspICI\nWeD222+/3YNYJEmqYWlpibm5OYC5lNJSU/fryoEkScoYDiRJUsZwIEmSMoYDSZKUMRxIkqSM4UCS\nJGUMB5IkKWM4kCRJGcOBJEnKGA4kSVLGcCBJkjLDXLJZ0oTo9XpUVVW73/T0NHv37i1QkaStwHAg\n7VC9Xo+ZmRYrKwdr9929ew8HDnQNCNKEMhxIO1RVVYNgsA9o1ejZZWVlnqqqDAfShDIcSDteC/By\n6ZL+jgckSpKkjOFAkiRlDAeSJCljOJAkSRnDgSRJyhgOJElSxnAgSZIyhgNJkpQxHEiSpIzhQJIk\nZQwHkiQpUzscRMS5EXFdRNwTEYci4oKjbPv+wTa/tLkyJUnSqAyzcnACcAdwKZCOtFFEXAi8ALhn\nuNIkSdI41L4qY0rpeuB6gIiI9baJiH8MXA2cD+zfTIGSJGm0Gj/mYBAYPgS8I6XUbfr+JUlSWSUO\nSLwceCyl9N4C9y1JkgqrvVvhaCJiDvgl4Pl1+y4sLDA1NZW1tdtt2u12Q9VJkrR9dTodOp1O1ra8\nvFxkrEbDAfDPgacB31x1OMJTgHdFxC+nlE49UsfFxUVmZ2cbLkeSpMmw3gfmpaUl5ubmGh+r6XDw\nIeCP1rR9ZtD+wYbHkiRJBdQOBxFxAnAacHhp4NSIOB14IKX0TeA7a7Z/HPhWSunrmy1WkiSVN8zK\nwZnAjfTPcZCAqwbt1wKvX2f7I54LQZIkbT3DnOfgJmp8y+FoxxlIkqStx2srSJKkjOFAkiRlDAeS\nJCljOJAkSRnDgSRJyhgOJElSxnAgSZIyhgNJkpQxHEiSpIzhQJIkZQwHkiQpYziQJEkZw4EkScoY\nDiRJUsZwIEmSMoYDSZKUMRxIkqSM4UCSJGUMB5IkKWM4kCRJGcOBJEnKGA4kSVLGcCBJkjKGA0mS\nlDEcSJKkjOFAkiRlDAeSJClTOxxExLkRcV1E3BMRhyLiglW/Oy4i3h4RX4mIhwfbXBsRz2i2bEmS\nVMowKwcnAHcAlwJpze/2AGcAvw48H7gQmAE+uYkaJUnSCB1Xt0NK6XrgeoCIiDW/ewg4f3VbRFwG\n3BYRz0wp3b2JWiVJ0giM4piDk+ivMDw4grEkSdImFQ0HEXE88DbgIymlh0uOJUmSmlEsHETEccDH\n6a8aXFpqHEmS1KzaxxxsxKpg8CzgxzeyarCwsMDU1FTW1m63abfbJUqUJGlb6XQ6dDqdrG15ebnI\nWI2Hg1XB4FTgRSml72yk3+LiIrOzs02XI0nSRFjvA/PS0hJzc3ONj1U7HETECcBpwOFvKpwaEacD\nDwD3Af+N/tcZXw48NSJOHmz3QErp8c2XLEmSShpm5eBM4Eb6xxIk4KpB+7X0z2/wikH7HYP2GPz8\nIuB/bqZYSZJU3jDnObiJox/I6CmZJUnaxnwjlyRJGcOBJEnKGA4kSVLGcCBJkjKGA0mSlDEcSJKk\njOFAkiRlDAeSJCljOJAkSRnDgSRJyhgOJElSxnAgSZIyhgNJkpQxHEiSpIzhQJIkZQwHkiQpYziQ\nJEkZw4EkScoYDiRJUsZwIEmSMoYDSZKUMRxIkqSM4UCSJGUMB5IkKWM4kCRJGcOBJEnKGA4kSVLG\ncCBJkjK1w0FEnBsR10XEPRFxKCIuWGeb34iIeyPiYET8UUSc1ky5kiSptGFWDk4A7gAuBdLaX0bE\nG4HLgF8AzgIeAW6IiO/ZRJ2SJGlEjqvbIaV0PXA9QETEOpu8AbgypfSHg21eA9wPvBL42PClSpKk\nUWj0mIOIeDbwdOCzh9tSSg8BtwFnNzmWJEkqo+kDEp9Of1fD/Wva7x/8TpIkbXG1dytI0jB6vR5V\nVdXqMz09zd69ewtVJOlImg4H3wICOJl89eBk4H8frePCwgJTU1NZW7vdpt1uN1yipFHr9XrMzLRY\nWTlYq9/u3Xs4cKBrQJCATqdDp9PJ2paXl4uM1Wg4SCl9IyK+BbwY+ApARJwIvAB439H6Li4uMjs7\n22Q5kraIqqoGwWAf0Npgry4rK/NUVWU4kFj/A/PS0hJzc3ONj1U7HETECcBp9FcIAE6NiNOBB1JK\n3wTeDbwlIv4KuBO4Ergb+GQjFUvaxlqAHwKkrW6YlYMzgRvpH3iYgKsG7dcCr08pvSMi9gAfAE4C\nbgZellJ6rIF6JUlSYcOc5+AmjvEth5TSW4G3DleSJEkaJ6+tIEmSMoYDSZKUMRxIkqSM4UCSJGUM\nB5IkKWM4kCRJGcOBJEnKGA4kSVLGcCBJkjKGA0mSlDEcSJKkjOFAkiRlDAeSJCkzzCWbJWni9Ho9\nqqqq1Wd6epq9e/cWqkgaH8OBpB2v1+sxM9NiZeVgrX67d+/hwIGuAUETx3AgacerqmoQDPYBrQ32\n6rKyMk9VVYYDTRzDgST9fy1gdtxFSGPnAYmSJCljOJAkSRnDgSRJyhgOJElSxnAgSZIyhgNJkpQx\nHEiSpIzhQJIkZQwHkiQpYziQJEkZw4EkSco0Hg4iYldEXBkRfxMRByPiryLiLU2PI0mSyihx4aXL\ngX8DvAb4KnAmcE1EPJhSem+B8SRJUoNKhIOzgU+mlK4f/NyLiIuAswqMJUmSGlbimINbgRdHxHMA\nIuJ04Bxgf4GxJElSw0qsHLwNOBH4WkQ8ST+AvDml9NECY0kTp9frUVVV7X7T09Ps3bu3QEXr63a7\nRbaVNH4lwsHPARcBr6J/zMEZwNURcW9K6cNH6rSwsMDU1FTW1m63abfbBUqUtqZer8fMTIuVlYO1\n++7evYcDB7ojCAj3AbuYn58vPI6k1TqdDp1OJ2tbXl4uMlaJcPAO4LdSSh8f/PwXEXEK8CbgiOFg\ncXGR2dnZAuVI20dVVYNgsA9o1ejZZWVlnqqqRhAOHgQOUa/G/cAVxSqSdoL1PjAvLS0xNzfX+Fgl\nwsEe4Mk1bYfwnApSDS1gq4flOjW6W0HaTkqEg08Bb4mIu4G/oP/qsQD8boGxJElSw0qEg8uAK4H3\nAd8P3Av850GbJEna4hoPBymlR4BfGdwkSdI243EAkiQpYziQJEkZw4EkScoYDiRJUsZwIEmSMoYD\nSZKUMRxIkqSM4UCSJGUMB5IkKWM4kCRJGcOBJEnKGA4kSVLGcCBJkjIlLtksaaDX61FV1Ya373a7\nBavZOZx3aXMMB1IhvV6PmZkWKysHx13KjuK8S5tnOJAKqapq8Aa1D2htsNd+4IpyRe0Azru0eYYD\nqbgWMLvBbV3ebo7zLg3LAxIlSVLGcCBJkjKGA0mSlDEcSJKkjOFAkiRlDAeSJCljOJAkSRnDgSRJ\nyhgOJElSxnAgSZIyhgNJkpQpEg4i4gci4sMRUUXEwYj4ckRs9CTnkiRpjBq/8FJEnATcAnwWOB+o\ngOcA32l6LEmS1LwSV2W8HOillC5e1XZXgXEkSVIBJXYrvAL4UkR8LCLuj4iliLj4mL0kSdKWUGLl\n4FTg3wJXAb8JnAW8JyIeTSl9uMB40kj0ej2qqtrw9t1ut2A1klROiXCwC/hiSumKwc9fjogfAi4B\njhgOFhYWmJqaytra7TbtdrtAiVI9vV6PmZkWKysHx12KpB2q0+nQ6XSytuXl5SJjlQgH9wFrPzJ1\ngZ86WqfFxUVmZ/1Cg7amqqoGwWAf0Npgr/3AFcfcSpI2Yr0PzEtLS8zNzTU+VolwcAsws6ZtBg9K\n1ERoARsNse5WkLQ9lTggcRH40Yh4U0T8YERcBFwMvLfAWJIkqWGNh4OU0peAC4E28GfAm4E3pJQ+\n2vRYkiSpeSV2K5BS2k9/h6skSdpmvLaCJEnKGA4kSVLGcCBJkjKGA0mSlDEcSJKkjOFAkiRlDAeS\nJCljOJAkSRnDgSRJyhgOJElSxnAgSZIyhgNJkpQxHEiSpIzhQJIkZQwHkiQpYziQJEkZw4EkScoY\nDiRJUsZwIEmSMoYDSZKUMRxIkqSM4UCSJGUMB5IkKWM4kCRJGcOBJEnKGA4kSVLGcCBJkjKGA0mS\nlCkeDiLi8og4FBHvKj2WJEnavKLhICJ+BPgF4Mslx5EkSc0pFg4i4nuBfcDFwIOlxpEkSc0quXLw\nPuBTKaU/KTiGJElq2HEl7jQiXgWcAZxZ4v6lzej1elRVVatPt9stVI2Ope7cj/rfapjxpqen2bt3\nb60+wzxvhx1LajwcRMQzgXcDL0kpPb7RfgsLC0xNTWVt7XabdrvdcIXayXq9HjMzLVZWDo67FB3T\nfcAu5ufnx13IEQxf3+7dezhwoLvhN+3NPG/rjqWtq9Pp0Ol0srbl5eUiY5VYOZgDngYsRUQM2p4C\nnBcRlwHHp5TS2k6Li4vMzs4WKEf6O1VVDV5g9wGtGj33A1eUKUpH8CBwiK37bzVsfV1WVuapqmrD\nb9jDP2/rj6Wta70PzEtLS8zNzTU+Volw8MfAD69puwboAm9bLxhIo9cC6oRRdyuMz1b/t6pb33YZ\nSztZ4+EgpfQI8NXVbRHxCPDtlJKvsJIkbXGjOkOiqwWSJG0TRb6tsFZK6cdHMY4kSdo8r60gSZIy\nhgNJkpQxHEiSpIzhQJIkZQwHkiQpYziQJEkZw4EkScoYDiRJUsZwIEmSMoYDSZKUMRxIkqSM4UCS\nJGVGcuElqZRer0dVVRvevtud7KuG13l8kz4XkoZnONC21ev1mJlpsbJycNylbAH3AbuYn58fdyGS\nJoDhQNtWVVWDYLAPaG2w137ginJFjc2DwCGcC0lNMBxoArSA2Q1uO+lL6c6FpM3zgERJkpQxHEiS\npIzhQJIkZQwHkiQpYziQJEkZw4EkScoYDiRJUsZwIEmSMoYDSZKUMRxIkqSM4UCSJGUaDwcR8aaI\n+GJEPBQR90fE/4iI5zY9jiRJKqPEysG5wH8CXgC8BHgq8JmI+AcFxpIkSQ1r/KqMKaWfWP1zRLwO\n+FtgDvhc0+NJkqRmjeKYg5OABDwwgrEkSdImFQ0HERHAu4HPpZS+WnIsSZLUjMZ3K6zx28DzgHMK\nj6MCXv7yC/n0pz9Ru9/rXncxH/zg79Tu1+v1qKpqw9t3u93aY0iSjq1YOIiI9wI/AZybUrrvWNsv\nLCwwNTWVtbXbbdrtdqEKdSyf+9zNwIX0/xk36uPcdNPNtcfq9XrMzLRYWTlYu68k7QSdTodOp5O1\nLS8vFxmrSDgYBIOfBF6YUuptpM/i4iKzs7MlytGmnAVcXGP7A8BdtUepqmoQDPYBrQ322g9cUXss\nSdqO1vvAvLS0xNzcXONjNR4OIuK3gTZwAfBIRJw8+NVySmml6fE0aVrARkOiuxUkqYQSByReApwI\n/Clw76rbzxYYS5IkNazEeQ48JbMkSduYb+SSJCljOJAkSRnDgSRJyhgOJElSxnAgSZIyhgNJkpQx\nHEiSpIzhQJIkZQwHkiQpYziQJEkZw4EkScoYDiRJUsZwIEmSMo1flVF64onHWVpaqtWn2+0Wqkba\nmuo85zf791G3/6OPPsrxxx9fq8/09DR79+6t1Udbl+FADXuIu+++i7m5uXEXIm1R9wG7mJ+f38Jj\nPQV4slaP3bv3cOBA14AwIQwHath3SelJYB/QqtFvP3BFmZKkLeVB4BD1/kaG/fvYzFh1+nRZWZmn\nqirDwYQwHKiQFjBbY3t3K2inqfM3stm/j2HGqvs3rEniAYmSJCljOJAkSRnDgSRJyhgOJElSxnAg\nSZIyhgNJkpQxHEiSpIzhQJIkZQwHkiQpYziQJEkZw8EO1el0xl3CDuScj55zPmq+tkyGYuEgIn4x\nIr4REd+NiC9ExI+UGkv1+Qc8Ds756Dnno+Zry2QoEg4i4ueAq4BfA54PfBm4ISKmS4wnSZKaU2rl\nYAH4QErpQymlrwGXAAeB1xcaT5IkNaTxcBARTwXmgM8ebkspJeCPgbObHk+SJDXruAL3OQ08Bbh/\nTfv9wMw62+8G6HY3e71y1bG8vMzS0tJRt3nyySeALwL/pcY9/+Xgv/updw36W4boN6o+TY11N/D7\nW7i+rTZWE/VtZM6HHWurz98ox/oG0H8d38hri5qz6r1zd5P3G/0P9Q3eYcQzgHuAs1NKt61qfztw\nXkrp7DXbX8TG/nolSdL6Xp1S+khTd1Zi5aACngROXtN+MvCtdba/AXg1cCewUqAeSZIm1W7gFPrv\npY1pfOUAICK+ANyWUnrD4OcAesB7UkrvbHxASZLUmBIrBwDvAq6JiNvp77ReAPYA1xQaT5IkNaRI\nOEgpfWxwToPfoL874Q7g/JTS/ykxniRJak6R3QqSJGn78toKkiQpYziQJEmZkYSDuhdhiogfi4jb\nI2IlIv4yIl47ijonSZ05j4gLI+IzEfG3EbEcEbdGxEtHWe8kGPZiYxFxTkQ8HhGeOaamIV5bvici\nfjMi7hy8vvxNRLxuROVOhCHm/NURcUdEPBIR90bE70XEPxxVvdtdRJwbEddFxD0RcSgiLthAn02/\nhxYPB3UvwhQRpwB/SP/0y6cDVwO/GxH/onStk2KIC1+dB3wGeBkwC9wIfCoiTh9BuRNh2IuNRcQU\ncC3904urhiHn/OPAi4CfB54LtIEDhUudGEO8np9D//n9O8DzgJ8GzqLeaVd3uhPoH9R/KXDMgwQb\new9NKRW9AV8Arl71c9A/p+mvHmH7twNfWdPWAfaXrnVSbnXn/Aj38efAW8b9WLbLbdg5Hzy3f53+\ni+3SuB/HdroN8dryL4EHgJPGXft2vQ0x5/8O+PqatsuA3rgfy3a8AYeAC46xTSPvoUVXDoa8CNOP\n8vc/Rd1wlO21ShMXvhqctOr76L+Q6hiGnfOI+Hng2fTDgWoYcs5fAXwJeGNE3B0RByLinRHR6Dnp\nJ9WQc/554FkR8bLBfZwM/Azw6bLV7miNvIeW3q1wtIswPf0IfZ5+hO1PjIjjmy1vIg0z52v9e/pL\nWR9rsK5JVnvOI+I5wH+kfz70Q2XLm0jDPM9PBc4F/inwSuAN9Je531eoxklTe85TSrcC88AfRMRj\nwH3Ad+ivHqiMRt5D/baCMoMLYV0B/ExKqRp3PZMoInbRv9jYr6WU/vpw8xhL2il20V+WvSil9KWU\n0vXArwCv9YNHGRHxPPr7vN9K/3im8+mvln1gjGVpA0qdPvmwuhdhYtC+3vYPpZQebba8iTTMnAMQ\nEa+if6DQT6eUbixT3kSqO+ffB5wJnBERhz+17qK/R+cx4KUppT8tVOukGOZ5fh9wT0rp4VVtXfrB\n7JnAX6/bS4cNM+eXA7eklN41+PnPI+JS4OaIeHNKae0nXG1eI++hRVcOUkqPA7cDLz7cNtif/WLg\n1iN0+/zq7QdeOmjXMQw550REG/g94FWDT1TaoCHm/CHgh4Az6B9NfDrwfuBrg/+/bZ0+WmXI5/kt\nwA9ExJ5VbTP0VxPuLlTqxBhyzvcAT6xpO0T/qHtXy8po5j10BEdX/ixwEHgN8E/oLyd9G3ja4Pe/\nBVy7avtTgP9L/4jLGfpf33gMeMm4jxTdLrch5vyiwRxfQj9hHr6dOO7Hsl1uded8nf5+W6HwnNM/\njuYu4A+AFv2v8B4A3j/ux7JdbkPM+WuBRwevLc8GzqF/Mb5bx/1Ytstt8Lw9nf6HiUPALw9+ftYR\n5ryR99BRPbhLgTuB79JPL2eu+t0HgT9Zs/159BPqd4GvA/963P9A2+1WZ87pn9fgyXVu/3Xcj2M7\n3eo+z9f0NRyMYM7pn9vgBuDhQVB4B3D8uB/HdroNMee/CPzZYM7vpn/eg2eM+3FslxvwwkEoWPf1\nudR7qBdekiRJGb+tIEmSMoYDSZKUMRxIkqSM4UCSJGUMB5IkKWM4kCRJGcOBJEnKGA4kSVLGcCBJ\nkjKGA0mSlDEcSJKkzP8Dz+MCsOjh3s8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110d996d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAFkCAYAAAC0KZhSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAF5FJREFUeJzt3X2QnVd9H/Dvz4hYyA1iWlFDChrwAOqSdAwrB8pQkxII\nLkxxQycvLGwgMJRShwlROm1gwBMIQwNkQJBCAi1pMChsgOl0gNRjkxCgxA4wsDV5YVFJgiNsjNPF\nWBSbNcY6/eOuGh1V8r7oeXa1q89n5s5oz577nN89evbe7z33uc9TrbUAABx33mYXAACcXYQDAKAj\nHAAAHeEAAOgIBwBARzgAADrCAQDQEQ4AgI5wAAB0hAMAoLPmcFBVl1bVh6vq5qo6VlWX30vfdyz3\n+fkzKxMA2CjrWTm4IMkNSa5IctoLM1TVs5I8PsnN6ysNANgMO9Z6h9baNUmuSZKqqlP1qap/kOSt\nSS5LcvWZFAgAbKzBjzlYDgzvSfLG1trC0NsHAMa15pWDVXh5ku+21t62ms5V9fcyWWG4McnSCPUA\nwHa1M8nDklzbWvvGUBsdNBxU1f4kP5/ksWu422VJfmfIOgDgHPPcJO8bamNDrxz8kyQPTPLVEw5H\nuE+SN1fVL7TWLjrFfW5MkkOHDmVqamrgcjidAwcO5ODBg5tdxjllO835wsJCZmdnk7w2ycNX6P2V\nJFduyt/4dprzrcKcb6y//VucvJYOZehw8J4kv39S20eX23/7NPdZSpKpqalMT08PXA6ns3v3bvO9\nwbbnnD8jyUqPaT7JlZvyN7495/zsZs43zaAfy685HFTVBUkekeT40sBFVXVxkttaa19N8s2T+t+d\n5OuttS+fabEAwPjWs3JwSZKPZ3KOg5bkTcvtVyV54Sn6n/ZcCADA2Wc95zn4ZNbwFcjTHGcAAJyl\nXFvhHDUzM7PZJZxzzPnGM+cbz5xvD8LBOcof8MYz5xvPnG88c749CAcAQEc4AAA6wgEA0BEOAICO\ncAAAdIQDAKAjHAAAHeEAAOgIBwBARzgAADrCAQDQEQ4AgI5wAAB0hAMAoCMcAAAd4QAA6AgHAEBH\nOAAAOsIBANARDgCAjnAAAHSEAwCgIxwAAB3hAADoCAcAQEc4AAA6wgEA0BEOAICOcAAAdIQDAKAj\nHAAAnTWHg6q6tKo+XFU3V9Wxqrr8hN/tqKo3VNWfVNW3l/tcVVUPHrZsAGAs61k5uCDJDUmuSNJO\n+t2uJI9J8pokj03yrCT7knzoDGoEADbQjrXeobV2TZJrkqSq6qTffSvJZSe2VdVLk3ymqh7SWrvp\nDGoFADbARhxz8IBMVhhu34CxAIAzNGo4qKrzk7w+yftaa98ecywAYBijhYOq2pHkg5msGlwx1jgA\nwLDWfMzBapwQDB6a5EdXs2pw4MCB7N69u2ubmZnJzMzMGCUCwJYyNzeXubm5ru3o0aOjjDV4ODgh\nGFyU5MmttW+u5n4HDx7M9PT00OUAwLZwqjfM8/Pz2b9//+BjrTkcVNUFSR6R5Pg3FS6qqouT3Jbk\nliT/NZOvM/7zJPetqguX+93WWrv7zEsGAMa0npWDS5J8PJNjCVqSNy23X5XJ+Q2eudx+w3J7Lf/8\n5CT/40yKBQDGt57zHHwy934go1MyA8AW5oUcAOgIBwBARzgAADrCAQDQEQ4AgI5wAAB0hAMAoCMc\nAAAd4QAA6AgHAEBHOAAAOsIBANARDgCAjnAAAHSEAwCgIxwAAB3hAADoCAcAQEc4AAA6wgEA0BEO\nAICOcAAAdIQDAKAjHAAAHeEAAOgIBwBARzgAADrCAQDQEQ4AgI5wAAB0hAMAoCMcAAAd4QAA6AgH\nAEBnzeGgqi6tqg9X1c1VdayqLj9Fn1+pqq9V1Z1V9ftV9YhhygUAxraelYMLktyQ5Iok7eRfVtUv\nJXlpkhcneVySO5JcW1XfdwZ1AgAbZMda79BauybJNUlSVXWKLi9L8trW2u8t93lekluT/HiSD6y/\nVABgIwx6zEFVPTzJg5J87Hhba+1bST6T5AlDjgUAjGPoAxIflMlHDbee1H7r8u8AgLPcmj9WGMuB\nAweye/furm1mZiYzMzObVBEAnD3m5uYyNzfXtR09enSUsYYOB19PUkkuTL96cGGS/3lvdzx48GCm\np6cHLgcAtodTvWGen5/P/v37Bx9r0I8VWmtfySQgPOV4W1XdP8njk1w/5FgAwDjWvHJQVRckeUQm\nKwRJclFVXZzkttbaV5O8JcmrquovktyY5LVJbkryoUEqBgBGtZ6PFS5J8vFMDjxsSd603H5Vkhe2\n1t5YVbuSvDPJA5J8KsnTW2vfHaBeAGBk6znPwSezwscRrbVXJ3n1+koCADaTaysAAB3hAADoCAcA\nQEc4AAA6wgEA0BEOAICOcAAAdIQDAKAjHAAAHeEAAOgIBwBARzgAADrCAQDQEQ4AgI5wAAB0dmx2\nAbCdHTlyJIuLi6vqu2fPnuzdu3db1gBsLcIBjOTIkSPZt28qS0t3rqr/zp27cvjwwqAvzmdDDcDW\nIxzASBYXF5dflA8lmVqh90KWlmazuLg46Avz2VADsPUIBzC6qSTTagC2DAckAgAd4QAA6AgHAEBH\nOAAAOsIBANARDgCAjnAAAHSEAwCgIxwAAB3hAADoCAcAQEc4AAA6wgEA0BEOAIDO4OGgqs6rqtdW\n1V9V1Z1V9RdV9aqhxwEAxrFjhG2+PMm/TvK8JF9MckmSd1fV7a21t40wHgAwoDHCwROSfKi1ds3y\nz0eq6jlJHjfCWADAwMY45uD6JE+pqkcmSVVdnOSJSa4eYSwAYGBjrBy8Psn9k3ypqu7JJIC8srX2\nuyOMBWwRCwsLq+q3Z8+e7N27d+RqgHszRjj46STPSfLsTI45eEySt1bV11pr7z3dnQ4cOJDdu3d3\nbTMzM5mZmRmhRGDj3JLkvMzOzq6q986du3L48IKAACeZm5vL3Nxc13b06NFRxhojHLwxya+21j64\n/POfV9XDkrwiyWnDwcGDBzM9PT1COcDmuj3JsSSHkkyt0HchS0uzWVxcFA7gJKd6wzw/P5/9+/cP\nPtYY4WBXkntOajsW51SAc9xUEm8AYCsYIxx8JMmrquqmJH+eybPBgSTvGmEsAGBgY4SDlyZ5bZK3\nJ/n7Sb6W5DeX2wCAs9zg4aC1dkeSX1y+AQBbjOMAAICOcAAAdIQDAKAjHAAAHeEAAOgIBwBARzgA\nADrCAQDQEQ4AgI5wAAB0hAMAoCMcAAAd4QAA6AgHAEBHOAAAOjs2uwDgby0sLKyq3549e7J3796R\nqwHOVcIBnBVuSXJeZmdnV9V7585dOXx4QUAARiEcwFnh9iTHkhxKMrVC34UsLc1mcXFROABGIRzA\nWWUqyfRmFwGc4xyQCAB0hAMAoCMcAAAd4QAA6AgHAEBHOAAAOsIBANARDgCAjnAAAHSEAwCgIxwA\nAB3hAADoCAcAQEc4AAA6o4SDqvqBqnpvVS1W1Z1V9YWqch1aANgCdgy9wap6QJLrknwsyWVJFpM8\nMsk3hx4LABje4OEgycuTHGmtveiEtr8eYRwAYARjfKzwzCSfq6oPVNWtVTVfVS9a8V4AwFlhjHBw\nUZJ/k+Rwkqcl+c0kv15VPzPCWADAwMb4WOG8JJ9trV25/PMXquqHkrwkyXtHGA/OSQsLC4P02ajt\nDl1Dktx11105//zzV9V3z5492bt375mUBeeMMcLBLUlO/steSPIv7+1OBw4cyO7du7u2mZmZzMzM\nDFsdbHm3JDkvs7OzW2S7Y9ZwnyT3rKrnzp27cvjwgoDAljU3N5e5ubmu7ejRo6OMNUY4uC7JvpPa\n9mWFgxIPHjyY6WnfdoSV3Z7kWJJDSaZW6Ht1kitX6DP2dtdiPTWspu9ClpZms7i4KBywZZ3qDfP8\n/Hz2798/+FhjhIODSa6rqlck+UCSxyd5UZJ/NcJYcA6bSrJSoF7P8v9Y2x2rhtX0BdZi8AMSW2uf\nS/KsJDNJ/jTJK5O8rLX2u0OPBQAMb4yVg7TWrs5kzQ8A2GJcWwEA6AgHAEBHOAAAOsIBANARDgCA\njnAAAHSEAwCgIxwAAB3hAADoCAcAQEc4AAA6wgEA0BEOAICOcAAAdIQDAKAjHAAAHeEAAOgIBwBA\nRzgAADrCAQDQEQ4AgI5wAAB0hAMAoCMcAAAd4QAA6AgHAEBHOAAAOsIBANARDgCAjnAAAHSEAwCg\nIxwAAB3hAADojB4OqurlVXWsqt489lgAwJkbNRxU1Q8neXGSL4w5DgAwnNHCQVX9nSSHkrwoye1j\njQMADGvMlYO3J/lIa+0PRxwDABjYjjE2WlXPTvKYJJeMsX0AYDyDh4OqekiStyR5amvt7qG3z9mt\ntZYXv/iKfOpT162q//3utzPvf/+hPOpRjxq8liNHjmRxcXFVfffs2ZO9e/cOXgPnBvsa280YKwf7\nkzwwyXxV1XLbfZI8qapemuT81lo7+U4HDhzI7t27u7aZmZnMzMyMUCJjOXbsWN71rnckeUqSH1yx\nf9Xbc8011wweDo4cOZJ9+6aytHTnqvrv3Lkrhw8veNJmzexrbJS5ubnMzc11bUePHh1lrDHCwR8k\n+Ucntb07yUKS158qGCTJwYMHMz09PUI5bI7nJnnBir2q3jXK6IuLi8tP1oeSTK3QeyFLS7NZXFz0\nhM2a2dfYKKd6wzw/P5/9+/cPPtbg4aC1dkeSL57YVlV3JPlGa21h6PHg3k0lETrZCPY1to+NOkPi\nKVcLAICzzyjfVjhZa+1HN2IcAODMubYCANARDgCAjnAAAHSEAwCgIxwAAB3hAADoCAcAQEc4AAA6\nwgEA0BEOAICOcAAAdIQDAKAjHAAAHeEAAOgIBwBARzgAADo7NrsAOFssLCysqt9dd92V888/f7Dt\nsXHW8n+yZ8+e7N27d8Rq4OwlHEBuSXJeZmdnV9n/PknuGbEehrfW/+Nk585dOXx4QUDgnCQcQG5P\ncizJoSRTK/S9OsmVa+zL5lvL/3GSLGRpaTaLi4vCAeck4QD+n6kk0yv0Ob4svZa+nD1W8/8GOCAR\nAOgIBwBARzgAADrCAQDQEQ4AgI5wAAB0hAMAoCMcAAAd4QAA6AgHAEBHOAAAOsIBANARDgCAjnAA\nAHQGDwdV9Yqq+mxVfauqbq2q/1ZVjxp6HABgHGOsHFya5D8meXySpya5b5KPVtX9RhgLABjYjqE3\n2Fp7xok/V9XPJvmbJPuT/NHQ4wEAw9qIYw4ekKQluW0DxgIAztDgKwcnqqpK8pYkf9Ra++KYY3Fu\nOHLkSBYXF1fst7CwsAHVsN2tZj8ac19b7f6eJHv27MnevXtHq4Vzy6jhIMlvJHl0kieu1PHAgQPZ\nvXt31zYzM5OZmZmRSmOrOXLkSPbtm8rS0p2bXQrb3i1Jzsvs7OymVbDW/X3nzl05fHhBQNjG5ubm\nMjc317UdPXp0lLFGCwdV9bYkz0hyaWvtlpX6Hzx4MNPT02OVwzawuLi4/ER5KMnUCr2vTnLl+EWx\nTd2e5Fg2c19b2/6+kKWl2SwuLgoH29ip3jDPz89n//79g481SjhYDgb/IsmPtNaOjDEG57KpJCsF\nSR8rMISzYV9bTQ0wrMHDQVX9RpKZJJcnuaOqLlz+1dHW2tLQ4wEAwxrj2wovSXL/JJ9I8rUTbj81\nwlgAwMDGOM+BUzIDwBbmhRwA6AgHAEBHOAAAOsIBANARDgCAjnAAAHSEAwCgIxwAAB3hAADoCAcA\nQEc4AAA6wgEA0BEOAICOcAAAdIQDAKCzY7MLgNtuuy3z8/Mr9ltYWNiAamB8q9mX7e9sJuGATdXa\nsbzuda/Pa17zms0uBTbALUnOy+zs7GYXAvdKOGDTfe97dyU5lGRqhZ5XJ7ly/IJgNLcnORb7O2c7\n4YCzxFSS6RX6WGZlu7C/c3ZzQCIA0BEOAICOcAAAdIQDAKAjHAAAHeEAAOgIBwBARzgAADrCAQDQ\nEQ4AgI5wAAB0hAMAoCMcAAAd4eAcNTc3t9klnIPM+cYz5xvNc8v2MFo4qKqfq6qvVNV3qurTVfXD\nY43F2vkD3gzmfOOZ843muWV7GCUcVNVPJ3lTkl9O8tgkX0hybVXtGWM8AGA4Y60cHEjyztbae1pr\nX0rykiR3JnnhSOMBAAMZPBxU1X2T7E/yseNtrbWW5A+SPGHo8QCAYe0YYZt7ktwnya0ntd+aZN8p\n+u9MkoWFhRFK4XSOHj2a+fn5wbd7zz33LP/rE0nuXrF/a99b/tfVSVbaB67b4n1vSvI7m1zDudb3\n+JyvZbtj1DF2368kOTueR8d6buHUTvg/3znkdmvypn7ADVY9OMnNSZ7QWvvMCe1vSPKk1toTTur/\nnJz+GRMAWNlzW2vvG2pjY6wcLCa5J8mFJ7VfmOTrp+h/bZLnJrkxydII9QDAdrUzycMyeS0dzOAr\nB0lSVZ9O8pnW2suWf64kR5L8emvt1wYfEAAYzBgrB0ny5iTvrqrPJ/lsJt9e2JXk3SONBwAMZJRw\n0Fr7wPI5DX4lk48TbkhyWWvtf48xHgAwnFE+VgAAti7XVgAAOsIBANDZkHCw1oswVdU/rarPV9VS\nVf2vqnr+RtS5naxlzqvqWVX10ar6m6o6WlXXV9XTNrLe7WC9FxurqidW1d1V5cwxa7SO55bvq6rX\nVdWNy88vf1VVP7tB5W4L65jz51bVDVV1R1V9rap+q6r+7kbVu9VV1aVV9eGqurmqjlXV5au4zxm/\nho4eDtZ6EaaqeliS38vk9MsXJ3lrkndV1Y+NXet2sY4LXz0pyUeTPD3JdJKPJ/lIVV28AeVuC+u9\n2FhV7U5yVSanF2cN1jnnH0zy5CQvSPKoJDNJDo9c6raxjufzJ2ayf//nJI9O8hNJHpfkP21IwdvD\nBZkc1H9FkhUPEhzsNbS1NuotyaeTvPWEnyuTc5r++9P0f0OSPzmpbS7J1WPXul1ua53z02zjz5K8\narMfy1a5rXfOl/ft12TyZDu/2Y9jK93W8dzyz5LcluQBm137Vr2tY87/bZIvn9T20iRHNvuxbMVb\nkmNJLl+hzyCvoaOuHKzzIkz/OP//u6hr76U/JxjiwlfLJ636/kyeSFnBeue8ql6Q5OGZhAPWYJ1z\n/swkn0vyS1V1U1Udrqpfq6pBz0m/Xa1zzv84yUOr6unL27gwyU8m+e/jVntOG+Q1dOyPFe7tIkwP\nOs19HnSa/vevqvOHLW9bWs+cn+zfZbKU9YEB69rO1jznVfXIJP8hk/OhHxu3vG1pPfv5RUkuTfKD\nSX48ycsyWeZ++0g1bjdrnvPW2vVJZpO8v6q+m+SWJN/MZPWAcQzyGurbCnSWL4R1ZZKfbK0tbnY9\n21FVnZfJxcZ+ubX2l8ebN7Gkc8V5mSzLPqe19rnW2jVJfjHJ873xGEdVPTqTz7xfncnxTJdlslr2\nzk0si1UY6/TJx631IkxZbj9V/2+11u4atrxtaT1zniSpqmdncqDQT7TWPj5OedvSWuf8+5NckuQx\nVXX8Xet5mXyi890kT2utfWKkWreL9ezntyS5ubX27RPaFjIJZg9J8penvBfHrWfOX57kutbam5d/\n/rOquiLJp6rqla21k9/hcuYGeQ0ddeWgtXZ3ks8necrxtuXPs5+S5PrT3O2PT+y/7GnL7axgnXOe\nqppJ8ltJnr38jopVWsecfyvJDyV5TCZHE1+c5B1JvrT878+c4j6cYJ37+XVJfqCqdp3Qti+T1YSb\nRip121jnnO9K8r2T2o5lctS91bJxDPMaugFHV/5UkjuTPC/JP8xkOekbSR64/PtfTXLVCf0fluT/\nZHLE5b5Mvr7x3SRP3ewjRbfKbR1z/pzlOX5JJgnz+O3+m/1YtsptrXN+ivv7tsLIc57JcTR/neT9\nSaYy+Qrv4STv2OzHslVu65jz5ye5a/m55eFJnpjJxfiu3+zHslVuy/vtxZm8mTiW5BeWf37oaeZ8\nkNfQjXpwVyS5Mcl3Mkkvl5zwu99O8ocn9X9SJgn1O0m+nORnNvs/aKvd1jLnmZzX4J5T3P7LZj+O\nrXRb635+0n2Fgw2Y80zObXBtkm8vB4U3Jjl/sx/HVrqtY85/LsmfLs/5TZmc9+DBm/04tsotyY8s\nh4JTPj+P9RrqwksAQMe3FQCAjnAAAHSEAwCgIxwAAB3hAADoCAcAQEc4AAA6wgEA0BEOAICOcAAA\ndIQDAKDzfwHE2mcKykpKCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110d07240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "S_1 = 2\n",
    "S_2 = 10\n",
    "M = S_1*S_2\n",
    "\n",
    "a = np.ones(M)\n",
    "\n",
    "N = 100\n",
    "P = np.random.dirichlet(a, size=N)\n",
    "A = np.zeros((N,S_1))\n",
    "B = np.zeros((N*S_1,S_2))\n",
    "\n",
    "for n in range(N):\n",
    "    temp = P[n,:].reshape((S_1,S_2))\n",
    "    A[n,:] = np.sum(temp, axis=1)\n",
    "    for i in range(S_1):\n",
    "        B[(n*S_1+i),:] = temp[i,:]/A[n,i]\n",
    "\n",
    "import pylab as plt\n",
    "\n",
    "plt.hist(A[:,0],bins=20)\n",
    "plt.gca().set_xlim([0,1])\n",
    "#plt.plot(B[:,0],B[:,1],'.')\n",
    "plt.show()\n",
    "\n",
    "P2 = np.random.dirichlet(S_2*np.ones(S_1), size=N)\n",
    "plt.hist(P2[:,0],bins=20)\n",
    "plt.gca().set_xlim([0,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.980829253012\n",
      "model, consistent, inconsistent\n",
      "i  j -6.99291308732 -7.37775890823\n",
      "i->j -7.02108396429 -7.15461535691\n",
      "i<-j -7.02108396429 -7.27239839257\n",
      "i--j -7.02108396429\n",
      "[[2 1]\n",
      " [0 1]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from notes_utilities import randgen, log_sum_exp, normalize_exp, normalize\n",
    "\n",
    "import scipy as sc\n",
    "from scipy.special import gammaln\n",
    "\n",
    "# Log of multivariate Beta function\n",
    "def log_mbeta(a):\n",
    "    return sum(gammaln(a.flatten())) - gammaln(sum(a.flatten()))\n",
    "\n",
    "a = 1\n",
    "b = 1\n",
    "\n",
    "\n",
    "S = np.array([[2,1],[0,1]])\n",
    "#S = np.array([[2,1,0],[1,0,1]])\n",
    "\n",
    "S_ip = np.sum(S, axis=1)\n",
    "S_pj = np.sum(S, axis=0)\n",
    "S_pp = np.sum(S)\n",
    "\n",
    "gamma = 1\n",
    "alpha = gamma*np.ones_like(S)\n",
    "\n",
    "alpha_ip = np.sum(alpha, axis=1)\n",
    "alpha_pj = np.sum(alpha, axis=0)\n",
    "alpha_pp = np.sum(alpha)\n",
    "\n",
    "#\n",
    "log_const = a*np.log(b) -  (a+S_pp)*np.log(b+1) + gammaln(a + S_pp) - gammaln(a) - np.sum(gammaln(S+1))\n",
    "print(log_const)\n",
    "\n",
    "# inconsistent independent \n",
    "alpha_wrong_ip = gamma*np.ones_like(alpha_ip)\n",
    "alpha_wrong_pj = gamma*np.ones_like(alpha_pj)\n",
    "\n",
    "# i    j\n",
    "marglik_1_wrong =  log_mbeta(S_ip+alpha_wrong_ip) - log_mbeta(alpha_wrong_ip) \n",
    "marglik_1_wrong += log_mbeta(S_pj+alpha_wrong_pj) - log_mbeta(alpha_wrong_pj)\n",
    "# i -> j\n",
    "marglik_2_wrong  =  log_mbeta(S_ip+alpha_wrong_ip) - log_mbeta(alpha_wrong_ip) \n",
    "for i in range(S.shape[0]):\n",
    "    marglik_2_wrong +=  log_mbeta(S[i,:]+alpha[i,:]) - log_mbeta(alpha[i,:])\n",
    "# i <- j\n",
    "marglik_3_wrong  =  log_mbeta(S_pj+alpha_wrong_pj) - log_mbeta(alpha_wrong_pj) \n",
    "for j in range(S.shape[1]):\n",
    "    marglik_3_wrong +=  log_mbeta(S[:,j]+alpha[:,j]) - log_mbeta(alpha[:,j])\n",
    "\n",
    "# independent\n",
    "marglik_1  = log_mbeta(S_ip+alpha_ip) - log_mbeta(alpha_ip) \n",
    "marglik_1 += log_mbeta(S_pj+alpha_pj) - log_mbeta(alpha_pj)\n",
    "\n",
    "# i -> j\n",
    "marglik_2  =  log_mbeta(S_ip+alpha_ip) - log_mbeta(alpha_ip) \n",
    "for i in range(S.shape[0]):\n",
    "    marglik_2 +=  log_mbeta(S[i,:]+alpha[i,:]) - log_mbeta(alpha[i,:])\n",
    "# i <- j\n",
    "marglik_3  =  log_mbeta(S_pj+alpha_pj) - log_mbeta(alpha_pj) \n",
    "for j in range(S.shape[1]):\n",
    "    marglik_3 +=  log_mbeta(S[:,j]+alpha[:,j]) - log_mbeta(alpha[:,j])\n",
    "\n",
    "# dependent\n",
    "marglik_4 =  log_mbeta(S+alpha) - log_mbeta(alpha)\n",
    "\n",
    "# \n",
    "print('model, consistent, inconsistent')\n",
    "print('i  j', marglik_1+log_const, marglik_1_wrong+log_const)\n",
    "print('i->j', marglik_2+log_const, marglik_2_wrong+log_const)\n",
    "print('i<-j', marglik_3+log_const, marglik_3_wrong+log_const)\n",
    "print('i--j', marglik_4+log_const)\n",
    "\n",
    "print(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Are two given histograms drawn from the same distribution?\n",
    "\n",
    "$[3,5,12,4]$\n",
    "\n",
    "$[8, 14, 31, 14]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the Dirichlet Distribution\n",
    "\n",
    "[http://blog.bogatron.net/blog/2014/02/02/visualizing-dirichlet-distributions/]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.tri as tri\n",
    "from functools import reduce\n",
    "from scipy.special import gammaln\n",
    "\n",
    "corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])\n",
    "triangle = tri.Triangulation(corners[:, 0], corners[:, 1])\n",
    "\n",
    "refiner = tri.UniformTriRefiner(triangle)\n",
    "trimesh = refiner.refine_triangulation(subdiv=4)\n",
    "\n",
    "\n",
    "# Mid-points of triangle sides opposite of each corner\n",
    "midpoints = [(corners[(i + 1) % 3] + corners[(i + 2) % 3]) / 2.0 \\\n",
    "             for i in range(3)]\n",
    "\n",
    "def xy2bc(xy, tol=1.e-3):\n",
    "    '''Converts 2D Cartesian coordinates to barycentric.'''\n",
    "    s = [(corners[i] - midpoints[i]).dot(xy - midpoints[i]) / 0.75 \\\n",
    "         for i in range(3)]\n",
    "    return np.clip(s, tol, 1.0 - tol)\n",
    "\n",
    "class Dirichlet(object):\n",
    "    def __init__(self, alpha):\n",
    "        self._alpha = np.array(alpha)\n",
    "        self._coef = gammaln(np.sum(self._alpha)) - np.sum(gammaln(self._alpha))\n",
    "    def log_pdf(self, x):\n",
    "        return self._coef + np.sum(np.log(x)*(self._alpha - 1)) \n",
    "    def pdf(self, x):\n",
    "        '''Returns pdf value for `x`.'''\n",
    "        return np.exp(self.log_pdf(x))\n",
    "        \n",
    "def draw_pdf_contours(dist, nlevels=200, subdiv=8, **kwargs):\n",
    "\n",
    "    refiner = tri.UniformTriRefiner(triangle)\n",
    "    trimesh = refiner.refine_triangulation(subdiv=subdiv)\n",
    "    pvals = [dist.pdf(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]\n",
    "\n",
    "    plt.tricontourf(trimesh, pvals, nlevels, **kwargs)\n",
    "    plt.axis('equal')\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(0, 0.75**0.5)\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAFdCAYAAAB1gNVOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3U+obl9dx/Hvvvdc7/01CInCQRBSRhA0MGwoZBCFgyCo\nBgpCYZMGgtEgEgQLQxKchBkEIpFhBA5skKhEBEISOBAsaVJIkk5+gSA/z++e421w7n6e9exn/d9r\nrf1da71fEN7z7P2sc+616/5+1/P9rLu8ePFCAADAXB4d/QMAAID2KAAAAJgQBQAAABOiAAAAYEIU\nAAAATIgCAACACVEAAAAwIQoAAAAmRAEAIMmyfOjFsnyIE8SAzlEAAMhCEQD0jQIAQDQe+sA4KAAA\nZKMgAPpFAQAgCg97YCwUAAB2oTAA+kQBACCIhzwwHgoAAF4xD38KBKA/FAAAiqAIAPpCAQDAiYc6\nMC4KAADFUDAA/aAAAGDFwxwYGwUAgCt7Hv4UDkAfKAAAAJgQBQCACyU6eHYBAP0oAABUQREA6EYB\nAOCEhzYwDwoAANVQUAB6UQAAEBEe1sBsKAAAVH34U1gAOlEAAAAwIQoAYHItOnR2AQB9KAAANHAj\ny/JhigBAEQoAYGJtOvOb+t8CQDIKAADNsAsA6EEBAEyK7h+YGwUAMKEjh/LYBQB0oAAAUIm7+6cI\nAI5HAQBMhkgeABEKAABVhD/7ZxcAOBYFADARBv8ArCgAABzkFVmWj7ELAByEAgCYBN0/ABMFADAB\nfYN/r5x+xS4AcAwKAACF5Hf/FAFAexQAwOA0d/8AjkMBAKCA2O7f/fBnFwBoiwIAGBiDfwBcKACA\nQfW39f9EluXPlf3MwLgoAADsQPcP9IoCABhQj93/il0AoA0KAACZ6nX/FAFAfRQAwGD0Df7Fd/8A\n2qEAAFBR/sOfXQCgLgoAYCD6uv8evg8wJwoAYBA9D/5de3j4L8snlP2egHFQAABI0KIrp/MHWqAA\nAAYwVvd/iV0AoA4KAACRjuv+KQKA8igAgM7pG/yrEfvjXxAESqMAAFDQ/sE/15rL8il2AYCCKACA\njunr/gH0ggIA6NRYg3/+7n/FLgBQDgUAAI8jY3/2goIiACiDvT2gQ2N1/yn4dwOAUtgBAOAQ2x/s\nmdBP6/5Xy/JpZQUQ0B8KAKAz/Q3+0f0DGlEAANih/uCfC7sAwD4UAEBH+uv+S7EVEhp/TqAfFABA\nJ8Ya/Evp/t0P/2X5O2V/JkA/KAAAGPTF/kIoAoA8FABAB8bq/lOw9Q/UQgEA4KVSsb96g3+u77cs\nn1VWIAH6UQAAyjH4Z9L4cwJ9ogAAFBtr67/M4J/rPnYBgDQUAMD0juyq9239A8hHAQAoNUf3H7tO\nzPtfkWX5R2V/ZoBeFADA1PqL/T1wFxsUAUAcCgBAIX2Df9pjf3uKCWBOFAAAAvTF/kKW5UvsAgAB\nFACAMvq6/1ZKdf8MBwIxKAAAReYY/CsT+7OvZcYC2QUAfCgAgOmMN/gHIB0FAKDEWN1/ilpb/zey\nLP+s7M8U0IMCAJhKqfP+c77HMd0/RQBgRwEAKNDf4J/+7h+AHwUAgA3tsb/QWtuf4RVZlq+wCwBs\nUAAAB+uv+y+l1L/2x3AgkIMCADjQWIN/x8b+7Guc72cXALhEAQAMb9bY33WxsCxfpQgAXqIAAA4y\nVvef4ojBPz4mALYoAIChlYr9HT34l8r9vdkFAB5QAAAHYPDP1Lr71/jnArRHAQBM78juP2fwL7SG\n7/6He5fla+wCYHoUAEBjdP8t1+Gzf8CFAgBoaKzBP92xP//3eyLL8h/K/rsA2qIAAIZD7O96PfvP\nSxGAmVEAAI2M1f2n0DD4F3sNmAcFADAUYn/XfH8mr8iy/LeywgxogwIAaIDBP5OW7n9PwQL0jwIA\nqGysrf8jBv9Ca8QP/rksy/8o++8IqI8CABjCkd1/7a3/GrE/82fWuHMC1EcBAFQ0R/cfu06Nrf/9\n3f+KXQDMhgIA6B6xv+v1/IN/9vueyLJ8hyIA06AAACrRN/hH7I/BP+CMAgCYArE//33nn4ddAMyC\n6RegAn3dfyu9dP+u+QAOCcI82AEACptj8K+P8/5zLcuryv47BMqjAAC6NOvgX+56qd3/K7Isr1EE\nYGgUAEBBY3X/KfqL/bkxKIg5UAAA3Sl13n/O99A0+Lc/9mf/9dmyPFdW0AHlUAAAhfQ3+Ndb9x/6\nnrZrMSf+pd4PjIECABiS9thfrfP+Y8UXP+wCYFQUAEAB/XX/peR2/zHr5Ny7J/YXcz8wDgoAYKex\nBv/mjP2F1lwWUfbfMbAfBQDQhVljf6FiIeZ9qd2/HUUARkMBAOwwVvefQsPgn+vBnlqwsPWPOVEA\nAOqViv0dPfiXqt55/+GfzRULZBcA46AAADIx+GfSHvtLHfwDxkcBAGQYa+v/iMG/0BqlY38h8Vv/\n7AJgFJS8gFoj//XUFvuLtKS/BdCKHQAgEd1/6P39x/5ClkfsAqB/FACASsT+rtcred5/bIFjMLv/\nxyLLE4oA9I0CAEigb/Bvptif61qp2F+kx+lvATSiAAC6RezPf9/+2N+J7bP/V0SWH2YXAP2iAAAi\n6ev+W+ml+28Q+6P7x0AoAIAIDP6F3j/w4J+j+1+XWX6EXQD0iQIAUGPWwb/c9RrH/tbuf88fFaAI\nBQAQMFb3n6LV1n+L7n/nef+2rX9jmeXH2AVAfygAABVm7f5bxf5sdmz9m14WBxQB6A0FAOChb/Bv\ntO4/9D1t12IKlsqDfzeWa3w0gM5QAADd0B7725PH39ob+wsp1P1vri0/wS4A+kEBADjo6/5bye3+\nQ+sMEvuzdf87vgVwFAoAwGKswT9if1mxP9v2vu2PYnNt+Ul2AdAHCgDgMLMO/uWud9C/9ucZ/HNd\nW36aIgD6UQAAG2N1/ymI/Z2kbu/7PhYAlKIAAA4R2/3v6dQ1dv/jxP68f4xPRJafZRcAulEAAIb+\nBv966/5D39N2rZPBP7p/dIYCAFDryNhfzuBfaI0Ssb9c+877T7q2fqtnIsvPswsAvSgAgJf66/5L\nqbmLkHOvgu7ft/Xve1/p+UigIgoAQEYb/CP2tzv2Z1tm5ev+V8/Oay6/wC4AdKIAAJqZNfbXavDv\nuNjfxY/y7PoWigBopHE/EmhqrO4/RS+DfzEUxv74WADKsQMANFEq9nf0ef+pap73ryv2JyIXW/9X\n3+6X2AWALhQAmBqDfyaN3X8Hsb+Y35blYwHgaBQAmNZYW/9HDP6F1mgd+8sY/DPlxv5W2+7f8vBf\nfpldAOhBAQBUpbH7jxGz9a8t9hfJfPiXiv3FdPgcFARlKAAwJbr/mPf71uo49rcqGfvbrumaCXgi\nsryTXQDoQAEAVEPs73q92rG/BoN/oYd8BIoAaNDr/iSQTd/gH7G/crG/SDXP+/d0/xdrAQdjBwA4\nFLE//30HxP5iroUe8hGW32QXAMeiAMBU9HX/rfTS/Xce+1uFuv+nEWsAlVEAYBoM/oXeP8F5/6ba\nsb/QxwVPRJZ3sQuA41AAAEXNOviXu97Asb8V3T+UogDAFMbq/lOMFPtrdN7/tvvPKQxCW//G68t7\n2AXAMSgAgGJKnfef8z00Df7Viv3Z6DrvP9fyOxQBaI8CAMPrb/Cvt+4/9D1t12IG+RQN/vns6P4v\n/hNojAIAaEp77G9PHn8r58lW6nhhKTv4VyD2Z7V+r2ciy++xC4C2KAAwtP66/1Jyu//QOh3H/vYO\n/vnkdv/AgSgAMKyxBv+I/WXF/kqe9x8b+0vd+n92/np5H7sAaIcCANiF2F/aeo1jf75lcs77r8Uo\nCpbfpwhAGxQAGNJY3X8KzbG/Uuf9F479ba9Vjv1d/GftwgLwoAAAshH7u16vk/P+D479XaxrKQqW\nP2IXAPVRAGA4/Q3+9db9h76n7dpA5/3XjP1tiwKgIgoAoKojY385g3+hNUrH/lIoPO8/lmfwz/Uz\nLB9kFwB1UQBgKP11/6XU3EXIubdx91/jvP8Ysd2/T0JRAJREAYBhjDX4R+yv6Xn/Ng1if04vry1/\nzC4A6qEAAJLMGvtrdd7/fLG/i++3/aO8EVn+lCIAdVAAYAhjdf8pehn8izF57G97L1FBVEYBAEQr\nFfs7+rz/VCUG/zqO/ZWq1Tyxvyuba8vH2AVAeYyZoHsM/pk0dv+dx/5WLWJ/tns1/r8ehsAOALo2\n1tb/EYN/oTVax/4yzvs3dRz78+4M3Igsf8EuAMqiAACCjmzBam/9a4v9RaoR+4t5yNeM/Zk4EAgN\nUACgW3N0/7HraDnvf4+DY3/bNVvE/lyDf7Z7noksn2QXAOVQAABexP6u16sd+2sw+Bd6yJeWMfh3\n8cdn/pPBf00RgDIoANAlfYN/xP7iChZFg38+R8T+Yq8BhVAAALsR+/PfV2rOQOK7/5hrtWJ/ocG/\nmO7fZFln+Vt2AbAfBQC6o6/7b6WX7n/y2F/Mmqnd//a1Wh9VYCoUAOjKHIN/nPdvtXb/5sOvduwv\n9KCNHfzLiP1dvbZ6+TMtn2MXAPtQAABXZh38y11vgvP+Xd1/jL2xP1sR8fL3RxGAPTTucwJWY3X/\nKUaK/U103n+B2N/Va48trwGZ2AEALsT+r+qeTl1j998q9mej4Lz/UvYM/jlifye239+NyPIFdgGQ\nhwIAXehv8K+37j/0PW3XYgb5FA3++ZTq/m1yB/+210KvAYkoAIBk2mN/R5/33yj2lzr4d1TsL+Y9\npoTuf71/+Rd2AZCOAgDq9df9l5Lb/cesk3Ovgtjf3vP+fUrG/mKKAt81z+Bf8KMCIBIFAFQba/CP\n2N/u2J9tmVXJ2F/JwT/ZXCsx+GcpFpavsAuANBQAwLSxv1CxEPO+QWN/OVKKgpikgMmx9X+yng1A\nEYAEbBxBrbG6/xQaBv9cD/bUh/lEsT/fz1Mr9ldjkwbTYAcAkysV+zt68C9VzfP+if3VjP1drW/+\nnp6JLP/OLgDiUABAJQb/TNpjfx2e968x9rdn8C/22AbAQAEABB3Z/ecM/oXWKB37C8kY/DPVPu8/\nVugz/pgHvS/2Zwp99GErhIx1lv9kFwBhFABQh+6/5Totuv9I5sO/RezP9R5X9++TM/i3p/sXsf9e\nT+vw/EcYBQBUGWvwj9hf1cG/yWN/Fz/fxToPf4WWb94p+7sEbSgAMCFif9fr+f5MSnT/DQb/eov9\n+e5PHPw73bd2/jf3D7d953sUAXCiAIAaY3X/KXoZ/Iuh8Lx/LbG/mB0BV+wvZfDv5cP/0dNbzw8K\nUABgOsT+rnUS+4u5Vuu8/9We2J9vHRF/7E8s923XNj73Xx/+N0/u5emr31VWWEMLCgCowOCfSWP3\n33nsbzV67O9l9y/y8PAXEXlsvAaYKABwuLG2/o8Y/AutoXjw74jYX+gwoNjBv9TYn29HYGfsz9b9\ni5wf/m98/dvK/o5BAwoATOLI7r/21n/Bf373pNPYX8rgn6v7j5E7+Ld9f6nYn4i1+xcRecMzZgFg\nRwGAQ83R/ceuU2PrX3H3v8qN/R0x+BfzMG8V+7Nk/m3dv4jIzc29/OgPvqXs7xqORgGACRD7u15v\ngNjf9lvVOu9/lTP4F/osP/O8//O9l7E/kevu/2aNBD66lzfJNykCcEIBgMPoG/wj9hdXsCga/PM5\nIvaXes31Wmbszxz8M7f+Hz96+bowEIgzCgBARIj9he4rOGfQQ+wvNPjXIva3c/BPRE7d/8MyD79+\ns3yDXQCICAUADqKv+2+ll+5/8thfzJq1Y38iRQb/1u5fROQNwkAgzigA0Nwcg3+c92+1dv/mA75V\n7G/7o6UO/tWK/aUM/p2u53X/D0vfy1vk68r+DuIIFAAY1KyDf7nrNY79+Zapfd5/zv9r1Iz92dYQ\ncZ/3f1ozvfu/kfMswM/I1ygCJkcBgKbG6v5TjBT7iykWPGJjf9trM8f+Lu6N6/7NwT9z63/dDeDj\nAFAAYECxLd6eTl1j998q9mczUOxvz+Bfw9ifed7/+vA3Y38Py1xu/a/W3YC3yr8qK8jREgUAmulv\n8K+37j/0PW3XOhv88ynV/dvkDv5tr7le23nev0jc4B8xQJgoADCpI2N/OYN/oTVKxP5yNTzvP/SQ\nz5kJML9XzIPe9R6Totjfw9LuWYC3yZfZBZgUBQCa6K/7L6XmLkLOvQq6/73n/ad8rz2xv5iiwHdN\nWezPNQvArsC8KABQ3ViDf8T+dsf+bMusasT+Sgz+yeZai8G/nd3/+rC/kXvnbsB639vli8r+jqIF\nCgAMYtbYX6hYiHnf4LG/FClFQUxSwJQz+JfR/fu2/s/f6vq+d8jnKQImQwGAqsbq/lNoGPxzPdhT\nH+YTxf58P0/MNn9s929b3/z5LtZ+IeZ5/yvzvH8z9ifiHvzb7gYwIDg3CgAMoFTsj/P+7b+2IfYX\nXEekXOxPxBr7E0kf/LPdt379TvmssoIdNVEAoBoG/0zE/q6W2Xvev8bYX+rgX4PYn6vL3xYFHAw0\nHwoAVDHW1v8Rg3+hNVrH/jIG/0y1z/uPFfqMP+ZBn3Pe//Z+89eVYn+hwT/X+35N/l7Z313UQgGA\njmns/kvRFvuLZD78Ncb+fEoN/h0c+/MN/l3+uO6zATAHCgAUR/cfev/Asb8VsT//GhVjf+vX52/l\nv8+WDvgN+bSyv8OogQIAnSL2d72e78+kRPffYPCvt9if7/7c2J9x3v/p7Zvuf33474n9ud4nIvJU\nbuXd8kmKgMFRAKAofYN/M8X+XNdKxf4i1TjvX2vsz/Y+V+wvY/Dv0dPb08PfjP2JxJ3374v92d73\nWOwzBRgTBQAGRuzPf98Bsb+Ya6GH/F5Hxf7Ect927Qrn/dvu8339VG5PX/+2/CW7AAOjAEAx+rr/\nVnrp/juP/a2I/RWL/W2LhKccDDQVCgAUweBf6P0DD/4dEfsLHQYUO/g3eewv9PVjuZf3yseV/d1G\nKRQA6Misg3+563Ua+0sZ/HN1/zEmiP2lfESw7f7X9z7lgKBhUQBgt7G6/xSttv5bdP+NzvvfPiCP\nGPyLeZi3iv2d7o3r/teHf2zsTyT8EcG247fNA7xPPqrs7zhKoABAJ2bt/lvF/mwGOu9/lTP4F/os\nv9R5/zf31vP+zdjfwzJ1Yn/bQmC7I/B++QhFwGAoALCLvsG/0br/0Pe0XYsZ5FM0+OdzROwv9Zrr\ntYaDf+clr2N/sUXCdq3tPXwUMB4KAAxEe+xvTx5/K2dHpNTxwlJ28K9W7C80+Nci9ld58G/li/3Z\n3ueK/a33bh/26z0fkA+yCzAQCgBk09f9t5Lb/YfWGST2l3KtZuwvZs3asT+RKoN/56XLxv588wDb\n9dE/CgBkGWvwj9hfVuzPfMC3iv3tHfyrFftLGfw7XW8X+0uJAW53BLZff0j+UNnffeSiAIBisw7+\n5a7XOPbnW2ak8/5d64jYf38x5/2f1ozr/tcHcG7sb/t1qPv33fNUbuUj8n6KgAFQACDZWN1/Cs2x\nv1Ln/ReO/W2v9XDe//aeUrG/i3vTYn8ieef958b+THwUMC4KACgV2/3v6dQ1dv+5sT8Tsb9dg38d\nxf5KnPe/vnfb7W8HAbevfVTep6wRQCqNE1ZQrL/Bv966/9D3tF1TMPinLfZnkzv4t73mek3pef+m\nlNif+XVoOBB9YgcAHTsy9pcz+Bdao3TsL4Wi2F/qTEDO4N/2vZ3E/nyDf6ZSsT/fa0/lVj4u72UX\noGMUAIjWX/dfSs1dhJx7G3f/Nc77jxHb/fvEFAW+a0pif77BP999JWJ/vuFAdgL6RgGAKGMN/hH7\nqzr4pyn2Z9Ny8G9n9x9z3n9M7C9l8G/9OvTa+vVfyXuU/W8DYlEAQJFZY3+tzvsn9nclJvZnyhn8\n2xH7e3j7/tifyPVRvjmDf+s620Lgk/JuioAOadxvhTJjdf8pehn8i0Hsz/qfoXV83b9tffPnu1j7\nxenh/+jp7enh//jm/vTwPyL2l3uP+Rr/RkC/2AGAEqVif0ef95+qxOAfsb/DzvtPif1JvfP+98b+\nUu+xFQOfkV9X1igghAIAXgz+mTR2/x3E/mJ+W5pif6mDf8pif6F0QO6OQMxHCOgLBQCcxtr6P2Lw\nL7RG69hfwtZ/ydjf6sjYX8w2f+i8/+395q8PiP3Fnu8fG+kL3bOu6/so4HPyK8r+NwM+FAA42JHd\nf+2tf22xv0g1Yn8xD/masT9T7nn/ymJ/D98mP/aXe4/5GkcE940CAFZzdP+x69TY+p8s9rddk9jf\nwy0KYn++e0KnAq6vm2t+Xt6h7H874EIBgAMR+7ter3bsr8HgX+ghX9re2J/v/tzYn3He/+ntm+5/\nffjvif35YoC1Yn/ma/ZC4k6+KG+nCOgABQCu6Bv8I/ZXLvYXqYfz/kvF/mzvc8X+Mgb/XLE/Ef/g\n33lJd+zP9r7WsT/7wUGvC/SjAIByxP789x0Q+4u5FnrI5woN/tWO/Zm/HjT2F/ooYL3XvoNwJ0/l\ndXks9/JleRu7AMpRAOCCvu6/lV66/85jf6vc7j9mzdqxP5PC2F/p8/63r7ne/1jurl6HbhQAOJlj\n8I/z/q2OiP2FPi6IHfwj9hf8em/sL6ZAeFjj9Yv7/01+Ttn/psBEAYDGZh38y11PYeyv9Hn/ru4/\nxtGxvx3n/Z9eD8T+Uj4iyI392YYD0wqEu6vP/W2vQRcKAIjIaN1/ipFifxOd9x/zMG8R+7u4N677\nXx/+odifa/CvRuwv9TXbmub95scBX5e3KPvfFqwoANBQbPe/p1PX2P23iv3ZcN7/xdeu10qd939z\nf3r43zy5Pz38zdjfwzLh2J/tvj2DfzVjf+vW/+V7Xj9d/4a8mSJAIQoAdDj411v3H/qetmsxg3yK\nBv98NJ73v73meq3hef/nJa+7/9giwfV17CBg6OMBEXcxcH7P9dY/HwXoRAEAZbTH/vbk8bdyiqJS\nxwtL2cG/o2J/Me8xdRr7s72vRuxvfT0n9re13R34pryJXQBlKAAm11/3X0pu9x9ap+PY397z/n1a\nx/581/YO/m3XESH2J9cP/C1igfpQAExsrME/Yn9Zsb+S5/3Hxv72Dv7lbvMrOe//vHRa7C8UA/Tt\nCJSK/bk/57c98O+urn9b3qjsf3PmRgGAymaN/YWKhZj3NYj9+ZapFfvLkZIGiEkKmGy/v5jz/k9r\nxnX/60MwN/a3/drVtW+/3hv78w3+Xb5+dzH4t7Ve/z/5IYoAJSgAJjVW95/iiME/3/ercd7/gLE/\n389zZOzPOO9/5Yr9ieSd93907M9cM/bjgIf3vH5xfS0YHt/zUYAWFACoiNjf9XqxfybE/rqI/Yns\njv3VOu8/NfZnWzu0tX+5zutXD/ur997fy+O7O/ne7WNlDcicKAAm1N/gX8/dv/LBP1NK7C/mt0Xs\nr9jgnyk15mf7XN891V8m9nf5fV8WDEbn//T2B1f3oz0KABzsyNhfzuBfaI3Ssb8UCs/7j1Vi8K+T\n2F/K4J/rfa7P7Ldfpwz+2V7zxf58swGuif/Hd3fy9PYH8vhO5MWrC7sAB6MAmEx/3X8pNXcRcu5t\n3P2nnvdvk5Ptj+3+fXIG/xTG/nyDf6ZeYn8P1yLPAXi59W9aGAU4HAXARMYa/CP2V3XwT1Psz6aj\n2J+vq/fdt2fwL/Zz/T2xv+tzAO6810Xk1P0v9yJyJyL/yy7AkSgAUBixv+v1fH8mJWJ/DQb/eov9\n+e7PGfzbEft7ePv+2J/I9Tb/3sE/22uxsb+H9/gH/8zuf334X/i+iPwXRcBRNO7VooKxuv8UvQz+\nxag4+GfTU+wvZkfA1v3b7jN/vou1L2N/68P/8c396eF/ROwv9x7ztdjBP1usL+W6iNH9f1/EMSqA\nRtgBQEGlYn9Hn/efqsTgX8exv1K12p7Yn28dkSqxP5H4wb9Vzdhf6j05g39bocG/bfdv/dz/TkS+\nwS7AEdgBmACDfyaN3X/nsb+Vpthf6uCfsthf6B8Fyt0RiPl4wLw/57x/W+zvgtn934lYvj0aYQdg\ncGNt/R8x+BdaQ/Hg3xGxv9jDgErH/nw7AqGPPhTG/rZfx0T6asT+HtZIO+//4vq2+18f/lt3IvIV\ndgFaowBAAUd2/7W3/gv+87snjc/7LxX7Sxn8ax37s71fcewv5SMCX+wv9R7bxwO5sT/vdWPr38rs\n/u+M/0NTFAADm6P7j12nxta/4u5/lRv7O2LwL+Zh3ir2Zw7+vVQy9idS/7z/0Gf929dSY38x10Xk\nuvtff23zT+wCtEQBgJ2I/V2vN0Dsb/utap33v8oZ/KsZ+xM5P/w93f/68N8T+/MVCcPF/kzb7n/9\n+gsUAa1QAAxK3+Afsb+4gkXR4J/PEbG/1Guu2F/G4J8r9ifiH/w7L+mO/dneN3TsLzT4RzSwGQoA\nNEDsz39fwTmD2O4/5lqt2F9o8K927M/89aCxv5SPAqrH/ly7ALbuf73/H9gFaIECYED6uv9Weun+\nJ4/9xaxZO/YnUmXw77z0vtiflvP+U2N/2/P+T7bdf2jwj2hgExQAg5lj8I/z/q3W7t98wLeK/W1/\ntNTBv1qxv5TBv9P142N/tq/3xv5yBv+2QrE/Ect5/7bYn2nb/a9ff4ZdgNooAJBh1sG/3PUax/58\ny9Q+7z/n/zVqxv5sa4jYt/47j/3ZhgNjCoTzGtef41eP/a1f29yJyN9QBNSkcR8Xmcbq/lOMFPuL\nKRY8ejrvX0vs7+LeuO5/ffiHYn+uwb8asb/U1/zFQKPY38pWIPAxQHXsACBRbM24p1PX2P23iv3Z\nDBT72zP4F5rkL3Xe/8396eF/8+T+9PA3Y38Py7Qf/Ksd+3N/HFAx9udLBtyJyMfZBaiFAmAQ/Q3+\n9db9h76n7VrMIJ+iwT+fUt2/Te7g3/aa67WG5/2bQoN/rntdX9t2CGyvhT4eEHEXA+f3xMX+bOf9\nO2N/MSf9bYuDtUBAFRQAqEB77K/Wef+xGsX+Ugf/jor9xbzH1Gnsz/a+vbE/12u5sb+U3YHTPbGx\nPxF/9+/b+v8zdgFqoAAYQH/dfym53X/MOjn3Koj97T3v36dk7C+mKPBdI/YXHA40X4uN/T1cOyD2\nt2V2/+ssTQGkAAAJ3klEQVRaKI4CoHNjDf4R+9sd+7Mts6oR+ysx+Cebay0G/3Z2/+vDPjX2F/rI\nwLcj0FXsz/fADnX/tof/nYj8CbsApVEAIMKssb9QsRDzvsFjfylSioKYpIApZ/Avo/v3bf2fv5X/\nfH9f97++tv16b+zPt7V/+Xqh2J+IP/aX09F/X0Q+QBFQksZ9XUQaq/tPoWHwz/VgT32YTxT78/08\nJbt/23vNn+/i+70Q87z/lXnevxn7E8k77/+I2F/Oef9b29ifdTDQFfsLxfhSun8+BqiCHQAElIr9\nHT34l6rmef/E/lTF/kSssT+R/YN/pWJ/scOBMQXCw+t5sT8RsZ/37xKK/cXOBJjr/QG7AKWwA9Ap\nBv9M2mN/HZ73rzH2lzr41yD2lzL4Z4qN+cXck3IuQI3Yn4i4Y3+x0b+Y7p8DgopjBwAF9Db4F1qj\ndOwvJGPwz1T7vP9Yoc/4Yx70Oef9b+83f10p9pcy+Od6X8xQ39GxP2vBkBL7E4mL/aUkBEREfpdd\ngBIoADo0b/dfirbYXyTz4d8i9ud6j6v79yk1+Hdw7M83+Hf54+bH/nLvMV/Ljf3FDP5Z+WJ/KR17\nTPefWizAiQKgM2MN/hH7qzr4R+yveOxv/fr8rdJif3sH/2xb+Xtifznn/SfH/kTyYn+hf0fgXewC\n7EUBAAtif9fr+f5MSnT/DQb/eov9+e7Pjf0Z5/2f3l4h9ueLAYY+s0/9XN/2Wmzs7+E9hc77Lx37\n8629/vq3KAL2GHmfdzj6tv5niv25rpWK/UWqcd6/lthfzI6AK/aXMfj36Ont6eFfOvZne1/ukJ/v\nHvO1nNhfzOCfTbPYn+3eUkUF2AFALmJ//vsOiP3FXKt13v9qT+zPt46IP/Ynlvu2azc4739v7M93\nT+ijgPXe0ODf1nbw7+p66dhfyc/v70TkV9kFyEUB0Al93X8rvXT/ncf+VsT+uj3vf8v1OX7wgW9e\nrxH7EynT/eeeJYATCoAOMPgXev/Ag39HxP5ChwHFDv5NHvuL+bpF7O9hjfzz/q2xv+9f3XaWEvsr\ntY3/i+wC5KAAwEtHdv9Hbv3nrtdp7C9l8M/V/ceYIPaX8hFBydhf6OMB8/1Vz/svEfsT2d/9EwvM\nRgGg3Bzdf+w6Nbb+FXf/q9zY3xGDfzEP81axv9O9dWJ/IvXP+0854tdVDJj358T+RMR93r/L3thf\n6sP/VkR+gV2AVBQAkDbdfw+Df4PF/rbfqtZ5/6ucwb/QZ/mlzvu/ub8473/1hme3p4c/sT/jekzs\nT8Qf+2vdld+JyFspAlJonPrCS/oG/4j9xRUsigb/fI6I/aVec722c/DPjP2J+Af/zksS+7MO/vm0\n6v7NNRGNHQBEIvbnv6/U8cLSR+wvNPjXIvZXefBvRexPwtv9MYN/ta1FxE+xCxCLAkApfd1/K710\n/5PH/mLWrB37E6ky+HdeWnfsz9fpX78eH/vLOu8/Ru3u3/w+iEIBoNAcg3/E/qzW7t98wLeK/e0d\n/KsV+0sZ/DtdHzP2Fzv4tyf2J+I4798nJfYX828H+Li2/s0i4sfZBYhBATCtWQf/ctdrHPvzLaPp\nvP9VzdifbQ0R+9Z/Rve/PgQ1xP5sw4ExBcJ5jYNifzldd81Df25Fnr+RIiCEAkCZsbr/FCPF/mKK\nBY/Y2N/22syxv4t747r/vYN/NWJ/qa+pif2JtI/9ubr/W5HndyKuTzJwRgEwpdjuf0+nrrH7bxX7\nsxko9rdn8O+g2N/68Ddjfw/LtB/8qx37c38cEB/7c573H4r9tRj2c9k8/F+7FXn1CbsAPhQAivQ3\n+Ndb9x/6nrZrnQ3++ZTq/m1yB/+211yvNTzv3xQa/HPd6/ratkNgey308YCIuxg4vycu9pd83n/I\nUd2/xXN2AbwoAOBwZOwvZ/AvtEaJ2F+uhuf9hx7yuTMBKYN/rveYOo392d63N/bnei039peyO3C6\np0bsr+bDN7D1/9r6axH51sIugAsFgBL9df+l1NxFyLlXQfe/97z/lO+1J/YXUxT4rimJ/a0P+1Fi\nfw/XFMT+ROp1/wHrw/+5iHw34cedDQWAAmMN/hH72x37sy2z0hT7s2k5+Lez+4857992X+gjA9+O\ngG3bvnTs7/ocgDvvdRFH7M8X1cuJ/cXECW3fI3Hw7zXLHMLX2AWwogCYxqyxv1CxEPM+Yn8nKUVB\nTFLAlDP4tyP29/D29Njf9mtX1779OvbfBMgZ/Nvyxf5O12PO+w99Bp8jtvtPYNv6fy4ir738zzsR\n+SpFwBWNe8JTGav7T6Fh8M/1YE99mE8U+/P9PDHb/LHdv2198+e7WPvF6eH/6Omt9bz/XmN/Jc77\n33b/yef9+yiJ/bmsD//XAr+NWbEDMIVSsT/O+7f/2obYX3AdkXKxP9F73r+ts895uG8VP+/fJRT7\nO2LSPjD4t3b/puci8iV2AS5QAByIwT8Tsb+rZfae968x9pc6+Ncg9pcy+GeKjfnF3JNyLkDWA38z\n+Gfliv3FRv9adv8e261/89fPI34rM6EAOMhYW/9HDP6F1mgd+8sY/DPVPu8/Vugz/pgHfc55/9v7\nzV9X6v5TBv9c74uJ9LWO/UVdt8X+fAN6sbG/GrsBkYN/a/fvWmItDD7LLsAJBcDQNHb/pWiL/UUy\nH/4aY38+pQb/lMT+1q/P38o96FfyvH/fPeZrubG/mPP+nXxdeKzS3X8E1+Df+uuLexN+K6OjADgA\n3X/o/QPH/lbE/vxrKI797R38s23ll4z9RV1Pjf2J5Mf+XAN7hc77t8X+RK63/s2vP8UugIhQAAyM\n2N/1er4/kxLdf4PBv95if777c2N/xnn/p7dvuv/14b8n9ueLAcZG+nz3rOuUjv3FnvdvFXoQl9ji\nL/BRQUzs7+o9cnnfJygCKABa0zf4N1Psz3WtVOwvUo3z/rXG/mzvc8X+Mgb/XLE/Ef/g33lJd+zP\n9r7cIT/fPeZrJWJ/MddFMmN/It3E/rYP/K0jwgvaUABMjdif/74DYn8x10IP+Vyhwb/asT+x3Ldd\nu4PYn++e0EcB672hwb+t5rG/lk/PiME/10N+uyuwfr3++sOT7wJQADSkr/tvpZfuv/PY3yq3+49Z\nc7LYn9bz/g+P/Ym06/49fIN/tt/Ga3JdCMxsefFi6gIIAIApsQMAAMCEKAAAAJgQBQAAABOiAAAA\nYEIUAAAATIgCAACACVEAAAAwIQoAAAAmRAEAAMCEKAAAAJgQBQAAABOiAAAAYEIUAAAATIgCAACA\nCVEAAAAwIQoAAAAmRAEAAMCEKAAAAJgQBQAAABOiAAAAYEIUAAAATIgCAACACVEAAAAwIQoAAAAm\nRAEAAMCEKAAAAJgQBQAAABOiAAAAYEIUAAAATOj/AQchMybioTQRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106d1dba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_pdf_contours(Dirichlet([1, 3, 1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAFdCAYAAAB1gNVOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJztvTuQG1mWpvmjEYFHkBEMPiqZycqsyuzKqs20qq6t7uld\na7NVVxm1lVFHHWHUUcdabXWEVlttpZRZs1bGbKWxUXpHmBmzrrHtsmR1dpNFbpIZZJCBAAgYVrh+\n3I9fXHc4XhF4fJ8ZDA53hwMIEjj//c8597am06kAAADgsPiD234DAAAAcPMgAAAAAA4QBAAAAMAB\nggAAAAA4QBAAAAAABwgCAAAA4ABBAAAAABwgCAAAAIADBAEAAAvRav3FtNX6C2YQA9hxEAAAsBSI\nAIDdBgEAAI0h6APsDwgAAFgaBAHA7oIAAIBGEOwB9gsEAACsBMIAYDdBAADAXAjyAPsHAgAAamkS\n/BEIALsHAgAA1gIiAGC3QAAAQCUEdYD9BQEAAGsDwQCwOyAAACAJwRxgv2lNp3zHAaBMs+DfT+wb\nSJKm03/fWu87AoB1gwMAAEuQCv51+wFg28ABAIAS80f/cZD3jwf5PS4AwHaDAwAAC9CPtqvEQJ8a\nAoAt5+i23wAAbA/1QTsO/kb8M0IaAGAXwAEAgAWxAH+k9Bii2IcLALC94AAAgKSmo38f/P3jmEHN\nMQDYBigCBIAVg391O6B0qVAQ+O8oCATYMnAAAKCGOLj74G/HjhPPO5b0YVNvCgDWAA4AwIHTfPRf\nFfzr0gAfJL2WNNZ0+m9xAQC2CIoAAaCClMVvwf9YRfC37TN3848fSDpSq/WXjDYAtghSAAAHTPMq\n/SPN9v1b8K8SClYHcJbdj5d9mwCwAXAAACBByvo3/MjfRMGpirZAu50qjP7NCThVq/UfcAEAtgQc\nAIADpXr0P8/6j4N/yh0wBtH2IHEOANwGFAECHCCLFf7VBf9Td/5x9HwL9m+z7ZeSXkl6ren031AQ\nCHDL4AAAgGOZ4O/3VV3rbelIq/WXU+YGALhdEAAAB8Zi1r9nXvBPzQegxP4Pkj5q/H4BYDNQBAgA\nCVKj/zjnf6yiyO9Yamm2DvBIYb/OJD3Ozn8s6YyCQIBbBgcA4IBoNvr3Pwux9X+kEMSPi/uWpHbi\nqT1J19n2paTJsTR9nO2wSYIA4LbAAQCABKnZ/uKRvwv+p5LuSjrPbg8l3cnuzyX9IDveOlZwAD6S\n9FCt1l/hAgDcEjgAAAdC89F/lfXfV7Dys2M++PcUAn5fUje6/FAh+Peyy1+YEzBQXBwIADcHAgDg\nAJg/419c+Jey/p0I8MH/PNs+z07vZZewFMBpdt/NzpOkizNp+rmkS7VafzOdTv8VHQEANwwCAOCg\n8YE/ntAnlffPjvngf1chuJsTcKTCBbiT3Q8VZgL2Mwd/fybpx5IGarX+asrcAAA3CwIAYM9ZrPCv\nyvp3eX+b5deCv9n7JgTMBTAnQAoOwHX2Mn7/958qrBZIKgDgpkEAAECGFwQp6z/L+/dVBHwr/DvS\nrBDw3FFwAXrRsbGky68VUgF/PZ1O/zUuAMANgQAA2GOWK/zz+53131IR9O8oBHIr/jtWEAJ33FNT\nxYDnki7cywwkDY6l8ZfCBQC4WRAAAHvK4oV/ti+e3S8TAW2FX4w7CoHcWv/s3uoAzBHoRZe+Vhjx\nx+JgIumbx5J+SEEgwA2CAAA4OKoK/45Va/1bcLeZgO+qnAaw0b+lAtrRSw0Ugv1dSf9f9JauJT0P\nqQAAuBkQAAB7yPJtf3HhX3aez+/Hgd+Cvy8MNAcgLgSUpHfZft9wcCnp3bF0+WO1Wn87nU7/JS4A\nwIZBAAAcFE3b/pwIsDn9Le9/7u4faTb4W42AXTYWAXeycyylYLMCDyX9j0+l6WsKAgFuAAQAwJ6x\n+Gp/8XK+Nvo/K6r+vfV/x9374G+3fNa/qdQbqdMb5q80uu5K1x2p3ypcAGOo4A5887Wkl4t+bABY\nEAQAwEHiR/ne7j9VmO73QTitrWK0bqN7X/yXCv53p+qcX6p9NFGnN1S3Myq98nDU0eDdiUYXp1I7\nG+SPVQiA746ly6/Uav16Op3+OS4AwIZAAADsEctN+uP3Jwr/epod3cdB3+4fDdW/e6XTs0u1NVFH\nQx1pUn4nnbZOHgx0dfdKb3r3VWoJeKdQD/B3nwoXAGCzIAAA9oT5hX+e1KQ/NYV/pyoK+8z6N0fg\nkWaCf19XOtFAXQ3VzgRARyON1JEkXelEJ50rdZ+MdNE716h/Fl5zoCAAXkn65ksKAgE2CAIAYO9p\nOulPTeGf7wAw678kAkLwPz+70ImudKpL9Z0A6MjVAairU13qUqfqaKT2g4leSRoNzoIAeJ/dvjuT\nLr/c2F8F4NBBAADsAYuN/o2a0b9N9+sL/86jWy4CpqXgf64LnSqkAOy+6wSAJA3VVV9XGuhEXQ3V\neTDU859KI52F4H+hMFfA3/1UrdZ/nk6n/wcuAMCaQQAA7DUrjv79qD+2/s8l3Zc655c6PbssBX9L\nAZzoSh2NkgLgRFe60kBtTUKa4IGCCHh/FtoFv5NLBTBDIMC6QQAA7DgbGf1b258t3hO7AflkQEOd\nPwhB327nCk5APxcAQ3UVOgHammiitiRpqI4GOsnPkSQ9kH73VVe67gYBcCHp28fS+OvFPyIA1IIA\nANhb6kb/cyb9idv+bPQfWf/3Hn2vfpbzP9dFKfif63t1NVJfV+pqlBcDStJEbQ3V0YkG6ucuQRAJ\noyddPX/3RSgG/E5ZKuBrCgIB1gwCAGCHWX7BH7t3o38/6c+5itG/L/xrq2T9n3QGWRAfZLZ/SANY\nCuBUl3kaoB21A07U1kgdXekkD/75sZ+19fL6R1kxoKR/Ppae/1HzPwwAzAUBALCj1Af/JRb8sal5\n29kpFvht2l+/AuDdUPV/In8rhEBwAgoBEHcCGAOdqO/qAKxWYKK2Lj8/1eDd/eAAXEh6/ikFgQBr\nBAEAsLcsMPq3XXHu388BkE/4E6r+TzqDUt4/3L7Piv8GeTqgo2E+J4AU5gOQpFGWArhSX10NS4WC\nQ3U0POvq7z8/kb7rSv8s6UtJ//BHFAQCrAkEAMAOsvbRf1Xu34/+TQT0RnnVv7f+zQXwtQDmANic\nAB1n9ftJgS51mrkE4XgoG+zq6klfv/vqq5AG+E7SizPpkoJAgHWAAADYK6oW/ImPrTb677vgbla/\nv/mWQKsBONGVJM2IABMQQSAMS8cmamvwsxO9/O5H0j9Jeirpv/ySgkCANYAAANgxmrf9rWH0bw5A\nJgKKnv+Q7zf7v6+rPJDHQiCIhXdqaxy6ASauG6AdugFO9U4djUrdAJOsMuBSp7r46lyj786CC/BU\nFAQCrAEEAMDesMLov6fy6L+n2cr/u1O1jyaR5V9Y/2b/+7kA8o6AyUDd4Uida6lVbgbQtD3Q5b2O\n2u2x2hrn3QIjdTTMugQuHpzrt1/9PLgAfyrp/6IgEGBVEAAAO8Tyo3+pcsW/eM7/U/fYV/7nuf+B\nOhrl8/37QkAvDnzwP30zUuu9wpK/4+gzHUln70c6vTNS90HhAAzzKYS6utSprn52oufffSH9XsEF\n+O+4AACrgAAA2BGaFf5Vjf4rZv1ra3aWP//YRv895bn/ePTf1iTP4ftWQAv+Z69H0huFwG/3MXel\n1nvp/mSg9oOJsskCJYVUgL3axS/ONfj2flYPcKZW679Op9M/wQUAWAIEAMDescDo/1Rl2/9YhSNg\nbkBb0t2hTs8udZRV6sejfwv6Rf7/XTn4v5f0ToUDEE8J8D57zaF0dj1S9+Fz6SQE/4FOsmmD+7o4\nO9ff//F96YWCC/B//4larV9Pp9M/RwQALAgCAGAHWM/oPzrfRv828c9R9ti6AUwE3A+jf1vdLzX6\nL+oAimOnb1zwf6UQ+N9m99fR2+xlx87Cse5QevIkiIBRlgoIswyc69XPHurlP/1I+gcFEfANqQCA\nZUAAAOwVVaN/lffZtL8W8M9VdAAcucvckXQ0Vac3zPv5bYnfjoYJJyCz/q8GIedvI3/bvs7upbIL\n0M1e09yBcdj18ItXGraDABhkywtd6L4ufnWu0e/PggB4+lO1Wv9pOp3+n7gAAAuAAADYctYz+p8z\n7a+N/uPWv34x578t61s4AYMZJ8Cs/64P/q+zbXMDfCHgdfa671UIAOcOnGkkffFck/ZR3mh4oXO9\nfPCRfvuLn4eOgKeiIBBgCRAAAFvMYkv9Nhn9z2n9i4v/sta/MJ1vUQRo+f5OZtDb6L+tcWH9W8B/\nH21fqywCzBE4UhAMd7PjmUNw1h3pyWfPdKnTzAEIqYDLX57q5dMfZS7AYwoCARYEAQCws9SN/qVC\nBDRo/asq/ota/7pusp4i/z/IugBGOpkMQp+/Bfu3Kmx/EwJj1acB4lbBI+l+b6DPfvBtPvHwhc71\nUo/18hefSf/Ukn4j6b/gAgAsAgIAYEtZbfTfj24VrX+p4r+e8uK/Tm+Yr9LnC/zMAejnU/iOwj4b\n/VvQt2CfcgGiCYH0PntvYxUCYKxcpDxpv9Llg2e5ALjQub77w4f63a++CgLg6bFarf82nU5/iQsA\n0AAEAMBOssLo39v+5yrm/bfiv7vZ8aMw738ny/v7e3MAuhplpXnBASiN/ocKQf+tZoO/HZeKYH+U\nvZeJys5A9ivV6kpP7j3XoB0EwAt9pC/1W738xWMN/ux+SAX8R9YJAGgKAgBgC1l99C+3z5Eq/vPz\n/vfc8d5I3c4oX+zHpwHMFbD7jkY6uRqE4H6tENSt+v+dyk6AFwCpdsCu229CoCvpjnR2NNInXz7L\nHYAL3dezsyf6+1/cDwWBfy/pH0gFADQBAQCwczQZ/Vvlv3vsR/+p4r+o99/b/z4N4O1/X/zX8bb/\nWxWFfEMVgd/OeaMw0o8FgO2z5z1QURugsP3xvTca/OCprnQS5gXQQ7365cNQEPgbSU8/pS0QoAEI\nAIAtY7nRvzQ7+rd71xXg2/188Z/v/V/A/rf7vPjPgr3l8W30b/t98H+TvUVf8GftgPcSH7Ur6bmk\nu9LDB6/0pP1ML/RYL/WtfqLfhoLAP2tlSwb/783/hAAHCgIAYIuYH/yrRv+Jmf6qiv+qev8T9r9N\n/RtX/8f2f3foiv/M2jcXwNcEDFUE/0wYfIjWBjj28wGYQ2DvrSvpTDq7U04FvNJDPfvDT/S7X3wl\n/UrSb1gnAGAeCACAnSX++vr+/5rivyr7v1c8Nvv/VJelYJ8K/rn9b6P+2AUwQWB1AFnw//BeGgyl\nQZQG6Pek/lg6jlsEu5KeKW9X/PjeG1384FmeBnihx3r5p481eJoVBP7tnyz/pwU4ABAAAFvCcqP/\n1MQ/8T6Vi/+q7H+r/neT/7Q1yVv9LAVgk/9Y8M+r/22E74O+v3dCwAf/saQr94nG12F/vxeWBpAU\nZhOUe59ZKuDJg+e6aN/XKz3VCz0OBYE/zwTAb6RW6/+dTqc/xQUASIAAANhJUl/d48SxaOY/O22O\n/d/pDUv2f13+vz2ZlKv/fcAfK1j+NvrP9g+G0uV1CPyD6FMMsrc7zpyBM3/Qv+cz6ezBSJ999q1e\n6aE+1zchFfDHn+jN04+lf6FsnQDaAgFSIAAAtoDVR/81xX8W+ubZ/+7p3c4on/tf0vz8vx/5+/7+\n2AkYh5z/wAX/sWZFQL7Pi4CewqqCd7Lbs3D/6OErPTkpCgK/7PxW/88ffxxcgP8p6b+TCgBI8Qe3\n/QYAYFGqRv81xX/z7H8rEOzNtv81zv9LhRBIBf/s+GBYBHi7j2+vJV0qxPvLa+mtLSn8Jrt/rZAG\neC51/1H6RM/0hb7RZ/pWn+upPvnZN9IfKcwNcPRYrdZ/XqCzAuAwwAEAuGXWN/qXym2B7tSU/S+V\nnYCs/c/n/8Mp9fn/nHiVv5hhefTvg/+H7BRzAmZmOLiW+u+l49cqihXvKFgDD6SPH4aCwM/0rV7o\nIz3TEz3/1RPpaZd1AgAqQAAA3CKL9fzX4Yv/7HGN/S93X5H/l1TZ/1/6DFYA6F0AqZz/nxTHveVv\nwd+LAGXbx9E+vZE+knRsNQB3FOYLeCnpmfT4By/1mUI9wEs91ssnH+m3v/p5EAC/OWOdAIAIUgAA\nW008+o9H+HHrnzTX/pfK+X//WEX+P+7/98wUABrxFL5RCmDgWvvi4B+nAWzfWxXpgNe2poClAl4p\n1AK8kO5/O9ATPdMnepanAjpfvQ3zAvyZJP1SrdavSQUAZCAAAG6J1Ub/UZtf1fz/sf1v93H+Xyrl\n/6XC+jfiGoASFuiVuFf5sbf/pSLQS7Mi4K27txmG375WqAF4rRD8TQQ8kx5OXumJnutzPdXneqqf\nPPit9HOF2xdSKAwAAAkBALDFxJnwJsV/bn9s/9sl4vY/JxLaR7PBvyrolzoAjHiFv4h40p9BtD12\nt1gEWGHgpa0nYDcTAy+ls2cjfaZv9UTP9ETP9Lm+0b0//n0oBvwXkhTaAtPvDuCwQAAA3AKLjf6b\nFv9JSftfKoK9TwP4p7YlHU3V6Q2TBYCG7wAoURHwm2CBPr6EFwIfVAiBV8q6At4oBP5X7vaP0qOr\nV3qoV/pM3+oLPdUXnafS1wqD/z+SJNoCASQEAMCW0mT0b+fVLP9rc//bdjxhoBcDvVH2SkUBoNGZ\n8fOldlXU97t9K6AjLvrz++NbLALGCi7AB6sHyEb/eh62u8+UuwCf6Jk+11N99Mt/DLUAvxBtgQAZ\nCACAG2Z9o/+4BsAl+Vsq6wLL+0vp/n9XACgpH+H70b+1AHo63tKf1QgrEacH/BwBr5QVFL7SbFHg\ni+AChDTAcz3JigL1i6n0vyoIAWoBABAAADdJs+DfdPRv5yam/jX731f9231Xs2QFgFIR/K3yfybo\nx/b/mrFA77e9C/BBzgXwdQDvVXIBPtLLkgvwyR8+DXH/f5N0GlYL3OgHAdhyEAAAW0tqid+644lz\nrP1PKgf+RAFgOH0iTywGFsbXHkTvtimxCJhxASwV8ErSO+VOwKOrV3qsF7kL8ETPpF8NQxrgTyXp\nT2gLhIMGAQBwQyw3+o+PVdn/RrT0rxH3/8dkMwAaM21+UlQMWFP113B6seP5p8yQcgGqOgK8C5BP\nEfzkmfSVggj4RCIVAIcMAgBgK/F9e3XnJKr/41N8/7/f72cAdFjwLzoBGo7+U6mF+O3EHQgNSaUF\nBgpzCuQdASYCXFeAuQAP9SprC3wqfTUMdQC/kmgLhEMGAQBwA6w2+o/t/TjYR8+L8//SbAFgAmsB\nbMJk3jD/SGVB4NIAJ9Gpi6QDvO+QCwFbjMiLgEwUdF8FF+CxXuSdAT9+8k1wAX4u6UsJFwAOFQQA\nwIZZfMa/qtF/apQfl/pHxNP/GnEHQG92qt8mjOyafr4BT5T/73eL0+MGxqbEaYCxXDGgbw3MbqeT\nSz3Rs9wFeKLnYYrgrxVcgNantAXCQYIAANgK6kb/Vfu9/Z/h8//e2veFgA2s+pVIvYablCgla2x/\nRSljJZYGGAwV2hBNBNjtlXT2eqRzXZRcgE8ePAsOwNcK9QDBCgA4KBAAABtk+fn+49a+4+jeH5Mq\nw+ax0g5ATS6+6eg/TwN4cZFyGZw7cHxUvFvvWVT5GnWUagIsDeCLAs0FeCWd60IP9arsAvz0bVEQ\nePSYtkA4OBAAAFvFIuPfRPD3+f+4bCC+fOKljmqC/ySzFEbqzB70gT+edtgLg24oBDxROg1wFD2O\nFzmOMREwlgoHwO7fKW8LPH890Lm+17m+10d6UbgAP1EQAdnkQLQFwiGBAADYEM1H/4va/9Ks/V/T\nVBd77v6yPeWTANksgE0wETDsdorg7mcY9HgRcBTqAHzwP1bha8QioOpt+1kCrTDww1hF8I/aA1uv\nlaUBXma3F4UL8L8o1AHePxapADgkEAAAW0Oq+C9l/1cIgzj/n7Lkl2DoRvwjl9wvOQFx3t+EQNc9\nzmoSjo/SLkBKBMQOgH2E1F9gYA0MFvyHCosHZGmA+1dvslTAd3kq4JMHz6TPFUTAV5L0S9oC4WBA\nAABsgM2O/u34AjX0cc5/wX58LwKGLtJP2m1N42JDf7PA79MB3bQLENcEeDPBb6c+7ZVtjCVNFIL/\nWCEFkAmCblYLcD8rCCy5AJ8rCIBPJFwAOBQQAABrZrnCv6rRf4qKc9ozJwbikfmSjCIREGYN6BTX\ntRF/Kv9v23fCeSkXwDsAqVG/VP6L1MoemxfAUgLZ7XRymRcEWmfArAvwU7Va/wkXAPYeBADArdGk\n4M9X/acmBPL3KhcA+g6AhozV1sTdYlLHRlGOf+bWVbk4MeECnLq3fKa0GPC3yooHW364+EAhDZCt\nGXD6ZqS+rnSqy7wr4JFelV2AL5RtAOw3CACANbJ8258nJQzikBed02pw2bg3P3IDhqPZ6v6ROhqp\nq5G6JetfCm2AEx0VhYCx/W+3Oyofj1yAM80Gdy8GYkcgrohImhpDzaYC3kut9yENUNxCV0DuAvxY\nIfYzORAcAAgAgFshZWY3sf+rnjvnKZ4V0gBm/Y/cLa8DqEoD2L4s8Ksn6V54fHanSAWYCxCnA46j\nW+pjxtMLSwpOgM0NYELgvXQyGeTB/7Fe6pF3Ab5WEAI/kaQv1Wr9DSIA9hYEAMCaWM/o35MyuueI\nhCrLf04qYHSdnh5wEqUEJmrnboCvA6hMA9xR2QW4pyIl4FIBp04EnEp6oCIVkHII4hTBXFwq4PTN\nSCe60v2SE3Chhw9eSY8VXICvJR09lvTTJlcH2EkQAAA3Tl0pW+rxnPa/mDh1X5Uwr1jR1wf3edSm\nASzYexegrSAGMgdAd6TjO2URkGoDPHM3v993CfRTUxAP3S1LBbTeS31dqa+BTnSVFwM+1gvps2Fw\nAOymL2kLhL1lBTMQAIzVR/9N7X8jKgCMg74PvPO4lnQ3bI7U1YmbZmdmpK/OTCGgTwO0/Mh/6Lav\ns3upyM33iu3jO/lL6ug6CIGr7OkpneJbAk+k+csMD939MKQBTtuXUT3AhT568kIvP/+R9ELBBXhx\nJl3+cM7FAXYTHACAFVm/9S/Vzuy30HTBzfFpACv+Kx2PrH+rAZhJA9jI32/fVTkt4BwA3Qv7zQno\n98JpJ5IeqkgL+MyCOQDWRihl6wzEyyBLhYKwmoCx1B2GNMBJ1hGQpwH0KrgAj+VcgK9pC4S9BAEA\ncKMsav974kly+7OHm1Bh/Ru+3z/k/AsxkKoD8GmAqRX6WaD3dQA+6FsNgBcBWTrg7I704F5ICRwp\n2P4mBux2ku0/UhAMM/Z/qqTBpQM611JHI/U10Kku89bAc13o3qPvpc8k/VBBAJweS/qxWq2/RgTA\nXkEKAGAF1jP6r/saNsj/N2kBrMLs/+uOdPcqt/W9xZ+0/LPHA/XV1TBPA4x6UteC/VghSo9V2P82\nCr9X856OpONx1ibYLab4Hbj+frP8+1mgz0f/VmcgVRc+XkutidTWOE9unKjoDHjUeaU3Hz+WPm9J\n/6wwQdDf/VTS7+r/lgA7BgIA4MZoYt03tfcrzms68Y/L+2sQtifjEDnHWQQ1q/9EV6U6gLYm+WO/\nfaUTtbsTde6M1Bpm1x+qqAEYu9esop19hjfh4XG2hPCHcRHsPccW9H37oa996Cr9K3ctdTNv4yRz\nAU51mXcG9B9daPCD+8EF+ETSfUnff6VW62+m0+m/WkVyAWwNCACAJVlt9F+3vE1d/r/J8YZkgV+S\nNA4xbTjqqN8JAb+vQW712/ZIXXU1ytMAtj3Qiboa6ardV7c3UtccgDtqFvjdcsF6o+AQWP9+5ghU\n4oO/ay9MTo1s3QBDqT2Z6KR9pa6GOtFV3hlwqkudn11o8Ol96fcKaYDnkv7uU0mfzvkgALsDNQAA\nS7B48F+0cG/BRYKWrQt0troVAVq+32oBrA7AbH8rALzSSb5tLkBeC2Dtfr3s/kxFIWDX7bf8f1tF\njUBcIOgfx3MM2Dk9dw0rBJzTBdEdhuWP+xqoryuduLbAc12o8+it9LHKLoB+rFbr19QCwF6AAwBw\n4zQN7vG8dytU/19rtobQmBSbo+uuTjoDjdRRV8M84J+oPWP7dzTSSF1d6UQdjUougO5JZ5NRkfuX\nwui7Lvc/VDFqN6Fg7XtVeX2/36cCvBPgCwSjNEIn9zhGeRogFwEPLvTy0Zn0A4WB/0+UuQAvaz4E\nwO6AAABYkPW1/dUtcOupmvG+YSrgg8qBc6DZdrlrSUdFIaCvA7CA39YkSgF0coFwqVO1NcnFwKTd\n1vBOFm9tOl4vBoyjaNsCdE+FaJHKYiDGnpOqA/CfMw7+11L7JOQWrMzRTxB0qku9fDSUPu4GJ+CT\n7Pb8h2q1/nY6nf5LagFgp0EAAGycRUbua8rvN8WCc5ajn4zbpTqAtiY60VWpE8AXA17pRG1N1NUo\n7wi4tIudSJ3rgVqpwG/ESxVbrt/cAAveXgxUXcfSAPGMhF5YJOiUmhpHefA/1aXuPfpebx5+HByA\n3yvUAvz+sTT9oVqtv55Op/8aEQA7CwIAYAFWH/2vsxOgAWPVBr+cgaReS6Prrvp3rzRSV0ealNIA\n3gkwu99qAToa5oJAChX2bY11eW+iU42qOxWtBqEr6b2KQP9esymAOyqlK0rEqQC/8FA3sZ28RJA4\nnbwocKCTzkBv7kt6lN0+UXADnn+u0CMIsLsgAAA2yoptfWs7P8Ps/75CsLVV/Fw6YHTdVbczmkkD\nWMDvaliy/S3/L0ndSAhYYV6lCDhK3FugdhX7uYhpmgaI72sD/1gdjfL3H8TLJJ8psPPorUZ3z0Lg\n/1SZC3CGCwA7DwIAoCHrnfJ3Xv5/DamAa5UtdAusdecfdTTpDTUcddTuTFzwH2YuQMjxmyC4ymx/\nGz3ngd9TJwJ8zt8e2xoCvg7AUgPz0gC+DsCLgNgJkMLUxQ5b7cBETUfDMEvg3SuNHp2FKQhnXAAK\nAmF3QQAANGA9wX/eqD0O+puZ819SuRDQ6gDaytMA7aNJqRvAuwCW6/d1AB4LoCN18tH1pN3W5b1O\nWQT0JL3Lti3ov1O5A8C7E03qAOy6dg2bhtg7Af7cBNbUeKJBcAI6A725O5VOW0Uq4IfKXICPmBwI\ndhYEAMDzcAP7AAAf8klEQVTGWEcAT12j5rpWbT9UoSfi7b4a1QVYN8CsCzDRpU4lhVy/7etk3QFX\nFe/vpH1VFgHe9recv1Qe+fs6gCOFYsWqSYF8x4AP+nfdY+8Q5J+z+Bk0kWOf19oCO+eXGj08C4WA\nDxVaAz+W9PxrSZd1f0aArQUBADCH9a/2F3/t5gmFFb6mfspfvy2VR/4ly72cBjjKg3/hAkwyYWBC\noCnt9rhIB8TteWMVdv9blUf+vg5AStcCpNoBYycg2zdtS8NupxT8u1kdgKS8qNGETf/ulUb9s8IB\nOM/uf3+cuQC/nk6nf44LADsFAgDgRmg6kp83SdCStQF+LgBfCCiFyvqKNIAVA6ZcgPhxN7f9q0r1\nsyCb1QR0h9mcg77y3wSA3afSAFL6lyuuJ4iDvzkB3ZD/n7TLcwUP3SqIhqUDup2RdHcq3W2FtYkf\nKRQEfqe8FqDV+qvpdPpvEAGwMyAAAGpYfvS/ifx9zTUHUj4Yt/hrQd+3Atr2WLMTAlmhYI0LsMiI\nPyUEbMQ9aR9pcnKlYXei0zeRG2BB34SA7wCwX6x8DYPs3v+S+Tx/nAa4I03vhNG/TWNsN5UuYYsg\nj/IJgjrnlxrdOStG/89VcgGkx43/NgDbAAIA4FCwAG91AFJ55H+qShfACv1OdZlMBVSN+m0JYWOY\nRWcbbVtKoDscqdOWWuYCeCFg79Pm/69bGEjuOXEHQLY2gY3+JzrKJjWqL4joaKQjTdQ+moS/210F\nEfADSa9U6gigLRB2CQQAQAXrG/37x1VfueOK7QXwgXGgEPR80D9WUQcQzwfQc8+fqOQCXL491enZ\nZTIVENPVqNYl6GbJ++I+pAQmJ20NuxN1eyN1rqWWvQc/F8BdFSmAeasDSrPrA/SkYTT6N0cidP93\nZwSLfdaORur0hhpYGuCugqCw1sDfH0vTHysULwDsBggAgATrL/yL6Uf3TflQ/zzf+x8H/TglYGkA\nP/K3XwTnAnTmpALi0f+pLkvHfVC17eI+jMJP9U7t9rgQAuYIWErCOhvmFQLKnePSB9M7YeQ/7HZ0\n1e7nr+tXOZy9TPkFup2R1BtJvW74O55mN3MEvn8sXADYJRAAADvJILuf4xbEQT9OA/h0gC8GdLUA\nI6kyFeDxQiAWAcYoTwGELHtfg2y64aNszYGjGSGgO2HhnlZcCGiiwOPfUi9U+0uzwf9SdzMR0J1Z\n9tiwz2cuQFsTdXpDje52i6B/V8EBeCjpoqgFaLX+w3Q6/beIANhqEAAAEauN/jc4eY8Gza4fr7xX\nlwY4VZEOkGZrAdSSVBQE+q4A5U9J5/87Ljr7YGrHvBtggbiTrSHQ1Uid9khXJ2F72A3P84JAUnqF\nwQwL+lLI+Q8zy9+C/6VO89f2tQC22FHq87SPJkV7oRcB5yq5AKQCYBdAAABslCb5/xWZanaRnImK\nOgApnQbwb6/KBXCLBEnhGu1OdZtfFSYGbHGhE7WzNfg6OsldgBCMzREYqpMVGh5J7TBn/+QkiIar\nE6k9qX4f1uI3zEf3R9nkvp1S8L/Uqa7Uz/P/VQWBMyLHlhy+o3JBYO4CvKAtELYeBACAY/O5/w2S\nMgj8LIA+DSDVuwC+ODBKBUiSHyCn2/3iXH95tG/Pm2RCwFIKKSEQlhweZ9MQn+QL95goSGHFfWbp\nW+AfZSN9C/72+lc6yesBRllBYIpOb6jB0VTqtYq/k3UnPFSYF+ASFwB2AwQAwNpY1f6ve35D+18q\nFwJKRd2gpQFs2xcD2sv7kb8XBFkqwCoPJEkd6ahm0h+PBdSRujNugAmBvq5KwX+gvlt4aJgHfxME\nhp/BTypP6OOFgAV+vz10RYCWCjCxMnLbJXojqd0tgv+5QuC3CYLemQvwGhcAthoEAEDGzY3+N1kn\nkOEDeTfa510A+wUwFyBuC4xSAboO9QC2ZLBfCKhq1G+5fSmd+49H/WH63WG2zkBZDHiL3qr0B6kV\nCDPsfF/gFzsAVn9go/9hFPiTXQIW/K0e4FzShYKrclfOBWByINheEAAA2lTwv4FA3wQ/6o8fx7UA\nUnUq4HtJ91saXZyqc36Zzw9wpRONEyNlG/XbCF9SHmTteEex5V84AnY8LDlcBP8i8PeLdEAFPujb\n6w+zOf7CwkUnuSMwLImA2TkBZrA6AJtr4I6cADiWZC4AbYGwnSAAAG6NVQWCawW0QkA/ipdmuwGk\n8uJA9tieZykBP/I3rB6gNysC4pz5SJ28FdB3A5gb4C3/OiFwpRN1Ncyv7wWBpCwdMH8WQksL+GI/\nnw6wQsBCBBSioHamQBNJXYVaCysMPJI0LmoBaAuEbQQBAAfPekb/G1zRL58H19OwJsBEgVQuBjR6\n0WN7map6gPx5ZRHQ6Q3zwsBy8B+VRv/F6L4IyBb8/bYFdhMDUtlJGESB34sMu275z2BdAek6gFFJ\nBHRz29/PW1CJzajYU9ESeCrp+8IFoCAQthEEAMCts+jUvx+WeI6KYkBpthbA0gOpVEBP6SXvnQiw\nNsFJJ4zs4+AvlYN9MR/AKGn3Wx1ALAa87e8FwFVNHYBUrkeI0wEDnTgRkEoJzM4JkP9d7O/l1xw4\nVhABl5LGZwoiABcAtg8EABw0u9f2VzPyr0oDpKx8aXYWPZ8KiJkjAtRzufhOCM7jUvDvzAR7qVyc\n18nG37ErEI4V4sDb/vF0vSn86N2CfnhPsQNQtALGtQC2NmASXwdgaQBLDVyeSXog6Uy4ALBtIADg\nYFlf8E9FzHUVAM6z+htOCeyJR/227VMBfoIgf34sAgaS7ofugIGkybit0VG3VBdQDv6zQsBG/Knj\nUrlzwO+TVKoHqMK3BY6ULgb0nQBXmSMQj/4nahcTIqXIFhxST8EBeC/pnaSpiYC3tAXCVoEAALh1\nzNJvEuz7NY8zmrgARlwQ6EVBv+I6XgRYd0C/JY1DKO5kKwh2ekNNOu3cDfDtgJLy3v44NRC7AhO1\nSzn/uBAwrgeY/XPUFwOGxx1d6aQkDIbunFKXw3VFSsCInYA8DUAtAGwXCAA4SHbP+o+J6wCyx1OF\nOXuakppIr5fYlxIBY/c4Fxddja47Um+UcAOGskrBVJD3QkBSUgzYc6V04F+kGNAvAjRbEFjUApj9\nPxx1NBl7IRD9fawQ0JwAW1fheyn8W+ECwHaBAADYOlJFfr4ToKET0NQFSAV8vy/uDIivZd0CFouz\nCYOauAFW2DfP+i93APRLRYSeecWAPuhLmgn8sxMDdTS2moDrbkgBjFvpv2OMCYIjSeMQ/EMtwAe1\nWn85nU7/HSIAbhUEABwc6x393/bqf4lzrBgwnscm1dtvzBMB8XTBdjwWGflbylICkRuQEgKSdwE6\npZF+Vf+/CQCbJ6ApvgMgPC6K++Lgf6nT7D0mRv/+s3q8/W9OQF/ZxEBFLYASSyUD3DQIAICtxwd5\nv12RBjDiWoBUkI9bA2PifUPNThbkxca1QgHcOzUWAh0NNY5cAEkzYiDsKxcB+mNN8DUA/rk+8Nu2\nBf+rUb8Y/V93wmecqCx6YkwItLPtlqRpXyHwmwtAWyDcLggAOChuJvcfj9o38TWLRUH0unEtgAVs\nqewCxEXt82oC+tF2XBfQV5gT36rh7XUsLVAhBCQlXYFiJsByHYAf+c8rAvSUCwKLAsC4IDAO/pNx\nOwR/s//j4F+VErA0QF4M+FaFC8DPL9wu/A+Eg2H3C/9SMwIu4AIY81IBq4gAu2bsEJgb0FNSCEgq\nuQLhEtVOgO035q0JkHqOH/3bMQv8E7U1HHXy4J/n/m30b58x5QT0FBwQf99X1j1xptANgAsAtw8C\nAODWqMvx17UGVu2TKl2AVCrATxPsWVQESOWAH9cGVAmBd0VqoNMbavDuJHcFJM04A1Lo+7cliH3Q\nb1oL4CcFMjEwdumFXBTEwd+sf1/42KQQsKsgAkppgAfZHwwXAG4X/vfBQbA7q/1ViQLbn+oG8KP+\naGKguCAwTgX4p3jmiQBfC5DqEpB73dT+/GO6GgEpdwUkzTgDRuwQGJVT9kaUZyiMagHiwC+F4P8u\ns/5txG9/G7/tscAf3nBxPz5W+Lc5Ei4A3DYIAICdZJ4zUJMKiEVAW9Wj2VRhYLzfXjIlCswGb7v9\ncfFgyhWQZpwBI3YIjKXqAEZZUaCtZWDiw0b9Zvv74B8/ts8+UCEGYjOiNIFS392CGKAtEG4DBADs\nPbuR+69b4Mcid13Qj12ARCogrgdYRgT4l4nFQVwb8E7ltIAJgapjeRYjLQYk5YJAUkkU2LGm+Cl9\nS0FfKgK/lA7+ce4//rtVdQfkdQAW/E+zJ39Q2kYA2CwIAICtpUoU2FA6JQwqXICUCLDRuJQWAcea\nTRfEWKrApwHGKtyAONhL5dSAHbNrpcSAJI2zQj0/Da9bgMgLg6aU5vW366YCvz224D+IHg/d/Txa\nkqaWBvAuwCkuANw4CADYa7Zj9J+q3q+jyQRAKXzhoDTjDKSKAmMRIC2WDvCjeKkc+FNCwP8ZfG1C\nnRjw76nv4uPYFfTNm5+/Ch/wjarA7495MeDfX2r031eiW8D+bfqyWoDNTioFMEtrOt2C30eADbD5\n4F/1g53af1Rz3B4fV5xzHO3zi9Cn9h9X7MtoqRx4fb7eAnnbPdXWuO+5844aHvNv1eJ13WumPl4v\nuvfn+vOXxYK3D/r+3gJ9nAYYqBj9DxRW/7tWsPmHCuLFb1/IraHwQdILhZbAy2z7raSXwgWAmwIH\nAGApVok6VSP8ujqAFHEqIPUac6YKrnICfKogZW3HbkD81qscgVgHVY36fYD3Ffep10/V/tVNDpg6\nPxX0pbI4qAv+w4rn12JCbazCESAVADcHAgD2ku2w/ldhwaBeOh7PIRC1BkppESDNBtl5ve6LCgGp\nEAJ+X0oM2H0sCOxcqSxcqt5TFfHzUkHfn1cX/H1nRSyY4r9hXgdg/1aXCn80LwgANg8CAGDrSTkD\n8bwAsThYQATEl7O8vi/cW1QIpI6nHIA6MeBbBv114v56fz25403xn827A3Hgt3Orgv+1ymsr1K0V\nUOJIxb/ZB0kDXAC4ERAAsHfs3ui/qn+/6SyA8bEFRUDcHeBHs8aiboAF9JQjUBX068SAlC5Y9Ndd\nNPB7riu248Bv+/xjzzw3YAY/6h8IFwBuEgQAwE7QxAVIHVtCBEjllEBTNyDuCEi9fX+dpmJAKi9a\nFIsCe21PkzRAaoTu4659zmF0zI/87Trx6F9q4Ab4fw9cALh5EACwV+zO6L/JOgDxean1AapSAanr\nxSJAmmkRjOsClnEDqoJ80/O8GJDKA+IqUWDnGaukAaTZoG/n+P0+RfBBs+kAY+4Ehal/E1wA2Dy0\nAcLecLPBv64LoOpYHAWr2gGl+S2Bfr8/XtUa6J9b0yIoFXMFpNry4pa8+BwLyvFsgXZ/VHFe3bn+\nfP9aSpwXE4sEo8qaj2NulQvgg3wq+Nu+uH3wMns8tScOFNr/XmcXeJ09fpsdu6QtEDYGDgDAVjLP\nBfD7/fG6WQJTTkA8jbAWdwPsrTSlbrQvpQv/mjgA0qw4aDqIrnIA4mvEgd/OTf1tYmrfi80TXE4D\nSH21Wn8xnU7/PSIA1g4CAPaC7bL+m87kF5/XpMBPqq8HkNIiQKquCagQAVK5S8AuITUXAlXFgP45\nKTHgz48FgT3XB1T7Jasa2cdioa44Lw7U8RwBPvDb+al0QDz6r8T+PY7cY0sDLDsrJMB8EAAAW0uT\niYGq1gFI1QTMKww0KtwAqbkQiHv3PamCQHuurxeoEgT2uk1WKjSaOAGpIO33pQK/nROnAxZxRErE\nogwXADYHAgB2nu0a/ddRV61f9XiRVEC83VQESLOBRyoVCEqzaQF7mpQOnk3SA1UaJ3YHpFnrP3qb\nOYt0BNSNzH0VfyolEI/67by62QPz/H8d/t8HFwA2B0WAsNPcXvCf94PctBAwdW78OI5wTYoFU4WB\nfn9dcaDfTkTnuEjQn+5fqm4Of/9S3p5PFQXG2/HreKqK/hYhTg/ENQDx9rzgb9eYEQA2wh9k26/d\niVYcWN5PQSCsExwAgBtlVRcgPr6IE+D3V9UFSLVugDRbH5ByBKR6V8CnCPxbrasD8M+Pr2f4YsFF\nqUoTVE0SFKcE4uDvr7tUSsD+/lYgCLBecABgZ7ld67+JJbuKC5DaVzXarzs2zwmIt+vcgPh1MvyY\nNDXCT7kC8Xbd6n7xqL4uZVBH1XLGdcTHU2mBeWsH2Dkz9n/sAIS2v9l91hKICwDrBQEAO8t+CYDU\n+YukAuLjVQE8fv15tv+KQiB+WhMxUPd8I2X1L7KQYhPimf2qFvmpW0bYzrNzKu1/OzElAGbFAQWB\nsA5IAcBOsjuFfymqpsaLWSQVEB+P0wFyj1MpAb/tbf84LaDEORn+X8TPfuft76oUQaojwEh1FKTq\n4urmB1iUVJtgbOPHgd9v1wb/FBb0bUrg1IvGxwBWAwcAdpLbFwBNf4hXdQFS+xZxAuLjVcPxuvOW\ndASMJs5A6u3My/fH14tZNlZW5evjKX3nBX7/nKT1byd6N8AuknIBykWCuACwKjgAsHPcfvCXVm/N\nqnIBUtddxQmIj/tRfFwxV3XeIo6AneeI2wj92/LEH7vpcrpV+f+qZYOrmDdnf/x+UqsGSg2Cf3yi\np27CAloCYb3gAMBOsR3B31ilDsBoWhCY2p8adS9SONi0NiB+3NRhmJOUj8evcaBO/RlSf64mBYB1\n5zcVGvF5qaAvlYVEMvjHo3/btkI//ziuAyiKAXEBYFUQALBT7J4AmHfeIqmA1P5FRUDq+CaEwLxz\nI1JhLDVyXzSjsk5Sg/NGgV9KB3/b7/fVCQC/HxEAq0MKAHaG7Qr+62KRVEBqf2o6vCYpAblz5hUJ\nVp0b2/5xtV9d83tFmsCTsuSrLnkTrnjqteP3mPwfWhf8m77Qcc35AMuBAwA7w/YJgEWizjKpgLrn\npfbPG5E3OaeJ775qu2ITF6OCqrFu0zz/KlTVCFT+r0zl++PgH+8fR/s+aNYtmHUKcAFgGXAAYCfY\nvuAv3UxBVlMnQJo/0rdzpHK7oD8n9ZxUsWDdc6pew85pMi1ehSio+l8wr4BvXTT6X1hV6Fc18vfB\nP0X8N7PHTdc6BkjzB7f9BgAOg3lBr+7HvOq5qf0pm9i3mFWdF58zSOwba7byzZ8TP/4Q3aquO2jw\n3Dn29/SGbpVUfU7/GaVqcRDTRCQZ/S0VyLDt4ADA1rM/P27zHIO6CYIWdQKk2VF00/oBReelagSM\nea6A31dXM+DPrdpfJQLWPQVgE5rm7/3jupRAPIFA1TX83+dIzR0VgFmoAYCtZjeC/zprAaR6XV73\n/CY5/0XPbTIxkbRaO2PVa9ddo+nxpq/jWaTYbl7Ajq+XcgVS/YSpLoGUu1CuGaAWABYBAQBbzf4J\ngKbnr1MESKsLgapzm4qBqnPr9tcF63WLrjrmjbCrjtfZ/U2Df93jVNugEAHQGFIAsLXsRvCXFi8G\nbHL+vHSAKq5RdawuLRCfXzWrX116wO+P6xniNEHdNfyx1Eg8VViYou76yzLvOlX1F1WPq2YSWhbS\nAbAYOACwteyOAJCWG2Wu6gTMu0bdsUUt90XPX2bGnmU/i+em6gHq0gRN0gJVwb+uWLBp6yAuADQD\nBwC2kt0K/tJyLYGrOgF2DVVcp+7YIo6APz9+TpNRfaql0PCfrar4r+pY6vq3MVnOItX88edPBf+m\n1/b4dQJwAaAZOACwleyeAJCWzzU3fd4qbsC848vk25fN0S87KVLT56/7eYsE1Kpz6wJ//LiqdqDO\nAZjdjwsA88ABgK1jN4O/tJwLsMjzVnED5h2vcgTqnlezAmBlXn/eMSk9J8I8pyCmzhFZB02uNW/x\ngNS+ecG/DlwAWAwcANgqdjf4G6tUnK/LCWh6vXnHV6nCn5eLX2d73zaMY+bNyrdou2B8vG4OAf/6\nuADQHAQAbBW7LwCkmxEBUrPA1+R6mw7mTQrzmn7uVVv6lhULi067u+gsf4sG/3g7LQAkRABUsw3S\nGUDSvgR/aflUwKLPnZcSsOtpzjXnnVNn9fvnG3WpgqbXqHo/i7T/pdjE/PlN7PYmcwUscr15kAaA\n+eAAwNawPwLAuCknQFpfWmCR89Y5sl+0fW/TizCtyjrmDIjPW3ROAVwAqAcHALaC/Qv+q9Jk5O6J\nV+ybd915125y3jxnIL5O02vNu2bqujE3JRAWGWUvOnfAqi2NuABQDw4A3Dr7HfzXEYiWucYi2n4T\n0+ouMqJf5vPdxgJAi9IkgDcN/HWjf2n+xEK4ADALDgDARlmlHsBfQwtep6kj4K/f5DWantsk75+6\nZtP30SS4bsOsgCnqRuWbqgkAmAUHAG6V/R79e9ZlSS97nWW0/qYX3Fk2QG97/j+mSQBvurxwal/z\nSYZwAcCDAIBb5XAEgLTewLWtQmCV50jrGbXflkBYdKQ+zzmout7yAkBCBEABKQC4NQ4r+EvLWfnr\nvpYPFk2//k0L+VZ9jlQdFBcRBttsmS9bF9DkWBUUA0KaP7jtNwBweKx7Stplrzd2t2Vec9E58uPb\nInyouW0ri77XeX+X1LHF5zU4POENVeAAwK3Aj9A6igPj62mFay7jDPjXNRZ5/apgt+hnWEYErJpq\nWJfwWGUSIYDVoAYAbhyCf8wmctbruuaqY4RtqHvYNtaxuqBUPfpvtqwwtQCAAIAbBwGQYpPBbVvE\ngLELn3UdrDJyn/fcOuu/afEgIuDQIQUANwrBv4p1Fghu6tpx0Fn252Nd1v8i194Fmr73TaxnAIcI\nRYAAW8UmA9iyBXhVjLV8IWGKVKHgLgf0Jiz6OZdZdrgaBPlhgwMANwY/Nk3ZpBsQv8Y6X2ddDkHM\ntsz7vwqrChlG/bB+qAGAGwMBsCw3GeA2/VrbMubYRAfGJlgk8M97H9XHqQU4TLbl2wh7DsF/FW7C\nEYhfy1j3a1YFtJv+Kdrm1MIyo/1t/jywreAAwMYh+K+b27K8b+t1D2GcsorF3zT415+HC3B4UAQI\nsHPcVoHcbRXoxcWG6y4+vGnW+RnW92+AUD88DkFawy3Cj8qmucn0QN3rG7fxPpYNoOv++btJMYLl\nD6tDCgA2CgLgNti2qvhtez+7zLKBv/nzSAUcDjgAsDEI/rfFJlr8VmGTE/8cCoz4Yf3gAMBGIPhv\nI7sUcHfpvW6Cda8YuRi4AIcBDgDAwbAN+fqm7MPkP01hdA+3Aw4ArB1G/7vIPgXUFLf5+W4zwC//\n2rgA+w8CANYOAmAf2HdBcAisLjwQAfsNKQBYKwT/fSEVPBAFuwNpBZgPEwEBQEMObaW+XWS9/y4I\n+v0GBwDWBj8WhwhOwXaAGIPFoQYA1gLBH+aDMFg/NxP4qQXYT3AAAOCGqAtWiINmMNKH9YEDACvD\n6B9uhkMUCdsT8HEB9g8EAKwMAgC2j10RC9sT4JuACNgvSAHAShD8YTvZrcAKcBvQBggAAI1A8O8X\nCABYGn4MAAB2FwQALAXBH+Aw4bu/PyAAAABgIRAB+wECABaGLz8AwO6DAAAAgIVhILD7IABgIfjS\nAwDsBwgAAABYCgYEuw0CABrDlx0AYH9AAEAjCP4AkILfht0FAQAAACuBCNhNEAAwF77cAAD7BwIA\nAABWhoHC7oEAgFr4UgMA7CcIAKiE4A8Ai8Bvxm6BAAAAADhAEACQBCUPAMvAb8fugAAAAIC1ggjY\nDRAAMANfXgCA/QcBAAAAa4eBxPaDAIASfGkBAA4DBADkEPwBYJ3wm7LdIAAAAAAOEAQASEKpA8Bm\n4Ldle0EAAADARkEEbCcIAODLCQBwgCAAAABg4zDQ2D4QAAcOX0oAgMMEAXDAEPwB4CbhN2e7QAAA\nAMCNgQjYHhAABwpfQgCAwwYBAAAANwoDkO0AAXCA8OUDAAAEAAAA3DgMRG4fBMCBwZcOAAAkBMBB\nQfAHgG2C36TbBQEAAAC3BiLg9kAAHAh8yQAAwIMAAACAW4UByu2AADgA+HIBAEAMAmDPIfgDwC7A\nb9XNgwAAAAA4QBAAewyKGgB2CX6zbhYEAAAAbA2IgJsDAbCn8CUCAIA6EAAAALBVMIC5GRAAewhf\nHgAAmAcCYM8g+APAPsBv2eZBAAAAwFaCCNgsCIA9gi8LAAA0BQEAAABbCwObzYEA2BP4kgAAwCIg\nAAAAYKthgLMZEAB7AF8OAABYFATAjkPwB4BDgN+69YMAAACAnQARsF4QADsMXwYAAFgWBAAAAOwM\nDHzWBwJgR+FLAAAAq4AAAACAnYIB0HpAAOwg/OcHAIBVQQDsGAR/AAB+C9cBAgAAAHYSRMBqIAB2\nCP6zAwDAukAAAADAzsLAaHkQADsC/8kBAGCdIAB2AII/AEA1/EYuBwIAAADgAEEAbDkoWwCA+fBb\nuTgIAAAA2AsQAYuBANhi+M8MAACbAgEAAAB7AwOn5iAAthT+EwMAwCZBAGwhBH8AgOXhN7QZCAAA\nANg7EAHzQQBsGfynBQCAmwABAAAAewkDqnoQAFsE/1kBAOCmQAAAAMDewsCqGgTAlsB/UgAAuEkQ\nAFsAwR8AYHPwG5sGAQAAAHsPImAWBMAtw39KAAC4DRAAAABwEDDgKoMAuEX4zwgAALcFAuCWIPgD\nANw8/PYWIAAAAAAOEATALYACBQC4PfgNDiAAAADg4EAEIABuHP7TAQDANoAAAACAg+TQB2QIgBvk\n0P+zAQDA9tCaTolJAAAAhwYOAAAAwAGCAAAAADhAEAAAAAAHCAIAAADgAEEAAAAAHCAIAAAAgAME\nAQAAAHCAIAAAAAAOEAQAAADAAYIAAAAAOEAQAAAAAAcIAgAAAOAAQQAAAAAcIAgAAACAAwQBAAAA\ncIAgAAAAAA4QBAAAAMABggAAAAA4QBAAAAAABwgCAAAA4ABBAAAAABwgCAAAAIADBAEAAABwgCAA\nAAAADhAEAAAAwAGCAAAAADhAEAAAAAAHCAIAAADgAEEAAAAAHCD/P2JAUvKBMevUAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110ce0438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_pdf_contours(Dirichlet([1.99, 3.99, 10.99]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6-faced die with repeated labels\n",
    "\n",
    "Consider a die where the numbers on each face are labeled, possibly with repetitions, from the set $1\\dots 6$. \n",
    "A 'normal' die has labels $1,2,3,4,5,6$ but we allow other labelings, for example as $1,1,3,5,5,5$ or $1,1,1,1,1,6$.\n",
    "\n",
    "Can we construct a method to find how the die has been labeled from a sequence of outcomes?  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Does my data have a single cluster or are there two clusters?\n",
    "\n",
    "We observe a dataset of $N$ points and want to decide if there are one or two clusters. For example, the below dataset, when visualized seems to suggest two clusters; however the separation is not very clear; perhaps a single component might have been also sufficient. How can we derive a procedure that leads to a resonable answer in this ambigious situation?\n",
    "\n",
    "<img src=\"clusters.png\" width='180' align='center'>\n",
    "\n",
    "One principled approach is based on Bayesian model selection. Our approach will be describing two alternative generative models for data. Each generative model will reflect our assumption what it means to have clusters. In other words, we should describe two different procedures: how to generate a dataset with a single cluster, or a datset that has two clusters. Once we have a description of each generative procedure, we may hope to convert our qualitative question (how many clusters?) into a well defined computational procedure. \n",
    "\n",
    "Each generative procedure will be a different probability model and we will compute the marginal posterior distribution conditioned on an observed dataset.\n",
    "\n",
    "The single cluster model will have a single cluster center, denoed as $\\mu$. Once $\\mu$ is generated, each observation is generated by a Gaussian distribution with variance $R$, centered around $\\mu$.\n",
    "\n",
    "Model $M =1$: Single Cluster\n",
    "\\begin{eqnarray}\n",
    "\\mu & \\sim & {\\mathcal N}(\\mu; 0, P) \\\\\n",
    "x_i | \\mu & \\sim & {\\mathcal N}(x; \\mu, R)\n",
    "\\end{eqnarray}\n",
    "\n",
    "The parameter $P$ denotes a natural range for the mean, $R$ denotes the variance of data, the amount of spread around the mean. To start simple, we will assume that these parameters are known; we will able to relax this assumption later easily.\n",
    "\n",
    "Below, we show an example where each data point is a scalar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-ae8470b7f8a0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mN\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0mmu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m \u001b[0mx\u001b[0m  \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mR\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "# Parameters\n",
    "P = 100\n",
    "R = 10\n",
    "\n",
    "# Number of datapoints\n",
    "N = 5\n",
    "\n",
    "mu = np.random.normal(0, np.sqrt(P))\n",
    "x  = np.random.normal(mu, np.sqrt(R), size=(N))\n",
    "\n",
    "plt.figure(figsize=(10,1))\n",
    "plt.plot(mu, 0, 'r.')\n",
    "plt.plot(x, np.zeros_like(x), 'x')\n",
    "ax = plt.gca()\n",
    "ax.set_xlim(3*np.sqrt(P)*np.array([-1,1]))\n",
    "ax.set_ylim([-0.1,0.1])\n",
    "ax.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model $M = 2$: Two Clusters\n",
    "\\begin{eqnarray}\n",
    "\\mu_0 & \\sim & {\\mathcal{N}}(\\mu; 0, P) \\\\\n",
    "\\mu_1 & \\sim & {\\mathcal{N}}(\\mu; 0, P) \\\\\n",
    "c_i & \\sim & {\\mathcal{BE}}(r; 0.5) \\\\\n",
    "x_i | \\mu_0, \\mu_1, c_i & \\sim & {\\mathcal N}(x; \\mu_0, R)^{1-c_i} {\\mathcal N}(x; \\mu_1, R)^{c_i}\n",
    "\\end{eqnarray}\n",
    "\n",
    "The parameter $P$ denotes a natural range for both means, $R$ denotes the variance of data in each cluster. The variables $r_i$ are the indicators that show the assignment of each datapoint to one of the clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1QAAAB8CAYAAABudvGbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAABJFJREFUeJzt3T1uE1EUhuEzbCTZSdhDOppI2UI6CpAooGIDFJGgYB1k\nJ2QlQ5FRsKx4xv58x5mMn0eK8jPo3jvNsV5sJ13f9wUAAMDh3r32AQAAAN4qQQUAABASVAAAACFB\nBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUA\nABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQ\nElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJU\nAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAA\nACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBRyk6+qm6+pyx7XLrqubY9fd3mNz\n3Zf2aHGmue4LWJc9ZsWPFrNkYp+7rqu7Y/doxfzk3Akq4FAPVXW//eA5fH8/XD923eevN9cd2aPF\nmea6L2BdpmbFr4nr+86SsX2uq+p6QfPK/OSsdX3fv/YZgDdm40Hytu/r7/b3LdYdfvR7+Pxh+Lxz\njxZnmuu+gHWZmhWtZsnYOsM/Wcy8Mj85Z56hAg42PDje1tP/SF5VowfNzXWr6mLj0sXUHi3ONNd9\nAesyNStazZKxdZY2r5Z2Hjglz1ABseFB809Vve/7di/p2Fx3+NHee7Q401z3BazL1KxoNUvG1lna\nvFraeeAUPEMFRIaXc3yup+j5vOsNyUeu+2342GuPFmea676AdZmaFa1mydg6S5tXSzsPnIqgAg62\n9dr4h/r/Mo/LVutW1ePGpcepPVqcaa77AtZlala0miVj6yxtXi3tPHBKXvIHHGTXG42PfQPyrjdb\n7/p6n70POdNc9wWsyx6z4lNVfRm5vtcsmdjn+Rf2LGFemZ+cO0EFHKR7+nsiDyO/ae+q7+vnMetu\n77G57kt7tDjTXPcFrMses+JjVX09dpZM7HNXVdX39f2YPVoxPzl3ggoAACDkPVQAAAAhQQUAABAS\nVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQA\nAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAA\nIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFB\nBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUA\nABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQ\nElQAAAAhQQUAABASVAAAACFBBQAAEBJUAAAAIUEFAAAQElQAAAAhQQUAABD6B5FAoeZRGmAqAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11109cac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters\n",
    "P = 100\n",
    "R = 2\n",
    "\n",
    "# Number of datapoints\n",
    "N = 10\n",
    "\n",
    "# Number of clusters\n",
    "M = 2\n",
    "\n",
    "mu = np.random.normal(0, np.sqrt(P), size=(M))\n",
    "c = np.random.binomial(1, 0.5, size=N)\n",
    "x = np.zeros(N)\n",
    "for i in range(N):\n",
    "    x[i]  = np.random.normal(mu[c[i]], np.sqrt(R))\n",
    "\n",
    "plt.figure(figsize=(10,1))\n",
    "#plt.plot(mu, np.zeros_like(mu), 'r.')\n",
    "plt.plot(x, np.zeros_like(x), 'x')\n",
    "ax = plt.gca()\n",
    "ax.set_xlim(3*np.sqrt(P)*np.array([-1,1]))\n",
    "ax.set_ylim([-0.1,0.1])\n",
    "ax.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extension\n",
    "We dont know the component variances\n",
    "\n",
    "\n",
    "#### Combining into a single model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Capture\n",
    "\n",
    "Recapture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Change point\n",
    "Coin switch\n",
    "\n",
    "Coal Mining Data\n",
    " Single Change Point\n",
    " Multiple Change Point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bivariate Gaussian model selection\n",
    "\n",
    "\n",
    "Suppose we are given a dataset $X = \\{x_1, x_2, \\dots, x_N \\}$ where $x_n \\in \\mathbb{R}^K$ for $n=1 \\dots N$ and consider two competing models:\n",
    "\n",
    "- ### Model $m=1$:\n",
    "$\\newcommand{\\diag}{\\text{diag}}$\n",
    "\n",
    "- Observation\n",
    "\\begin{eqnarray}\n",
    "\\text{for}\\; n=1\\dots N&& \\\\\n",
    "x_n| s_{1:K} & \\sim &  \\mathcal{N}\\left(x; 0, \\diag\\{s_1, \\dots, s_K\\}\\right) = \\mathcal{N}\\left(x; 0, \\left(\\begin{array}{ccc} s_1 & 0 & 0\\\\0 & \\ddots & 0 \\\\ 0 & \\dots & s_K \\end{array} \\right) \\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "$$\n",
    "p(x_n| s_{1:K})  = \n",
    "\\prod_{k=1}^K \\mathcal{N}(x_{k,n}; 0, s_k) \n",
    "$$\n",
    "\n",
    "$$\n",
    "p(X|s_{1:K} ) = \\prod_{n=1}^N p(x_n| s_{1:K})  = \n",
    "\\prod_{n=1}^N \\prod_{k=1}^K \\mathcal{N}(x_{k,n}; 0, s_k) \n",
    "$$\n",
    "\n",
    "\n",
    "- Prior\n",
    "\\begin{eqnarray}\n",
    "\\text{for}\\; k=1\\dots K&& \\\\\n",
    "s_k & \\sim & \\mathcal{IG}(s_k; \\alpha, \\beta) \n",
    "\\end{eqnarray}\n",
    "\n",
    "\n",
    "\n",
    "- ### Model $m=2$:\n",
    "\n",
    "-- Observation\n",
    "\n",
    "\\begin{eqnarray}\n",
    "x_n \\sim \\mathcal{N}(x_n; 0, \\Sigma)=\\left|{ 2\\pi \\Sigma } \\right|^{-1/2} \\exp\\left(-\\frac12 {x_n}^\\top {\\Sigma}^{-1} {x_n} \\right)=\\exp\\left( -\\frac{1}{2}\\trace {\\Sigma}^{-1} {x_n}{x_n}^\\top -\\frac{1}{2}\\log \\left|2{\\pi}\\Sigma\\right|\\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "$$\n",
    "{\\cal IW}(\\Sigma; 2a, 2B)  =  \\exp( - (a + (k+1)/2) \\log |\\Sigma| - \\trace B\\Sigma^{-1}  - \\log\\Gamma_k(a) + a\\log |B|) \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\mathcal{N}(x_n; 0, \\Sigma)=\\left|{ 2\\pi \\Sigma } \\right|^{-1/2} \\exp\\left(-\\frac12 {x_n}^\\top {\\Sigma}^{-1} {x_n} \\right)=\\exp\\left(  -\\frac{1}{2}\\trace {\\Sigma}^{-1} {x_n}{x_n}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}\\Sigma\\right|\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "\\mathcal{N}(x_n; 0, \\Sigma)=\\left|{ 2\\pi \\Sigma } \\right|^{-1/2} \\exp\\left(-\\frac12 {x_n}^\\top {\\Sigma}^{-1} {x_n} \\right)=\\exp\\left(  -\\frac{1}{2}\\trace {\\Sigma}^{-1} {x_n}{x_n}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}\\Sigma\\right|\\right)"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table width=100%><tr><td align=\"center\">\\begin{eqnarray}\\mathcal{N}(x_n; 0, \\Sigma)\\end{eqnarray}</td><td>\\mathcal{N}(x_n; 0, \\Sigma)</td></tr><tr><td align=\"center\">\\begin{eqnarray}\\left|{ 2\\pi \\Sigma } \\right|^{-1/2} \\exp\\left(-\\frac12 {x_n}^\\top {\\Sigma}^{-1} {x_n} \\right)\\end{eqnarray}</td><td>\\left|{ 2\\pi \\Sigma } \\right|^{-1/2} \\exp\\left(-\\frac12 {x_n}^\\top {\\Sigma}^{-1} {x_n} \\right)</td></tr><tr><td align=\"center\">\\begin{eqnarray}\\exp\\left(  -\\frac{1}{2}\\trace {\\Sigma}^{-1} {x_n}{x_n}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}\\Sigma\\right|\\right)\\end{eqnarray}</td><td>\\exp\\left(  -\\frac{1}{2}\\trace {\\Sigma}^{-1} {x_n}{x_n}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}\\Sigma\\right|\\right)</td></tr><tr><td align=\"center\">\\begin{eqnarray}  -\\frac{1}{2}\\trace {\\Sigma}^{-1} {x_n}{x_n}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}\\Sigma\\right|\\end{eqnarray}</td><td>  -\\frac{1}{2}\\trace {\\Sigma}^{-1} {x_n}{x_n}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}\\Sigma\\right|</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# %load template_equations.py\n",
    "from IPython.display import display, Math, Latex, HTML\n",
    "import notes_utilities as nut\n",
    "from importlib import reload\n",
    "reload(nut)\n",
    "Latex('$\\DeclareMathOperator{\\trace}{Tr}$')\n",
    "\n",
    "#L = nut.pdf2latex_gauss(x=r's', m=r'\\mu',v=r'v')\n",
    "#L = nut.pdf2latex_mvnormal(x=r's', m=r'\\mu',v=r'\\Sigma')\n",
    "L = nut.pdf2latex_mvnormal(x=r'x_n', m=0,v=r'\\Sigma')\n",
    "#L = nut.pdf2latex_gamma(x=r'x', a=r'a',b=r'b')\n",
    "#L = nut.pdf2latex_invgamma(x=r'x', a=r'a',b=r'b')\n",
    "#L = nut.pdf2latex_beta(x=r'\\pi', a=r'\\alpha',b=r'\\beta')\n",
    "\n",
    "eq = L[0]+'='+L[1]+'='+L[2]\n",
    "display(Math(eq))\n",
    "display(Latex(eq))\n",
    "display(HTML(nut.eqs2html_table(L)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the marginal likelihood\n",
    "\n",
    "* #### Model 1\n",
    "\n",
    "\\begin{eqnarray}\n",
    "p(X| m=1) & = & \\int d{s_{1:K}} p(X|s_{1:K}) p(s_{1:K}) = \n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\mathcal{N}(x_t; (Ax_{t-1}), Q)=\\left|{ 2\\pi Q } \\right|^{-1/2} \\exp\\left(-\\frac12 ({x_t} - {(Ax_{t-1})} )^\\top {Q}^{-1} ({x_t} - {(Ax_{t-1})} ) \\right)=\\exp\\left(  -\\frac{1}{2}\\trace {Q}^{-1} {x_t}{x_t}^\\top + \\trace {Q}^{-1} {x_t}{(Ax_{t-1})}^\\top -\\frac{1}{2}\\trace {Q}^{-1} {(Ax_{t-1})}{(Ax_{t-1})}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}Q\\right|\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\mathcal{N}(y_t; (Cx_{t}), R)=\\left|{ 2\\pi R } \\right|^{-1/2} \\exp\\left(-\\frac12 ({y_t} - {(Cx_{t})} )^\\top {R}^{-1} ({y_t} - {(Cx_{t})} ) \\right)=\\exp\\left(  -\\frac{1}{2}\\trace {R}^{-1} {y_t}{y_t}^\\top + \\trace {R}^{-1} {y_t}{(Cx_{t})}^\\top -\\frac{1}{2}\\trace {R}^{-1} {(Cx_{t})}{(Cx_{t})}^\\top  -\\frac{1}{2}\\log \\left|2{\\pi}R\\right|\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# %load template_equations.py\n",
    "from IPython.display import display, Math, Latex, HTML\n",
    "import notes_utilities as nut\n",
    "from importlib import reload\n",
    "reload(nut)\n",
    "Latex('$\\DeclareMathOperator{\\trace}{Tr}$')\n",
    "\n",
    "#L = nut.pdf2latex_gauss(x=r's', m=r'\\mu',v=r'v')\n",
    "L = nut.pdf2latex_mvnormal(x=r'x_t', m=r'(Ax_{t-1})',v=r'Q')\n",
    "#L = nut.pdf2latex_mvnormal(x=r's', m=0,v=r'I')\n",
    "#L = nut.pdf2latex_gamma(x=r'x', a=r'a',b=r'b')\n",
    "#L = nut.pdf2latex_invgamma(x=r'x', a=r'a',b=r'b')\n",
    "#L = nut.pdf2latex_beta(x=r'\\pi', a=r'\\alpha',b=r'\\beta')\n",
    "\n",
    "eq = L[0]+'='+L[1]+'='+L[2]\n",
    "display(Math(eq))\n",
    "\n",
    "L = nut.pdf2latex_mvnormal(x=r'y_t', m=r'(Cx_{t})',v=r'R')\n",
    "eq = L[0]+'='+L[1]+'='+L[2]\n",
    "display(Math(eq))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\n",
      "  \"key\": \"5fe2a052-2599-465d-96fa-793a58b02ea2\",\n",
      "  \"transport\": \"tcp\",\n",
      "  \"signature_scheme\": \"hmac-sha256\",\n",
      "  \"shell_port\": 65197,\n",
      "  \"stdin_port\": 65199,\n",
      "  \"ip\": \"127.0.0.1\",\n",
      "  \"hb_port\": 65201,\n",
      "  \"control_port\": 65200,\n",
      "  \"iopub_port\": 65198\n",
      "}\n",
      "\n",
      "Paste the above JSON into a file, and connect with:\n",
      "    $> ipython <app> --existing <file>\n",
      "or, if you are local, you can connect with just:\n",
      "    $> ipython <app> --existing kernel-457030b1-3d25-4222-9094-06fb4bbdbefe.json \n",
      "or even just:\n",
      "    $> ipython <app> --existing \n",
      "if this is the most recent IPython session you have started.\n"
     ]
    }
   ],
   "source": [
    "%connect_info"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.1"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}