Computational Proba and Inference - Part 2 - Inference in Graphical Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Inference in Graphical Models\n",
    "\n",
    "As we saw earlier, inference about hidden random variables from observations requires an exponential amount of computation. In this section, we will see that exploiting structure and probability distributions can lead to much more efficient inference.\n",
    "\n",
    "Interestingly, the class of probability distributions that have that structure is quite broad and applicable to a wide range of real world problems. We will formally describe and characterize these probability distributions. Surprisingly, the complexity of the inference algorithms for this class of probability distributions becomes linear in the number of variables involved.\n",
    "\n",
    "In this section, the core concept that we will introduce and use is *graphical models*. Graphical models combine graph theory with probability to capture structure in probability distributions. We will use graphical models to describe the structure, to reason about conditional independence relationships in these probability distributions, and to do derive efficient inference algorithms.\n",
    "\n",
    "The algorithms you will see in the section are closely related to two famous algorithms. The first one is the *Kalman filter*. It was used to guide the Apollo mission, and is used today in many important engineering applications.\n",
    "\n",
    "The second algorithm is the *Viterbi algorithm*. It was originally designed to decode messages in noisy communication networks. And today it is used in almost every wireless communication protocol."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I— Big O notation\n",
    "\n",
    "## I.1— Definition\n",
    "\n",
    "We say that a function $f$ that depends on a variable $n$ is $\\mathcal{O}(g(n))$ (read as \"big O of g of n\") if there exists some minimum $n_0$ such that for all $n \\ge n_0$,\n",
    "$$f(n) \\le c \\cdot g(n)$$\n",
    "\n",
    "for some constant $c > 0$. Notationally, we write either $f(n) = \\mathcal{O}(g(n))$ or $f(n) \\in \\mathcal{O}(g(n))$. Note that the $\\mathcal{O}$ stands for “order\" so you could think of $f$ being $\\mathcal{O}(g(n))$ as $f$ growing on the order of $g$ as a function of $n$.\n",
    "\n",
    "**Example**: Storing the probability table for a random variable with alphabet size $k$ takes $\\mathcal{O}(k)$ amount of space. Even if we wanted to be efficient and store $k-1$ numbers, $k-1$ is still $\\mathcal{O}(k)$, and in fact $k-1 \\le k$ for all $k$, satisfying the definition of big O notation (with $n_0$ being set to whatever we want, and $c$ set to 1).\n",
    "\n",
    "**Example**: With joint probability table $p_{X,Y}$ where $X$ and $Y$ each take on $k$ values, then marginalization to compute $p_ X$ costs $\\mathcal{O}(k^2)$ basic operations. For example, using the way we counted the number of basic operations earlier, $1 + 3k + 3k^2 \\le 6k^2$ for all $k \\ge 2$, satisfying the definition of big O notation with $n_0 = 2$ and $c=6$.\n",
    "\n",
    "The basic idea is that when $n$ is large enough (specifically larger than some $n_0$), then we'll always have that $f(n)$ grows at most as fast as $g(n)$ scaled by a constant that does not depend on $n$.\n",
    "\n",
    "## I.2— Big O Notation with Multiple Variables\n",
    "\n",
    "Big O notation can be used with functions that depend on multiple variables.\n",
    "\n",
    "*Example:* In our material coverage for Bayes' theorem for random variables, we saw in an exercise where we have $n$ random variables that we want a posterior distribution over, where each of these random variables has alphabet size $k$. Then computing the denominator of Bayes' theorem involves summing over $k^n$ table entries, a computation that takes running time $\\mathcal{O}(k^ n)$.\n",
    "\n",
    "## I.3— Important Remarks Regarding Big O Notation\n",
    "\n",
    "Big O notation is defined as an upper bound on how fast a function grows. For example, a function that grows linearly in $n$\n",
    "is $\\mathcal{O}(n^2)$ **and also** $\\mathcal{O}(2^ n)$, both of which grow much faster than linear.\n",
    "\n",
    "Typically, when using big O notation such as $f(n) = \\mathcal{O}(g(n))$, we will pick $g$ that grows at a rate that *matches or closely matches* the actual growth rate of $f$, but not always !\n",
    "\n",
    "An example of a mathematical operation for which people often use a function $g$ that doesn't actually match $f$ in order of growth is multiplying two n-by-n matrices, for which the fastest known algorithm takes time $\\mathcal{O}(n^{2.373})$ but since the exponent 2.373 is peculiar and the algorithm that achieves it is quite complicated and not what's typically used in practice, often times for convenience people just say that multiplying two n-by-n matrices takes time $\\mathcal{O}(n^3)$. Note that when they say such a statement, they are not explicitly saying which matrix multiplication algorithm is used either.\n",
    "\n",
    "Finally, some terminology:\n",
    "\n",
    "- If $f(n) = \\mathcal{O}(1)$, then we say that $f$ is *constant* with respect to input $n$.\n",
    "\n",
    "- If $f(n) = \\mathcal{O}(\\log n)$, then we say that $f$ *grows logarithmically* respect to input $n$. For example, if $f$ is the running time of a computer program, then we would say that $f$ has logarithmic running time in $n$.\n",
    "\n",
    "- If $f(n) = \\mathcal{O}(n)$, then we say that $f$ *grows linearly* in $n$. For example, if $f$ is the running time of a computer program, then we would say that $f$ has linear running time in $n$.\n",
    "\n",
    "- If $f(n) = \\mathcal{O}(n^2)$, then we say that $f$ *grows quadratically* in $n$. For example, if $f$ is the running time of a computer program, then we would say that $f$ has quadratic running time in $n$.\n",
    "\n",
    "- If $f(n) = \\mathcal{O}(b^ n)$ for some constant $b > 1$, then we say that $f$ *grows exponentially* in n. For example, if $f$ is the running time of a computer program, then we would say that $f$ has exponential running time in $n$. \n",
    "\n",
    "The above are of course just a few examples. There are many other “categories\" of functions such as $\\mathcal{O}(n^3)$ corresponding to cubic growth in $n$.\n",
    "\n",
    "In much of our discussion to follow in 6.008.1x, we will analyze either the “time complexity\" (i.e., how long it takes code to run in terms of number of basic operations, in big O) or “space complexity\" (i.e., how much space a code uses for storing variables, also in big O). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# II— Graphical Models\n",
    "\n",
    "To introduce graphical models, we start with three simple examples.\n",
    "\n",
    "Note that previously we have been using random variables $X$\n",
    "and $Y$, but here we will use $X_1$, $X_2$, up to $X_ n$, as it will provide a clean way to write out the general case.\n",
    "\n",
    "## II.1— Graphical Model of Two Independent Random Variables\n",
    "\n",
    "If $X_1$ and $X_2$ are independent random variables, then we know that $p_{X_1,X_2}(x_1,x_2)=p_{X_1}(x_1)p_{X_2}(x_2)$. Graphically, we represent this distribution as two circles. Because they are independent, we do not “connect\" these two circles:\n",
    "![indep var](./images/images_sec-graphical-models-2-rv-indep.png)\n",
    "On a computer, we store a separate table for each of the two circles, one for $p_{X_1}$ and one for $p_{X_2}$. Later we will see that tables that we store that depend only on a single random variable need not be a marginal distribution (such as in this example). Thus, we will give the tables new names. We let $\\phi _1 = p_{X_1}$ and $\\phi _2 = p_{X_2}$. (Later, $\\phi _ i$ is the table we store that is associated with a single random variable $X_ i$. In this example, $\\phi _1$ and $\\phi _2$ are just set to be the same as the marginal distributions $p_{X_1}$ and $p_{X_2}$.)\n",
    "\n",
    "To summarize, the graphical model here includes the picture above (which is called a graph) along with the tables $\\phi _1$ and $\\phi _2$.\n",
    "\n",
    "## II.2— Two Possibly Dependent Random Variables\n",
    "\n",
    "Now suppose that we do not know whether $X_1$ and $X_2$ are independent. Then without further assumptions, we can work directly with the joint probability table $p_{X_1, X_2}$, or we can use the product rule (which importantly always holds and does not require $X_1$ and $X_2$ to be independent).\n",
    "\n",
    "Here are three different ways we can store the distribution:\n",
    "\n",
    "- (a) Store $p_{X_1, X_2}$. If we do this, then there is a single table $p_{X_1, X_2}$ that we are storing that depends on two random variables $X_1$ and $X_2$.\n",
    "\n",
    " We introduce new notation here: a table we store that depends on *exactly two random variables* $X_ i$\n",
    "and $X_ j$ will be called $\\psi _{i, j}$, so in this case, we store $\\psi _{i, j} = p_{X_1, X_2}$. Later, we will see examples where $\\psi _{i, j}$ is not a joint probability table for two random variables.\n",
    "\n",
    " There are no tables here that depend on exactly 1 random variable.\n",
    "\n",
    "- (b) Store $p_{X_1}$ and $p_{X_2 \\mid X_1}$.\n",
    "Here, there is one table that depends on *exactly 1 random variable*: $p_{X_1}$. We call this $\\phi _{X_1} = p_{X_1}$.\n",
    "\n",
    " There is one table that depends on *exactly 2 random variables*: $\\psi _{1, 2} = p_{X_2 \\mid X_1}$. (This $\\psi _{1, 2}$ is different from the one in (a).)\n",
    "\n",
    "- (c) Store $p_{X_2}$ and $p_{X_1 \\mid X_2}$.\n",
    "\n",
    " Here, there is one table that depends on exactly 1 random variable: $p_{X_2}$. We call this $\\phi _{X_2} = p_{X_2}$.\n",
    "\n",
    " There is one table that depends on exactly 2 random variables: $\\psi _{1, 2} = p_{X_1 \\mid X_2}$. (This $\\psi _{1, 2}$ is different from the ones in (a) and (b).) \n",
    "\n",
    "The tables being stored in each of the above three cases are different, but the joint probability distribution they describe is the same.\n",
    "\n",
    "A common feature of all three different ways: **there is a table that depends on both $X_1$ and $X_2$**. For this reason, when we draw out a graphical representation in this case, we still have two circles, one for $X_1$ and one for $X_2$, but now we connect the two with a line:\n",
    "![connected](./images/images_sec-graphical-models-2-rv-possibly-dependent.png)\n",
    "Again, the line is there between $X_1$ and $X_2$ precisely because to store the associated joint probability table, regardless of which of the different ways we store the tables, we have to use a table that depends on both $X_1$ and $X_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II.3— Markov Chain\n",
    "\n",
    "Next, suppose we have three random variables $X_1$, $X_2$, and $X_3$, where we make an assumption here that the joint probability table has factorization\n",
    "$$p_{X_1,X_2,X_3}(x_1,x_2,x_3) = p_{X_1}(x_1) p_{X_2|X_1}(x_2|x_1) p_{X_3|X_2}(x_3|x_2).$$\n",
    "\t(3.1)\n",
    "\n",
    "Note that this factorization is **in general not true** for three random variables $X_1$, $X_2$, and $X_3$.\n",
    "\n",
    "To see, this, we can apply the product rule (which is true for any three random variables $X_1$, $X_2$, and $X_3$):\n",
    "$$p_{X_1,X_2,X_3}(x_1,x_2,x_3) = p_{X_1}(x_1) p_{X_2|X_1}(x_2|x_1) p_{X_3|X_1, X_2}(x_3|x_1, x_2).$$\n",
    "\t(3.2)\n",
    "\n",
    "Thus, in the case of a Markov Chain:\n",
    "$$p_{X_3|X_1,X_2}(x_3|x_1,x_2)=p_{X_3|X_2}(x_3|x_2).$$\n",
    "This means that given $X_2$, knowing $X_1$ does not tell us anything new about $X_3$: $X_1 \\perp X_3 \\mid X_2$.\n",
    "\n",
    "So assuming equation (3.1) holds means that we are making an assumption about conditional independence structure!\n",
    "\n",
    "Next, let's see how to store the distribution when it has factorization given by equation (3.1). We see that it suffices to keep track of three tables $\\phi _1$, $\\psi _{1, 2}$, and $\\psi _{2, 3}$:\n",
    "$$p_{X_1,X_2,X_3}(x_1,x_2,x_3) = \\underbrace{p_{X_1}(x_1)}_{\\phi _1(x_1)} \\underbrace{p_{X_2|X_1}(x_2|x_1)}_{\\psi _{1, 2}(x_1, x_2)} \\underbrace{p_{X_3| X_2}(x_3|x_2)}_{\\psi _{2, 3}(x_2, x_3)}.$$\n",
    "\n",
    "The graph associated with this representation has three circles, one for each of the random variables $X_1$, $X_2$, and $X_3$. We have a table $\\psi _{1, 2}$ that depends on $X_1$ and $X_2$ so we draw a line between the circles for $X_1$ and $X_2$. Next we have a table $\\psi _{2, 3}$ that depends on $X_2$ and $X_3$ so we draw a line between the circles for $X_2$ and $X_3$. This yields the following:\n",
    "![markov chain](./images/images_sec-graphical-models-3-rv-markov-chain.png)\n",
    "This line-shaped graph is called a **Markov chain**. We will encounter Markov chains more later on. Notationally, when $X_1$, $X_2$, and $X_3$ form a Markov chain, we write $X_1 \\leftrightarrow X_2 \\leftrightarrow X_3$.\n",
    "\n",
    "## II.4— The General Case\n",
    "\n",
    "We are almost ready to mathematically define what a graphical model is. As we saw from the above examples, each time we had a graph (a picture with circles and possibly lines) along with tables that we store.\n",
    "\n",
    "For a graph, the circles are formally called nodes or vertices, and the lines are formally called edges. The graph is undirected because the edges do not have directionality associated with them. Furthermore:\n",
    "\n",
    "- Each node corresponds to a random variable\n",
    "\n",
    "- Each edge indicates there being some **possible** dependence between the two nodes that it connects. Specifically: an edge being present between two nodes $X_ i$ and $X_ j$ means that the equation for the probability distribution has a factor that depends on both $X_ i$ and $X_ j$; this factor is stored as table $\\psi _{i, j}$. \n",
    "\n",
    "From these definitions, an important concept emerges:\n",
    "\n",
    "More edges implies that the model can encode more probability distributions but we have to store more tables. (Think about why the second example we presented for two random variables where we don't know whether they're independent or not (and had a factor $\\psi _{1, 2}$) can actually encode a distribution in which $X_1$ and $X_2$ are independent!)\n",
    "\n",
    "Meanwhile, in the graphs that we will consider, we will assume that there are **no loops**. We do this for simplicity: probabilistic graphical models with loops are beyond the scope of this course.\n",
    "\n",
    "Note that a graph is specified by saying what the nodes (circles) are, and what the lines (edges) are. To give specific examples:\n",
    "\n",
    "- When $X_1$ and $X_2$ were independent, the graph we had consisted of the set of nodes $V=\\{ 1, 2\\}$ and the set of edges $E=\\emptyset$.\n",
    "\n",
    "- When $X_1$ and $X_2$ were not known as to whether they are independent or not, the graph we had consisted of the set of nodes $V=\\{ 1, 2\\}$ and the set of edges $E=\\{ (1, 2)\\}$. Note that in general we use $(i,j)$ to mean an edge between nodes $i$ and $j$, where $(i,j)$ and $(j,i)$ mean the same thing, which is why the set of the edges here includes only $(1,2)$ and not both $(1,2)$ and $(2,1)$.\n",
    "\n",
    "- When we have $X_1 \\leftrightarrow X_2 \\leftrightarrow X_3$, the set of nodes is $V=\\{ 1, 2, 3\\}$ and the set of edges is $E=\\{ (1, 2), (2, 3)\\}$. \n",
    "\n",
    "**Definition of an undirected pairwise graphical model** (we will just call this a graphical model):\n",
    "\n",
    "An undirected pairwise graphical model for random variables $X_1, \\dots , X_ n$ consists of an undirected graph with vertices $V=\\{ 1,\\dots ,n\\}$ and edges $E$, and tables $\\phi _ i$'s and $\\psi _{i,j}$'s that have nonnegative entries. The joint probability table of $X_1, \\dots , X_ n$ is given by\n",
    "$$p_{X_1, \\dots , X_ n}(x_1, \\dots , x_ n) =\\frac{1}{Z}\\prod _{i\\in V}\\phi _{i}(x_{i})\\prod _{(i,j)\\in E}\\psi _{ij}(x_{i},x_{j}),$$\n",
    "\t \n",
    "where $Z$ is the normalization constant that ensures that the probability distribution actually sums to 1 (we'll give a concrete example shortly).\n",
    "\n",
    "It can be easily shown that the normalization constant $Z$ can be computed as:\n",
    "$$Z=\\sum _{x_{1}}\\cdots \\sum _{x_{n}}\\prod _{i\\in V}\\phi _{i}(x_{i})\\prod _{(i,j)\\in E}\\psi _{ij}(x_{i},x_{j}).$$\n",
    "\n",
    "Note that in earlier examples, we didn't always specify that a random variable $X_ i$ needed to have a table $\\phi _ i$. This is not a problem: we can just set $\\phi _ i(x_ i) = 1$ for all values of $x_ i$ in this case.\n",
    "\n",
    "**Terminology:**\n",
    "\n",
    "-  Each table $\\phi _ i$ depends only on random variable $X_ i$ and is called the **node potential function** or **node potential** of node $i$.\n",
    "\n",
    "- Each table $\\psi _{i, j}$ depends only on random variables $X_ i$ and $X_ j$ and is called the **pairwise potential function** or **pairwise potential** or **edge potential** of nodes $i$ and $j$. \n",
    "\n",
    "**Important**: We require that the potential tables consist of *nonnegative entries* but each potential table *does not have to sum to 1*. The constant $Z$ will ensure that the joint probability table actually sums to 1. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I.5— Trees\n",
    "\n",
    "A **tree** is a graph for which\n",
    "- there are no loops, and\n",
    "- we can reach from any node to any other node (moving along edges in the graph).\n",
    "\n",
    "We'll be seeing trees quite a bit so here are some basics of trees.\n",
    "\n",
    "Theorem: For any graph that has $n$ nodes, if the graph is a tree, then it will always have exactly $n-1$ edges. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I.6— Incorporating Observations in Graphical Models\n",
    "\n",
    "Let's start from the 3-node Markov chain we had earlier $X_1 \\leftrightarrow X_2 \\leftrightarrow X_3$.\n",
    "\n",
    "We had the factorization:\n",
    "$$p_{X_1,X_2,X_3}(x_1,x_2,x_3) = \\frac{1}{Z} \\phi _1(x_1) \\phi _2(x_2) \\phi _3(x_3) \\psi _{1,2}(x_1, x_2) \\psi _{2,3}(x_2, x_3).$$\n",
    "\n",
    "Suppose we condition on $X_2 = v$ for some fixed value $v$ in the alphabet of $X_2.$ We want to figure out the distribution $p_{X_1,X_3\\mid X_2}(\\cdot ,\\cdot \\mid v)$.\n",
    "By the definition of conditional probability,\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_1,X_3\\mid X_2}(x_1,x_3\\mid v)\n",
    "& =& \\frac{p_{X_1,X_2,X_3}(x_1,v,x_3)}{p_{X_2}(v)}\\\\\n",
    "& =& \\frac{\\frac{1}{Z} \\phi _1(x_1) \\phi _2(v) \\phi _3(x_3) \\psi _{1,2}(x_1, v) \\psi _{2,3}(v, x_3)}{p_{X_2}(v)}\\\\\n",
    "& =& \\frac{1}{Z \\frac{p_{X_2}(v)}{\\phi _2(v)}} \\phi _1(x_1)\\psi _{1,2}(x_1, v) \\phi _3(x_3)\\psi _{2,3}(v, x_3).\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Let's define the following:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    " Z' & \\triangleq & Z\\frac{p_{X_2}(v)}{\\phi _2(v)},\\\\\n",
    " \\phi _1'(x_1) & \\triangleq & \\phi _1(x_1)\\psi _{1,2}(x_1, v),\\\\\n",
    " \\phi _3'(x_3) & \\triangleq & \\phi _3(x_3)\\psi _{2,3}(v, x_3).\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "We can then write:\n",
    "$$p_{X_1,X_3\\mid X_2}(x_1,x_3\\mid v) = \\frac{1}{Z'}\\phi _1'(x_1)\\phi _3'(x_3)$$ \n",
    "\n",
    "This corresponds to a new graphical model with only 2 nodes and no edges (minimum) ! This new graphical model is the same as if we remove the node corresponding to $X_2$ in the original graph.\n",
    "\n",
    "We can always view conditioning (and thus incorporating observations) as fixing the value(s) of whichever random variable(s) we observe.\n",
    "This always has the effect that we just saw: the pairwise potentials involving the observed random variables become part of new node potentials involving the **unobserved** (also called *hidden* or *latent*) random variables. Since these pairwise potentials corresponded to edges that were present, and now they have become part of node potentials instead, the *effect on the graph is that we deleted the nodes that we made observations for*. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# II— Fundamental Inference Tasks in Graphical Models\n",
    "\n",
    "Now that we have introduced graphical models for representing probability distributions, we proceed to solving problems of inference. As we have seen, observations can be taken care of quite easily and once we account for them, we are left with a new graphical model that is only over the random variables that we don't get to observe, which are called “hidden\" or “latent\" random variables. Thus, our discussion to follow will be about graphical models where we don't explicitly mention conditioning, with the idea that the conditioning has already happened!\n",
    "\n",
    "Given a graphical model with graph $G=(V, E)$ and its associated node and edge potentials, the two fundamental inference tasks we focus on are as follows:\n",
    "\n",
    "- **Marginalization**: Compute marginal probability table $p_{X_ i}$ for every $i \\in V$.\n",
    "\n",
    "- **Most probable configuration**: Compute the most probable configuration ($\\widehat{x}_1, \\widehat{x}_2, \\dots , \\widehat{x}_ n$) such that \n",
    "$$(\\widehat{x}_1, \\widehat{x}_2, \\dots , \\widehat{x}_ n) = \\arg \\max _{x_1, x_2, \\dots , x_ n} p_{X_1, X_2, \\dots , X_ n}(x_1, x_2, \\dots , x_ n).$$\n",
    "\n",
    "Again, we are assuming conditioning has already happened.\n",
    "\n",
    "This means that if we observed random variable $Y=y$ and already accounted for it, then the graphical model that remains is only over random variables $X_1, \\dots , X_ n$ and so in the first task above, the table $p_{X_ i}$ actually corresponds to $p_{X_ i \\mid Y}(\\cdot \\mid y)$ and the second task would compute the most probable configuration in the posterior distribution $p_{X_1, \\dots , X_ n \\mid Y}(\\cdots \\mid y)$. (Note that there was nothing special about us calling the random variables in a graphical model $X_1, \\dots , X_ n$. Doing so was just to make the notation simple, but the random variables could be called whatever you want, including having $X_ i$'s and $Y_ i$'s, and also observing multiple random variables rather than just observing one random variable is handled the same way: fix the values for what has been observed, which corresponds to removing those nodes from the original graphical model and changing potential tables to account for the observations.)\n",
    "\n",
    "We now focus on the inference task of marginalization in graphical models.\n",
    "\n",
    "## II.1— The Sum-Product Algorithm\n",
    "\n",
    "The Sum-product algorithm provides a way to find all the marginal probability distributions in any tree-structured undirected graphical model.\n",
    "\n",
    "What about graphical models that aren't tree-structured and that don't have loops? Think about why in such cases, the different disconnected pieces can be treated separately (and each of these pieces of the graph will be a tree)! For example, when computing the marginal distribution $p_{X_i}$ for a specific node $i$, any node that isn't reachable from $i$ doesn't matter in the computation. So we can always turn the problem into just looking at one tree at a time.\n",
    "\n",
    "### II.1.1— Marginalization: A Worked Example\n",
    "\n",
    "Let's work through an example of how to marginalize in an undirected graphical model. This example will illustrate a general algorithm for marginalization, which we can write a computer program for.\n",
    "\n",
    "Suppose random variables $X_{1},X_{2},X_{3},X_{4},X_{5}$ each take on values in the set $\\mathcal{X}$ and have the following joint distribution:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_1, \\dots, X_n}(x_1, \\dots, x_n)\n",
    "&=& \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3})\\phi_{4}(x_{4})\\phi_{5}(x_{5}) \\\\\n",
    "&& \\;\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})\\psi_{24}(x_{2},x_{4})\\psi_{25}(x_{2},x_{5}).\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "This joint distribution corresponds to the following undirected graphical model:\n",
    "![five-node-example](./images/images_sec-graphical-models-five-node-example.png)\n",
    "\n",
    "Let's compute the marginal distribution $p_{X_{1}}$. This distribution is a table assigning a probability to each value of $\\mathcal{X}$. We obtain the marginal distribution for $X_{1}$ by summing out the other random variables:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{1}}(x_{1}) & = &\\sum_{x_{2}}\\sum_{x_{3}}\\sum_{x_{4}}\\sum_{x_{5}}p_{X_1,\\dots,X_n}(x_1,\\dots,x_n)\\\\\n",
    "& = & \\sum_{x_{2}}\\sum_{x_{3}}\\sum_{x_{4}}\\sum_{x_{5}}\n",
    "\\bigg\\{\\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3})\\phi_{4}(x_{4})\\phi_{5}(x_{5}) \\\\\n",
    "& & \\qquad\\qquad\\qquad\\;\\;\\;\\;\n",
    "\\cdot\n",
    "\\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})\\psi_{24}(x_{2},x_{4})\\psi_{25}(x_{2},x_{5})\\bigg\\}.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "To compute the summation on the right-hand side efficiently on a computer, notice that we can push the summations around a bit taking advantage of the distributive property of arithmetic: $ab+ac=a(b+c)$ for any choice of constants $a,b,c$. For example, only two factors depend on $x_{5}$ (namely $\\phi _{5}$ and $\\psi _{25}$) so we can push the summation over $x_{5}$ inward as follows:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{1}}(x_{1}) & = & \\sum_{x_{2}}\\sum_{x_{3}}\\sum_{x_{4}}\\sum_{x_{5}}\n",
    "\\bigg\\{\\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3})\\phi_{4}(x_{4})\\phi_{5}(x_{5}) \\\\\n",
    "& & \\qquad\\qquad\\qquad\\;\\;\\;\\;\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})\\psi_{24}(x_{2},x_{4})\\psi_{25}(x_{2},x_{5}) \\bigg\\} \\\\\n",
    "& = & \\sum_{x_{2}}\\sum_{x_{3}}\\sum_{x_{4}}\n",
    "\\bigg\\{\\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3})\\phi_{4}(x_{4}) \\\\\n",
    "& & \\qquad\\qquad\\quad\\;\\;\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})\\psi_{24}(x_{2},x_{4})\\underbrace{\\sum_{x_{5}}\\phi_{5}(x_{5})\\psi_{25}(x_{2},x_{5})}_{\\triangleq m_{5 \\rightarrow 2}(x_{2})} \\bigg\\}.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Here we have introduced some notation: $m_{5\\rightarrow 2}(x_{2})$ is a value we get from summing out (and thus eliminating) $x_{5}$ where the result depends on $x_{2}$. Note that $m_{5\\rightarrow 2}$ is a table: there is one entry in the table per value of $x_{2}\\in \\mathcal{X}$. We use the letter “m\" since we can interpret $m_{5\\rightarrow 2}$ as a message that node 5 sends to node 2.\n",
    "\n",
    "We can keep pushing sums around:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{1}}(x_{1}) & = & \\sum_{x_{2}}\\sum_{x_{3}}\n",
    "\\bigg\\{ \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3}) \\\\\n",
    "& & \\qquad\\quad\\;\\;\\;\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})m_{5\\rightarrow2}(x_{2})\\underbrace{\\sum_{x_{4}}\\phi_{4}(x_{4})\\psi_{24}(x_{2},x_{4})}_{\\triangleq m_{4\\rightarrow2}(x_{2})} \\bigg\\} \\\\\n",
    " & = & \\sum_{x_{2}}\\sum_{x_{3}}\n",
    "\\bigg\\{ \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3}) \\\\\n",
    "& & \\qquad\\quad\\;\\;\\;\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})m_{4\\rightarrow2}(x_{2})m_{5\\rightarrow2}(x_{2}) \\bigg\\} \\\\\n",
    " & = & \\sum_{x_{2}}\\bigg\\{ \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\psi_{12}(x_{1},x_{2}) \\\\\n",
    "& & \\qquad\\;\n",
    "\\cdot m_{4\\rightarrow2}(x_{2})m_{5\\rightarrow2}(x_{2})\\underbrace{\\sum_{x_{3}}\\phi_{3}(x_{3})\\psi_{13}(x_{1},x_{3})}_{\\triangleq m_{3\\rightarrow1}(x_{1})} \\bigg\\} \\\\\n",
    " & = & \\frac{1}{Z}\\phi_{1}(x_{1})m_{3\\rightarrow1}(x_{1})\\underbrace{\\sum_{x_{2}}\\phi_{2}(x_{2})\\psi_{12}(x_{1},x_{2})m_{4\\rightarrow2}(x_{2})m_{5\\rightarrow2}(x_{2})}_{\\triangleq m_{2\\rightarrow1}(x_{1})}\\\\\n",
    " & = & \\frac{1}{Z}\\underbrace{\\phi_{1}(x_{1})m_{2\\rightarrow1}(x_{1})m_{3\\rightarrow1}(x_{1})}_{\\triangleq\\widetilde{p}_{X_{1}}(x_{1})},\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "where in the last line, we define $\\widetilde{p}_{X_{1}}(\\cdot )$\n",
    "to be a table that just corresponds to an unnormalized version of the marginal distribution of $X_{1}$. We can readily compute the normalization constant:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "Z & = & \\sum_{x_{1}}\\widetilde{p}_{X_{1}}(x_{1})\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Once we have computed this, we know the marginal distribution $p_{X_1}$ for $X_1$.\n",
    "\n",
    "### II.1.2—Computing Another Marginal Distribution in the Graph\n",
    "\n",
    "Now suppose that we want the marginal distribution $p_{X_{3}}$ for $X_{3}$. We can approach solving this problem the same way as how we computed the marginal distribution $p_{X_{1}}$ for $X_{1}$. But something nice happens, assuming that we've already computed the intermediate tables that we computed to obtain the marginal for $X_{1}$. In particular, we have\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{3}}(x_{3}) & = & \\sum_{x_{1}}\\sum_{x_{2}}\\sum_{x_{4}}\\sum_{x_{5}}\n",
    "\\bigg\\{ \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3})\\phi_{4}(x_{4})\\phi_{5}(x_{5}) \\\\\n",
    "&&\\qquad\\qquad\\qquad\\;\\;\\;\\;\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})\\psi_{24}(x_{2},x_{4})\\psi_{25}(x_{2},x_{5}) \\bigg\\} \\\\\n",
    " & = & \\sum_{x_{1}}\\sum_{x_{2}}\\sum_{x_{4}}\n",
    "\\bigg\\{ \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3})\\phi_{4}(x_{4})\\\\\n",
    "&&\\qquad\\qquad\\quad\\;\\;\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})\\psi_{24}(x_{2},x_{4})\\underbrace{\\sum_{x_{5}}\\phi_{5}(x_{5})\\psi_{25}(x_{2},x_{5})}_{m_{5\\rightarrow2}(x_{2})\\text{ - we already computed this}} \\bigg\\} \\\\\n",
    " & = & \\sum_{x_{1}}\\sum_{x_{2}}\n",
    "\\bigg\\{ \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{2}(x_{2})\\phi_{3}(x_{3}) \\\\\n",
    "&&\\qquad\\qquad\n",
    "\\cdot \\psi_{12}(x_{1},x_{2})\\psi_{13}(x_{1},x_{3})m_{5\\rightarrow2}(x_{2})\\underbrace{\\sum_{x_{4}}\\phi_{4}(x_{4})\\psi_{24}(x_{2},x_{4})}_{m_{4\\rightarrow2}(x_{2})\\text{ - also already computed}} \\bigg\\} \\\\\n",
    " & = & \\sum_{x_{1}}\n",
    "\\bigg\\{ \\frac{1}{Z}\\phi_{1}(x_{1})\\phi_{3}(x_{3}) \\\\\n",
    "&&\\qquad\\;\n",
    "\\cdot \\psi_{13}(x_{1},x_{3})\\underbrace{\\sum_{x_{2}}\\phi_{2}(x_{2})\\psi_{12}(x_{1},x_{2})m_{4\\rightarrow2}(x_{2})m_{5\\rightarrow2}(x_{2})}_{m_{2\\rightarrow1}(x_{1})\\text{ - already computed}} \\bigg\\} \\\\\n",
    " & = & \\frac{1}{Z}\\phi_{3}(x_{3})\\underbrace{\\sum_{x_{1}}\\phi_{1}(x_{1})\\psi_{13}(x_{1},x_{3})m_{2\\rightarrow1}(x_{1})}_{\\triangleq m_{1\\rightarrow3}(x_{3})}\\\\\n",
    " & = & \\frac{1}{Z}\\underbrace{\\phi_{3}(x_{3})m_{1\\rightarrow3}(x_{3})}_{\\triangleq\\widetilde{p}_{X_{3}}(x_{3})}.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "### II.1.3—The General Case: The Sum-Product Algorithm\n",
    "\n",
    "We are now ready to state the general algorithm for computing marginal distributions for every node in an undirected graphical model. This algorithm specifically works when the corresponding graph is a tree, which means that it has no loops and we can reach from any node in the graph to any other node by traversing along edges. The resulting general algorithm for marginalization is called the **sum-product algorithm** (and popularly also goes by the name **belief propagation**):\n",
    "\n",
    " 1. Choose a root note (arbitrarily) and identify corresponding leaf nodes.\n",
    " (In our earlier example, we choose node 1 as the root node; the leaf nodes then are nodes 3, 4, and 5.)\n",
    "\n",
    " 2. Proceed from the leaf nodes to the root node computing required messages along the way.\n",
    " (In our earlier example, we computed messages in the order $m_{5\\rightarrow 2}$, $m_{4\\rightarrow 2}$, $m_{3\\rightarrow 1}$, and $m_{2\\rightarrow 1}$.)\n",
    " \n",
    " When the root node is reached, reverse direction and calculate messages that go back to the leaf nodes\n",
    "(In our earlier example, we then computed messages $m_{1\\rightarrow 3}$ followed by $m_{1\\rightarrow 2}$; to further get the marginals at $X_{4}$ and $X_{5}$ compute $m_{2\\rightarrow 4}$ and $m_{2\\rightarrow 5}$.)\n",
    "\n",
    " The equation for computing the table of messages is given by\n",
    "$$\n",
    "m_{i\\rightarrow j}(x_{j})=\\sum _{x_{i}}\\bigg[\\phi _{i}(x_{i})\\psi _{i,j}(x_{i},x_{j})\\prod _{k\\in \\mathcal{N}(i)\\text { such that }k\\ne j}m_{k\\rightarrow i}(x_{i})\\bigg],\n",
    "$$\n",
    "where $\\mathcal{N}(i)$ denotes the neighboring nodes of node $i$ in the graph.\n",
    "\n",
    " 3. For each node $i$, use the incoming messages to compute the marginal distribution $p_{X_{i}}$:\n",
    " $$\n",
    "p_{X_{i}}(x_{i})=\\frac{1}{Z}\\underbrace{\\phi _{i}(x_{i})\\prod _{j\\in \\mathcal{N}(i)}m_{j\\rightarrow i}(x_{i})}_{\\widetilde{p}_{X_{i}}(x_{i})}.\n",
    "$$\n",
    "\n",
    "Typically, this marginal is computed by first computing the unnormalized table $\\widetilde{p}_{X_{i}}$ and then normalizing it; the normalization constant $Z$ is not saved (although it could be). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## II.1.4 — The Sum-Product Algorithm: Computational Complexity\n",
    "\n",
    "We study the case of running the sum-product algorithm on a tree.\n",
    "\n",
    "### II.1.4.1 — Rough calcul of Sum-Product algorithm complexity\n",
    "\n",
    "#### Computing messages tables\n",
    "- Compute a message table $m_{i\\rightarrow j}$ for a specific $i$ and $j$ (We considere that the messages that $m_{i\\rightarrow j}$ depends on have already been computed): $\\mathcal{O}(n k^2)$\n",
    "- How many messages do we compute? $2\\cdot (n-1)$\n",
    "- => how many operations does computing all the messages take? $\\mathcal{O}(n^2 k^2)$\n",
    "\n",
    "#### Compute the marginal distribution\n",
    "After the message tables have been computed, it takes $\\mathcal{O}(n k)$ to compute the marginal distribution $p_{X_ i}$ for a specific $i$.\n",
    "\n",
    "#### Running the sum-product algorithm\n",
    "The cost for running the sum-product algorithm is thus:\n",
    "$$\\mathcal{O}(n^2 k^2)$$\n",
    "\n",
    "### II.1.4.2 — Book-keeping Sum-Product algorithm complexity\n",
    "\n",
    "We can improve the running time of the algorithm by some careful bookkeeping!\n",
    "\n",
    "#### computing messages from leaves to the root\n",
    "\n",
    "For a non-root node $i$, we break the computation of $m_{i\\rightarrow \\pi (i)}$ into two steps.\n",
    "- First we compute the table $\\xi _ i$:\n",
    " $$\\xi _ i(x_ i) = \\phi _{i}(x_{i})\\prod _{\\ell \\in \\mathcal{N}(i)\\text { such that }\\ell \\ne \\pi (i)}m_{\\ell \\rightarrow i}(x_{i}).$$\n",
    " Time complexity: $\\mathcal{O}(d_ i k)$\n",
    "-  After computing the table $\\xi _ i$, we compute\n",
    " $$m_{i\\rightarrow \\pi (i)}(x_{\\pi (i)}) = \\sum _{x_{i}} \\psi _{i,\\pi (i)}(x_{i},x_{\\pi (i)})\\xi _ i(x_ i).$$\n",
    " Time complexity: $\\mathcal{O}(k^2)$\n",
    "\n",
    "For a tree with $n$ nodes, we call $d_ i$ is the number of edges that node $i$ participates in, and so $\\sum _{i=1}^ n d_ i$ is the sum of the number of edges that every node participates in, where we count each edge multiple times since multiple nodes can participate in the same edge.\n",
    "\n",
    "The total number of operations for computing messages from leaves to the root can be decomposed into the number of operations for computing all the $\\xi _ i$'s plus the number of operations for computing every message $m_{i\\rightarrow \\pi (i)}$ where $\\xi _ i$ is already computed.\n",
    "- number of operations for computing all the $\\xi _ i$'s: $\\mathcal{O}(nk)$\n",
    "- number of operations for computing all the messages of the form $m_{i\\rightarrow \\pi (i)}$ (where when computing the message from $i$ to $\\pi (i)$, we have already computed $\\xi _ i$): $\\mathcal{O}(nk^2)$\n",
    "- total number of operations for passing messages from leaves to the root: $\\boxed {\\mathcal{O}(nk^2)}$\n",
    "\n",
    "**Note**: computing messages from the leaves to the root does not require any modification to the sum-product (or the later detailed max-product algorithm) to produce a running time of $\\mathcal{O}(nk^2)$. The issue with the star graph only was an issue when we are passing messages from the root back to the leaves (hereafter).\n",
    "\n",
    "#### Computing the Marginal Distribution at the Root\n",
    "\n",
    "Let's denote the unnormalized marginal distribution for $X_ r$ where $r$ is the root node as\n",
    "$$\\widetilde{p}_{X_ r}(x_ r) \\triangleq \\phi _ r(x_ r) \\prod _{j \\in \\mathcal{N}(r)} m_{j\\rightarrow r}(x_ r).$$\n",
    "\n",
    "The number of operations for computing unnormalized marginal $\\widetilde{p}_{X_ r}$ is then $\\mathcal{O}(d_ r k)$.\n",
    "\n",
    "Of course, computing the normalized marginal $p_{X_ r}$ is straightforward as we just divide the entries of unnormalized marginal $\\widetilde{p}_{X_ r}$ by the sum of its entries, which costs $\\mathcal{O}(k)$ operations.\n",
    "\n",
    "#### Message Passing from the Root to the Leaves\n",
    "\n",
    "The basic idea is that we write the message passing equation now in terms of the unnormalized marginal of whichever node we are passing the message from.\n",
    "\n",
    "We start by passing messages from the root node to its neighboring nodes:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "m_{r\\rightarrow j}(x_{j})\n",
    "&=&\n",
    "\\sum_{x_{r}}\n",
    "  \\bigg[\\phi_{r}(x_{r})\\psi_{r,j}(x_{r},x_{j})\\prod_{\\ell \\in\\mathcal{N}(r)\\text{ such that }\\ell \\ne j}m_{\\ell \\rightarrow r}(x_{r})\\bigg] \\\\\n",
    "&=&\n",
    "\\sum_{x_{r}}\n",
    "  \\frac{\\psi_{r,j}(x_r,x_j) \\widetilde{p}_{X_r}(x_r)}\n",
    "       {m_{j\\rightarrow r}(x_r)},\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "where here is where we make the assumption that all the message table entries are strictly positive, so as to not run into a division by 0 issue. Importantly, we have already computed all the tables inside the summation!\n",
    "\n",
    "The total number of operations for computing table $m_{r\\rightarrow j}$ given that we already have tables $\\psi _{r,j}$, $\\widetilde{p}_{X_ r}$, and $m_{j\\rightarrow r}$ available, is then $\\mathcal{O}(k^2)$.\n",
    "\n",
    "Right after we compute $m_{r\\rightarrow j}$ for each neighbor $j$ of root node $r$, we can compute the unnormalized marginals for each neighbor node $j$:\n",
    "$$\n",
    "\\widetilde{p}_{X_ j}(x_ j) \\triangleq \\xi _ j(x_ j) m_{\\pi (j) \\rightarrow j}(x_ j).\n",
    "$$\t \n",
    "This is where the $\\xi _ i$'s that we computed earlier come back into play!\n",
    "\n",
    "The number of operations for computing unnormalized marginal $\\widetilde{p}_{X_ j}$ for a specific node $j$ is thus $\\mathcal{O}(k)$.\n",
    "\n",
    "At this point, we can actually just repeat the same idea since in general for any node $i$ and its neighbor $j$,\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "m_{i\\rightarrow j}(x_{j})\n",
    "&=&\n",
    "\\sum_{x_{i}}\n",
    "  \\frac{\\psi_{i,j}(x_i,x_j) \\widetilde{p}_{X_i}(x_i)}\n",
    "       {m_{j\\rightarrow i}(x_i)}.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "In particular, we keep passing messages, and after passing each message, we compute an unnormalized marginal (and also normalizing so that we get all the marginal distributions).\n",
    "\n",
    "Using this strategy outlined above, computing all the messages from the root back to the leaves and computing all the marginals take $\\mathcal{O}(nk^2)$.\n",
    "\n",
    "Thus, the total number of operations for passing messages from the root back to the leaves and computing the marginals is\n",
    "$$\\mathcal{O}(n k) + (n - 1) (\\mathcal{O}(k^2) + \\mathcal{O}(k) + \\mathcal{O}(k)) = \\boxed {\\mathcal{O}(n k^2)}.$$\n",
    "\n",
    "Overall, the running time is now linear in $n$ with careful storage of some extra tables and reordering some of the computations.\n",
    "\n",
    "In practice, how we presented the sum-product algorithm previously is cleaner to code up and is typically fast enough as we aren't dealing with bad graph cases such as the star graph. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II.2—The Max-Product Algorithm\n",
    "\n",
    "### II.2.1— The Key Idea Behind Sum-Product and Max-Product\n",
    "\n",
    "In deriving sum-product, the key idea was to push in summations. As a simple example, this amounts to the following:\n",
    "$$\\sum _{x_2}f(x_1)g(x_1,x_2) = f(x_1) \\sum _{x_2}g(x_1,x_2),$$\n",
    "for arbitrary functions $f$ and $g$.\n",
    "\n",
    "A similar idea applies when we want the maximum:\n",
    "$$\n",
    "\\max _{x_1, x_2}\\{ f(x_1)g(x_1,x_2)\\}  = \\max _{x_1}\\Big[ \\max _{x_2} f(x_1)g(x_1,x_2) \\Big] = \\max _{x_1}\\Big[ f(x_1) \\max _{x_2} g(x_1,x_2) \\Big],\n",
    "$$\n",
    "\n",
    "where the second equality holds provided that $f(x_1)$ is nonnegative for all $x_1$.\n",
    "\n",
    "For the equation above, think of the right-hand side function $f(x_1) \\max _{x_2} g(x_1,x_2)$ as some 1-dimensional table with a number associated with each value of $x_1$. (In particular, suppose we've precomputed what $\\max _{x_2} g(x_1, x_2)$ is for every value of $x_1$.) Then the argument $\\hat{x}_1$ achieving the maximum is whichever $x_1$ has the highest value in the table. In other words, the highest value is\n",
    "$$f(\\hat{x}_1) \\max _{x_2} g(\\hat{x}_1,x_2),$$\n",
    "\n",
    "from which we see that to get the optimal value of $x_2$, we just need to figure out which $x_2$ maximizes the function $g(\\hat{x}_1,x_2)$.\n",
    "\n",
    "So we first found the optimal value $\\hat{x}_1$ for $x_1$. Then we traced backward and found the best value for $\\hat{x}_2$ given $\\hat{x}_1$, where we could readily answer this query if when we first compute $\\max _{x_2} g(x_1, x_2)$, we also store what value of $x_2$ achieves the maximum for each value of $x_1$.\n",
    "\n",
    "The intuition for this two-variable example carries over to the case when we have many variables represented by a tree-structured undirected graphical model and we want the most probable configuration. The resulting algorithm is the max-product algorithm.\n",
    "\n",
    "### II.2.2— Max-Product\n",
    "\n",
    "As a reminder of notation, for random variables $X_1,\\dots ,X_ n$ corresponding to a tree $G=(V, E)$,\n",
    "$$\n",
    "p_{X_1,\\dots ,X_ n}(x_1,\\dots ,x_ n) = \\frac{1}{Z} \\Bigg( \\prod _{i\\in V} \\phi _ i(x_ i) \\Bigg) \\Bigg( \\prod _{(i,j)\\in E} \\psi _{i,j}(x_ i, x_ j) \\Bigg).\n",
    "$$\n",
    "\n",
    "The max-product algorithm is as follows. First, we pick a root node $r$ and compute *messages going from the leaves to the root*:\n",
    "$$\n",
    "m_{i \\to j}(x_ j) = \\max _{x_ i} \\left[ \\phi _ i(x_ i) \\psi _{i,j}(x_ i,x_ j) \\prod _{k \\in \\mathcal{N}(i) \\setminus \\{ j\\} } m_{k \\to i}(x_ i) \\right].\n",
    "$$\n",
    "\n",
    "Note that this is formula is the same as that of sum-product except that **the sum has been replaced by a max**.\n",
    "\n",
    "As with our simple two-variable example, we want to **keep track of which argument achieves the maximum**. Otherwise, all we'd be computing is a *scaled* version of the probability of the most probable configuration. (Remember that in general we don't know the normalization constant Z, in which case we won't actually find out the probability of the most probable configuration! Luckily, the argument achieving the maximum doesn't care what the value of Z is.)\n",
    "\n",
    "To keep track of which argument achieves the maximum, we compute *traceback messages*:\n",
    "$$\n",
    "t_{j \\to i}(x_ j) = \\arg \\max _{x_ i} \\left[ \\phi _ i(x_ i) \\psi _{ij}(x_ i,x_ j) \\prod _{k \\in \\mathcal{N}(i) \\setminus \\{ j\\} } m_{k \\to i}(x_ i) \\right].\n",
    "$$\n",
    "\n",
    "After computing the messages and traceback messages going from the leaves to the root, we're ready to find the most probable configuration, starting from the root and going to the leaves. This is like how we figured out the optimal value for $x_1$ in our two-variable example first (i.e., node 1 was the root).\n",
    "\n",
    "For root node $r$, we compute the optimal value:\n",
    "$$\n",
    "\\hat{x}_ r = \\arg \\max _{x_ r} \\left[\\phi _ r(x_ r) \\prod _{k \\in \\mathcal{N}(r)} m_{k \\to r}(x_ r) \\right].\n",
    "$$\n",
    "\n",
    "We then work backwards to the leaves, using the traceback messages to assign values:\n",
    "$$\n",
    "\\hat{x}_ i = t_{j \\to i}(\\hat{x}_ j).\n",
    "$$\n",
    "\n",
    "### II.2.3— Max-Product Complexity\n",
    "\n",
    "Computing max-product messages from leaves to the root has the same time complexity as computing sum-product messages from leaves to the root, which takes time $\\mathcal{O}(nk^2)$ as we saw previsouly (note that computing messages from the leaves to the root does not require any modification to the sum-product or the max-product algorithm to produce a running time of $\\mathcal{O}(nk^2)$, i.e., the issue we encountered with the star graph only was an issue when we were passing messages from the root back to the leaves).\n",
    "\n",
    "The reason why the running time is the same as for sum-product specifically for passing messages from the leaves to the root is that we are still iterating over the same number of elements doing maximization instead of summation, and separately we are also storing the argmax for each maximization, which can be done within the same for loop as the maximization itself.\n",
    "\n",
    "The main change in going from sum-product to max-product is precisely doing maximization instead of summation and keeping track of the argmax and storing the argmax in a table. The extra storage needed is roughly twice the amount we used before where now we keep track of, for each max-product message, also a traceback table.\n",
    "\n",
    "Following backpointers just involves going from the root to the leaves each time looking up a traceback table entry, so the time complexity is linear in the number of edges, which for a tree is $n-1$. Thus, following backpointers has running time $\\mathcal{O}(n).$\n",
    "\n",
    "Putting together the pieces above, max-product has running time $$\\mathcal{O}(nk^2) + \\mathcal{O}(d_ r k) + \\mathcal{O}(n) = \\boxed {\\mathcal{O}(nk^2)}. $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### II.2.4— Node max-marginal\n",
    "\n",
    "#### II.2.4.1— The Most Probable Values for Node Marginals\n",
    "\n",
    "Suppose we run the max-product algorithm to produce\n",
    "$$(\\widehat{x}_{1},\\widehat{x}_{2},\\dots ,\\widehat{x}_{n})=\\arg \\max _{x_{1},x_{2},\\dots ,x_{n}}p_{X_{1},X_{2},\\dots ,X_{n}}(x_{1},x_{2},\\dots ,x_{n}).\n",
    "$$\n",
    "\n",
    "Suppose we separately run the sum-product algorithm to produce $p_{X_{i}}$ for every $i$, and then we compute\n",
    "$$\n",
    "\\widehat{x}_{i}^{(\\text {sum-product})}=\\arg \\max _{x_{i}}p_{X_{i}}(x_{i}).\n",
    "$$\n",
    "\n",
    "In general, we don't always have $\\widehat{x}_{i}=\\widehat{x}_{i}^{(\\text {sum-product})}$ for every $i$.\n",
    "\n",
    "#### II.2.4.2— The Max-Marginals Variant of the Max-Product Algorithm\n",
    "\n",
    "The Max-Marginals is a variant of the max-product algorithm that more closely resembles the sum-product algorithm, where we do not keep track of traceback messages.\n",
    "\n",
    "In particular we pass max-product messages from leaves to the arbitrarily chosen root and then pass max-product messages from the root back to the leaves. Finally, for every node $i$, we compute what is called a node “max-marginal\" table $\\overline{p}_{X_{i}}$, which is a table where\n",
    "$$\n",
    "\\overline{p}_{X_{i}}(x_{i})=\\phi _{i}(x_{i})\\prod _{\\ell \\in \\mathcal{N}(i)}m_{\\ell \\rightarrow i}(x_{i}),\n",
    "$$\n",
    "and $m_{\\ell \\rightarrow i}$ here is a max-product message (not a sum-product message).\n",
    "\n",
    "Note that\n",
    "$$\n",
    "\\overline{p}_{X_{i}}(x_{i})\\propto \\underbrace{\\max _{x_{1},x_{2},\\dots ,x_{i-1},x_{i+1},\\dots ,x_{n}}}_{\\text {maximize over everything except }x_{i}}p_{X_{1},X_{2},\\dots ,X_{n}}(x_{1},x_{2},\\dots ,x_{n}).\n",
    "$$\n",
    "\n",
    "However, it is not always possible to readily tell which value of $x_{2}$ goes with the optimal value of $x_{1}$. We run into a problem when the argmax is not unique (for instance in a XOR table).\n",
    "\n",
    "The max-marginal variante is usefull when each Node Max-Marginal has a unique Most Probable Value: then traceback tables aren't needed.\n",
    "\n",
    "**Theorem:** If every node max-marginal $\\overline{p}_{X_ i}$ had a unique maximum probable value $x_ i^*$, i.e.,\n",
    "$$x_ i^* = \\arg \\max _{x_ i} \\overline{p}_{X_ i}(x_ i), \\qquad \\text {and} \\qquad \\overline{p}_{X_ i}(x_ i) < \\overline{p}_{X_ i}(x_ i^*) \\quad \\text {for all }x_ i \\ne x_ i^*.$$\n",
    "\n",
    "Then the joint probability table $p_{X_1,X_2,\\dots ,X_ n}$ has a unique most probable configuration given by $x\\triangleq (x_1^*,x_2^*,\\dots ,x_ n^*)$ (i.e., we put together the unique maximum probable values for each of the node max-marginals).\n",
    "\n",
    "**Proof**:\n",
    "Recall that node max-marginal $\\overline{p}_{X_ i}$ satisfies\n",
    "$$\\overline{p}_{X_{i}}(x_{i})\\propto \\underbrace{\\max _{x_{1},x_{2},\\dots ,x_{i-1},x_{i+1},\\dots ,x_{n}}}_{\\text {maximize over everything except }x_{i}}p_{X_{1},X_{2},\\dots ,X_{n}}(x_{1},x_{2},\\dots ,x_{n}).\n",
    "$$\n",
    "\n",
    "We start by supposing that there is a configuration $\\widehat{x}\\triangleq (\\widehat{x}_1, \\widehat{x}_2, \\dots , \\widehat{x}_ n)$ that is a most probable configuration, i.e.,\n",
    "$$p_{X_1,\\dots ,X_ n}(\\widehat{x}_1, \\dots , \\widehat{x}_ n) \\ge p_{X_1,\\dots ,X_ n}(x_1, \\dots , x_ n)\\qquad \\text {for all }x_1,\\dots ,x_ n,$$\n",
    "which also means that\n",
    "$$\n",
    "p_{X_1,\\dots ,X_ n}(\\widehat{x}_1, \\dots , \\widehat{x}_ n) \\ge p_{X_1,\\dots ,X_ n}(x_1^*, \\dots , x_ n^*).\n",
    "$$\n",
    "\n",
    "Then\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "&&\\overline{p}_{X_i}(\\widehat{x}_i) \\\\\n",
    "(\\text{since }\\widehat{x}\\text{ is a most probable configuration})\n",
    "&&= p_{X_1, \\dots, X_n}(\\widehat{x}_1, \\dots, \\widehat{x}_n) \\\\\n",
    "(\\text{since }\\widehat{x}\\text{ is a most probable configuration})\n",
    "&&= \\max_{x_1, \\dots, x_n}\n",
    "      p_{X_1, \\dots, X_n}(x_1, \\dots, x_n) \\\\\n",
    "&&= \\max_{x_i} \\max_{x_1, \\dots, x_{i-1}, x_{i+1}, x_n}\n",
    "      p_{X_1, \\dots, X_n}(x_1, \\dots, x_n) \\\\\n",
    "&&= \\max_{x_i} \\overline{p}_{X_i}(x_i).\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "Look at the very left-most side of this equation and the very right-most side: what the equation says is that $\\widehat{x}_ i$ achieves the highest possible value in the node max-marginal $p_{X_ i}$. And by assumption there is only one possible value for $X_ i$ that achieves this highest node max-marginal value: $x_ i^*$. So we conclude that indeed $\\widehat{x}_ i = x_ i^*$. \n",
    "\n",
    "**Practical implication:** If you run the max-marginal variant of the max-product algorithm that more closely resembles the sum product algorithm, and you check the node max-marginals to find that the most probable value for each of these max-marginals is unique, then you *don't need traceback messages and you also don't need to follow backpointers*.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II.2.5— Numerical Stability Issues: Max-Product to Min-Sum\n",
    "\n",
    "Messages are computed by multiplying many potential function values and then propagating these products onwards to compute more messages. Meanwhile, potential functions can often consist of small numbers. Multiplying many small numbers may lead to underflow issues on computers, where multiplication just yields 0 when the product isn't actually 0.\n",
    "\n",
    "A popular fix to this issue is to work with negative logs. Log turns small values into large values since values in $(0,1]$ get mapped to $(-\\infty ,0]$, and negating then turns these numbers positive. Large potential values aren't affected as much.\n",
    "\n",
    "So taking the negative log of our messages, we get the following equations:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{i \\to j}'(x_j)\n",
    "&=& \\min_{x_i}\n",
    "      \\left[\n",
    "        -\\log\\phi_i(x_i) - \\log\\psi_{ij}(x_i,x_j) +\n",
    "        \\sum_{k \\in \\mathcal{N}(i) \\setminus \\{j\\}} m_{k \\to i}'(x_i)\n",
    "      \\right], \\\\\n",
    "t_{j \\to i}(x_j)\n",
    "&=& \\arg\\min_{x_i}\n",
    "      \\left[\n",
    "        -\\log\\phi_i(x_i) - \\log\\psi_{ij}(x_i,x_j) +\n",
    "        \\sum_{k \\in \\mathcal{N}(i) \\setminus \\{j\\}} m_{k \\to i}'(x_i)\n",
    "      \\right], \\\\\n",
    "\\hat{x}_r\n",
    "&=& \\arg\\min_{x_r}\n",
    "      \\left[\n",
    "        -\\log\\phi_r(x_r) +\n",
    "        \\sum_{k \\in \\mathcal{N}(r)} m_{k \\to r}'(x_r)\n",
    "      \\right], \\\\\n",
    "\\hat{x}_i\n",
    "&=& t_{j \\to i}(\\hat{x}_j).\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "Following the usual lack of imagination in naming these algorithms, this algorithm is called the *min-sum algorithm *since it takes minimums of sums. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II.3— The Viterbi Algorithm: Max-product / min-sum specialized to HMMs\n",
    "\n",
    "Max-product or min-sum specialized to HMMs results in the Viterbi algorithm, proposed by MIT alum Andrew Viterbi back in 1967, which pre-dates both the forward-backward algorithm (1980) and the sum-product algorithm (1982). The Viterbi algorithm produces an MAP estimate for the sequence of hidden states given the observations of an HMM.\n",
    "\n",
    "Let's work through an example. I have two coins, one fair and one biased, and I keep picking a coin and flipping it. You get to observe whether each flip was heads or tails, and from that information, you have to guess the sequence of coins that I used.\n",
    "\n",
    "At first, I choose either the fair coin or the biased coin uniformly at random. I prefer not to switch between flips: that is, after each flip, my next flip will use the same coin with probability 3/4, and switch coins with probability 1/4. The fair coin is equally likely to be heads or tails, and the biased coin has probability 3/4 of coming up tails and probability 1/4 of coming up heads.\n",
    "\n",
    "If you observe the sequence $HHTTT$ (where $H$ is heads and $T$ is tails), what's the most likely sequence of coins that I used?\n",
    "\n",
    "We reuse notation from how we formulated HMM's when we presented the forward-backward algorithm. In particular, for this HMM with 5 hidden states, we define the following potentials:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "\\widetilde{\\phi}_i(x_i)\n",
    "&=&\\begin{cases}\n",
    "   p_{X_1}(x_1)p_{Y_1|X_1}(y_1|x_1) & \\text{for }i=1, \\\\\n",
    "   p_{Y_i|X_i}(y_i|x_i) & \\text{for }i=2,\\dots,5, \\\\\n",
    "   \\end{cases} \\\\\n",
    "\\psi(x_{i-1}, x_i)\n",
    "&=&p_{X_i \\mid X_{i-1}}(x_i \\mid x_{i-1})\n",
    "\\quad\\text{for }n=2,\\dots,5.\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "Thus, we have (importantly, observations HHTTT have been folded into the node potentials!):\n",
    "![images_sec-viterbi-example1](./images/images_sec-viterbi-example1.png)\n",
    "\n",
    "Defining $\\gamma \\triangleq \\log _2(3) \\approx 1.6$, we have:\n",
    "![images_sec-viterbi-example2](./images/images_sec-viterbi-example2.png)\n",
    "\n",
    "Below, we show what is called the *trellis diagram* for this problem: note that while this diagram has nodes and edges, it is not a probabilistic graphical model in that the nodes do not correspond to random variables, and the nodes/edges aren't associated with potential tables, etc. Instead, the i-th column consists of all the possible states that $X_ i$ can be, and the edges denote all possible transitions, so that a path going from left to right (a path in a trellis diagram can't go rightward and then leftward; it can only go rightward) corresponds to a specific possible sequence of states that $X_1,\\dots ,X_ n$ can take on.\n",
    "![images_sec-viterbi-trellis1](./images/images_sec-viterbi-trellis1.png)\n",
    "The trellis diagram for our example\n",
    "\n",
    "We will see that the Viterbi algorithm just finds a path along this trellis.\n",
    "\n",
    "For a length-$n$ HMM, the min-sum messages are\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{1 \\rightarrow 2}'(x_2)\n",
    "&=& \\min_{x_1}\n",
    "      \\left[\n",
    "        -\\log_2\\widetilde{\\phi}_1(x_1) - \\log_2\\psi(x_1,x_2)\n",
    "      \\right], \\\\\n",
    "t_{2 \\rightarrow 1}(x_2)\n",
    "&=& \\arg\\min_{x_1}\n",
    "      \\left[\n",
    "        -\\log_2\\widetilde{\\phi}_1(x_1) - \\log_2\\psi(x_1,x_2)\n",
    "      \\right], \\\\\n",
    "m_{(i-1) \\rightarrow i}'(x_i)\n",
    "&=& \\min_{x_{i-1}}\n",
    "      \\left[\n",
    "        -\\log_2\\widetilde{\\phi}_{i-1}(x_{i-1})\n",
    "        - \\log_2\\psi(x_{i-1},x_i)\n",
    "        + m_{(i-2) \\rightarrow (i-1)}'(x_{i-1})\n",
    "      \\right]\n",
    "& & \\text{for }i=3,\\dots,n, \\\\\n",
    "t_{i \\rightarrow (i-1)}(x_i)\n",
    "&=& \\arg\\min_{x_{i-1}}\n",
    "      \\left[\n",
    "        -\\log_2\\widetilde{\\phi}_{i-1}(x_{i-1})\n",
    "        - \\log_2\\psi(x_{i-1},x_i)\n",
    "        + m_{(i-2) \\rightarrow (i-1)}'(x_{i-1})\n",
    "      \\right]\n",
    "& & \\text{for }i=3,\\dots,n.\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "Note that $m_{1 \\rightarrow 2}'(x_2)$ could be re-written in the same way as $m_{(i-1) \\rightarrow i}'(x_i)$ by including the prior distribution in a previous message $m_{0 \\rightarrow 1}'(x_1)$ instead of $\\widetilde{\\phi}_1(x_1)$.\n",
    "\n",
    "Let's determine the messages for our length $n=5$ HMM.\n",
    "\n",
    "We compute the first message:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{1\\rightarrow 2}(\\text{fair})\n",
    "&=& \\min_{x_1}\n",
    "      \\big[ -\\log_2\\widetilde\\phi_1(x_1) - \\log_2(\\psi(x_1, \\text{fair})) \\big] \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\overbrace{-\\log_2\\widetilde\\phi_1(\\text{fair})\n",
    "          - \\log_2(\\psi(\\text{fair}, \\text{fair}))}^{x_1 = \\text{fair}}, \\\\\n",
    "&&\\qquad\\;\\;\n",
    "        \\underbrace{-\\log_2\\widetilde\\phi_1(\\text{biased})\n",
    "          - \\log_2(\\psi(\\text{biased}, \\text{fair}))}_{x_1 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{2 + (2 - \\gamma)}_{x_1 = \\text{fair}},\n",
    "        \\underbrace{3 + 2}_{x_1 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{4 - \\gamma}_{x_1 = \\text{fair}},\n",
    "        \\underbrace{5}_{x_1 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& 4 - \\gamma,\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "where the $\\arg \\min is x_1 = \\text {fair}$, so $t_{2\\rightarrow 1}(\\text {fair}) = \\text {fair}$.\n",
    "\n",
    "Next:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{1\\rightarrow 2}(\\text{biased})\n",
    "&=& \\min_{x_1}\n",
    "      \\big[ -\\log_2\\widetilde\\phi_1(x_1) - \\log_2(\\psi(x_1, \\text{biased})) \\big] \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\overbrace{-\\log_2\\widetilde\\phi_1(\\text{fair})\n",
    "          - \\log_2(\\psi(\\text{fair}, \\text{biased}))}^{x_1 = \\text{fair}}, \\\\\n",
    "&&\\qquad\\;\\;\n",
    "        \\underbrace{-\\log_2\\widetilde\\phi_1(\\text{biased})\n",
    "          - \\log_2(\\psi(\\text{biased}, \\text{biased}))}_{x_1 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{2 + 2}_{x_1 = \\text{fair}},\n",
    "        \\underbrace{3 + (2 - \\gamma)}_{x_1 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{4}_{x_1 = \\text{fair}},\n",
    "        \\underbrace{5 - \\gamma}_{x_1 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& 5 - \\gamma,\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "where the $\\arg \\min is x_1 = \\text {biased}$, so $t_{2\\rightarrow 1}(\\text {biased}) = \\text {biased}$.\n",
    "\n",
    "Our message and traceback tables are then:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{1 \\rightarrow 2}'(x_2)\n",
    "&=& \\begin{cases}\n",
    "      4-\\gamma & \\text{if }x_2 = \\text{fair} \\\\\n",
    "      5-\\gamma & \\text{if }x_2 = \\text{biased}\n",
    "    \\end{cases} \\\\\n",
    "t_{2 \\rightarrow 1}(x_2)\n",
    "&=& \\begin{cases}\n",
    "      \\text{fair} & \\text{if }x_2 = \\text{fair} \\\\\n",
    "      \\text{biased} & \\text{if }x_2 = \\text{biased}\n",
    "    \\end{cases}\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "This traceback table tells us how to select $x_1$ once we decide on $x_2$. This is illustrated in the first-step trellis diagram below:\n",
    "![images_sec-viterbi-trellis2](./images/images_sec-viterbi-trellis2.png)\n",
    "The first traceback\n",
    "\n",
    "The second message is computed similarly:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{2\\rightarrow 3}'(\\text{fair})\n",
    "&=& \\min_{x_2}\n",
    "      \\big[ -\\log_2\\widetilde\\phi_2(x_2) - \\log_2(\\psi(x_2, \\text{fair}))\n",
    "            + m_{1\\rightarrow 2}'(x_2) \\big] \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\overbrace{-\\log_2\\widetilde\\phi_2(\\text{fair})\n",
    "          - \\log_2(\\psi(\\text{fair}, \\text{fair}))\n",
    "          + m_{1\\rightarrow 2}'(\\text{fair})}^{x_2 = \\text{fair}}, \\\\\n",
    "&&\\qquad\\;\\;\n",
    "        \\underbrace{-\\log_2\\widetilde\\phi_2(\\text{biased})\n",
    "          - \\log_2(\\psi(\\text{biased}, \\text{fair}))\n",
    "          + m_{1\\rightarrow 2}'(\\text{biased})}_{x_2 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{1 + (2 - \\gamma) + (4 - \\gamma)}_{x_2 = \\text{fair}},\n",
    "        \\underbrace{2 + 2 + (5 - \\gamma)}_{x_2 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{7 - 2\\gamma}_{x_2 = \\text{fair}},\n",
    "        \\underbrace{9 - \\gamma}_{x_2 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& 7 - 2\\gamma,\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "where the $\\arg \\min is x_2 = \\text {fair}$, so $t_{3\\rightarrow 2}(\\text {fair}) = \\text {fair}$.\n",
    "\n",
    "Next:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{2\\rightarrow 3}'(\\text{biased})\n",
    "&=& \\min_{x_2}\n",
    "      \\big[ -\\log_2\\widetilde\\phi_2(x_2) - \\log_2(\\psi(x_2, \\text{biased}))\n",
    "            + m_{1\\rightarrow 2}'(x_2) \\big] \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\overbrace{-\\log_2\\widetilde\\phi_2(\\text{fair})\n",
    "          - \\log_2(\\psi(\\text{fair}, \\text{biased}))\n",
    "          + m_{1\\rightarrow 2}'(\\text{fair})}^{x_2 = \\text{fair}}, \\\\\n",
    "&&\\qquad\\;\\;\n",
    "        \\underbrace{-\\log_2\\widetilde\\phi_2(\\text{biased})\n",
    "          - \\log_2(\\psi(\\text{biased}, \\text{biased}))\n",
    "          + m_{1\\rightarrow 2}'(\\text{biased})}_{x_2 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{1 + 2 + (4 - \\gamma)}_{x_2 = \\text{fair}},\n",
    "        \\underbrace{2 + (2 - \\gamma) + (5 - \\gamma)}_{x_2 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& \\min\n",
    "      \\big\\{\n",
    "        \\underbrace{7 - \\gamma}_{x_2 = \\text{fair}},\n",
    "        \\underbrace{9 - 2\\gamma}_{x_2 = \\text{biased}}\n",
    "      \\big\\} \\\\\n",
    "&=& 7 - \\gamma,\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "where the $\\arg \\min is x_2 = \\text {fair}$, so $t_{3\\rightarrow 2}(\\text {biased}) = \\text {fair}$.\n",
    "\n",
    "Our message and traceback tables are then:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{2 \\rightarrow 3}'(x_3)\n",
    "&=& \\begin{cases}\n",
    "      7-2\\gamma & \\text{if }x_3 = \\text{fair} \\\\\n",
    "      7-\\gamma & \\text{if }x_3 = \\text{biased}\n",
    "    \\end{cases} \\\\\n",
    "t_{3 \\rightarrow 2}(x_3)\n",
    "&=& \\begin{cases}\n",
    "      \\text{fair} & \\text{if }x_3 = \\text{fair} \\\\\n",
    "      \\text{fair} & \\text{if }x_3 = \\text{biased}\n",
    "    \\end{cases}\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "The traceback is illustrated below:\n",
    "![images_sec-viterbi-trellis3](./images/images_sec-viterbi-trellis3.png)\n",
    "The second traceback\n",
    "\n",
    "We see that after observing wo heads in a row, we will eventually pick $x_2 = fair$ no matter what $x_3$ is.\n",
    "\n",
    "Make sure that you can also compute these messages! The rest of the message and traceback tables are as follows:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "m_{3 \\rightarrow 4}'(x_4)\n",
    "&=& \\begin{cases}\n",
    "      10-3\\gamma & \\text{if }x_4 = \\text{fair} \\\\\n",
    "      11-3\\gamma & \\text{if }x_4 = \\text{biased}\n",
    "    \\end{cases} \\\\\n",
    "t_{4 \\rightarrow 3}(x_4)\n",
    "&=& \\begin{cases}\n",
    "      \\text{fair} & \\text{if }x_4 = \\text{fair} \\\\\n",
    "      \\text{biased} & \\text{if }x_4 = \\text{biased}\n",
    "    \\end{cases} \\\\\n",
    "m_{4 \\rightarrow 5}'(x_5)\n",
    "&=& \\begin{cases}\n",
    "      13-4\\gamma & \\text{if }x_5 = \\text{fair} \\\\\n",
    "      15-5\\gamma & \\text{if }x_5 = \\text{biased}\n",
    "    \\end{cases} \\\\\n",
    "t_{5 \\rightarrow 4}(x_5)\n",
    "&=& \\begin{cases}\n",
    "      \\text{fair} & \\text{if }x_5 = \\text{fair} \\\\\n",
    "      \\text{biased} & \\text{if }x_5 = \\text{biased}\n",
    "    \\end{cases} \\\\\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "The traceback is illustrated below:\n",
    "![images_sec-viterbi-trellis5](./images/images_sec-viterbi-trellis5.png)\n",
    "The final traceback\n",
    "\n",
    "Finally, we can compute the most likely value for $X_5$:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "\\hat{x}_5 &=&\n",
    "\\arg\\min \\Big( \\underbrace{1 + (13 - 4\\gamma) }_{x_5 = \\text{fair}} , \\underbrace{(2-\\gamma) + (15-5\\gamma)}_{x_5 = \\text{biased}} \\Big)\n",
    "    = \\text{biased}\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "We can see from the trellis diagram that as soon as we decide on a value for $X_5$, we just have to follow our “trail of breadcrumbs\" from the traceback messages to get the path highlighted below:\n",
    "![images_sec-viterbi-trellis6](./images/images_sec-viterbi-trellis6.png)\n",
    "The result of the Viterbi algorithm\n",
    "\n",
    "\n",
    "We can do the same thing with the traceback messages:\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "\\hat{x}_4 &=& t_{5 \\rightarrow 4}(\\hat{x}_5)\n",
    "          = \\text{biased} \\\\\n",
    "\\hat{x}_3 &=& t_{4 \\rightarrow 3}(\\hat{x}_4)\n",
    "          = \\text{biased} \\\\\n",
    "\\hat{x}_2 &=& t_{3 \\rightarrow 2}(\\hat{x}_3)\n",
    "          = \\text{fair} \\\\\n",
    "\\hat{x}_1 &=& t_{2 \\rightarrow 1}(\\hat{x}_2)\n",
    "          = \\text{fair} \\\\\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "We conclude that the MAP estimate for the sequence of hidden states given observed sequence $HHTTT$ is “fair, fair, biased, biased, biased\".\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II.4— The sum-product algorithm applied to Hidden Markov Models (HMM)\n",
    "\n",
    "### II.4.1— Introduction to Hidden Markov Models (HMM's)\n",
    "\n",
    "We now show the sum-product algorithm applied to the special case of a tree-structured graphical model called a hidden Markov model (HMM), for which we have hidden states $X_{1},X_{2},\\dots ,X_{n}$ that we infer given observations $Y_{1},Y_{2},\\dots ,Y_{n}$, and the graphical model looks as follows:\n",
    "![images_sec-graphical-models-hmm](./images/images_sec-graphical-models-hmm.png)\n",
    "\n",
    "For example, in a basic robot localization problem setup, $X_{i}$ is the robot's location at time point $i$, and $Y_{i}$ is the robot's sensor readings at time point $i$. Then $p_{X_{i}\\mid Y_{1},\\dots ,Y_{n}}(\\cdot \\mid y_{1},\\dots ,y_{n})$ is the posterior distribution of the robot's location at time point $i$.\n",
    "\n",
    "The graphical model is\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{1},\\dots ,X_{n},Y_{1},\\dots ,Y_{n}}(x_{1},\\dots ,x_{n},y_{1},\\dots ,y_{n}) \\\\\n",
    "\\propto \\prod _{i=1}^{n}\\big \\{ \\phi _{X_{i}}(x_{i})\\phi _{Y_{i}}(y_{i})\\psi _{X_{i},Y_{i}}(x_{i},y_{i})\\big \\} \\prod _{i=1}^{n-1}\\psi _{X_{i},X_{i+1}}(x_{i},x_{i+1}).\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Note that by incorporating observations, we are left with a new graphical model that is only over the $X_{i}$'s:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{1},\\dots ,X_{n}\\mid Y_{1},\\dots ,Y_{n}}(x_{1},\\dots ,x_{n}\\mid y_{1},\\dots ,y_{n}) \\\\\t \n",
    "\\propto \\prod _{i=1}^{n}\\big \\{ \\underbrace{\\phi _{X_{i}}(x_{i})\\phi _{Y_{i}}(y_{i})\\psi _{X_{i},Y_{i}}(x_{i},y_{i})}_{\\triangleq \\widetilde{\\phi }_{i}(x_{i})}\\big \\} \\prod _{i=1}^{n-1}\\underbrace{\\psi _{X_{i},X_{i+1}}(x_{i},x_{i+1})}_{\\triangleq \\psi _{i,i+1}(x_{i},x_{i+1})}.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "![images_sec-graphical-models-markov-chain](./images/images_sec-graphical-models-markov-chain.jpg)\n",
    "\n",
    "We can run the sum-product algorithm choosing the right-most node $X_{n}$ as the root, so the initial message passing happens from the single leaf node $X_{1}$ going rightward until reaching $X_{n}$ (this is called the forward pass), and then we pass messages from the root node $X_{n}$ all the way back to $X_{1}$ (this is called the backward pass).\n",
    "\n",
    "The resulting algorithm for this specific case of HMM's is called the **forward-backward algorithm.**\n",
    "\n",
    "Hidden Markov models are widely used in practice for more than just robot localization (note that a robot could be, for instance, a self-driving car)! But a few other examples include recognizing speech, labeling the part of speech of each word in a sentence, aligning DNA sequences, detecting protein folds, and decoding a certain code in digital communications called a convolutional code.\n",
    "\n",
    "### II.4.2— Hidden Markov Models: Three Ingredients\n",
    "\n",
    "#### II.4.2.1— Notations\n",
    "\n",
    "Often for HMM's, a problem is modeled with a prior distribution $p_{X_{1}}$ on the initial state $X_{1}$, a transition model $p_{X_{i+1}\\mid X_{i}}$ that is the same for all $i$ (so we can denote it as just $p_{X_{\\text {next}}\\mid X_{\\text {current}}}$), and an observation model $p_{Y_{i}\\mid X_{i}}$ that is the same for all $i$ (so we can denote it as just $p_{Y\\mid X}$), in which case\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{1},\\dots ,X_{n},Y_{1},\\dots ,Y_{n}}(x_{1},\\dots ,x_{n},y_{1},\\dots ,y_{n}) \\\\\t \t \n",
    "=p_{X_{1}}(x_{1})\\bigg\\{ \\prod _{i=1}^{n}p_{Y\\mid X}(y_{i}\\mid x_{i})\\bigg\\} \\bigg\\{ \\prod _{i=1}^{n-1}p_{X_{\\text {next}}\\mid X_{\\text {current}}}(x_{i+1}\\mid x_{i})\\bigg\\} \\\\\n",
    "\\end{eqnarray}\n",
    "$$. \t \t \n",
    "\n",
    "Then as a graphical model, incorporating observations $Y_{1}=y_{1},\\dots ,Y_{n}=y_{n}$,\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p_{X_{1},\\dots ,X_{n}\\mid Y_{1},\\dots ,Y_{n}}(x_{1},\\dots ,x_{n}\\mid y_{1},\\dots ,y_{n}) \\\\\t \t \n",
    "\\propto \\underbrace{p_{X_{1}}(x_{1})p_{Y\\mid X}(y_{1}\\mid x_{1})}_{\\triangleq \\widetilde{\\phi }_{1}(x_{1})}\\prod _{i=2}^{n}\\big \\{ \\underbrace{p_{Y\\mid X}(y_{i}\\mid x_{i})}_{\\triangleq \\widetilde{\\phi }_{i}(x_{i})}\\big \\} \\prod _{i=1}^{n-1}\\underbrace{p_{X_{\\text {next}}\\mid X_{\\text {current}}}(x_{i+1}\\mid x_{i})}_{\\triangleq \\psi _{i,i+1}(x_{i},x_{i+1})}\n",
    "\\end{eqnarray}\n",
    "$$. \t \t \n",
    "\n",
    "The messages are:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "m_{1\\rightarrow 2}(x_{2}) \t& = \t& \\sum _{x_{1}}\\widetilde{\\phi }_{1}(x_{1})\\psi _{1,2}(x_{1},x_{2}) \\\\ \t \t \n",
    "\\text {For }i=2,\\dots ,n-1:\\qquad m_{i\\rightarrow (i+1)}(x_{i+1}) \t& = \t& \\sum _{x_{i}}\\widetilde{\\phi }_{i}(x_{i})\\psi _{i,i+1}(x_{i},x_{i+1})m_{(i-1)\\rightarrow i}(x_{i}) \t\\\\ \t \n",
    "m_{n\\rightarrow n-1}(x_{n-1}) \t& = \t& \\sum _{x_{n}}\\widetilde{\\phi }_{n}(x_{n})\\psi _{n-1,n}(x_{n-1},x_{n})\\\\ \t \t \n",
    " \\text {For }i=2,\\dots ,n-1:\\qquad m_{i\\rightarrow (i-1)}(x_{i-1}) \t& = \t& \\sum _{x_{i}}\\widetilde{\\phi }_{i}(x_{i})\\psi _{i-1,i}(x_{i-1},x_{i})m_{(i+1)\\rightarrow i}(x_{i})\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "The posterior distributions for each $X_{i}$ given all the observations $Y_{1}=y_{1},\\dots ,Y_{n}=y_{n}$ are:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    " p_{X_{1}\\mid Y_{1},\\dots ,Y_{n}}(x_{1}\\mid y_{1},\\dots ,y_{n}) \t& \\propto \t& \\widetilde{\\phi }_{1}(x_{1})m_{2\\rightarrow 1}(x_{1}) \t \t \\\\\n",
    "text {For }i=2,\\dots ,n-1:\\qquad p_{X_{i}\\mid Y_{1},\\dots ,Y_{n}}(x_{i}\\mid y_{1},\\dots ,y_{n}) \t& \\propto \t& \\widetilde{\\phi }_{i}(x_{i})m_{(i-1)\\rightarrow i}(x_{i})m_{(i+1)\\rightarrow i}(x_{i}) \t \t \\\\\n",
    "p_{X_{n}\\mid Y_{1},\\dots ,Y_{n}}(x_{n}\\mid y_{1},\\dots ,y_{n}) \t& \\propto \t& \\widetilde{\\phi }_{n}(x_{n})m_{n-1\\rightarrow n}(x_{n}) \t \t \n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "We can define conveniant notations as follows:\n",
    "\n",
    "-  The forward messages are denoted by $\\alpha _{i\\rightarrow (i+1)}(x_{i+1})\\triangleq m_{i\\rightarrow (i+1)}(x_{i+1})$\n",
    "\n",
    "- The backward messages are denoted by $\\beta _{(i+1)\\rightarrow i}(x_{i})\\triangleq m_{(i+1)\\rightarrow i}(x_{i})$\n",
    "\n",
    "- Since $\\psi _{i,i+1}$ is actually the same table/function for all $i$, we just call it $\\psi$, i.e., $\\psi (x_{i},x_{i+1})=p_{X_{\\text {next}}\\mid X_{\\text {current}}}(x_{i+1}\\mid x_{i})$\n",
    "\n",
    "- Since $\\widetilde{\\phi }_{i}$ above is actually the same for all $i\\ge 2$, we just call it $\\widetilde{\\phi }$, i.e., $\\widetilde{\\phi }(x_{i})=p_{Y\\mid X}(y_{i}\\mid x_{i}$); we actually extend this to the case $i=1$ as well with the following notational convenience: we define $\\alpha _{0\\rightarrow 1}(x_{1})\\triangleq p_{X_{1}}(x_{1})$ (you can think of table $\\alpha _{0\\rightarrow 1}$ as a message being sent from a non-existent node 0 to node 1 that corresponds to the prior distribution $p_{X_1}$), so that the forward messages can now be expressed in 1 formula:\n",
    "    $$\\alpha _{i\\rightarrow (i+1)}(x_{i+1})=\\sum _{x_{i}}\\widetilde{\\phi }(x_{i})\\psi (x_{i},x_{i+1})\\alpha _{(i-1)\\rightarrow i}(x_{i})\\qquad \\text {for all }i=1,2,\\dots ,n-1$$\n",
    "    \t \n",
    "- As another notational convenience, we define $\\beta _{(n+1)\\rightarrow n}(x_{n})=1$ (you can think of table $\\beta _{(n+1)\\rightarrow n}$ as a message being sent from a non-existent node $n+1$ to node $n$ that consists of all 1's), so that the backward messages can now be expressed in 1 formula:\n",
    "$$\\beta _{i\\rightarrow (i-1)}(x_{i-1})=\\sum _{x_{i}}\\widetilde{\\phi }(x_{i})\\psi (x_{i-1},x_{i})\\beta _{(i+1)\\rightarrow i}(x_{i})\\qquad \\text {for all }i=2,3,\\dots ,n$$\n",
    "\n",
    "Clarification: In terms of edges in the graphical model, any particular state variable $X_ k$ has edges only to $X_{k-1}$ (what came before it), $X_{k+1}$ (what comes after it), and $Y_ k$ (the observation associated with it). When we condition on all three of these, $X_ k$ becomes independent of everything else in the graph. However, when we don't condition on anything, then from any node we can reach any other node, so at least from looking at the graph, it could very well be that nodes far apart depend on each other. Even if we condition on all the $Y_ i$'s and nothing else, in general, all the $X_ i$'s can be related to each other!\n",
    "\n",
    "#### II.4.2.2— Study\n",
    "\n",
    "As a reminder, here is the graphical model for an HMM:\n",
    "![images_sec-graphical-models-hmm](./images/images_sec-graphical-models-hmm.png)\n",
    "\n",
    "Here's a quick recap of the important facts:\n",
    "\n",
    "-  We observe $Y_1$ through $Y_ n$, which we model as coming from hidden states $X_1$ through $X_ n$\n",
    "\n",
    "- The goal of the forward-backward algorithm is to find the conditional distribution over hidden states given the data.\n",
    "\n",
    "- In order to specify an HMM, we need three quantities:\n",
    "\n",
    "    1. A *transition distribution*, $p_{X_{k+1}|X_ k}(x_{k+1}|x_ k)$, which describes the distribution for the next state given the current state. This is often represented as a matrix that we'll call $A$. Rows of $A$ correspond to the current state, columns correspond to the next state, and each entry corresponds to the transition probability. That is, the entry at row $i$ and column $j$, $A_{ij}$, is $p_{X_{k+1}|X_ k}(j|i)$.\n",
    "    \n",
    "    2. An *observation distribution* (also called an “emission distribution\") $p_{Y_ k|X_ k}(y_ k|x_ k)$, which describes the distribution for the output given the current state. We'll represent this with matrix $B$. Here, rows correspond to the current state, and columns correspond to the observation. So, $B_{ij} = p_{Y_ k|X_ k}(j|i)$: the probability of observing output $j$ from state $i$ is $B_{ij}$. Since the number of possible observations isn't necessarily the same as the number of possible states, $B$ won't necessarily be square.\n",
    "\n",
    "    3. An *initial state distribution* $p_{X_1}(x_1)$, which describes the starting distribution over states. We'll represent this with a vector called $\\pi _0$, where item $i$ in the vector represents $p_{X_1}(i)$. \n",
    "\n",
    "- We compute forward and backwards messages as follows:\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\alpha _{(k-1)\\to k}(x_{k}) \t& = \t& \\sum _{x_{k-1}} \\overbrace{\\alpha _{(k-2)\\to (k-1)}(x_{k-1})}^{\\text {previous message}} \\overbrace{\\tilde{\\phi }(x_{k-1})}^{\\text {obs.}} \\overbrace{\\psi (x_{k-1}, x_{k})}^{\\text {transition}} \t\\\\ \t \n",
    "\\beta _{(k+1)\\to k}(x_{k}) \t& = \t& \\sum _{x_{k+1}} \\underbrace{\\beta _{(k+2)\\to (k+1)}(x_{k+1})}_{\\text {previous message}} \\underbrace{\\tilde{\\phi }(x_{k+1})}_{\\text {obs.}} \\underbrace{\\psi (x_{k}, x_{k+1})}_{\\text {transition}} \t\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "    These messages are illustrated below:\n",
    "![images_sec-graphical-models-hmm-messages](./images/images_sec-graphical-models-hmm-messages.png)\n",
    "\n",
    "    Here, $\\tilde{\\phi }(x_ k) = p_{Y_ k|X_ k}(y_ k|x_ k)$ and $\\psi (x_ k, x_{k+1}) = p_{X_{k+1}|X_ k}(x_{k+1}|x_ k)$. The first forward message $\\alpha _{0 \\to 1}(x_1)$ is initialized to $\\pi _0(x_1) = p_{X_1}(x_1)$. The first backward message $\\beta _{(n+1)\\to n}(x_ n)$ is initialized to uniform (this is equivalent to not including it at all).\n",
    "\n",
    "    The picture below illustrates the computation of one forward message $\\alpha _{2 \\to 3}(x_3)$.\n",
    "![images_sec-graphical-models-hmm-forward-message](./images/images_sec-graphical-models-hmm-forward-message.png)\n",
    "    In order for node 2 to summarize its belief about $X_3$, it must incorporate the previous message $\\alpha _{1 \\to 2}(x_2)$, its observation $\\tilde{\\phi }(x_2)$, and the relationship $\\psi (x_2,x_3)$ between $X_2$ and $X_3$.\n",
    "\n",
    "- To obtain a marginal distribution for a particular state given all the observations, $p_{X_ k|Y_1,\\ldots ,Y_ n}$, we simply multiply the incoming messages together with the observation term, and then normalize:\n",
    "$$ p_{X_ k | Y_1, \\ldots , Y_ n}(x_ k|y_1,\\ldots ,y_ n) \t \\propto \t \\alpha _{(k-1)\\to k}(x_{k})\\beta _{(k+1)\\to k}(x_{k}) \\tilde{\\phi }(x_ k) \t \t \n",
    "$$\n",
    "\n",
    "    Here, the symbol $\\propto$ means “is proportional to\", and indicates that we must normalize at the end so that the answer sums to 1. \n",
    "\n",
    "**Technical note:** The original version of the forward-backward algorithm was developed before the sum-product algorithm as well as the general framework of probabilistic graphical models. In the original version of the forward-backward algorithm, the backward messages $\\beta$ are computed as described here, but the forward messages $\\alpha$ are computed differently as $\\alpha ^\\prime _{i \\to (i+1)}(x_{i+1}) = \\alpha _{i \\to (i+1)}(x_{i+1}) \\tilde{\\phi }(x_{i+1})$. In general, this means that the original forward message to $X_ i$ would also include the observation $Y_ i$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}