Chapter 18 - Loglinear Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 18. Loglinear Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18.1 The Loglinear Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $X = (X_1, \\dots, X_m)$ be a random vector with probability\n",
    "\n",
    "$$ f(x) = \\mathbb{P}(X = x) = \\mathbb{P}(X_1 = x_1, \\dots, X_m = x_m) $$\n",
    "\n",
    "Let $r_j$ be the number of values that $X_j$ takes; without loss of generality, assume $X_j \\in \\{ 0, 1, \\dots, r_j - 1 \\}$.  Suppose we have $n$ such vectors.\n",
    "\n",
    "We can think of the data as a sample from a Multinomial with $N = r_1 \\times r_2 \\times \\dots \\times r_m$ categories.  The data can be represented as counts in a $r_1 \\times r_2 \\times \\dots \\times r_m$ table.  Let $p = \\{ p_1, \\dots, p_N \\}$ denote the multinomial parameter.\n",
    "\n",
    "Let $S = \\{ 1, \\dots, m \\}$.  Given a vector $x = (x_1, \\dots, x_m)$ and a subset $A \\subset S$, let $x_A = (x_j : j \\in A)$.  For example, if $A = \\{1, 3\\}$ then $x_A = (x_1, x_3)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 18.1**.  The joint probability function $f(x)$ of a single random vector $X = (X_1, \\dots, X_m)$ can be written as \n",
    "\n",
    "$$ \\log f(x) = \\sum_{A \\subset S} \\psi_A(x) $$\n",
    "\n",
    "where the sum is over all subsets $A$ of $S = \\{1, \\dots, m \\}$ and the $\\psi$'s satisfy the following conditions:\n",
    "\n",
    "1. $\\psi_\\varnothing(x)$ is a constant;\n",
    "2. For every $A \\subset S$, $\\psi_A(x)$ is only a function of $x_A$ and not the rest of the $x_j$'s.\n",
    "3. If $i \\in A$ and $x_i = 0$, then $\\psi_A(x) = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The formula in this theorem is known as the **log-linear expansion** of $f$.  Note that this is the probability function for a single draw.  Each $\\psi_A(x)$ will depend on some unknown parameters $\\beta_A$.  Let $\\beta = (\\beta_A : A \\subset S)$ be the set of all these parameters.  We will write $f(x) = f(x; \\beta)$ when we want to estimate the dependence on the unknown parameters $\\beta$.\n",
    "\n",
    "In terms of the multinomial, the parameter space is\n",
    "\n",
    "$$ \\mathcal{P} = \\left\\{ p = (p_1, \\dots, p_N) : p_j \\geq 0, \\sum_{j=1}^N p_j = 1 \\right\\} $$\n",
    "\n",
    "This is an $N - 1$ dimensional space.  In the log-linear representation, the parameter space is\n",
    "\n",
    "$$ \\Theta = \\Bigg\\{ \\beta = (\\beta_1, \\dots, \\beta_N) : \\beta = \\beta(p), p \\in \\mathcal{P} \\Bigg\\} $$\n",
    "\n",
    "where $\\beta(p)$ is the set of $\\beta$ values associated with $p$.  The set $\\Theta$ is a $N - 1$ dimensional surface in $\\mathbb{R}^N$.  We can always go back and forth between the two parametrizations by writing $\\beta = \\beta(p)$ and $p = p(\\beta)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 18.14**.  Let $(X_a, X_b, X_c)$ be a partition of vectors $(X_1, \\dots, X_m)$.  Then $X_b \\text{ ⫫ } X_c \\; | \\; X_a$ if and only if all the $\\psi$-terms in the log-linear expansion that have at least one coordinate in $b$ and one coordinate in $c$ are 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To prove this Theorem, we will use the following Lemma whose proof follows from the definition of conditional independence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Lemma 18.5**.  A partition $(X_a, X_b, X_c)$ satisfies $X_b \\text{ ⫫ } X_c \\; | \\; X_a$ if and only if $f(x_a, x_b, x_c) = g(x_a, x_b) h(x_a, x_c)$ for some functions $g$ and $h$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof of Theorem 18.14**.  Suppose that $\\psi_t$ is 0 whenever $t$ has coordinates in $b$ and $c$.  Hence, $\\psi_t$ is 0 if $t$ is not a subset of $a \\cup b$ or $t$ is not a subset of $a \\cup c$.  Therefore,\n",
    "\n",
    "$$ \\log f(x) = \\sum_{t \\subset a \\cup b} \\psi_t(x) + \\sum_{t \\subset a \\cup c} \\psi_t(x) - \\sum_{t \\subset a} \\psi_t(x) $$\n",
    "\n",
    "Exponentiating, we see that the joint density is of the form $g(x_a, x_b) h(x_a, x_c)$.  By Lemma 18.5, $X_b \\text{ ⫫ } X_c \\; | \\; X_a$.  The converse follows by reversing the argument."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18.2 Graphical Log-Linear Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\log f(x) = \\sum_{A \\subset S} \\psi_A(x)$ be a log-linear model.  Then $f$ is **graphical** if all $\\psi$-terms are non-zero except for any pair of coordinates not in the edge set for some graph $\\mathcal{G}$.  In other words, $\\psi_A(x) = 0$ if and only if $\\{i, j\\} \\subset A$ and $(i, j)$ is not an edge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a way to think about this definition: if you can add a term to the model and the graph does not change, then the model is not graphical."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18.3 Hierarchical Log-Linear Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a set of log-linear models that is larger than the set of graphical models and that are used quite a bit.  These are hierarchical log-linear models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A log-linear model is **hierarchical** if $\\psi_a = 0$ and $a \\subset t$ implies that $\\psi_t = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Lemma 18.9**.  A graphical model is hierarchical but the reverse need not be true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18.4 Model Generators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hierarchical models can be written succintly using **generators**.  This is most easily explained by example.  Suppose that $X = (X_1, X_2, X_3)$.  Then, $M = 1.2 + 1.3$ stands for:\n",
    "\n",
    "$$ \\log f = \\psi_\\varnothing + \\psi_1 + \\psi_2 + \\psi_3 + \\psi_{12} + \\psi_{13}$$\n",
    "\n",
    "The formula $M = 1.2 + 1.3$ says: \"include $\\psi_{12}$ and $\\psi_{13}$\".  We have to also include the lower terms or it won't be hierarchical.  The generator $M = 1.2.3$ is the **saturated** model\n",
    "\n",
    "$$ \\log f = \\psi_\\varnothing + \\psi_1 + \\psi_2 + \\psi_3 + \\psi_{12} + \\psi_{13} + \\psi_{23} + \\psi_{123}$$\n",
    "\n",
    "The saturated models corresponds to fitting and unconstrained multinomial.  Consider $M = 1 + 2 + 3$ which means\n",
    "\n",
    "$$ \\log f = \\psi_\\varnothing + \\psi_1 + \\psi_2 + \\psi_3 $$\n",
    "\n",
    "This is the mutual independence model.  Finally, consider $M = 1.2$ which has log-linear expansion\n",
    "\n",
    "$$ \\log f = \\psi_\\varnothing + \\psi_1 + \\psi_2 + \\psi_{12} $$\n",
    "\n",
    "This models makes $X_3 | X_2 = x_2, X_1 = x_1$ a uniform distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18.5 Lattices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hierarchical models can be organized into something called a **lattice**.  This is the set of all hierarchical models partially ordered by inclusion.  The set of all hierarchical models for two variables can be illustrated as in the figure below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"134pt\" height=\"260pt\"\r\n",
       " viewBox=\"0.00 0.00 134.00 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-256 130,-256 130,4 -4,4\"/>\r\n",
       "<!-- 1.2 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>1.2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.2</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>1 + 2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63\" cy=\"-162\" rx=\"29.795\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1 + 2</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2&#45;&gt;1 + 2 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>1.2&#45;&gt;1 + 2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M63,-215.697C63,-207.983 63,-198.712 63,-190.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"66.5001,-190.104 63,-180.104 59.5001,-190.104 66.5001,-190.104\"/>\r\n",
       "</g>\r\n",
       "<!-- 1 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2&#45;&gt;1 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>1 + 2&#45;&gt;1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.4685,-144.411C50.1645,-136.042 44.8508,-125.71 40.0473,-116.37\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"43.1563,-114.762 35.4703,-107.47 36.9313,-117.964 43.1563,-114.762\"/>\r\n",
       "</g>\r\n",
       "<!-- 2 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">2</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2&#45;&gt;2 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>1 + 2&#45;&gt;2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M71.5315,-144.411C75.8355,-136.042 81.1492,-125.71 85.9527,-116.37\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"89.0687,-117.964 90.5297,-107.47 82.8437,-114.762 89.0687,-117.964\"/>\r\n",
       "</g>\r\n",
       "<!-- 0 -->\r\n",
       "<g id=\"node5\" class=\"node\"><title>0</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">0</text>\r\n",
       "</g>\r\n",
       "<!-- 1&#45;&gt;0 -->\r\n",
       "<g id=\"edge4\" class=\"edge\"><title>1&#45;&gt;0</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M35.3496,-72.7646C39.7115,-64.2831 45.1469,-53.7144 50.0413,-44.1974\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"53.2346,-45.6409 54.6957,-35.1473 47.0096,-42.4395 53.2346,-45.6409\"/>\r\n",
       "</g>\r\n",
       "<!-- 2&#45;&gt;0 -->\r\n",
       "<g id=\"edge5\" class=\"edge\"><title>2&#45;&gt;0</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M90.6504,-72.7646C86.2885,-64.2831 80.8531,-53.7144 75.9587,-44.1974\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"78.9904,-42.4395 71.3043,-35.1473 72.7654,-45.6409 78.9904,-42.4395\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x10d35869be0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Digraph\n",
    "\n",
    "d = Digraph()\n",
    "\n",
    "d.edge('1.2', '1 + 2')\n",
    "d.edge('1 + 2', '1')\n",
    "d.edge('1 + 2', '2')\n",
    "d.edge('1', '0')\n",
    "d.edge('2', '0')\n",
    "\n",
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$M = 1.2$ is the saturated model, $M = 1 + 2$ is the independence model, $M = 1$ is the independence model plus $X_2 | X_1$ is uniform, $M = 2$ is the independence model plus $X_1 | X_2$ is uniform, $M = 0$ is the uniform distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lattice of trivariate models is shown in the figure below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"320pt\" height=\"548pt\"\r\n",
       " viewBox=\"0.00 0.00 319.55 548.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 544)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-544 315.546,-544 315.546,4 -4,4\"/>\r\n",
       "<!-- 1.2.3 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>1.2.3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"163\" cy=\"-522\" rx=\"28.6953\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"163\" y=\"-518.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.2.3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 1.3 + 2.3 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>1.2 + 1.3 + 2.3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"163\" cy=\"-450\" rx=\"64.9885\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"163\" y=\"-446.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.2 + 1.3 + 2.3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2.3&#45;&gt;1.2 + 1.3 + 2.3 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>1.2.3&#45;&gt;1.2 + 1.3 + 2.3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M163,-503.697C163,-495.983 163,-486.712 163,-478.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"166.5,-478.104 163,-468.104 159.5,-478.104 166.5,-478.104\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 1.3 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>1.2 + 1.3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"58\" cy=\"-378\" rx=\"43.5923\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"58\" y=\"-374.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.2 + 1.3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 1.3 + 2.3&#45;&gt;1.2 + 1.3 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>1.2 + 1.3 + 2.3&#45;&gt;1.2 + 1.3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M139.174,-433.116C124.164,-423.11 104.645,-410.096 88.508,-399.339\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"90.2035,-396.263 79.9415,-393.628 86.3206,-402.087 90.2035,-396.263\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 2.3 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>1.2 + 2.3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"163\" cy=\"-378\" rx=\"43.5923\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"163\" y=\"-374.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.2 + 2.3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 1.3 + 2.3&#45;&gt;1.2 + 2.3 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>1.2 + 1.3 + 2.3&#45;&gt;1.2 + 2.3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M163,-431.697C163,-423.983 163,-414.712 163,-406.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"166.5,-406.104 163,-396.104 159.5,-406.104 166.5,-406.104\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.3 + 2.3 -->\r\n",
       "<g id=\"node5\" class=\"node\"><title>1.3 + 2.3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"268\" cy=\"-378\" rx=\"43.5923\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"268\" y=\"-374.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.3 + 2.3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 1.3 + 2.3&#45;&gt;1.3 + 2.3 -->\r\n",
       "<g id=\"edge4\" class=\"edge\"><title>1.2 + 1.3 + 2.3&#45;&gt;1.3 + 2.3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M186.826,-433.116C201.836,-423.11 221.355,-410.096 237.492,-399.339\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"239.679,-402.087 246.058,-393.628 235.797,-396.263 239.679,-402.087\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 3 -->\r\n",
       "<g id=\"node6\" class=\"node\"><title>1.2 + 3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"61\" cy=\"-306\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"61\" y=\"-302.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.2 + 3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 1.3&#45;&gt;1.2 + 3 -->\r\n",
       "<g id=\"edge5\" class=\"edge\"><title>1.2 + 1.3&#45;&gt;1.2 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M58.7416,-359.697C59.0722,-351.983 59.4695,-342.712 59.838,-334.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"63.3355,-334.245 60.267,-324.104 56.3419,-333.945 63.3355,-334.245\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.3 + 2 -->\r\n",
       "<g id=\"node7\" class=\"node\"><title>1.3 + 2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"166\" cy=\"-306\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"166\" y=\"-302.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.3 + 2</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 1.3&#45;&gt;1.3 + 2 -->\r\n",
       "<g id=\"edge7\" class=\"edge\"><title>1.2 + 1.3&#45;&gt;1.3 + 2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M80.3777,-362.496C96.4775,-352.061 118.368,-337.873 135.984,-326.455\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"138.265,-329.147 144.753,-320.771 134.458,-323.273 138.265,-329.147\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 2.3&#45;&gt;1.2 + 3 -->\r\n",
       "<g id=\"edge6\" class=\"edge\"><title>1.2 + 2.3&#45;&gt;1.2 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M141.368,-362.155C126.467,-351.929 106.492,-338.22 90.184,-327.028\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"91.7825,-323.88 81.5568,-321.108 87.8215,-329.652 91.7825,-323.88\"/>\r\n",
       "</g>\r\n",
       "<!-- 2.3 + 1 -->\r\n",
       "<g id=\"node8\" class=\"node\"><title>2.3 + 1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"264\" cy=\"-306\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"264\" y=\"-302.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">2.3 + 1</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 2.3&#45;&gt;2.3 + 1 -->\r\n",
       "<g id=\"edge9\" class=\"edge\"><title>1.2 + 2.3&#45;&gt;2.3 + 1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M184.42,-362.155C199.175,-351.929 218.954,-338.22 235.102,-327.028\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"237.419,-329.681 243.645,-321.108 233.432,-323.927 237.419,-329.681\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.3 + 2.3&#45;&gt;1.3 + 2 -->\r\n",
       "<g id=\"edge8\" class=\"edge\"><title>1.3 + 2.3&#45;&gt;1.3 + 2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M246.368,-362.155C231.467,-351.929 211.492,-338.22 195.184,-327.028\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"196.782,-323.88 186.557,-321.108 192.822,-329.652 196.782,-323.88\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.3 + 2.3&#45;&gt;2.3 + 1 -->\r\n",
       "<g id=\"edge10\" class=\"edge\"><title>1.3 + 2.3&#45;&gt;2.3 + 1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M267.011,-359.697C266.57,-351.983 266.041,-342.712 265.549,-334.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"269.042,-333.888 264.977,-324.104 262.054,-334.288 269.042,-333.888\"/>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2 + 3 -->\r\n",
       "<g id=\"node9\" class=\"node\"><title>1 + 2 + 3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"116\" cy=\"-234\" rx=\"44.393\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"116\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1 + 2 + 3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 3&#45;&gt;1 + 2 + 3 -->\r\n",
       "<g id=\"edge11\" class=\"edge\"><title>1.2 + 3&#45;&gt;1 + 2 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M73.7563,-288.765C80.5305,-280.143 88.9992,-269.365 96.5719,-259.727\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"99.5048,-261.659 102.931,-251.633 94.0005,-257.334 99.5048,-261.659\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.2 -->\r\n",
       "<g id=\"node10\" class=\"node\"><title>1.2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.2</text>\r\n",
       "</g>\r\n",
       "<!-- 1.2 + 3&#45;&gt;1.2 -->\r\n",
       "<g id=\"edge14\" class=\"edge\"><title>1.2 + 3&#45;&gt;1.2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M52.9425,-288.411C48.9199,-280.129 43.9635,-269.925 39.4644,-260.662\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"42.5171,-258.936 34.9997,-251.47 36.2205,-261.994 42.5171,-258.936\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.3 + 2&#45;&gt;1 + 2 + 3 -->\r\n",
       "<g id=\"edge12\" class=\"edge\"><title>1.3 + 2&#45;&gt;1 + 2 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M154.403,-288.765C148.307,-280.23 140.702,-269.582 133.871,-260.019\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"136.541,-257.736 127.881,-251.633 130.845,-261.805 136.541,-257.736\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.3 -->\r\n",
       "<g id=\"node11\" class=\"node\"><title>1.3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"205\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"205\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1.3</text>\r\n",
       "</g>\r\n",
       "<!-- 1.3 + 2&#45;&gt;1.3 -->\r\n",
       "<g id=\"edge15\" class=\"edge\"><title>1.3 + 2&#45;&gt;1.3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M175.242,-288.411C179.905,-280.042 185.662,-269.71 190.865,-260.37\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"194.014,-261.909 195.824,-251.47 187.899,-258.502 194.014,-261.909\"/>\r\n",
       "</g>\r\n",
       "<!-- 2.3 + 1&#45;&gt;1 + 2 + 3 -->\r\n",
       "<g id=\"edge13\" class=\"edge\"><title>2.3 + 1&#45;&gt;1 + 2 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M238.558,-292.967C214.9,-281.777 179.387,-264.981 152.67,-252.344\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"154.157,-249.176 143.621,-248.064 151.164,-255.503 154.157,-249.176\"/>\r\n",
       "</g>\r\n",
       "<!-- 2.3 -->\r\n",
       "<g id=\"node12\" class=\"node\"><title>2.3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"277\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"277\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">2.3</text>\r\n",
       "</g>\r\n",
       "<!-- 2.3 + 1&#45;&gt;2.3 -->\r\n",
       "<g id=\"edge16\" class=\"edge\"><title>2.3 + 1&#45;&gt;2.3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M267.147,-288.055C268.579,-280.346 270.309,-271.027 271.918,-262.364\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"275.407,-262.746 273.792,-252.275 268.525,-261.468 275.407,-262.746\"/>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2 -->\r\n",
       "<g id=\"node13\" class=\"node\"><title>1 + 2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"77\" cy=\"-162\" rx=\"29.795\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"77\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1 + 2</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2 + 3&#45;&gt;1 + 2 -->\r\n",
       "<g id=\"edge20\" class=\"edge\"><title>1 + 2 + 3&#45;&gt;1 + 2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M106.559,-216.055C101.893,-207.679 96.1679,-197.404 91.0033,-188.134\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"94.0099,-186.339 86.0853,-179.307 87.8949,-189.746 94.0099,-186.339\"/>\r\n",
       "</g>\r\n",
       "<!-- 1 + 3 -->\r\n",
       "<g id=\"node14\" class=\"node\"><title>1 + 3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"155\" cy=\"-162\" rx=\"29.795\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"155\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1 + 3</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2 + 3&#45;&gt;1 + 3 -->\r\n",
       "<g id=\"edge21\" class=\"edge\"><title>1 + 2 + 3&#45;&gt;1 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M125.441,-216.055C130.107,-207.679 135.832,-197.404 140.997,-188.134\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"144.105,-189.746 145.915,-179.307 137.99,-186.339 144.105,-189.746\"/>\r\n",
       "</g>\r\n",
       "<!-- 2 + 3 -->\r\n",
       "<g id=\"node15\" class=\"node\"><title>2 + 3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"233\" cy=\"-162\" rx=\"29.795\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"233\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">2 + 3</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2 + 3&#45;&gt;2 + 3 -->\r\n",
       "<g id=\"edge22\" class=\"edge\"><title>1 + 2 + 3&#45;&gt;2 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M139.959,-218.666C158.511,-207.566 184.285,-192.146 204.001,-180.35\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"205.802,-183.351 212.587,-175.213 202.208,-177.344 205.802,-183.351\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.2&#45;&gt;1 + 2 -->\r\n",
       "<g id=\"edge17\" class=\"edge\"><title>1.2&#45;&gt;1 + 2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M38.0966,-217.465C44.4162,-208.617 52.4737,-197.337 59.6148,-187.339\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"62.6234,-189.149 65.5877,-178.977 56.9273,-185.08 62.6234,-189.149\"/>\r\n",
       "</g>\r\n",
       "<!-- 1.3&#45;&gt;1 + 3 -->\r\n",
       "<g id=\"edge18\" class=\"edge\"><title>1.3&#45;&gt;1 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M193.903,-217.465C187.584,-208.617 179.526,-197.337 172.385,-187.339\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"175.073,-185.08 166.412,-178.977 169.377,-189.149 175.073,-185.08\"/>\r\n",
       "</g>\r\n",
       "<!-- 2.3&#45;&gt;2 + 3 -->\r\n",
       "<g id=\"edge19\" class=\"edge\"><title>2.3&#45;&gt;2 + 3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M267.016,-217.116C261.58,-208.469 254.733,-197.575 248.609,-187.833\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"251.535,-185.911 243.25,-179.307 245.609,-189.636 251.535,-185.911\"/>\r\n",
       "</g>\r\n",
       "<!-- 1 -->\r\n",
       "<g id=\"node16\" class=\"node\"><title>1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"83\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"83\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">1</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2&#45;&gt;1 -->\r\n",
       "<g id=\"edge23\" class=\"edge\"><title>1 + 2&#45;&gt;1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M78.4831,-143.697C79.1443,-135.983 79.9389,-126.712 80.6761,-118.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"84.167,-118.367 81.5339,-108.104 77.1926,-117.769 84.167,-118.367\"/>\r\n",
       "</g>\r\n",
       "<!-- 2 -->\r\n",
       "<g id=\"node17\" class=\"node\"><title>2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"155\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"155\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">2</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 2&#45;&gt;2 -->\r\n",
       "<g id=\"edge24\" class=\"edge\"><title>1 + 2&#45;&gt;2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M93.1617,-146.496C104.289,-136.51 119.247,-123.086 131.664,-111.942\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"134.287,-114.292 139.391,-105.008 129.611,-109.082 134.287,-114.292\"/>\r\n",
       "</g>\r\n",
       "<!-- 1 + 3&#45;&gt;1 -->\r\n",
       "<g id=\"edge25\" class=\"edge\"><title>1 + 3&#45;&gt;1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M139.731,-146.155C129.669,-136.372 116.33,-123.404 105.114,-112.5\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"107.365,-109.806 97.7549,-105.345 102.485,-114.825 107.365,-109.806\"/>\r\n",
       "</g>\r\n",
       "<!-- 3 -->\r\n",
       "<g id=\"node18\" class=\"node\"><title>3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"227\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"227\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">3</text>\r\n",
       "</g>\r\n",
       "<!-- 1 + 3&#45;&gt;3 -->\r\n",
       "<g id=\"edge26\" class=\"edge\"><title>1 + 3&#45;&gt;3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M170.269,-146.155C180.331,-136.372 193.67,-123.404 204.886,-112.5\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"207.515,-114.825 212.245,-105.345 202.635,-109.806 207.515,-114.825\"/>\r\n",
       "</g>\r\n",
       "<!-- 2 + 3&#45;&gt;2 -->\r\n",
       "<g id=\"edge27\" class=\"edge\"><title>2 + 3&#45;&gt;2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M216.838,-146.496C205.711,-136.51 190.753,-123.086 178.336,-111.942\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"180.389,-109.082 170.609,-105.008 175.713,-114.292 180.389,-109.082\"/>\r\n",
       "</g>\r\n",
       "<!-- 2 + 3&#45;&gt;3 -->\r\n",
       "<g id=\"edge28\" class=\"edge\"><title>2 + 3&#45;&gt;3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M231.517,-143.697C230.856,-135.983 230.061,-126.712 229.324,-118.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"232.807,-117.769 228.466,-108.104 225.833,-118.367 232.807,-117.769\"/>\r\n",
       "</g>\r\n",
       "<!-- 0 -->\r\n",
       "<g id=\"node19\" class=\"node\"><title>0</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"155\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"155\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">0</text>\r\n",
       "</g>\r\n",
       "<!-- 1&#45;&gt;0 -->\r\n",
       "<g id=\"edge29\" class=\"edge\"><title>1&#45;&gt;0</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M97.5703,-74.8345C107.75,-64.9376 121.524,-51.5462 133.031,-40.3591\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"135.474,-42.865 140.204,-33.3847 130.595,-37.8461 135.474,-42.865\"/>\r\n",
       "</g>\r\n",
       "<!-- 2&#45;&gt;0 -->\r\n",
       "<g id=\"edge30\" class=\"edge\"><title>2&#45;&gt;0</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M155,-71.6966C155,-63.9827 155,-54.7125 155,-46.1124\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"158.5,-46.1043 155,-36.1043 151.5,-46.1044 158.5,-46.1043\"/>\r\n",
       "</g>\r\n",
       "<!-- 3&#45;&gt;0 -->\r\n",
       "<g id=\"edge31\" class=\"edge\"><title>3&#45;&gt;0</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M212.43,-74.8345C202.25,-64.9376 188.476,-51.5462 176.969,-40.3591\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"179.405,-37.8461 169.796,-33.3847 174.526,-42.865 179.405,-37.8461\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x10d358bfb38>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Digraph\n",
    "\n",
    "d = Digraph()\n",
    "\n",
    "d.edge('1.2.3', '1.2 + 1.3 + 2.3')\n",
    "\n",
    "d.edge('1.2 + 1.3 + 2.3', '1.2 + 1.3')\n",
    "d.edge('1.2 + 1.3 + 2.3', '1.2 + 2.3')\n",
    "d.edge('1.2 + 1.3 + 2.3', '1.3 + 2.3')\n",
    "\n",
    "d.edge('1.2 + 1.3', '1.2 + 3')\n",
    "d.edge('1.2 + 2.3', '1.2 + 3')\n",
    "\n",
    "d.edge('1.2 + 1.3', '1.3 + 2')\n",
    "d.edge('1.3 + 2.3', '1.3 + 2')\n",
    "\n",
    "d.edge('1.2 + 2.3', '2.3 + 1')\n",
    "d.edge('1.3 + 2.3', '2.3 + 1')\n",
    "\n",
    "d.edge('1.2 + 3', '1 + 2 + 3')\n",
    "d.edge('1.3 + 2', '1 + 2 + 3')\n",
    "d.edge('2.3 + 1', '1 + 2 + 3')\n",
    "\n",
    "d.edge('1.2 + 3', '1.2')\n",
    "d.edge('1.3 + 2', '1.3')\n",
    "d.edge('2.3 + 1', '2.3')\n",
    "\n",
    "d.edge('1.2', '1 + 2')\n",
    "d.edge('1.3', '1 + 3')\n",
    "d.edge('2.3', '2 + 3')\n",
    "\n",
    "d.edge('1 + 2 + 3', '1 + 2')\n",
    "d.edge('1 + 2 + 3', '1 + 3')\n",
    "d.edge('1 + 2 + 3', '2 + 3')\n",
    "\n",
    "d.edge('1 + 2', '1')\n",
    "d.edge('1 + 2', '2')\n",
    "d.edge('1 + 3', '1')\n",
    "d.edge('1 + 3', '3')\n",
    "d.edge('2 + 3', '2')\n",
    "d.edge('2 + 3', '3')\n",
    "\n",
    "d.edge('1', '0')\n",
    "d.edge('2', '0')\n",
    "d.edge('3', '0')\n",
    "\n",
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18.6 Fitting Log-Linear Models to Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\beta$ denote all the parameters in a log-linear model $M$.  The log-likelihood for $\\beta$ is:\n",
    "\n",
    "$$ \\ell(\\beta) = \\sum_j x_j \\log p_j(\\beta) $$\n",
    "\n",
    "where the sum is over the cells and $p(\\beta)$ denotes the cell probabilities corresponding to $\\beta$.  The MLE $\\hat{\\beta}$ generally has to be found numerically.  The model with all possible $\\psi$-terms is called the  **saturated models**.  We can also fit any **sub-model** which corresponds to setting some subset of $\\psi$ terms to 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For any submodel $M$, define the **deviance** $\\text{dev}(M)$ by\n",
    "\n",
    "$$ \\text{dev}(M) = 2 (\\hat{\\ell}_\\text{sat} - \\hat{\\ell}_M) $$\n",
    "\n",
    "where $\\ell_\\text{sat}$ is the log-likelihood of the saturated model evaluated at the MLE and $\\hat{\\ell}_M$ is the log-likelihood of the model $M$ evaluated at its MLE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 18.14**.  The deviance is the likelihood test statistic for\n",
    "\n",
    "$$\n",
    "H_0 : \\text{the model is } M\n",
    "\\quad \\text{versus} \\quad\n",
    "H_1 : \\text{the model is not } M\n",
    "$$\n",
    "\n",
    "Under $H_0$, $\\text{dev}(M) \\leadsto \\chi^2_\\nu$ with $\\nu$ degrees of freedom equal to the difference in the number of parameters between the saturated model and $M$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to find a good model is to use the deviance to test every sub-model.  Every model that is not rejected by this test is then considered a plausible model.  However, this is not a good strategy for two reasons.  First, we will end up doing many tests, which means there is ample opportunity for making Type I and Type II errors.  Second, we will end up using models where we failed to reject $H_0$.  But we might fail to reject $H_0$ due to low power.  The result is that we end up with a bad model just due to low power."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many model searching strategies.  A common approach is to use some form of *penalized likelihood*.  One version of penalized is the AIC that we used in regression.  For any model $M$ define\n",
    "\n",
    "$$ \\text{AIC}(M) = -2 \\left( \\hat{\\ell}(M) - |M|\\right) $$\n",
    "\n",
    "where $|M|$ is the number of parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a set of models $\\{ M_1, M_2, \\dots \\}$.  Let $\\hat{f}_j(x)$ denote the estimated probability function obtained by using the maximum likelihood estimator of model $M_j$.  Thus, $\\hat{f}_j(x) = \\hat{f}(x; \\hat{\\beta}_j)$ where $\\hat{\\beta}_j$ is the MLE of the set of parameters $\\beta_j$ for model $M_j$.  We will use the loss function $D(f, \\hat{f})$ where\n",
    "\n",
    "$$ D(f, g) = \\sum_x f(x) \\log \\frac{f(x)}{g(x)} $$\n",
    "\n",
    "is the Kullback-Leibler divergence between two probability density functions.  The corresponding risk function is $R(f, \\hat{f}) = \\mathbb{E}(D(f, \\hat{f}))$.\n",
    "\n",
    "Notice that $D(f, \\hat{f})  = c - A(f, \\hat{f})$ where $c = \\sum_x f(x) \\log f(x)$ does not depend on $\\hat{f}$ and \n",
    "\n",
    "$$ A(f, \\hat{f}) = \\sum_x f(x) \\log \\hat{f}(x) $$\n",
    "\n",
    "Thus minimizing the risk is equivalent to minimizing $a(f, \\hat{f}) = \\mathbb{E}(A(f, \\hat{f}))$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is tempting to estimate $a(f, \\hat{f})$ by $\\sum_x \\log \\hat{f}(x)$ but, just as the training error in regression is highly biased estimate of prediction risk, it is also the case that $\\sum_x \\log \\hat{f}(x)$ is a highly biased estimate of $a(f, \\hat{f})$.  In fact, the bias is approximately equal to $|M_j|$.  Thus:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 18.15**.  $\\text{AIC}(M_j)$ is an approximately unbiased estimate of $a(f, \\hat{f})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After finding a \"best model\" this way we can draw the corresponding graph.  We can also check the overall fit of the selected model using the deviance as described above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18.8 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 18.8.1**.  Solve for the $p_{ij}$'s in terms of the $\\beta$'s in Example 18.3:\n",
    "\n",
    "*Example*: Let $X = (X_1, X_2)$ where $X_1 \\in \\{0, 1\\}$ and $X_2 \\in \\{ 0, 1, 2 \\}$.  The joint distribution of $n$ such random vectors is a multinomial with 6 categories.  The multinomial parameters can be written as a 2-by-3 table as follows:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccccc}\n",
    "\\hline\n",
    "\\text{multinomial} & x_2 & 0 & 1 & 2 \\\\\n",
    "\\hline\n",
    "x_1 & 0 & p_{00} & p_{01} & p_{02} \\\\\n",
    "    & 1 & p_{10} & p_{11} & p_{12} \\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The $n$ data vectors can be summarized as follows:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccccc}\n",
    "\\hline\n",
    "\\text{multinomial} & x_2 & 0 & 1 & 2 \\\\\n",
    "\\hline\n",
    "x_1 & 0 & C_{00} & C_{01} & C_{02} \\\\\n",
    "    & 1 & C_{10} & C_{11} & C_{12} \\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "For $x = (x_1, x_2)$, the log-linear expansion takes the form\n",
    "\n",
    "$$ \\log f(x) = \\psi_\\varnothing(x) + \\psi_1(x) + \\psi_2(x) + \\psi_{12}(x) $$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& \\psi_\\varnothing(x) = \\log p_{00} \\\\\n",
    "& \\psi_1(x) = x_1 \\log \\left( \\frac{p_{10}}{p_{00}} \\right) \\\\\n",
    "& \\psi_2(x) = I(x_2 = 1) \\log \\left( \\frac{p_{01}}{p_{00}} \\right) \n",
    "            + I(x_2 = 2) \\log \\left( \\frac{p_{02}}{p_{00}} \\right) \\\\\n",
    "& \\psi_{12}(x) = I(x_1 = 1, x_2 = 1) \\log \\left( \\frac{p_{11}p_{00}}{p_{01}p_{10}} \\right)\n",
    "               + I(x_1 = 1, x_2 = 2) \\log \\left( \\frac{p_{12}p_{00}}{p_{02}p_{10}} \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The six parameters of this model are:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccc}\n",
    "\\beta_1 = \\log p_{00} &\n",
    "\\beta_2 = \\log \\left( \\frac{p_{10}}{p_{00}} \\right) &\n",
    "\\beta_3 = \\log \\left( \\frac{p_{01}}{p_{00}} \\right) \\\\\n",
    "\\beta_4 = \\log \\left( \\frac{p_{02}}{p_{00}} \\right) &\n",
    "\\beta_5 = \\log \\left( \\frac{p_{11}p_{00}}{p_{01}p_{10}} \\right) &\n",
    "\\beta_6 = \\log \\left( \\frac{p_{12}p_{00}}{p_{02}p_{10}} \\right)\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exponentiate the six definitions for $\\beta_k$ above and get:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccc}\n",
    "p_{00} = \\exp \\beta_1 &\n",
    "\\frac{p_{10}}{p_{00}} = \\exp \\beta_2 &\n",
    "\\frac{p_{01}}{p_{00}} = \\exp \\beta_3 \\\\\n",
    "\\frac{p_{02}}{p_{00}} = \\exp \\beta_4 &\n",
    "\\frac{p_{11}p_{00}}{p_{01}p_{10}} = \\exp \\beta_5 &\n",
    "\\frac{p_{12}p_{00}}{p_{02}p_{10}} = \\exp \\beta_6\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "and then, isolating the $p_{ij}$s,\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccc}\n",
    "p_{00} = \\exp \\beta_1 &\n",
    "p_{10} = \\exp \\left( \\beta_1 + \\beta_2 \\right) &\n",
    "p_{01} = \\exp \\left( \\beta_1 + \\beta_3 \\right) \\\\\n",
    "p_{02} = \\exp \\left( \\beta_1 + \\beta_4 \\right) &\n",
    "p_{11} = \\exp \\left( \\beta_1 + \\beta_2 + \\beta_3 + \\beta_5 \\right) &\n",
    "p_{12} = \\exp \\left( \\beta_1 + \\beta_2 + \\beta_4 + \\beta_6 \\right)\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 18.8.2**. Repeat example 18.17 using 7 covariates and $n = 1000$. To avoid numerical problems, replace any zero count with a one.\n",
    "\n",
    "*Example*: Here is a synthetic example.  We generate $n = 100$ random vectors $X = (X_1, \\dots, X_5)$ of length 5.  We generated the data as follows:\n",
    "\n",
    "$$ X_1 \\sim \\text{Bernoulli}(1/2)$$\n",
    "\n",
    "and\n",
    "\n",
    "$$ X_j | X_1, \\dots, X_{j-1} = \\begin{cases}\n",
    "1/4 & \\text{if } X_{j-1} = 0 \\\\\n",
    "3/4 & \\text{if } X_{j-1} = 1\n",
    "\\end{cases}$$\n",
    "\n",
    "It follows that\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f(x_1, \\dots, x_5) &= \n",
    "\\left( \\frac{1}{2} \\right) \\left( \\frac{3}{4} \\right)^{x_1} \\left( \\frac{1}{4} \\right)^{1 - x_1} \\left( \\frac{3}{4} \\right)^{x_2} \\left( \\frac{1}{4} \\right)^{1 - x_2} \\left( \\frac{3}{4} \\right)^{x_3} \\left( \\frac{1}{4} \\right)^{1 - x_3} \\left( \\frac{3}{4} \\right)^{x_4} \\left( \\frac{1}{4} \\right)^{1 - x_4} \\\\\n",
    "&= \\left( \\frac{1}{2} \\right) \\left( \\frac{3}{4} \\right)^{x_1 + x_2 + x_3 + x_4} \\left( \\frac{1}{4} \\right)^{4 - (x_1 + x_2 + x_3 + x_4)} \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We estimated $f$ using three methods: (i)  maximum likelihood treating this as a multinomial with 32 categories (ii) maximum likelihood from the best loglinear model using AIC and forward selection and (iii) maximum likelihood from the best loglinear model using BIC and forward selection. We estimated the risk by simulating the sample 100 times.  The average risks were:\n",
    "\n",
    "| Method | Risk |\n",
    "|--------|------|\n",
    "| MLE    | 0.63 |\n",
    "| AIC    | 0.54 |\n",
    "| BIC    | 0.53 |\n",
    "\n",
    "In this example, there is little difference between AIC and BIC.  Both are better than maximum likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate $n = 1000$ random vectors $X = (X_1, \\dots, X_k)$ of length $k = 7$ as follows:\n",
    "\n",
    "$$ X_1 \\sim \\text{Bernoulli}(1/2) \n",
    "\\quad \\text{and} \\quad\n",
    "X_j | X_1, \\dots, X_{j-1} = \\begin{cases}\n",
    "1/4 & \\text{if } X_{j-1} = 0 \\\\\n",
    "3/4 & \\text{if } X_{j-1} = 1\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "It follows that\n",
    "\n",
    "$$ f(x) = \\frac{1}{2} \\left( \\frac{3}{4} \\right)^{\\sum_{i=1}^{k - 1} x_i} \\left( \\frac{1}{4} \\right)^{(k - 1) - \\left(\\sum_{i=1}^{k - 1} x_i \\right)} $$\n",
    "\n",
    "We will use the KL divergence as our loss function:\n",
    "\n",
    "$$ D(f, g) = \\sum_{x \\in \\Omega} f(x) \\log \\frac{f(x)}{g(x)} $$\n",
    "\n",
    "and estimate the risk function $R(f, \\hat{f}) = \\mathbb{E}[D(f, \\hat{f})]$ by bootstraping the estimation process and calculating the average of the loss functions in each bootstrap step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### MLE estimate\n",
    "\n",
    "The MLE estimate for the full multinomial model is relatively simple: consider the adjusted counts $\\tilde{C}_{\\xi}$ to be the number of times an observation $\\xi$ appears, or 1 if that number of observations is 0.  There are $2^k = 128$ possible observations, so we get $2^k$ adjusted counts $\\tilde{C}_1, \\dots, \\tilde{C}_{2^k}$.  The MLE estimate is then computed as $\\hat{p} = (\\hat{p}_1, \\dots, \\hat{p}_{2^k})$, with $\\hat{p}_i = \\frac{\\tilde{C}_i}{\\sum_\\xi \\tilde{C}_\\xi}$.\n",
    "\n",
    "The function density estimate corresponding to the MLE estimate is a simple lookup table $\\hat{f}(\\xi) = \\hat{p}_\\xi$, as we already have a probability estimate associated with every single possible event.  The function density estimate can then be used to compute the loss function $D(f, \\hat{f})$, which is the KL divergence computed over the universe.   Repeating this process gives us the risk estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from itertools import chain, combinations\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "\n",
    "# Recipe from itertools documentation, https://docs.python.org/2.7/library/itertools.html#recipes\n",
    "def powerset(iterable):\n",
    "    \"powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)\"\n",
    "    s = list(iterable)\n",
    "    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))\n",
    "\n",
    "\n",
    "def generate_samples(n, k):\n",
    "    \"\"\"\n",
    "    Generates n samples of size k according to the synthetic distribution\n",
    "    \n",
    "    Args:\n",
    "       n:  number of samples\n",
    "       k:  sample size\n",
    "       \n",
    "    Returns:\n",
    "       X:  2D array of shape (n, k), representing n samples of size k\n",
    "    \"\"\"\n",
    "    \n",
    "    # Generate a random unifom value between 0 and 1 for each variable in each sample\n",
    "    random_seeds = np.random.uniform(low=0, high=1, size=(k, n))\n",
    "    \n",
    "    # Create a variable to store the generated samples\n",
    "    output = np.empty_like(random_seeds).astype(int)\n",
    "    \n",
    "    # Generate x_1's as Bernoulli(1/2)\n",
    "    output[0] = random_seeds[0] > 0.5\n",
    "    for j in range(1, k):\n",
    "        # Generate x_j's recursively\n",
    "        output[j] = random_seeds[j] > (output[j - 1] == 0) * (1/4) + (output[j - 1] == 1) * (3/4)\n",
    "        \n",
    "    return output.T\n",
    "\n",
    "\n",
    "def generate_universe(k):\n",
    "    \"\"\"\n",
    "    Generates 2**k samples of size k iterating through the valid universe\n",
    "    \n",
    "    Args:\n",
    "       k: sample size\n",
    "       \n",
    "    Returns:\n",
    "       X_universe:  2D array of shape (2**k, k), representing the universe of all samples of size k\n",
    "    \"\"\"\n",
    "    X_universe = np.zeros(shape=(2**k, k), dtype=int)\n",
    "    for i, line in enumerate(powerset(range(k))):\n",
    "        X_universe[i, line] = 1\n",
    "        \n",
    "    return X_universe\n",
    "\n",
    "\n",
    "def row_to_binary(x):\n",
    "    \"\"\"\n",
    "    Translates a single row into a binary low-endian representation\n",
    "    \n",
    "    Args:\n",
    "       x:  1D array of 0s and 1s\n",
    "       \n",
    "    Returns:\n",
    "       sum(x[i] * (2**i) for index i)\n",
    "    \"\"\"\n",
    "    \n",
    "    # Translate a row of 0s and 1s into a binary representation (low-endian)\n",
    "    return sum([x[i] * (2**i) for i in range(len(x))])\n",
    "\n",
    "\n",
    "def samples_to_count(X):\n",
    "    \"\"\"\n",
    "    Counts the number of occurrences of each sample in the dataset\n",
    "    \n",
    "    Args:\n",
    "       X:   2D array  n samples of size k\n",
    "       \n",
    "    Returns:\n",
    "       count:  1D array of size 2**k, where count[i] is the number of\n",
    "               occurrences of i as a row in binary (low-endian) in X,\n",
    "               or 1 if the count would be 0.\n",
    "    \"\"\"\n",
    "    \n",
    "    binary_samples = np.apply_along_axis(row_to_binary, 1, X)\n",
    "    k = X.shape[1]\n",
    "    count = np.zeros(2**k)\n",
    "    for sample in binary_samples:\n",
    "        count[sample] += 1\n",
    "    \n",
    "    # Replace all zeroes with ones\n",
    "    return np.where(count == 0, 1, count)\n",
    "\n",
    "\n",
    "def true_density(x):\n",
    "    \"\"\"\n",
    "    True density function for the synthetic distribution.\n",
    "    \n",
    "    Args:\n",
    "       x:  1D array  sample\n",
    "       \n",
    "    Returns:\n",
    "       value of PDF = (1/2) * (3/4)**(sum(x)) * (1/4)**(len(x) - sum(x))\n",
    "    \"\"\"\n",
    "    \n",
    "    s = sum(x)\n",
    "    k = len(x)\n",
    "    return (1/2) * ((3/4)**s) * ((1/4)**(k - s))\n",
    "\n",
    "\n",
    "def KL_divergence(f, g, X):\n",
    "    \"\"\"\n",
    "    Returns the KL divergence between f and g calculated on universe X\n",
    "    \n",
    "    Args:\n",
    "       f:  1D array => double  probability density function\n",
    "       g:  1D array => double  probability density function\n",
    "       X:  2D array (n, k) of n samples of size k\n",
    "       \n",
    "    Returns:  D(f, g) = sum over samples of f(x) * log(f(x) / g(x))\n",
    "    \"\"\"\n",
    "    \n",
    "    def term(x):\n",
    "        \"\"\" \n",
    "        Return f(x) * log(f(x) / g(x)) \n",
    "        \n",
    "        Args:\n",
    "           x:  1D array sample of values\n",
    "        \"\"\"\n",
    "        fx = f(x)\n",
    "        gx = g(x)\n",
    "        return fx * np.log(fx / gx)\n",
    "    \n",
    "    return sum(np.apply_along_axis(term, 1, X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can now calculate the expected value of the loss function by simulation\n",
    "\n",
    "def p_hat_mle(X):\n",
    "    \"\"\" MLE estimate for the multinomial (replacing zeros with ones) \"\"\"\n",
    "    count = samples_to_count(X)\n",
    "    return count / sum(count)\n",
    "\n",
    "\n",
    "def f_hat_mle(p_hat_mle):\n",
    "    \"\"\" Density function for the estimated multinomial \"\"\"\n",
    "    return lambda x : p_hat_mle[row_to_binary(x)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wrapper for parallel processing on progress bars\n",
    "# Source: https://stackoverflow.com/a/58936697\n",
    "\n",
    "import contextlib\n",
    "import joblib\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "@contextlib.contextmanager\n",
    "def tqdm_joblib(tqdm_object):\n",
    "    \"\"\"Context manager to patch joblib to report into tqdm progress bar given as argument\"\"\"\n",
    "    class TqdmBatchCompletionCallback:\n",
    "        def __init__(self, time, index, parallel):\n",
    "            self.index = index\n",
    "            self.parallel = parallel\n",
    "\n",
    "        def __call__(self, index):\n",
    "            tqdm_object.update()\n",
    "            if self.parallel._original_iterator is not None:\n",
    "                self.parallel.dispatch_next()\n",
    "\n",
    "    old_batch_callback = joblib.parallel.BatchCompletionCallBack\n",
    "    joblib.parallel.BatchCompletionCallBack = TqdmBatchCompletionCallback\n",
    "    try:\n",
    "        yield tqdm_object\n",
    "    finally:\n",
    "        joblib.parallel.BatchCompletionCallBack = old_batch_callback\n",
    "        tqdm_object.close()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8ec99281f688408dba7254a1f5177d4d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='MLE', max=500, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "MLE mean: 0.75\n"
     ]
    }
   ],
   "source": [
    "# Bootstrap MLE\n",
    "import multiprocessing\n",
    "\n",
    "B = 500\n",
    "n = 1000\n",
    "k = 7\n",
    "\n",
    "X_universe = generate_universe(k=k)\n",
    "\n",
    "n_jobs = max(multiprocessing.cpu_count() - 1, 1)\n",
    "\n",
    "def bootstrap_step_mle(i):\n",
    "    XX = generate_samples(n=n, k=k)\n",
    "    p_hat = p_hat_mle(XX)\n",
    "    f_hat = f_hat_mle(p_hat)\n",
    "    return KL_divergence(true_density, f_hat, X_universe)\n",
    "\n",
    "with tqdm_joblib(tqdm_notebook(desc=\"MLE\", total=B)) as progress_bar:\n",
    "    risk_mle = np.array(Parallel(n_jobs=n_jobs)(delayed(bootstrap_step_mle)(i) for i in range(B)))\n",
    "\n",
    "print('MLE mean: %.2f' % risk_mle.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEICAYAAAC+iFRkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAamElEQVR4nO3df5RdZX3v8fenMdxoEkCSIeb3oI1IQIh0muilpaGIJlMg3pariYIE7Yrxyl1wr/GW6m1LV7uuVK8/SmPlRuGGIE20ix+lECW5CIJVLENWCEGgxDg0w8QkBCFBSSXwvX+cZ+jJyTkzZ84+c86Zsz+vtc46+8fz7P09e/Z8Z59n9rMfRQRmZtbefq3ZAZiZ2chzsjczywEnezOzHHCyNzPLASd7M7MccLI3M8sBJ3treZKuk/QnVZTrlfTuRsRUtM9PS/p6I/dpVgsne2u6lKRfkvSipJ9JWitpwsD6iFgZEX/RzBgBJC2U1Fe8LCL+V0T84Qjtr+F/vKx9Odlbq7ggIiYA84B3AH/c5HjM2oqTvbWUiPgZcDeFpA9AutL/yzQ9WdKdkp6X9JykByQddR5Lepukn0paWm4/af3mtI0nJb2/aF23pB9LOijpGUmrJI0Hvg1MS99AXpQ0TdLVkr6R6nVKCkmXSdol6eeSVkr6TUnbUsyri/bzFknflbRf0rOSbpZ0fFp3EzAL+Me0r/+Rlr9T0g/Sth6RtLBoe8sl7Uxx/1TShzL8KKzNONlbS5E0A1gM7KhQ5JNAH9ABTAE+DRzxzA9JZwKbgP8aERvK7GM8sBn4O+BEYBnwt5JOTUWuBz4WEROB04DvRsQvUlz9ETEhvforxLgAmAN8APgy8Bng3cCpwPsl/c5AKMBngWnAKcBM4GqAiLgE+FfSN56I+Jyk6cBdwF8CJwCrgFskdaTPdC2wOMX9H4GtFeKzHHKyt1Zxu6SDwC5gL/BnFcq9DEwFZkfEyxHxQBz5gKffBu4ALo2IOyts43ygNyL+b0QcjogtwC3ARUX7mCvp2Ij4eVo/HH8REYciYhPwC2B9ROyNiGeAByg0UxEROyJic0T8W0TsA74I/E7lzXIxsDEiNkbEqxGxGegButP6V4HTJL0+InZHxGPDjNvamJO9tYr3pSvShcDbgMkVyn2ewlX/ptRkcVXJ+pXADyLi3kH2NRtYkJpCnpf0PPAh4E1p/R9QSKBPS/qepHcN87PsKZp+qcz8BABJJ0rakJqKDgDfoPLnHoj7P5fE/VvA1PTN4wMUPv9uSXdJetsw47Y25mRvLSUivgesBf53hfUHI+KTEfFm4ALgv0s6t6jISmCWpC8NsptdwPci4vii14SI+Hjax0MRsYRCE8/twLcGdp/pwx3ts2mbp0fEsRSu3FW0vnR/u4CbSuIeHxHXpLjvjojzKHzzeQL4Wp3jtVHMyd5a0ZeB8yTNK10h6XxJvy5JwAHglfQacBBYBJwt6ZoK278TeKukSySNTa/flHSKpGMkfUjScRHxctE+oHCFPknScXX6nBOBF4HnU3v8p0rW7wHeXDT/DeACSe+VNEbSuHQ76AxJUyRdmNru/y1t9xXMEid7azmp/XodUK4j1Rzg/1FIZj8E/jYi7iup/zxwHrBY0lH350fEQeA9wFKgH/gZ8FfAf0hFLgF6U9PKSgpX3ETEE8B6YGdqRpmW7ZPy58CZwAsU/vF6a8n6zwL/M+1rVUTsApZQ+Kf0PgpX+p+i8Hv8axT+ed0PPEeh7f+/ZIzP2og8eImZWfvzlb2ZWQ442ZuZ5YCTvZlZDjjZm5nlwOuaHUA5kydPjs7OzmaHYWY2ajz88MPPRkRHpfUtmew7Ozvp6elpdhhmZqOGpKcHW+9mHDOzHHCyNzPLASd7M7McaMk2+3Jefvll+vr6OHToULNDqbtx48YxY8YMxo4d2+xQzKxNjZpk39fXx8SJE+ns7KTwDKz2EBHs37+fvr4+TjrppGaHY2ZtatQ04xw6dIhJkya1VaIHkMSkSZPa8huLmbWOUZPsgbZL9APa9XOZWesYVcnezMxqM2ra7Et1dq7h6acP1G17s2cfS2/vikHLSOLiiy/mpptuAuDw4cNMnTqVBQsWcOedd7J27Vp6enpYvXp1SaydTJw4kTFjxgBw9tlnc+2119YtdjOzoYzaZP/00weIWFW37UllR8E7wvjx49m+fTsvvfQSr3/969m8eTPTp0+vavv33nsvkycPNryo2ZGKL2iquRipd31rL27GGabFixdz1113AbB+/XqWLVvW5IisXQ1c0ESsqulbbNb61l6c7Idp6dKlbNiwgUOHDrFt2zYWLFhQVb1zzjmHefPmMW/ePL70pcHGwjYzq79R24zTLKeffjq9vb2sX7+e7u7uquu5GcfMmsnJvgYXXnghq1at4r777mP//v3NDsfMbEhO9jX4yEc+wnHHHcfb3/527rvvvmaHY2Y2pFGb7GfPPraqO2iGs71qzZgxgyuuuKLsurVr13L77be/Nv/ggw8ChTb7gVsvTz/9dNatW5chWjOz4Rky2UuaCawD3gS8CqyJiL+WdALwTaAT6AXeHxE/L1N/EfDXwBjg6xFxTT0Cb8ZtZC+++OJRyxYuXMjChQsBWL58OcuXLz+qTG9v78gGZmY2hGruxjkMfDIiTgHeCXxC0lzgKuCeiJgD3JPmjyBpDPAVYDEwF1iW6pqZWQMNmewjYndEbEnTB4HHgenAEuDGVOxG4H1lqs8HdkTEzoj4FbAh1TMzswYaVpu9pE7gHcCPgCkRsRsKfxAknVimynRgV9F8H1D2xnRJK4AVALNmzSq7/4hoy4eGRUSzQ7AWUdrrtZoy7hlr1ag62UuaANwCXBkRB6pMuuUKlc1sEbEGWAPQ1dV1VJlx48axf//+tnvM8cDz7MeNG9fsUKwFVPMYkOIy9bxJwdpbVcle0lgKif7miLg1Ld4jaWq6qp8K7C1TtQ+YWTQ/A+ivJdAZM2bQ19fHvn37aqne0gZGqjIzGynV3I0j4Hrg8Yj4YtGqO4BLgWvS+z+Uqf4QMEfSScAzwFLgg7UEOnbsWI/kZGZWo2ruxjkLuAT4XUlb06ubQpI/T9JTwHlpHknTJG0EiIjDwOXA3RT+sfutiHhsBD6HmZkNYsgr+4j4PuXb3gHOLVO+H+gumt8IbKw1QDMzy85PvTQzywEnezOzHHCyNzPLASd7M7MccLI3M8sBJ3szsxxwsjczywEnezOzHHCyNzPLASd7M7MccLI3M8sBJ3szsxxwsjczywEnezOzHHCyNzPLgWpGqroBOB/YGxGnpWXfBE5ORY4Hno+IeWXq9gIHgVeAwxHRVae4zaxGxQOWF/Pg5e2tmjFo1wKrgXUDCyLiAwPTkr4AvDBI/XMi4tlaAzSz+qo0qLkHL29v1YxUdb+kznLr0vi07wd+t75hmZlZPWVts/9tYE9EPFVhfQCbJD0syd8PzcyapJpmnMEsA9YPsv6siOiXdCKwWdITEXF/uYLpj8EKgFmzZmUMy8zMitV8ZS/pdcDvA9+sVCYNPk5E7AVuA+YPUnZNRHRFRFdHR0etYZmZWRlZmnHeDTwREX3lVkoaL2niwDTwHmB7hv2ZmVmNhkz2ktYDPwROltQn6aNp1VJKmnAkTZO0Mc1OAb4v6RHgn4G7IuI79QvdzMyqVc3dOMsqLF9eZlk/0J2mdwJnZIzPzMzqwD1ozcxyIOvdOGbWJmbPPva1jlWVetMW9751j9vRxcnezACOSNyVetMW9751j9vRxc04ZmY54GRvZpYDTvZmZjngZG9mlgNO9mZmOeBkb2aWA072ZmY54GRvZpYD7lRlNoIq9Tgdbk/U0t6tw5W1vo1+TvZmI6hSj9Ph9kTN+lgCP9bA3IxjZpYDTvZmZjngZG9mlgPVjFR1g6S9krYXLbta0jOStqZXd4W6iyQ9KWmHpKvqGbiZmVWvmiv7tcCiMsu/FBHz0mtj6UpJY4CvAIuBucAySXOzBGtmZrUZMtlHxP3AczVsez6wIyJ2RsSvgA3Akhq2Y2ZmGWVps79c0rbUzPPGMuunA7uK5vvSsrIkrZDUI6ln3759GcIyM7NStSb7rwJvAeYBu4EvlCmjMsui0gYjYk1EdEVEV0dHR41hmZlZOTUl+4jYExGvRMSrwNcoNNmU6gNmFs3PAPpr2Z+ZmWVTUw9aSVMjYnea/U/A9jLFHgLmSDoJeAZYCnywpijNrCqlj2EYDTyIeWMMmewlrQcWApMl9QF/BiyUNI9Cs0wv8LFUdhrw9YjojojDki4H7gbGADdExGMj8inMDDjyMQyjhQcxb4whk31ELCuz+PoKZfuB7qL5jcBRt2WamVljuQetmVkOONmbmeWAk72ZWQ442ZuZ5YCTvZlZDjjZm5nlgJO9mVkOeAxasybzYODWCE72Zk3mxwNYI7gZx8wsB5zszcxywMnezCwHnOzNzHLAyd7MLAec7M3McsDJ3swsB4ZM9pJukLRX0vaiZZ+X9ISkbZJuk3R8hbq9kh6VtFVSTz0DNzOz6lVzZb8WWFSybDNwWkScDvwL8MeD1D8nIuZFRFdtIZqZWVbVDEt4v6TOkmWbimYfBC6qb1hmo0s7D5o9Ggcxt6PV43EJHwG+WWFdAJskBfB/ImJNpY1IWgGsAJg1a1YdwjJrnHYeNHs0DmJuR8v0D1pJnwEOAzdXKHJWRJwJLAY+IensStuKiDUR0RURXR0dHVnCMjOzEjUne0mXAucDH4qIKFcmIvrT+17gNmB+rfszM7Pa1ZTsJS0C/gi4MCJ+WaHMeEkTB6aB9wDby5U1M7ORVc2tl+uBHwInS+qT9FFgNTAR2Jxuq7wulZ0maWOqOgX4vqRHgH8G7oqI74zIpzAzs0FVczfOsjKLr69Qth/oTtM7gTMyRWdmZnXhHrRmZjngZG9mlgNO9mZmOeAxaM1q1M49Sz0IevtxsjerUTv3LG2nxz1YgZtxzMxywMnezCwHnOzNzHLAyd7MLAec7M3McsDJ3swsB5zszcxywMnezCwH3KnK2lY148K289ixI624l23pch/H1uNkb22rmnFh23ns2JFWKaH7OLamagYvuUHSXknbi5adIGmzpKfS+xsr1F0k6UlJOyRdVc/AzcysetW02a8FFpUsuwq4JyLmAPek+SNIGgN8hcJg43OBZZLmZorWzMxqMmSyj4j7gedKFi8BbkzTNwLvK1N1PrAjInZGxK+ADamemZk1WK1340yJiN0A6f3EMmWmA7uK5vvSsrIkrZDUI6ln3759NYZlZmbljOStlyqzLCoVjog1EdEVEV0dHR0jGJaZWf7Umuz3SJoKkN73linTB8wsmp8B9Ne4PzMzy6DWZH8HcGmavhT4hzJlHgLmSDpJ0jHA0lTPzMwarJpbL9cDPwROltQn6aPANcB5kp4CzkvzSJomaSNARBwGLgfuBh4HvhURj43MxzAzs8EM2akqIpZVWHVumbL9QHfR/EZgY83RmZlZXbgHrVkZWR6j4MG6rRU52ZuVkeUxCn4ujLUiP/XSzCwHnOzNzHLAyd7MLAec7M3McsDJ3swsB5zszcxywMnezCwHnOzNzHLAnaosd0p7xw6o1PN1sIG1hyOPPWs9oHvrcLK33CnuHVusUiKqV4LKY6LzgO6tw804ZmY54GRvZpYDTvZmZjngZG9mlgM1J3tJJ0vaWvQ6IOnKkjILJb1QVOZPs4dsZmbDVfPdOBHxJDAPQNIY4BngtjJFH4iI82vdj5mZZVevZpxzgZ9ExNN12p6ZmdVRvZL9UmB9hXXvkvSIpG9LOrXSBiStkNQjqWffvn11CsvMzKAOyV7SMcCFwN+XWb0FmB0RZwB/A9xeaTsRsSYiuiKiq6OjI2tYZmZWpB5X9ouBLRGxp3RFRByIiBfT9EZgrKTJddinmZkNQz2S/TIqNOFIepMkpen5aX/767BPMzMbhkzPxpH0BuA84GNFy1YCRMR1wEXAxyUdBl4ClkZEZNmnmZkNX6ZkHxG/BCaVLLuuaHo1sDrLPszMLDv3oDUzywEnezOzHHCyNzPLASd7M7MccLI3M8sBJ3szsxzwGLTWdNUMSl1cZrBylbTrYN+lg6G3wmfLcqxL6+Zx3N6R4mRvTVfNoNSlg4QPd/Dqdk0arfi5ssRUXNcDlNeXm3HMzHLAyd7MLAec7M3McsDJ3swsB5zszcxywMnezCwHnOzNzHIgU7KX1CvpUUlbJfWUWS9J10raIWmbpDOz7M/MzGpTj05V50TEsxXWLQbmpNcC4Kvp3axuSnvXDmiF3qRWXmmvaRt5I92DdgmwLg1F+KCk4yVNjYjdI7xfy5HS3rXW+vwza7ysbfYBbJL0sKRyfaSnA7uK5vvSMjMza6CsV/ZnRUS/pBOBzZKeiIj7i9arTJ2yA46nPxYrAGbNmpUxLDMzK5bpyj4i+tP7XuA2YH5JkT5gZtH8DKC/wrbWRERXRHR1dHRkCcvMzErUnOwljZc0cWAaeA+wvaTYHcCH01057wRecHu9mVnjZWnGmQLcJmlgO38XEd+RtBIgIq4DNgLdwA7gl8Bl2cI1M7Na1JzsI2IncEaZ5dcVTQfwiVr3YWZm9eEetGZmOeBkb2aWAx6W0Ealdh1Ttp1l+ZlVM05xM7Y1mjjZ26iUl1/QdpLlZ1bNOMXN2NZo4mYcM7MccLI3M8sBJ3szsxxwsjczywEnezOzHHCyNzPLASd7M7MccLI3M8sBJ3szsxxwD1qri3p1QS/tUu+esvlV6fEK1Z4jeX0sQiVO9lYX9eqCXvwLmaeu7Ha0Ssm52nMkr49FqMTNOGZmOZBlWMKZku6V9LikxyRdUabMQkkvSNqaXn+aLVwzM6tFlmacw8AnI2JLGov2YUmbI+LHJeUeiIjzM+zHzMwyqvnKPiJ2R8SWNH0QeByYXq/AzMysfurSZi+pE3gH8KMyq98l6RFJ35Z06iDbWCGpR1LPvn376hGWmZklmZO9pAnALcCVEXGgZPUWYHZEnAH8DXB7pe1ExJqI6IqIro6OjqxhmZlZkUzJXtJYCon+5oi4tXR9RByIiBfT9EZgrKTJWfZpZmbDl+VuHAHXA49HxBcrlHlTKoek+Wl/+2vdp5mZ1SbL3ThnAZcAj0rampZ9GpgFEBHXARcBH5d0GHgJWBoRkWGfVmeVehlm6X1Yr16wHlTchlJ8jpRbZ/+u5mQfEd8HNESZ1cDqWvdhI69SL8MsvQ/r1Qs2793bbWg+R6rnHrRmZjngZG9mlgNO9mZmOeBkb2aWA072ZmY54GRvZpYDTvZmZjngZG9mlgNtNyxhHsedHKnxX+uxr0q9aUu3YzaS6tUbeyTyS6NyVtsl+zyOOzkS47/Wa1+VetMWb8dspNUrgY5EfmlUznIzjplZDjjZm5nlgJO9mVkOONmbmeWAk72ZWQ442ZuZ5UDWMWgXSXpS0g5JV5VZL0nXpvXbJJ2ZZX9mZlabLGPQjgG+AiwG5gLLJM0tKbYYmJNeK4Cv1ro/MzOrXZYr+/nAjojYGRG/AjYAS0rKLAHWRcGDwPGSpmbYp5mZ1UC1jv8t6SJgUUT8YZq/BFgQEZcXlbkTuCaNV4uke4A/ioieMttbQeHqH+BkYD/wbE3BNcZkHF8Wji8bx5dNO8Y3OyI6Kq3M8riEcoONl/7lqKZMYWHEGmDNaxWlnojoqj28keX4snF82Ti+bPIYX5ZmnD5gZtH8DKC/hjJmZjbCsiT7h4A5kk6SdAywFLijpMwdwIfTXTnvBF6IiN0Z9mlmZjWouRknIg5Luhy4GxgD3BARj0lamdZfB2wEuoEdwC+By4axizVDF2kqx5eN48vG8WWTu/hq/getmZmNHu5Ba2aWA072ZmY50JBkP9RjFVKZhZK2SnpM0vfSspmS7pX0eFp+RVH5qyU9k+psldTd6PjS8l5Jj6Z1PUXLT5C0WdJT6f2NjY5P0slFx2erpAOSrkzrGnb8JH2qaD/bJb0i6YTB6jby+FWKr1XOvyGOX9PPv0GOX6ucf8dJ+kdJj6Sf42VD1W3w8SsbX93Pv4gY0ReFf97+BHgzcAzwCDC3pMzxwI+BWWn+xPQ+FTgzTU8E/mWgLnA1sKqZ8aXpXmByme1+DrgqTV8F/FUz4ivZzs8odLxo6PErKX8B8N2h6jby+A0SX0ucf5Xia5Xzb7D4WuH8Az498PmBDuC5VLYlzr9B4qvr+deIK/tqHqvwQeDWiPhXgIjYm953R8SWNH0QeByY3irxDWEJcGOavhF4X5PjOxf4SUQ8XWMcWeIrtgxYX0XdRh6/svG10PlXNr4hNP34lWjm+RfAREkCJlBIpoeHqNvI41c2vnqff41I9tOBXUXzfRwd8FuBN0q6T9LDkj5cuhFJncA7gB8VLb5chadp3pDha1bW+ALYlJYXj2o8JVKfgvR+YpPiG7CUo38JG3X8AJD0BmARcEsVdRt5/CrFV7yuk+adf4PF1wrn32DxDWjm+bcaOIVCh85HgSsi4tUh6jby+FWK7zX1OP8akeyreWTC64DfAH4PeC/wJ5Le+toGpAkUTqArI+JAWvxV4C3APGA38IUmxXdWRJxJ4Qmfn5B0do1xjFR8qNDp7ULg74vqNPL4DbgA+KeIeK6GurXKEl9hA80//waLrxXOv8Hia4Xz773AVmBa2t9qScdWWTerLPEVNlCn868Ryb7axyp8JyJ+ERHPAvcDZwBIGkvhg94cEbcOVIiIPRHxSvoL+DUKX5caHl9E9Kf3vcBtRXHsUXrCZ3qvpumn7vEli4EtEbFnYEGDj9+A0qu7weo28vhViq9Vzr+K8bXI+VcxvqTZ599lFJo5IyJ2AD8F3jZE3UYev0rx1ff8q+WfDsN5Ubjq3AmcxL//g+LUkjKnAPeksm8AtgOnUfiruA74cpntTi2a/m/AhibENx6YmMqMB35A4UmgAJ/nyH/wfK7R8RWt3wBc1qzjl8odR6Etcnw1dRt5/AaJryXOv0Hia4nzr1J8rXL+UbgKvjpNTwGeofBUyZY4/waJr67n37CDr/EDd1P4T/JPgM+kZSuBlUVlPkXhjpLtFL6uAPwWha882yh8zdkKdKd1N1Fo39pG4Rk8U5sQ35vTD+8R4LGBumndJAoJ+Kn0fkKj40vL30DhcdHHlWyz0cdvebkTslzdJh2/o+JrsfOvXHytdP5V+vk2/fyj0DyyKe1vO3BxK51/leKr9/nnxyWYmeWAe9CameWAk72ZWQ442ZuZ5YCTvZlZDjjZm5nlgJO9mVkOONmbmeXA/wf9LiNPWOxVEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.hist(risk_mle, bins=100, label='MLE', density=True, histtype='step', color='darkblue')\n",
    "plt.legend(loc='upper left')\n",
    "plt.title('Risk estimates')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Parametrizing and fitting log-linear models\n",
    "\n",
    "A log-linear model $M_S$ has density of the form:\n",
    "\n",
    "$$ \\log f(x) = \\sum_{A \\subset S} \\psi_A(x) $$\n",
    "\n",
    "where $\\psi_A$ depends on $x_i$ only if $i \\in A$.  Since x = $(x_1, \\dots, x_k)$ has its components with values $x_i \\in \\{ 0, 1 \\}$, we can express that as\n",
    "\n",
    "$$ \\psi_A(x) = \\sum_{B \\subset A} c_{B, A} g_A(x) = \\sum_{B \\subset A} c_{B, A} \\left( \\prod_{j \\in B} x_j \\right) $$\n",
    "\n",
    "for some constants $c_{B, A}$, and where $g_A(x) = \\prod_{j \\in A} x_j$. We can then write the log-likelihood for a parameter $\\beta = \\left\\{ \\beta_A | A \\subset S \\right\\}$ as the sum of the log density evaluated at each observation:\n",
    "\n",
    "$$ \\ell(\\beta) = \\sum_{x \\in \\text{obs}} \\log f(x; \\beta) =  \\sum_{x \\in \\text{obs}} \\sum_{A \\subset S} \\sum_{B \\subset A} c_{B, A} g_A(x) = \\sum_{x \\in \\text{obs}} \\sum_{A \\subset S} \\beta_A g_A(x) = \\sum_{A \\subset S} \\beta_A \\sum_{x \\in \\text{obs}} g_A(x) $$\n",
    "\n",
    "We will want to find the MLE $\\hat{\\beta}$ subject to the constraint that $f(x)$ is a probability density function which adds up to one, that is, $\\sum_{x \\in \\Omega} f(x) = 1$ over the random variable universe $\\Omega = \\{ 0, 1 \\}^k$.\n",
    "\n",
    "We can now explicitly frame this as an optimization problem, with target being the maximization of the log likelihood, and a nonlinear equality constraint that the density function adds up to 1:\n",
    "\n",
    "$$\n",
    "\\hat{\\beta} = \\text{argmax}_\\beta \\left\\{ \\sum_{A \\subset S} \\beta_A \\sum_{x \\in \\text{obs}} g_A(x) \n",
    "\\;\\Bigg|\\;  \\sum_{x \\in \\Omega} \\exp \\left( \\sum_{A \\subset S} \\beta_A g_A(x) \\right) = 1\n",
    "\\right\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matrix representation\n",
    "\n",
    "We can use matrix multiplications to speed up the function and constraint evaluations during optimization.  If $X$ is a matrix of $n$ samples by $k$ dimensions, containing $n$ observations, we define the matrix $h_S(X)$, where\n",
    "\n",
    "$$ (h_S(X))_{i, j} = \\begin{cases}\n",
    "1 & \\text{if } X_{r}^j = 1 \\text{ for all } r \\in A \\text{, where the } j \\text{-th element of } X \\text{ is } X^j = (X^j_1, \\dots, X^j_k) \\text{ and } A \\text{ is the } i \\text{-th subset of } S \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}$$\n",
    "\n",
    "By definition, the log-likelihood then becomes the sum of the elements of the matrix product $\\beta \\cdot h_S(X_\\text{obs})$:\n",
    "\n",
    "$$ \\ell(\\beta) = \\sum_{A \\in S} \\sum_{x \\in \\text{obs}} \\beta_A g_A(x) = \\sum_i \\left( \\beta \\cdot h_S(X_\\text{obs}) \\right)_i $$\n",
    "\n",
    "The constraint can also be computed using a matrix multiplication, first defining a matrix $X_\\Omega$ of shape $(2^k, k)$ that contains one row for each possible observation:\n",
    "\n",
    "$$ (X_\\Omega)_{i, j} = \\begin{cases}\n",
    "1 & \\text{if the } j \\text{-th digit in the binary representation of } i \\text{ is } 1 \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}$$\n",
    "\n",
    "Then the sum of the density functions is the sum of the element-wise exponentials resulting from a matrix product $\\beta \\cdot h_S(X_\\Omega)$:\n",
    "\n",
    "$$ \\sum_{x \\in \\Omega} f(x) = \\sum_{x \\in \\Omega} \\exp \\left(\\sum_{A \\subset S} \\beta_A g_A(x) \\right) \n",
    "= \\sum_{i} \\exp \\left( \\beta \\cdot h_S(X_\\Omega) \\right)_i$$\n",
    "\n",
    "Now, after pre-computing $h_S(X_\\text{obs})$ and $h_S(X_\\Omega)$, an optimizer will only need to perform sums, matrix multiplications and exponentiations to compute the target function and constraint penalty at each step.\n",
    "\n",
    "Finally, note that $\\sum_{x \\in \\Omega} = 1$ if and only if $\\log \\sum_{x \\in \\Omega} f(x) = 0$, so we can use scipy's `logsumexp` function to represent the total probability density function constraint."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Model selection and bootstraping\n",
    "\n",
    "**Model selection** is to be performed using forward selection with AIC or BIC scores.  This means that, for each observation set $X$, given a scoring process to pick between models, we will:\n",
    "\n",
    "1. Fit the uniform model $M_\\varnothing$;\n",
    "2. Given the current best model $M_{S_0}$, fit all models $M_{S}$ where $S = S_0 \\cup \\{ i \\}$ for each $i$ not in $S_0$;\n",
    "3. If a new model was found as an improvement, select it as the new best model and return to step 2.  Otherwise, keep the current best model and stop.\n",
    "\n",
    "**AIC** will be computed as the log likelihood of the candidate model at its MLE minus its number of parameters:\n",
    "\n",
    "$$ \\text{AIC}(M_S) = \\hat{\\ell}_S - |M_S| $$\n",
    "\n",
    "**BIC** will be computed as the log likelihood of the candidate model at its MLE minus the BIC penalty on this formulation (half times number of parameters times log sample size):\n",
    "\n",
    "$$ \\text{BIC}(M_S) = \\hat{\\ell}_S - \\frac{1}{2} |M_S| \\log n $$\n",
    "\n",
    "Once a model is selected, we can compute its loss $D(f, \\hat{f})$, the KL divergence over the universe of observations for the fitted density function.\n",
    "\n",
    "**Bootstraping** means repeating this process multiple times, observing the loss at each point and reporting its average and distribution.  Note that the model selection itself is a part of the bootstraping process -- different model formulations could potentially be chosen for different observations, even though they come from the same synthetic data generation process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "from scipy.special import logsumexp\n",
    "\n",
    "def get_loglinear_mle(X, subsets):\n",
    "    \"\"\"\n",
    "    Estimates a loglinear model for observations X given variables S via maximizing the likelihood estimator.\n",
    "    \n",
    "    Args:\n",
    "       X:  2D array, shape (n, k), of samples and observations (0 or 1)\n",
    "       S:  iterable of variables between 0 and k-1 inclusive\n",
    "       \n",
    "    Returns:\n",
    "       beta_hat:  1D array of size 2**|S| with estimated parameters via MLE\n",
    "       log_likelihood:  Value of the log-likelihood of the model estimated at the MLE\n",
    "    \"\"\"\n",
    "    k = X.shape[1]\n",
    "    \n",
    "    n_subsets = len(subsets)\n",
    "    \n",
    "    # Speed up calculation of gA operations with vector operations\n",
    "    def get_h(XX):\n",
    "        \"\"\"\n",
    "        Calculate the matrix gA(XX) of shape (2**|S|, XX.shape[0]), where\n",
    "        gA(XX){i, j} = 1 if all elements of the i-th subset A of S are in the j-th sample of XX\n",
    "                       0 otherwise\n",
    "        \"\"\"\n",
    "        h = np.zeros(shape=(len(subsets), XX.shape[0]), dtype=int)\n",
    "        for i, A in enumerate(subsets):\n",
    "            h[i] = XX[:, A].all(axis = 1)\n",
    "            \n",
    "        return h\n",
    "    \n",
    "    h_obs = get_h(X)\n",
    "    h_universe = get_h(generate_universe(k))\n",
    "    \n",
    "    def neg_log_likelihood(beta):\n",
    "        return -np.sum(beta @ h_obs)\n",
    "    \n",
    "    def log_density_sum(beta):\n",
    "        # Use scipy's logsumexp to avoid overflows\n",
    "        # exp(sum(beta @ h_universe)) - 1 == 0 iff logsumexp(beta @ h_universe) == 0\n",
    "        return logsumexp(beta @ h_universe)\n",
    "    \n",
    "    # Constraint: PDF adds up to 1\n",
    "    pdf_constraint = {'type': 'eq', 'fun': log_density_sum}\n",
    "    \n",
    "    # Get initial guess: all values zero other than first\n",
    "    beta0 = np.zeros(len(subsets))\n",
    "    beta0[0] = -k * np.log(2)\n",
    "    \n",
    "    res = minimize(neg_log_likelihood, beta0, constraints=[pdf_constraint])\n",
    "    beta_hat = res.x\n",
    "    log_likelihood = -res.fun\n",
    "    \n",
    "    return beta_hat, log_likelihood\n",
    "\n",
    "\n",
    "def f_loglinear(subsets, beta):\n",
    "    \"\"\"\n",
    "    Computes the density function for a given set of variables S and corresponding parameters beta.\n",
    "    \n",
    "    f(x) = exp ( \\sum_{A in S} \\beta(A) * g_A(x) )\n",
    "    \"\"\"\n",
    "    def f(x):\n",
    "        return np.exp(np.sum([x[A].all() * beta[i] for i, A in enumerate(subsets)]))\n",
    "    \n",
    "    return f\n",
    "\n",
    "\n",
    "def get_AIC(X, subsets):\n",
    "    \"\"\" \n",
    "    Calculates AIC using the loglinear model log likelihood function.\n",
    "    \n",
    "    Args:\n",
    "       X:        2D array (n, k), observed data for log-likelihood function\n",
    "       subsets:  iterable, list of subsets for log-likelihood function\n",
    "       \n",
    "    Returns:\n",
    "       AIC score for the given submodel:  ll - |subsets|\n",
    "    \"\"\"\n",
    "    _, log_likelihood = get_loglinear_mle(X, subsets)\n",
    "    penalty = len(subsets)\n",
    "    \n",
    "    return log_likelihood - penalty\n",
    "\n",
    "\n",
    "def get_BIC(X, subsets):\n",
    "    \"\"\" \n",
    "    Calculates BIC using the loglinear model log likelihood function.\n",
    "    \n",
    "    Args:\n",
    "       X:        2D array (n, k), observed data for log-likelihood function\n",
    "       subsets:  iterable, list of subsets for log-likelihood function\n",
    "       \n",
    "    Returns:\n",
    "       BIC score for the given submodel:  ll - (|subsets| log n) / 2\n",
    "    \"\"\"\n",
    "    _, log_likelihood = get_loglinear_mle(X, subsets)\n",
    "    n = X.shape[0]\n",
    "    penalty = len(subsets) * np.log(n) / 2\n",
    "    \n",
    "    return log_likelihood - penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_selection(score_func, S):\n",
    "    \"\"\"\n",
    "    Uses forward selection to select a subset A of S, in a search to maximize score_func(A).\n",
    "\n",
    "    Args:\n",
    "       score_func:  (A) => score, a function to score subsets\n",
    "       S:           iterable to select a subset from\n",
    "       \n",
    "    Returns:\n",
    "       A:           a subset of S resulting from forward selection\n",
    "       \n",
    "    \"\"\"\n",
    "    all_subsets = [list(s) for s in powerset(S)]\n",
    "    current_subset = [[]]\n",
    "    current_score = score_func(current_subset)\n",
    "    \n",
    "    while True:\n",
    "        best_subset, best_score = current_subset, current_score\n",
    "        improved = False\n",
    "        for s in all_subsets:\n",
    "            if s not in current_subset:\n",
    "                candidate_subset = current_subset.copy()\n",
    "                candidate_subset.append(s)\n",
    "                candidate_subset.sort()\n",
    "                candidate_score = score_func(candidate_subset)\n",
    "                \n",
    "                if candidate_score > best_score:\n",
    "                    improved = True\n",
    "                    best_subset, best_score = candidate_subset, candidate_score\n",
    "        \n",
    "        if not improved:\n",
    "            break\n",
    "        current_subset, current_score = best_subset, best_score\n",
    "    \n",
    "    return current_subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "975e6be47d5047888e89d4ba20e3bb17",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='AIC', max=500, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "AIC mean: inf\n"
     ]
    }
   ],
   "source": [
    "# AIC bootstrapping\n",
    "\n",
    "B = 500\n",
    "\n",
    "n = 1000\n",
    "k = 7\n",
    "\n",
    "S_full = [i for i in range(k)]\n",
    "X_universe = generate_universe(k=k)\n",
    "\n",
    "n_jobs = max(multiprocessing.cpu_count() - 1, 1)\n",
    "\n",
    "def bootstrap_step_aic(i):\n",
    "    XX = generate_samples(n=n, k=k)    \n",
    "    score_func = lambda S: get_AIC(XX, S)\n",
    "    best_subset = forward_selection(score_func, S_full)\n",
    "    beta_hat, _ = get_loglinear_mle(XX, best_subset)\n",
    "    f_hat = f_loglinear(best_subset, beta_hat)\n",
    "    return KL_divergence(true_density, f_hat, X_universe)\n",
    "\n",
    "with tqdm_joblib(tqdm_notebook(desc=\"AIC\", total=B)) as progress_bar:\n",
    "    risk_aic = np.array(Parallel(n_jobs=n_jobs)(delayed(bootstrap_step_aic)(i) for i in range(B)))\n",
    "\n",
    "print('AIC mean: %.2f' % risk_aic.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "32d9245694784f25990dbe517f3c6ce0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='BIC', max=500, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "BIC mean: 1.18\n"
     ]
    }
   ],
   "source": [
    "# BIC bootstrapping\n",
    "\n",
    "B = 500\n",
    "\n",
    "n = 1000\n",
    "k = 7\n",
    "\n",
    "S_full = [i for i in range(k)]\n",
    "X_universe = generate_universe(k=k)\n",
    "\n",
    "n_jobs = max(multiprocessing.cpu_count() - 1, 1)\n",
    "\n",
    "def bootstrap_step_bic(i):\n",
    "    XX = generate_samples(n=n, k=k)    \n",
    "    score_func = lambda S: get_BIC(XX, S)\n",
    "    best_subset = forward_selection(score_func, S_full)\n",
    "    beta_hat, _ = get_loglinear_mle(XX, best_subset)\n",
    "    f_hat = f_loglinear(best_subset, beta_hat)\n",
    "    return KL_divergence(true_density, f_hat, X_universe)\n",
    "\n",
    "with tqdm_joblib(tqdm_notebook(desc=\"BIC\", total=B)) as progress_bar:\n",
    "    risk_bic = np.array(Parallel(n_jobs=n_jobs)(delayed(bootstrap_step_bic)(i) for i in range(B)))\n",
    "\n",
    "print('BIC mean: %.2f' % risk_bic.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MLE mean: 0.75\n",
      "AIC mean: 2.55\n",
      "BIC mean: 1.18\n"
     ]
    }
   ],
   "source": [
    "# Remove samples with unsuccessful optimization\n",
    "risk_aic = risk_aic[np.isfinite(risk_aic)]\n",
    "\n",
    "print('MLE mean: %.2f' % risk_mle.mean())\n",
    "print('AIC mean: %.2f' % risk_aic.mean())\n",
    "print('BIC mean: %.2f' % risk_bic.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this scenario, the best risk is obtained by the MLE, followed by the BIC and the AIC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArkAAAHiCAYAAADs/9QdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdE0lEQVR4nO3df7DldX3f8dfbhXT5JSpsqLLmXtJGQNCsZKPZaMkvTYQQSCIk60QkSRlrLQGbkExM2xGnGZO2mUg6Je1Qk4DB7k5HlEnVRukAlkQFdwkhq0tmMno3bPi1bCI/IkTET/+4h+2yP9hz7t5zz7mf+3jM3HHvPd97zvuMh73P+9nP+X6rtRYAAOjJCyY9AAAALDaRCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QusGJV1X+rqn83xHFzVfXGpZhpr8f8tar64FI+JkBPRC7QrUGcPllVT1TVg1V1XVUd++ztrbV3ttb+/SRnTJKq+v6q2rn311pr72+tXTqmx1vyaAdYaiIX6N2PtdaOTbIuyWuSvGfC8wCwBEQusCK01h5M8qnMx26SZLCy++uDP59YVR+vqq9W1d9W1e1Vtd/fkVV1WlV9pao2HuhxBrffPLiPv6yqn9rrtnOr6ktV9XhV/U1VXVlVxyT530leNlhxfqKqXlZVV1XVDYPvm62qVlU/V1X3VdXfVdU7q+q7q+qewcz/Za/H+SdVdUtV7a6qR6rqw1X1osFtf5jk25L8r8Fj/crg699TVZ8d3NefV9X373V/P1tVXx7M/ZWq+pnD+L8CYEmIXGBFqKq1Sc5J8lcHOeSXkuxMsibJSUl+LclzrnteVWcl+XSSX2itbT7AYxyT5OYk/yPJtyZ5a5LfraozBof8XpJ/0Vo7LsmZSW5prf39YK77W2vHDj7uP8iMr0vyHUl+OsnVSf5NkjcmOSPJT1XV9z07SpLfSPKyJKcneXmSq5KktXZxkr/OYIW7tfYfq+rkJJ9I8utJXpLkyiQ3VtWawXP6z0nOGcz9vUnuPsh8AFND5AK9u6mqHk9yX5KHk7z3IMc9neSlSWZaa0+31m5vre0duf8syR8luaS19vGD3Md5SeZaa3/QWvtGa+2uJDcmuXCvx3hlVb2wtfZ3g9tH8e9ba0+11j6d5O+TbGqtPdxa+5skt2d+O0Zaa3/VWru5tfYPrbVdSX47yfcd/G7ztiSfbK19srX2zdbazUm2JDl3cPs3k5xZVUe11h5orX1xxLkBlpzIBXr344MVyO9PclqSEw9y3H/K/Crvpwf/NP+r+9z+ziSfba3d+jyPNZPkdYN/8v9qVX01yc8k+ceD29+S+XDcUVWfqaoNIz6Xh/b685MH+PzYJKmqb62qzYMtEY8luSEHf97Pzn3RPnO/IclLByvNP5355/9AVX2iqk4bcW6AJSdygRWhtfaZJNcl+a2D3P54a+2XWmvfnuTHkvxiVf3QXoe8M8m3VdUHnudh7kvymdbai/b6OLa19i8Hj/GF1toFmd/KcFOS//nswx/Wk9vfbwzu89WttRdmfqW29rp938e7L8kf7jP3Ma213xzM/anW2psyv9J9b5L/vsjzAiw6kQusJFcneVNVrdv3hqo6r6r+aVVVkseSPDP4eNbjSd6c5Oyq+s2D3P/Hk7yiqi6uqiMHH99dVadX1bdU1c9U1fGttaf3eoxkfkX2hKo6fpGe53FJnkjy1cF+21/e5/aHknz7Xp/fkOTHqupHqmpVVa0enNZsbVWdVFXnD/bm/sPgfp8JwJQTucCKMdif+qEkB7oAxHck+T+Zj7jPJfnd1tpt+3z/V5O8Kck5VbXf+XVba48n+eEkG5Pcn+TBJP8hyT8aHHJxkrnBFoJ3Zn6FNa21e5NsSvLlwXaBlx3eM837kpyV5NHMv6Hso/vc/htJ/u3gsa5srd2X5ILMv9luV+ZXdn858z8jXpD5N+Xdn+RvM7+3912HOR/A2NVz31cBAADLn5VcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6c8Q47vTEE09ss7Oz47hrAABIkmzduvWR1tqaA902lsidnZ3Nli1bxnHXAACQJKmqHQe7zXYFAAC6I3IBAOiOyAUAoDtj2ZN7IE8//XR27tyZp556aqkecsmsXr06a9euzZFHHjnpUQAAyBJG7s6dO3PcccdldnY2VbVUDzt2rbXs3r07O3fuzCmnnDLpcQAAyBJuV3jqqadywgkndBW4SVJVOeGEE7pcoQYAWK6WdE9ub4H7rF6fFwDAcrWi3nhWVbn44ov3fP6Nb3wja9asyXnnnZckue6663LZZZft932zs7N51atelXXr1mXdunW5/PLLl2xmAABGt2R7cvc1O5vsOOjpe0c3M5PMzT3/Mcccc0y2bduWJ598MkcddVRuvvnmnHzyyUPd/6233poTTzzx8AcFAGDsJraSu2NH0trifQwbzOecc04+8YlPJEk2bdqUt771rWN8lgAATMKK2q6QJBs3bszmzZvz1FNP5Z577snrXve6ob7vB37gB/ZsV/jABz4w5ikBADgcE9uuMCmvfvWrMzc3l02bNuXcc88d+vtsVwAAWD5WXOQmyfnnn58rr7wyt912W3bv3j3pcQAAWGQrMnJ//ud/Pscff3xe9apX5bbbbpv0OAAALLIVGblr167NFVdcccDbrrvuutx00017Pv/85z+fZH5P7qpVq5LMb3n40Ic+NP5BAQBYkGqtLfqdrl+/vm3ZsuU5X9u+fXtOP/30PZ9P4hRi47Tv8wMAYLyqamtrbf2BbpvYSu4kgxQAgL6tuFOIAQDQP5ELAEB3VuQbzxbqnnuSr3/9wLc98kjyylfu//VJ7xUGAFiJRO4Ivv71ZP0BtzYn27fPX154X1XjnQkAgP3ZrgAAQHdWXOR+7GMfS1Xl3nvvTZLMzc3lzDPP3HP7nXfembPPPjunnnpqTjvttFx66aX52te+NqlxAQBYgIltV5i9ejY7Hl28E+XOHD+TuXfPHfK4TZs25Q1veEM2b96cq6666jm3PfTQQ7nooouyefPmbNiwIa213HjjjXn88cdz9NFHL9qsAACM18Qid8ejO9Leu3gXoqj3HXrz6xNPPJE//dM/za233przzz9/v8i95pprcskll2TDhg3z91mVCy+8cNFmBABgaayo7Qo33XRT3vzmN+cVr3hFXvKSl+Suu+56zu3btm3Ld33Xd01oOgAAFsuKitxNmzZl48aNSZKNGzdm06ZNE54IAIBxWDGnENu9e3duueWWbNu2LVWVZ555JlWVd73rXXuOOeOMM7J169ZccMEFE5wUAIDDtWJWcj/ykY/k7W9/e3bs2JG5ubncd999OeWUU7Jz5849x1x22WW5/vrrc8cdd+z52g033JAHH3xwEiMDALBAKyZyN23alJ/4iZ94ztfe8pa35P3vf/+ez0866aRs3rw5V155ZU499dScfvrpuf322/PCF75wqccFAOAwTGy7wszxM0OdEWGU+3s+t912235fu/zyy3P55Zc/52sbNmzI7bffvmhzAQCw9CYWucOc0xYAABZixWxXAABg5RC5AAB0R+QCANAdkQsAQHdELgAA3VlRkbtq1aqsW7cu3/md35mzzjorn/3sZ5Mkc3NzOfPMM/ccd+edd+bss8/OqaeemtNOOy2XXnppvva1r01qbAAARjS5y/rOziY7dize/c3MJHNzz3vIUUcdlbvvvjtJ8qlPfSrvec978pnPfOY5xzz00EO56KKLsnnz5mzYsCGttdx44415/PHHkxy9ePMCADA2k4vcHTuS1hbv/mq0C0s89thjefGLX7zf16+55ppccskl2bBhw+BuKxdeeGGS5L77Dn9MAADGb3KROwFPPvlk1q1bl6eeeioPPPBAbrnllv2O2bZtWy655JIJTAcAwGJZUZG793aFz33uc3n729+ebdu2TXgqAAAW24p649neNmzYkEceeSS7du16ztfPOOOMbN26dUJTAQCwGFZs5N5777155plncsIJJzzn65dddlmuv/763HHHHXu+dsMNN+TBBx9c6hEBAFigFbVd4dk9uUnSWsv111+fVatWPeeYk046KZs3b86VV16Zhx9+OC94wQty9tln5yd/8icnMTIAAAswucidmRn5jAiHvL9DeOaZZw749dnZ2efszd2wYUNuv/32RRsNAIClNbnIPcQ5bQEAYKFW7J5cAAD6JXIBAOjOkkZuW8wrnE2RXp8XAMBytWSRu3r16uzevbu7IGytZffu3Vm9evWkRwEAYGDJ3ni2du3a7Ny5c7+LLywnjzySbN++/9dXr16dtWvXLv1AAAAc0JJF7pFHHplTTjllqR5uLF75yqSzhWgAgC554xkAAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0ZKnKr6l9X1ReraltVbaqq1eMeDAAAFuqQkVtVJye5PMn61tqZSVYl2TjuwQAAYKGG3a5wRJKjquqIJEcnuX98IwEAwOE5ZOS21v4myW8l+eskDyR5tLX26XEPBgAACzXMdoUXJ7kgySlJXpbkmKp62wGOe0dVbamqLbt27Vr8SQEAYEjDbFd4Y5KvtNZ2tdaeTvLRJN+770GttWtba+tba+vXrFmz2HMCAMDQhoncv07yPVV1dFVVkh9Ksn28YwEAwMINsyf3jiQfSXJXkr8YfM+1Y54LAAAW7IhhDmqtvTfJe8c8CwAALApXPAMAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6I3IBAOiOyAUAoDsiFwCA7ohcAAC6M1TkVtWLquojVXVvVW2vqg3jHgwAABbqiCGP+50kf9xau7CqviXJ0WOcCQAADsshI7eqXpjk7CQ/myStta8n+fp4xwIAgIUbZrvCtyfZleQPqurPquqDVXXMvgdV1TuqaktVbdm1a9eiDwoAAMMaJnKPSHJWkv/aWntNkr9P8qv7HtRau7a1tr61tn7NmjWLPCYAAAxvmMjdmWRna+2OwecfyXz0AgDAVDpk5LbWHkxyX1WdOvjSDyX50linAgCAwzDs2RV+IcmHB2dW+HKSnxvfSAAAcHiGitzW2t1J1o95FgAAWBSueAYAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3Rk6cqtqVVX9WVV9fJwDAQDA4RplJfeKJNvHNQgAACyWoSK3qtYm+dEkHxzvOAAAcPiGXcm9OsmvJPnmGGcBAIBFccjIrarzkjzcWtt6iOPeUVVbqmrLrl27Fm1AAAAY1TArua9Pcn5VzSXZnOQHq+qGfQ9qrV3bWlvfWlu/Zs2aRR4TAACGd8jIba29p7W2trU2m2Rjkltaa28b+2QAALBAzpMLAEB3jhjl4NbabUluG8skAACwSKzkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0Rucvd7GxSNfzH7OykJwYAGLsjJj0Ah2nHjqS14Y+vGt8sAABTwkouAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2Ru9LMzCRVw3/Mzk56YgCAkR0x6QFYYnNzox1fNZYxAADGyUouAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdOWTkVtXLq+rWqtpeVV+sqiuWYjAAAFioI4Y45htJfqm1dldVHZdka1Xd3Fr70phnAwCABTnkSm5r7YHW2l2DPz+eZHuSk8c9WC9mZpKq4T9mZyc9MQDA8jfMSu4eVTWb5DVJ7hjHMD2amxvt+KqxjAEAsKIM/cazqjo2yY1J3t1ae+wAt7+jqrZU1ZZdu3Yt5owAADCSoSK3qo7MfOB+uLX20QMd01q7trW2vrW2fs2aNYs5IwAAjGSYsytUkt9Lsr219tvjHwkAAA7PMCu5r09ycZIfrKq7Bx/njnkuAABYsEO+8ay19idJvB0KAIBlwxXPAADojsgFAKA7IhcAgO6IXAAAuiNyp8xXMjvadYBnZiY9MgDA1Bnpsr6M32x2JK1NegwAgGXNSi4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdEfkAgDQHZELAEB3RC4AAN0RuTy/mZmkaviP2dlJTwwAkCMmPQBTbm5utOOrxjIGAMAorOQCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHdELgAA3RG5AAB0R+QCANAdkQsAQHeOmPQA7K/eV0MfO3P8TObePTe+YQAAliGRO4Xae9vQx44SxAAAK4XtCgAAdEfkAgDQHZELAEB3RC4AAN0RuQAAdMfZFVaY2atns+PRHUMf7xRlAMByJHKXuZnjZ0Y+r65TlAEAvRO5y5xVVgCA/dmTCwBAd6zkAvZqA9AdkQtkx6M77NUGoCu2KwAA0B2RCwBAd0QuAADdEbkAAHRH5LK4ZmaSqtE+ZmcnPTUA0BlnV2Bxzc2N/j3lnfoAwOKykgsAQHes5MIyMOrFGkY1c/zM2O4bACZB5MIyMOrFGlYiV20DYG8iF+iCq7YBsDeRO2ajri5N21rdzPEzI8WA1TEAYBqI3DEb+Z+Zr5qu1aVRg9XqGAAwDUQuMHbjfuNc4s1zADyXyGUq2BLRN2+cA2CpiVymgjcMLS8L2asNAEtJ5AIjs5IOwLRzxTMAALpjJZdFNeo/YyfTd9o0WAwuTgEwWSKXRbWgH9J/MJvU8GH8leOTvHf4uxcbTIKLUwBMlshl8ubmRjp8doQgTsQGAKxEIpdlyTv7AYDnI3JZlpxzlcPlktWHZqsPsJyJXGBF6uGS1eOOUFt9gOVM5I7ZVz6Q5Krh/+Kfy0xmxzYNsMfsbLJj+EC878WrUpmubTIi9NCsRsPKJXLHbPbRJG34H0KnlFNqwZLYsWOk/zbXVi37bTLL/Up1owZrMv8cpukXAdENS0fkwhJb6A9qlpkxrxQno78ulnssjbpyPY16WH0X6iwXIheWWA8/qJfEiJGYmZmRT0c3VitwpXglWkjwTZuFPIflHuqsDENFblW9OcnvJFmV5IOttd8c61RAX0YN1mQ+WkeIxFEuKLIgMzOjPcbMEsTMcv9FoAM9/NLaw3OAAzlk5FbVqiTXJHlTkp1JvlBVf9Ra+9K4h4MDWkhsjPEHew8rOWM34qrmgow7Qscdh6PO/+z3TNMvAiPy387KMO7T9S10C9g4H8MWjekwzErua5P8VWvty0lSVZuTXJBk2UfuQhZBmAKjxsaIP9gX8hfy1K2CrMQVPvMvvhFfR6NecnsaVxBX6pvzxmncp+tbyOto9urZsf49P21bNJbiF4FpNEzknpzkvr0+35nkdeMZZ2mNvLg0O5vUaC+SuePjlGCTNuIK2dzMTDI3XT94Rzbqi3t2dryrfH5DXJ5GfB3NVi3rQEym7815XfzSPaKl+EVj3P8/j/ocxm0hr4uF/CIwbf/9VDvEX2BVdVGSH2mtXTr4/OIkr22t/cI+x70jyTsGn56a5C8Xf1xGdGKSRyY9BFPJa4OD8drgYLw2eD6Ten3MtNbWHOiGYVZydyZ5+V6fr01y/74HtdauTXLtgsZjLKpqS2tt/aTnYPp4bXAwXhscjNcGz2caXx8vGOKYLyT5jqo6paq+JcnGJH803rEAAGDhDrmS21r7RlVdluRTmT+F2O+31r449skAAGCBhjpPbmvtk0k+OeZZWHy2j3AwXhscjNcGB+O1wfOZutfHId94BgAAy80we3IBAGBZEbkdqqrfr6qHq2rbpGdhelTVy6vq1qraXlVfrKorJj0T06OqVlfVnVX154PXx/smPRPTpapWVdWfVdXHJz0L06Oq5qrqL6rq7qraMul59ma7Qoeq6uwkTyT5UGvtzEnPw3SoqpcmeWlr7a6qOi7J1iQ/7hLdJElVVZJjWmtPVNWRSf4kyRWttc9PeDSmRFX9YpL1SV7YWjtv0vMwHapqLsn61trUnUPZSm6HWmv/N8nfTnoOpktr7YHW2l2DPz+eZHvmr2gIafOeGHx65ODDKghJkqpam+RHk3xw0rPAsEQurEBVNZvkNUnumOwkTJPBP0ffneThJDe31rw+eNbVSX4lyTcnPQhTpyX5dFVtHVz9dmqIXFhhqurYJDcmeXdr7bFJz8P0aK0901pbl/krW762qmx3IlV1XpKHW2tbJz0LU+n1rbWzkpyT5F8NtkxOBZELK8hgr+WNST7cWvvopOdhOrXWvprktiRvnvAoTIfXJzl/sPdyc5IfrKobJjsS06K1dv/gfx9O8rEkr53sRP+fyIUVYvDGot9Lsr219tuTnofpUlVrqupFgz8fleSNSe6d7FRMg9bae1pra1trs0k2Jrmltfa2CY/FFKiqYwZvZE5VHZPkh5NMzZmdRG6HqmpTks8lObWqdlbVP5/0TEyF1ye5OPOrMHcPPs6d9FBMjZcmubWq7knyhczvyXWqKOD5nJTkT6rqz5PcmeQTrbU/nvBMeziFGAAA3bGSCwBAd0QuAADdEbkAAHRH5AIA0B2RCwBAd0QuAADdEbkAAHRH5AIA0J3/B1U3BNZu4dZAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "r = (0.5, 5.0)\n",
    "plt.figure(figsize=(12, 8))\n",
    "plt.hist(risk_mle, bins=50, range=r, label='MLE', density=True, histtype='step', color='blue')\n",
    "plt.hist(risk_aic, bins=50, range=r, label='AIC', density=True, histtype='step', color='green')\n",
    "plt.hist(risk_bic, bins=50, range=r, label='BIC', density=True, histtype='step', color='red')\n",
    "plt.legend(loc='upper left')\n",
    "plt.title('Risk estimates')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 18.8.3**.  Prove Lemma 18.5.\n",
    "\n",
    "A partition $(X_a, X_b, X_c)$ satisfies $X_b \\text{ ⫫ } X_c \\; | \\; X_a$ if and only if $f(x_a, x_b, x_c) = g(x_a, x_b) h(x_a, x_c)$ for some functions $g$ and $h$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  We have $X_b \\text{ ⫫ } X_c \\; | \\; X_a$ if and only if $\\mathbb{P}(X_b = b, X_c = c | X_a = a) = \\mathbb{P}(X_b = b | X_a = a) \\mathbb{P}(X_c = c | X_a = a)$ for every $a$.  \n",
    "\n",
    "Let $\\mathbb{P}(X_b = b, X_c = c | X_a = a) = f_a(b, c)$, $\\mathbb{P}(X_b = b | X_a = a) = g_a(b)$ and $\\mathbb{P}(X_c = c | X_a = a) = h_a(c)$.  Then the statement is equivalent to $f_a(b, c) = g_a(b) h_a(c)$ for all $a, b, c$.  The result follows by defining $f(a, b, c) = f_a(b, c)$, $g_a(b) = g(a, b)$ and $h_a(c) = h(a, c)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 18.8.4**. Prove Lemma 18.9.\n",
    "\n",
    "A graphical model is hierarchical but the reverse need not be true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "A log-linear model $\\log f(x) = \\sum_{A \\subset S} \\psi_A(x)$ is graphical if $\\psi_A(x) \\neq 0$ except for any pair of coordinates not in the edge set for some graph $\\mathcal{G}$.\n",
    "\n",
    "A log-linear model $\\log f(x) = \\sum_{A \\subset S} \\psi_A(x)$ is hierarchical if $\\psi_A = 0$ and $A \\subset B$ implies $\\psi_B = 0$.\n",
    "\n",
    "If, on a graphical model with associated graph $\\mathcal{G}$, $\\psi_A = 0$, then by definition there exists $\\{ i, j \\} \\subset A$ such that $(i, j)$ is not an edge of $\\mathcal{G}$.  But then, for any $B$ containing $A$, we also have that $\\{ i, j \\} \\subset B$, and since the model is graphical we must have $\\psi_B = 0$.  Therefore a graphical model is hierarchical.\n",
    "\n",
    "We can prove that hierarchical models are not all graphical with a counterexample, namely the book's Example 18.11. This is the hierarchical model with generator $1.2 + 1.3 + 2.3$:\n",
    "\n",
    "$$\\log f(x) = \\psi_\\varnothing(x) + \\psi_1(x) + \\psi_2(x) + \\psi_3(x) + \\psi_{12}(x) + \\psi_{13}(x) + \\psi_{23}(x)$$\n",
    "\n",
    "This model is hierarchical, but it is not graphical since $\\psi_{123} = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 18.8.5**.  Consider random variables $(X_1, X_2, X_3, X_4)$.  Suppose the log-density is\n",
    "\n",
    "$$\\log f(x) = \\psi_\\varnothing(x) + \\psi_{12}(x) + \\psi_{13}(x) + \\psi_{24}(x) + \\psi_{34}(x) $$\n",
    "\n",
    "**(a)** Draw the graph $\\mathcal{G}$ for these variables.\n",
    "\n",
    "**(b)** Write down all the independence and conditional independence relations implied by the graph.\n",
    "\n",
    "**(c)** Is this model graphical?  Is it hierarchical?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**\n",
    "\n",
    "**(a)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"134pt\" height=\"188pt\"\r\n",
       " viewBox=\"0.00 0.00 134.00 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-184 130,-184 130,4 -4,4\"/>\r\n",
       "<!-- X1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>X1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X1</text>\r\n",
       "</g>\r\n",
       "<!-- X2 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>X2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X2</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X2 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>X1&#45;&#45;X2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.6504,-144.765C48.8345,-133.456 41.1103,-118.437 35.3043,-107.147\"/>\r\n",
       "</g>\r\n",
       "<!-- X3 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>X3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X3</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X3 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>X1&#45;&#45;X3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M71.3496,-144.765C77.1655,-133.456 84.8897,-118.437 90.6957,-107.147\"/>\r\n",
       "</g>\r\n",
       "<!-- X4 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>X4</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X4</text>\r\n",
       "</g>\r\n",
       "<!-- X2&#45;&#45;X4 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>X2&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M35.3496,-72.7646C41.1655,-61.456 48.8897,-46.4367 54.6957,-35.1473\"/>\r\n",
       "</g>\r\n",
       "<!-- X3&#45;&#45;X4 -->\r\n",
       "<g id=\"edge4\" class=\"edge\"><title>X3&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M90.6504,-72.7646C84.8345,-61.456 77.1103,-46.4367 71.3043,-35.1473\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Graph at 0x10d4faf7d30>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Graph\n",
    "\n",
    "g = Graph()\n",
    "\n",
    "g.edge('X1', 'X2')\n",
    "g.edge('X1', 'X3')\n",
    "g.edge('X2', 'X4')\n",
    "g.edge('X3', 'X4')\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**\n",
    "\n",
    "- $X_1 \\text{ ⫫ } X_4 |  X_2, X_3 $ \n",
    "- $X_2 \\text{ ⫫ } X_3 |  X_1, X_4 $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)** The model is not graphical, since $\\psi_1 = 0$ but the definition of a graphical model requires all $\\psi_A$ to be non-zero unless an edge from $\\mathcal{G}$ is contained in $A$ -- and $A = \\{ 1 \\}$ contains no edges.\n",
    "\n",
    "Since the model is not graphical, it is not hierarchical either."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 18.8.6**. Suppose that parameters $p(x_1, x_2, x_3)$ are proportional to the following values:\n",
    "\n",
    "$$\n",
    "\\begin{array}{cccccc}\n",
    "\\hline\n",
    "    & x_2 & 0  & 0   & 1  & 1 \\\\\n",
    "    & x_3 & 0  & 1   & 0  & 1 \\\\\n",
    "\\hline\n",
    "x_1 & 0   & 2  &   8 &  4 & 16 \\\\\n",
    "    & 1   & 16 & 128 & 32 & 256 \\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Find the $\\psi$-terms for the log-linear expansion.  Comment on the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  Rewriting the given parameters table, where $p'$ are the given values,\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccc | c | c}\n",
    "x_1 & x_2 & x_3 & p'   & \\log_2 p'\\\\\n",
    "\\hline\n",
    "0   &   0  &  0 & 2   & 1\\\\\n",
    "0   &   0  &  1 & 8   & 3\\\\\n",
    "0   &   1  &  0 & 4   & 2\\\\\n",
    "0   &   1  &  1 & 16  & 4\\\\\n",
    "1   &   0  &  0 & 16  & 4\\\\\n",
    "1   &   0  &  1 & 128 & 7\\\\\n",
    "1   &   1  &  0 & 32  & 5\\\\\n",
    "1   &   1  &  1 & 256 & 8\\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "and by inspection we can verify that\n",
    "\n",
    "$$ \\log_2 p' = 1 + 3 x_1 + x_2 + 2 x_3 + x_1 x_3 $$\n",
    "\n",
    "Now, since $f(x_1, x_2, x_3) \\propto p'(x_1, x_2, x_3)$, then the log density is the value above plus a constant:\n",
    "\n",
    "$$ \\log f(x) = c + \\frac{3}{\\log 2} x_1 + \\frac{1}{\\log 2} x_2 + \\frac{2}{\\log 2} x_3 + \\frac{1}{\\log 2} x_1 x_3  = \\psi_\\varnothing + \\psi_1(x) + \\psi_2(x) + \\psi_3(x) + \\psi_{13}(x)$$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \\psi_\\varnothing(x) = c\n",
    "\\quad\n",
    "\\psi_1(x) = \\frac{3}{\\log 2} x_1\n",
    "\\quad\n",
    "\\psi_2(x) = \\frac{1}{\\log 2} x_2\n",
    "\\quad\n",
    "\\psi_3(x) = \\frac{2}{\\log 2} x_3\n",
    "\\quad\n",
    "\\psi_{13}(x) = \\frac{1}{\\log 2} x_1 x_3\n",
    "$$\n",
    "\n",
    "All that remains is finding $c = \\psi_\\varnothing(x)$ so that the probability density adds up to 1.  But $f(0, 0, 0) = \\frac{p'(0, 0, 0)}{\\sum_x p'(x)} = \\frac{1}{231}$, so\n",
    "\n",
    "$$ \\psi_\\varnothing(x) = - \\log 231$$\n",
    "\n",
    "This is a graphical model with generator $1.3 + 2$; we have $X_1 \\text{ ⫫ } X_3 |  X_2$ and the following graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"62pt\" height=\"188pt\"\r\n",
       " viewBox=\"0.00 0.00 62.00 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-184 58,-184 58,4 -4,4\"/>\r\n",
       "<!-- X1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>X1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X1</text>\r\n",
       "</g>\r\n",
       "<!-- X2 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>X2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X2</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X2 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>X1&#45;&#45;X2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-143.697C27,-132.846 27,-118.917 27,-108.104\"/>\r\n",
       "</g>\r\n",
       "<!-- X3 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>X3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X3</text>\r\n",
       "</g>\r\n",
       "<!-- X2&#45;&#45;X3 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>X2&#45;&#45;X3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-71.6966C27,-60.8463 27,-46.9167 27,-36.1043\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Graph at 0x10d4ce96b38>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Graph\n",
    "\n",
    "g = Graph()\n",
    "\n",
    "g.edge('X1', 'X2')\n",
    "g.edge('X2', 'X3')\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 18.8.7**.  Let $X_1, \\dots, X_4$ be binary.  Draw the independence graphs corresponding to the following log-linear models.  Also, identify whether each is graphical and/or hierarchical (or neither).\n",
    "\n",
    "**(a)**  $\\log f = 7 + 11 x_1 + 2 x_2 + 1.5 x_3 + 17 x_4$\n",
    "\n",
    "**(b)**  $\\log f = 7 + 11 x_1 + 2 x_3 + 1.5 x_3 + 17 x_4 + 12 x_2 x_3 + 78 x_2 x_4 + 3 x_3 x_4 + 32 x_2 x_3 x_4$\n",
    "\n",
    "**(c)**  $\\log f = 7 + 11 x_1 + 2 x_3 + 1.5 x_3 + 17 x_4 + 12 x_2 x_3 + 3 x_3 x_4 + x_1 x_4 + 2 x_1 x_2$\n",
    "\n",
    "**(d)**  $\\log f = 7 + 5055 x_1 x_2 x_3 x_4$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  Let's assume these are valid density functions -- that is, that a constant is added to each of these so that the densities add up to 1, or that the log is applied in an applicable base.  In either scenario, we just need the zero or non-zero properties of each coefficient to determine the answers for this exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(a)** $\\log f = \\psi_\\varnothing + \\psi_1 + \\psi_2 + \\psi_3 + \\psi_4$\n",
    "\n",
    "This is a graphical model with mutual independence across all variables.  It is graphical and hierarchical, corresponding to model generator $1 + 2 + 3 + 4$.  The graph is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"278pt\" height=\"44pt\"\r\n",
       " viewBox=\"0.00 0.00 278.00 44.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 40)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-40 274,-40 274,4 -4,4\"/>\r\n",
       "<!-- X1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>X1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X1</text>\r\n",
       "</g>\r\n",
       "<!-- X2 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>X2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X2</text>\r\n",
       "</g>\r\n",
       "<!-- X3 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>X3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X3</text>\r\n",
       "</g>\r\n",
       "<!-- X4 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>X4</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"243\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X4</text>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Graph at 0x10d358bfeb8>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Graph\n",
    "\n",
    "g = Graph()\n",
    "\n",
    "g.node('X1')\n",
    "g.node('X2')\n",
    "g.node('X3')\n",
    "g.node('X4')\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** $\\log f = \\psi_\\varnothing + \\psi_1 + \\psi_2 + \\psi_3 + \\psi_4 + \\psi_{23} + \\psi_{24} + \\psi_{34} + \\psi_{234}$\n",
    "\n",
    "This is a graphical model with independence between 1 and all other variables.  It is graphical and hierarchical, corresponding the model generator $1 + 2.3.4$.  The graph is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"134pt\" height=\"188pt\"\r\n",
       " viewBox=\"0.00 0.00 134.00 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-184 130,-184 130,4 -4,4\"/>\r\n",
       "<!-- X1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>X1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X1</text>\r\n",
       "</g>\r\n",
       "<!-- X2 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>X2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X2</text>\r\n",
       "</g>\r\n",
       "<!-- X3 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>X3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"71\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"71\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X3</text>\r\n",
       "</g>\r\n",
       "<!-- X2&#45;&#45;X3 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>X2&#45;&#45;X3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M92.3644,-144.411C87.901,-133.252 82.0414,-118.604 77.588,-107.47\"/>\r\n",
       "</g>\r\n",
       "<!-- X4 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>X4</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"98\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"98\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X4</text>\r\n",
       "</g>\r\n",
       "<!-- X3&#45;&#45;X4 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>X3&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M77.3986,-72.411C81.6403,-61.4141 87.1894,-47.0274 91.4597,-35.9562\"/>\r\n",
       "</g>\r\n",
       "<!-- X4&#45;&#45;X2 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>X4&#45;&#45;X2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M101.654,-36.0925C103.676,-46.4315 105.981,-59.9098 107,-72 108.344,-87.9434 108.197,-92.0449 107,-108 106.095,-120.063 104.046,-133.542 102.248,-143.888\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Graph at 0x10d5040f438>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Graph\n",
    "\n",
    "g = Graph()\n",
    "\n",
    "g.node('X1')\n",
    "g.edge('X2', 'X3')\n",
    "g.edge('X3', 'X4')\n",
    "g.edge('X4', 'X2')\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)** $\\log f = \\psi_\\varnothing + \\psi_1 + \\psi_2 + \\psi_3 + \\psi_4 + \\psi_{23} + \\psi_{34} + \\psi_{14} + \\psi_{12}$\n",
    "\n",
    "This is a graphical and hierarchical model.  It corresponds to model generator $1.2 + 2.3 + 3.4 + 4.1$; the only conditional independences require the other two given variables.  The graph is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"89pt\" height=\"260pt\"\r\n",
       " viewBox=\"0.00 0.00 89.00 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-256 85,-256 85,4 -4,4\"/>\r\n",
       "<!-- X1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>X1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"54\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X1</text>\r\n",
       "</g>\r\n",
       "<!-- X2 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>X2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X2</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X2 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>X1&#45;&#45;X2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M47.6014,-216.411C43.3597,-205.414 37.8106,-191.027 33.5403,-179.956\"/>\r\n",
       "</g>\r\n",
       "<!-- X4 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>X4</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"54\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X4</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X4 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>X1&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M57.6538,-215.908C59.6758,-205.569 61.9808,-192.09 63,-180 67.0322,-132.17 67.0322,-119.83 63,-72 61.9808,-59.9098 59.6758,-46.4315 57.6538,-36.0925\"/>\r\n",
       "</g>\r\n",
       "<!-- X3 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>X3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X3</text>\r\n",
       "</g>\r\n",
       "<!-- X2&#45;&#45;X3 -->\r\n",
       "<g id=\"edge4\" class=\"edge\"><title>X2&#45;&#45;X3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-143.697C27,-132.846 27,-118.917 27,-108.104\"/>\r\n",
       "</g>\r\n",
       "<!-- X3&#45;&#45;X4 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>X3&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M33.3986,-72.411C37.6403,-61.4141 43.1894,-47.0274 47.4597,-35.9562\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Graph at 0x10d4ff2e240>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Graph\n",
    "\n",
    "g = Graph()\n",
    "\n",
    "g.edge('X1', 'X2')\n",
    "g.edge('X1', 'X4')\n",
    "g.edge('X3', 'X4')\n",
    "g.edge('X2', 'X3')\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(d)** $\\log f = \\psi_\\varnothing + \\psi_{1234}$\n",
    "\n",
    "This model has dependencies between all variables, but it is not graphical nor hierarchical; $\\psi_{12} = 0$ but $\\psi_{1234} \\neq 0$.  We can nonetheless draw a graph that would include this model without including all of its independence information by drawing a fully saturated graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"122pt\" height=\"260pt\"\r\n",
       " viewBox=\"0.00 0.00 122.47 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-256 118.473,-256 118.473,4 -4,4\"/>\r\n",
       "<!-- X1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>X1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"73\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"73\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X1</text>\r\n",
       "</g>\r\n",
       "<!-- X2 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>X2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"73\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"73\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X2</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X2 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>X1&#45;&#45;X2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M73,-215.697C73,-204.846 73,-190.917 73,-180.104\"/>\r\n",
       "</g>\r\n",
       "<!-- X3 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>X3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X3</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X3 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>X1&#45;&#45;X3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M59.7877,-217.887C51.7919,-207.846 42.1708,-193.98 37,-180 28.2078,-156.23 26.5578,-126.648 26.5193,-108.292\"/>\r\n",
       "</g>\r\n",
       "<!-- X4 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>X4</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"82\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"82\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">X4</text>\r\n",
       "</g>\r\n",
       "<!-- X1&#45;&#45;X4 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>X1&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M86.7504,-218.069C94.9608,-208.1 104.619,-194.247 109,-180 124.828,-128.529 102.288,-64.7748 89.6146,-35.3225\"/>\r\n",
       "</g>\r\n",
       "<!-- X2&#45;&#45;X3 -->\r\n",
       "<g id=\"edge4\" class=\"edge\"><title>X2&#45;&#45;X3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M62.5619,-145.116C55.0337,-133.66 44.9155,-118.263 37.3968,-106.821\"/>\r\n",
       "</g>\r\n",
       "<!-- X2&#45;&#45;X4 -->\r\n",
       "<g id=\"edge5\" class=\"edge\"><title>X2&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M74.0856,-143.871C75.8155,-116.578 79.1783,-63.5205 80.9105,-36.1893\"/>\r\n",
       "</g>\r\n",
       "<!-- X3&#45;&#45;X4 -->\r\n",
       "<g id=\"edge6\" class=\"edge\"><title>X3&#45;&#45;X4</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M38.9342,-73.811C48.1595,-62.0697 60.8696,-45.8932 70.0891,-34.1593\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Graph at 0x10d348b5748>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Graph\n",
    "\n",
    "g = Graph()\n",
    "\n",
    "g.edge('X1', 'X2')\n",
    "g.edge('X1', 'X3')\n",
    "g.edge('X1', 'X4')\n",
    "g.edge('X2', 'X3')\n",
    "g.edge('X2', 'X4')\n",
    "g.edge('X3', 'X4')\n",
    "\n",
    "g"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}