An Evening with EM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An Evening with… EM!\n",
    "===\n",
    "\n",
    "The EM algorithm ([Dempster/Laird/Rubin, 1977](#refDempster)) is one of the most widely used unsupervised learning methods in NLP. It is important to have a solid understanding of its properties and limitations in order to use it well. However, it is quite complex and a bit tricky, and can easily be confusing for the beginner. There are a lot of terms, a lot of implementations, and a lot of Greek letters to add to the confusion. This tutorial tries to avoid the confusion by providing the step-by-step implementation of a concrete case (with pseudocode and a [Python implementation](#secrunnable)), highlighting only the necessary terms and providing layman’s explanations of the underlying concepts.\n",
    "\n",
    "Table of contents\n",
    "--\n",
    "1. [Introduction](#secintro)\n",
    "2. [Preliminaries](#secpreliminaries)\n",
    "3. [Example Uses of EM](#secexample)\n",
    "4. [The Goal](#secgoal)\n",
    "5. [Implementation](#secimplementation)\n",
    "6. [Finally...](#secfinally)\n",
    "7. [Troubleshotting](#sectrouble)\n",
    "8. [Useful Reading](#secreading)\n",
    "9. [Advanced Topics](#secadvanced)\n",
    "10. [Going Further...](#secfurther)\n",
    "\n",
    "\n",
    "+ [Acknowledgements](#secacknowledgements)\n",
    "+ [References](#secref)\n",
    "+ [Answers](#secanswers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "from collections import defaultdict\n",
    "import time\n",
    "from functools import reduce\n",
    "\n",
    "sp.set_printoptions( precision = 3 )\n",
    "\n",
    "# CONSTANTS\n",
    "NINF = float('-1e303')\n",
    "LINF = sp.log(NINF)\n",
    "SMOOTHING = sp.log(0.1)\n",
    "HARDNESS = 1.0/1.0\n",
    "THRESHOLD = 0.00001\n",
    "MAX_NUM_ITERATIONS = 40\n",
    "logadd = sp.logaddexp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1 Introduction<a name=\"secintro\"></a>\n",
    "==\n",
    "\n",
    "Before we look at what EM is, let’s look at some of the common misconceptions about EM, and clear up what it is not: \n",
    "- a model (it’s a way to optimize one)\n",
    "- a silver bullet (you have to be careful what to use it for)\n",
    "- magic (really, it’s just code. . . )\n",
    "\n",
    "So, after we got that out of the way, let’s see what EM actually is. EM is a form of ***unsupervised learning***. It is not so much an algorithm, but rather a class of algorithms that use a 2-step procedure learning (E step, M step) to train a ***generative model***. Generative models are joint probabilities (written $P(x,y)$) that explain how the data was, well, generated (*discriminative* models, on the other hand, are simply weight vectors that explain the conditional probability $P(x|y)$). EM tries to find the parameters of a model which best explain the observed data. Or: \"If only I knew X, I could estimate Y. If only I knew Y, I could estimate X.\" We will see this sentence a couple of times.\n",
    "\n",
    "Throughout the tutorial, there will be little questions to see whether you are still awake. You can just ignore them and read on, or try to solve them and feel good when you look up the answer at the end.\n",
    "Don’t beat yourself up when you got one wrong, though!\n",
    "\n",
    "Also, do not get scared by complicated-looking formulas. We will translate them into plain English, and the mathematical notation will become just the shorthand it was meant to be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2 Preliminaries<a name=\"secpreliminaries\"></a>\n",
    "==\n",
    "\n",
    "In order to understand EM, we need to look at probabilities and graphical models. If you are not quite sure what they are, fear not, we will explain them in the next sections. If you are already familiar with the concepts, you can skip ahead to [Section 3](#section3).\n",
    "\n",
    "2.1 Probabilities\n",
    "--\n",
    "\n",
    "We will be talking about probabilities a lot, so we will review some basics. If you want to get deeper into it, look at the relevant sections in [Manning/Schütze (2000)](#refMS) and [Russell/Norvig (2003)](#refRN).\n",
    "Since we do NLP, let’s look at the probability of words. Say we have a corpus of 100 sentences, $x$ is \"unsupervised\" and occurs in 20 sentences, and y is \"learning\" and occurs in 50 sentences.\n",
    "Probabilities are basically counts that have been normalized. $P(x)$, or the probability of seeing a sentence with \"unsupervised\", is simply the count of sentences with \"unsupervised\" (20) divided by the number of all sentences in our corpus/text (100).\n",
    "\n",
    "$$\n",
    "P(x=\"unsupervised\") = \\frac{count(sentences\\ with\\ \"unsupervised\")}{number\\ of\\ all\\ sentences} = \\frac{20}{100}\n",
    "= \\frac{1}{5}\n",
    "= 0.2\n",
    "$$\n",
    "\n",
    "P(x,y) is a joint probability, i.e., how likely is it that we see x and y together, that is the words in one sentence (\"unsupervised learning\" or \"learning unsupervised\", or even separated by other words, order does not matter in this case). Say we have 10 sentences that contain both words, then\n",
    "\n",
    "$$\n",
    "P(unsupervised, learning) =\n",
    "\\frac{count(sentences\\ with\\ \"unsupervised\"\\ and\\ \"learning\")}{number\\ of\\ all\\ sentences}\n",
    "= \\frac{10}{100}\n",
    "= \\frac{1}{10}\n",
    "= 0.1\n",
    "$$\n",
    "\n",
    "$P(y|x)$ is a conditional probability, i.e., how likely is it to see $y$ after having seen $x$, or \"learning\" in a sentence that contains \"unsupervised\". We can compute this as\n",
    "\n",
    "$$\n",
    "P(y|x) =\n",
    "\\frac{P(x,y)}{P(x)} =\n",
    "\\frac{\\frac{count(sentences\\ with\\ \"unsupervised\"\\ and\\ \"learning\")}{number\\ of\\ all\\ sentences}}\n",
    "{\\frac{count(sentences\\ with\\ \"unsupervised\")}{number\\ of\\ all\\ sentences}}\n",
    "= \\frac{count(sentences\\ with\\ \"unsupervised\"\\ and\\ \"learning\")}{count(sentences\\ with\\ \"unsupervised\")}\n",
    "= \\frac{10}{20}\n",
    "= \\frac{1}{2}\n",
    "= 0.5\n",
    "$$\n",
    "\n",
    "Note that order ***does*** matter in this case!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q1\"></a>**QUESTION**: What would be the corresponding probability formulation for seeing \"unsupervised\" in a sentence that contains \"learning\"?[$^1$](#answer1)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secprob\"></a>2.1.1 Where Probabilities Come From\n",
    "--\n",
    "\n",
    "You might have wondered where the probabilities we are gonna use come from. That is a very good\n",
    "question, and there are two answers. The good, scientifically satisfying answer is “data\". We can use\n",
    "the statistics collected by statistical institutions (now you know what they are good for), or we can\n",
    "get them by counting and dividing from (preferably large amounts of) data. How likely is ist that it is\n",
    "sunny in LA? Get the weather data of the last 30 years, count the days marked *sunny* and divide by the\n",
    "total number of days. The larger your sample, the more accurate your probability.\n",
    "\n",
    "The other, less satisfying answer is “elaborate guessing\". Some things cannot be measured, or there\n",
    "is simply no data (some things are just too rare. The only exception to this rule is baseball. There are\n",
    "statistics for everything in baseball). What is the probability of being eaten by a tiger in Montana?\n",
    "Probably very small (luckily), but what number is “very small\"? We can make up a number, but it\n",
    "might not be accurate.\n",
    "\n",
    "The good news is that we can use EM to make up numbers and then adjust them until our model\n",
    "explains the data best.\n",
    "\n",
    "2.2 Graphical Models\n",
    "--\n",
    "\n",
    "Graphical Models are a nice way of visualizing probabilistic models, and also to express the dependencies\n",
    "that hold between the individual elements. There are several types of graphical models, but the\n",
    "ones we are interested in (and the ones we mean when we use the term) are ***Bayes Nets*** and ***Hidden\n",
    "Markov Models*** (HMMs). Graphical models consist of two elements, ***nodes*** and ***arcs***.\n",
    "\n",
    "The nodes are ***random variables***. Random variables are events that have some probability, and come\n",
    "in different flavors. If a random variable has exactly two values, like {`on`, `off`} or {`true`, `false`}, they\n",
    "are ***binary*** (in the latter case ***boolean***). If they have a list of values (something like {`red`, `green`, `blue`}\n",
    "or {`chocolate`, `vanilla`, `strawberry`, `pistacchio`}), they are ***discrete***. If they have numbers as values,\n",
    "they are called ***continuous***. \n",
    "The ***probabilities*** associated with each of the values of a random variable have to\n",
    "sum up to $1.0$, i.e., the variable `COLOR` with the values {`red`, `green`, `blue`} could have the respective probabilities {$0.2, 0.5, 0.3$} or\n",
    "{$0.33, 0.33, 0.33$} associated with the values, but not {$0.8, 0.4, 0.7$}.\n",
    "\n",
    "Arcs are the directed links between the random variables, and you can think of them as causal relations\n",
    "(there are other kinds, but it is easiest this way). They denote what influence the parent has on the\n",
    "child node. This influence is expressed by a conditional probability. \n",
    "\n",
    "<img src=\"pics/bn1.png\" width=\"200px\"/>\n",
    "<div align=\"center\">*Figure 1: A simple graphical model with three random variables*</div>\n",
    "\n",
    "E.g., in a network like the one in Figure 1, we can say how likely it is that traffic ($T$) is bad, given that the weather ($W$) is rain. A node\n",
    "$X$ can have several parents, which means that its value is influenced by several factors (traffic could\n",
    "also be influenced by a Lakers game, $G$). If there are no links between two variables, then they are\n",
    "***independent of one another***, i.e., whether the Lakers play or not luckily has no influence on the weather\n",
    "$W$ (the examples in this section are largely influenced by [Russell/Norvig 2003](#refRN)).\n",
    "\n",
    "2.2.1 Bayes Nets\n",
    "--\n",
    "<img src=\"pics/bn2.png\" width=\"450px\"/>\n",
    "<div align=\"center\">*Figure 2: A Bayes Net with three random variables and associated parameters*</div>\n",
    "\n",
    "If we combine several nodes in a network, we call it a ***Bayes Net***. Let’s look at a very simple example\n",
    "(Figure 2), inspired by [Russell/Norvig 2003](#refRN). Say we have three random variables, namely the weather\n",
    "($W$), with values {`sunny`, `rainy`}, traffic ($T$), which can be {`normal`, `bad`, `terrible`} and whether we are\n",
    "late for a meeting ($L$: {`true`, `false`})."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q2\"></a>**QUESTION**: What are the types of the random variables $W$, $T$, and $L$?[$^2$](#answer2)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From pop culture we know that it never rains in southern California, and from our meterological data\n",
    "(see [section 2.1.1](#secprob)) we know that *never* means 5%. So the probabilities for $W$ are $(0.95, 0.05)$. If we\n",
    "talk about the probability of a specific outcome of the variables values, we write $P(W =$ `sunny` $) = 0.95$\n",
    "or shorter $P($ `sunny` $) = 0.95$.\n",
    "\n",
    "If it rains, traffic tends to get worse, and if traffic is bad, we are more likely to be late for our meeting.\n",
    "If it is sunny, the traffic behaves different than when it is rainy, so we have to specify the probability of\n",
    "each value of $T$ for each value of $W$. We do that in a table, where each column is a value for a variable,\n",
    "$T$ and $W$. Notice that the rows with the same value forW have to sum up to $1.0$. You can imagine that\n",
    "each value for weather is a state you are in, and the different values for $T$ are options you can choose\n",
    "from. Some options are more likely than others, but all probability is distributed between them (thus\n",
    "summing to $1.0$). You cannot choose something that is not there.\n",
    "\n",
    "Whether I am late for a meeting ($L$) in turn depends on the state of the traffic ($T$), so we have to\n",
    "specify another table with probabilities for each value of $L$ given each value of $T$. Again you can see\n",
    "that with worse traffic, our chances of being late increase.\n",
    "Using the Bayes Net, we can now compute how likely we are to be late if the weather is bad but\n",
    "traffic is normal, and other interesting things."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secHMMs\"></a>2.2.2 Hidden Markov Models\n",
    "--\n",
    "\n",
    "<img src=\"pics/hmm-example.png\" width=\"300px\"/>\n",
    "<div align=\"center\">*Figure 3: A Hidden Markov Model, resulting from connecting several of the Bayes Net above in a sequence*</div>\n",
    "\n",
    "Things change over time, but they might be connected. Tomorrow’s weather does not just happen, it\n",
    "actually depends on the weather of today. If we want to capture this, we can include another kind of\n",
    "conditional probabilities, namely the ones expressing how a random variable changes over time. This\n",
    "is the Markov part of HMMs. To make things easier, we assume that each state depends only on the\n",
    "previous one, not all previous states. This is one of the so-called ***Markov properties***. In order to \n",
    "make it a hidden Markov model, we assume that the random variable we are actually interested in is\n",
    "unobservable, but related to something we can observe.\n",
    "\n",
    "Using our example from above, we have the following scenario: one year from now, we want to get\n",
    "the sequence of sunny and rainy days that occurred (see Figure 3). We do not remember the weather\n",
    "($W$ is hidden), but we do have our diary, in which we noted for each day whether we were late or not\n",
    "($L$ is our observed variable, and it is dependent on $W$). We just copy the Bayes net from above for\n",
    "each day, and add the new transition probabilities $P(W_t | W_{t-1})$ between each of the Bayes nets. The\n",
    "probability means “how likely is it to be {`rainy`, `sunny`} today if it was {`rainy`, `sunny`} yesterday\".\n",
    "\n",
    "In this case, we do have another hidden variable, $T$, but it is not necessary for HMMs in general. In most applications (such as the tagging we will see later), you only have one hidden layer. The important\n",
    "part is that whether I am late one day does not depend on whether I was late the day before, but on the\n",
    "*weather* on that day. Also, the traffic of today is independent of yesterday’s traffic. This is why there\n",
    "are no arcs between the $T$ and $L$ variables, only the $W$ nodes. This is another Markov property, that the\n",
    "observations (here, $L$) are independent of one another. We guesstimated the probabilities $P(T|W)$ and\n",
    "$P(L|T)$ based on intuitions or data (we will later say, we *initialized* them), and we could now use EM\n",
    "to adjust them to reflect observations, using our diary as data and reconstructing the weather one year\n",
    "ago."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secexample\"></a>3 Example Uses of EM\n",
    "==\n",
    "\n",
    "Let’s start with a practical NLP example of what EM is good for. Say you have a context free grammar\n",
    "(CFG), and would like to attach probabilities to each rule like $P(S\\ \\to\\ NP\\ VP|S)$ to say how likely it is\n",
    "that $S$ goes to $NP\\ VP$ (vs. that $S$ goes to $S\\ CC\\ S$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q3\"></a>**QUESTION**: How do you get the rule probabilites?[$^3$](#answer3)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That’s easy enough, but it presupposes that we do have some suitable treebank. What if we do not\n",
    "have a treebank, only plain text? This case is far more common... No problem, if you just had some\n",
    "rule probabilities, you could generate a treebank by applying the most likely rules to produce the text."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q4\"></a>**QUESTION**: What algorithm would you use to do that?[$^4$](#answer4)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oh wait, dangit! To get those rule probabilities, you’d have to have a treebank to collect them from.\n",
    "But that’s what we wanted in the first place. This is a pretty circular problem! “If only I knew the rule\n",
    "probabilities, I could make a treebank from text. If only I had a treebank, I could compute the rule\n",
    "probabilities...\"\n",
    "\n",
    "EM can help us solve this problem. But how?\n",
    "\n",
    "If you know $k$-means clustering, you already know how EM works! If not, don’t despair. \n",
    "<img src=\"pics/kmeans.png\" width=\"450px\"/>\n",
    "<div align=\"center\">*Figure 4: $k$-means clustering. Modified from [Bishop 2006](refBishop)*</div>\n",
    "\n",
    "Here is an example: We want to divide the green points in Figure 4 into a red and a blue cluster. If\n",
    "only we knew the cluster centroids, we could assign the data points to the closest clusters. If only we\n",
    "knew which clusters the data points belong to, we could compute their centroids... Sounds somehow\n",
    "familiar? Again, we have a circular problem. And here is how we gonna solve it:\n",
    "We start by randomly placing cluster centroids on the graph (a). Then, we assign each data point to\n",
    "a cluster (b)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q5\"></a>**QUESTION**: Based on the image, what do we use to assign the points to clusters?[$^5$](#answer5)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we compute the centers of those new clusters and move the centroids to that position (c).\n",
    "As we can see, we alternate between assigning the points to clusters and computing new centroids.\n",
    "Once the centroids stop moving around, we are done.\n",
    "\n",
    "EM works similarly, and in fact is a “soft\" version of $k$-means. Instead of assigning each point to\n",
    "just one cluster (hard clustering), EM will attach a probability to the membership of a point in each\n",
    "cluster ($P($ `cluster` $|$ `point` $)$). A data point can thus belong to several clusters (though with different\n",
    "probabilities).\n",
    "\n",
    "One the most well-known NLP tasks for EM is labeling unannotated text. We want to find labeling\n",
    "$argmax_t P(t|w)$, i.e., given some words $\\mathbf{w}$, what is the best tag sequence $\\mathbf{t}$? Or: “If only I knew the right\n",
    "tag sequence, I could compute their probabilities... If only I knew the tag probabilities, I could tag\n",
    "the words.\" This is going to be our running example, and I will use $\\mathbf{t}$ for a tag sequence and $\\mathbf{w}$ for a\n",
    "word sequence. The subscript $i$ denotes the position in the respective sequence, so $w_i$ is the $i$th word).\n",
    "\n",
    "As mentioned in the first section, there are several algorithms to do EM (especially the E step), and\n",
    "the one we will be using is the ***Forward-Backward*** (or Baum-Welch) algorithm for labeling (there are other\n",
    "algorithms for other tasks)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secgoal\"></a>4 The Goal\n",
    "==\n",
    "\n",
    "What we want from EM is the ***model parameters*** of our grpahical model (i.e., the conditional probability tables) that best explain the observed data. If we plot how \n",
    "well a specific model configuration explains the data over all possible parameter configurations, we get the\n",
    "graph in Figure 5.\n",
    "<img src=\"pics/localmaxima.png\" width=\"450px\"/>\n",
    "<div align=\"center\">*Figure 5: What we try to optimize with EM*</div>\n",
    "\n",
    "The configuration we want is the one at the highest peak. As you can see, there are several smaller\n",
    "peaks. These are local maxima (they are a bit of a hassle, but we come back to that later).\n",
    "\n",
    "Let’s first look at how we get those parameters. There will be a few formulas involved, but\n",
    "don’t be intimidated, they are easier than they look.\n",
    "\n",
    "Remember that EM models a joint probability distribution $P(x,y)$. We can write that as\n",
    "\n",
    "$$P(x,y) = P(y) \\times P(x|y) = P(x) \\times P(y|x)$$\n",
    "\n",
    "The last two parts are the same because we don’t care about the order of $x$ and $y$. To see why that is\n",
    "so, look at the diagram in Figure 6. \n",
    "<img src=\"pics/venn.png\" width=\"200px\"/>\n",
    "<div align=\"center\">*Figure 6: Venn diagram for $P(x,y)$*</div>\n",
    "\n",
    "$P(x,y)$ is the grey area. To get that, we can just look at $x$ and the\n",
    "part of it that overlaps with $y$. This is $P(x) \\times P(y|x)$. You could call this very $x$-centric. If you don’t like that (and I won't blame you), you can get the same result by looking at $y$ and the part of it that overlaps with $x$. That would be $P(y) \\times P(x|y)$, which is, as you can see, the same as $P(x) \\times P(y|x)$.\n",
    "\n",
    "In our case, however, we don’t want to model $x$s and $y$s, but the probability of a word and tag\n",
    "sequence occurring together, $P(\\mathbf{w}, \\mathbf{t})$. So let’s substitute $x$ and $y$ for $\\mathbf{w}$ and $\\mathbf{t}$ and see what we get.\n",
    "Something like\n",
    "\n",
    "$$P(\\mathbf{w},\\mathbf{t}) = P(\\mathbf{t}) \\times P(\\mathbf{w}|\\mathbf{t}) = P(\\mathbf{w}) \\times P(\\mathbf{t}|\\mathbf{w})$$\n",
    "\n",
    "If we translate this into English, it says that the probability of seeing the word and tag sequence together ($P(\\mathbf{w}, \\mathbf{t})$) is equal to the probability\n",
    "of the tag sequence $P(\\mathbf{t})$ times the probability that we generate the words from that tag sequence\n",
    "($P(\\mathbf{w}|\\mathbf{t})$). It is also equal to the probability of seeing those words ($P(\\mathbf{w})$) times the probability of\n",
    "turning those words into that tag sequence ($P(\\mathbf{t}|\\mathbf{w})$). Both of these formulations are equivalent to $P(\\mathbf{w}, \\mathbf{t})$, and that fact will come in handy.\n",
    "\n",
    "What we want to maximize is the last part, $P(\\mathbf{t}|\\mathbf{w})$ (i.e., what is the best tag sequence for the sentence\n",
    "we see), since we already have the words and want to know the tags. Let’s take the last two parts and move a\n",
    "few things around to get $P(\\mathbf{t}|\\mathbf{w})$ alone."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q6\"></a>**QUESTION**: What operations do we have to do to get $P(\\mathbf{t}|\\mathbf{w})$? (Don't peek ahead...)[$^6$](#answer6)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{P(\\mathbf{t}) \\times P(\\mathbf{w}|\\mathbf{t})}{P(\\mathbf{w})} = P(\\mathbf{t}|\\mathbf{w})$$\n",
    "\n",
    "Ok, that’s something. Take a deep breath and make sure you followed this.\n",
    "Since we observe $\\mathbf{w}$ (the sentence), we can say that $P(\\mathbf{w})$ is $1.0$. In that case we can forget\n",
    "about that denominator (anything divided by $1.0$ stays the same as it was)! So what we want to optimize ($P(\\mathbf{t}|\\mathbf{w})) becomes simply\n",
    "\n",
    "$$P(\\mathbf{t}|\\mathbf{w}) = P(\\mathbf{t}) \\times P(\\mathbf{w}|\\mathbf{t})$$\n",
    "\n",
    "Much cleaner, hmm?\n",
    "\n",
    "Let’s look at that P(\\mathbf{t}). It's a sequence of tags, where each depends on the previous one. The probability of seeing the whole tag sequence $t_1, t_2, \\cdots t_n$ (i.e., P(\\mathbf{t})) is therefore really\n",
    "just the product of seeing each of the tags following another tag. We can write that as\n",
    "\n",
    "$$P(\\mathbf{t}) = P(t_1) \\times P(t_2|t_1) \\times P(t_3|t_2) \\times \\cdots \\times P(t_n|t_{n-1})$$\n",
    "\n",
    "or for short\n",
    "\n",
    "$$P(\\mathbf{t}) = P(t_1) \\times \\prod^n_2 P(t_i|t_{i-1})$$\n",
    "\n",
    "Ok, so we played with the formula, made it cleaner, and dissected $P(\\mathbf{t})$. But how does that help us?\n",
    "Good question...\n",
    "\n",
    "We wanted to find the parameters for our model, and now we have them: $P(\\mathbf{t})$ and $P(\\mathbf{w}|{t})$. So we\n",
    "put it all together, and what we want to optimize is finally the product of our two parameters.\n",
    "\n",
    "$$ P(\\mathbf{t}|\\mathbf{w}) = P(t_1) \\times P(w_1|t_1) \\times \\prod^n_2 P(t_i|t_{i-1}) \\times P(w_i|t_i)$$\n",
    "\n",
    "By the way: parameters come in two flavors: ***free parameters*** (these are the ones we want EM to \n",
    "optimize) and ***fixed parameters*** (we already know these and don’t want EM to change them). You can \n",
    "make all parameters free and let EM handle them, or just some!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q7\"></a>**QUESTION**: What would happen if you fixed all parameters?[$^7$](#answer7)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secimplementation\"></a>5 Implementation Details\n",
    "==\n",
    "\n",
    "In order to compute the parameters, we have to develop a data structure that allows us to manipulate\n",
    "them. We will use a Hidden Markov Model to represent the model and ***lattice-based dynamic\n",
    "programming*** to compute and manipulate the probabilities. The following sections walk through the\n",
    "individual parts and explain them in detail using pseudocode. You can find a Python implementation\n",
    "(not optimized for performance) [further down](#secrunnable).\n",
    "\n",
    "5.1 HMM as Lattice\n",
    "--\n",
    "\n",
    "The idea behind our HMM is something like this: what we can see (the words) was generated by something we\n",
    "cannot see (the tag sequence). This is our ***generative story*** (see Figure 7). Sounds strange? Just wait... \n",
    "<img src=\"pics/hmm.png\" width=\"200px\"/>\n",
    "<div align=\"center\">*Figure 7: The Hidden Markov Model of the POS tagging task*</div>\n",
    "\n",
    "\n",
    "Tags are connected by ***transition probabilities*** $P(t_i|t_{i-1})$, and emit words with ***emission probabilities*** $P(w_i|t_i)$. These look strangely familiar... \n",
    "We could also say “Assume the first word was generated by the first tag, how likely is it the next word was generated by the next tag, and how likely is it that that tag followed the first?\" This translates to a sequence of conditional probabilities that we already know from earlier (aren’t you glad now we went through all those equations?):\n",
    "\n",
    "$$ P(\\mathbf{t}|\\mathbf{w}) = P(t_1) \\times P(w_1|t_1) \\times \\prod^n_2 P(t_i|t_{i-1}) \\times P(w_i|t_i)$$\n",
    "\n",
    "Again, this just means “the most likely tag sequence given a sentence is computed by concatenating\n",
    "the most likely tags that can emit those words\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q8\"></a>**QUESTION**: Why do we need $P(t_i|t_{i-1})$ in there? Why don’t we just take the best tag for each word ($P(t|w)$) and be done?[$^8$](#answer8)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can model HMMs as a matrix/***lattice*** (or automaton) of tags and words, as in Figure 8, by \n",
    "replacing each random variable in the HMM with all possible values and drawing all possible arcs\n",
    "between them. This can be tricky, and translating from a graphical model to a lattice takes some\n",
    "getting used to.\n",
    "<img src=\"pics/lattice.png\" width=\"600px\"/>\n",
    "<div align=\"center\">*Figure 8: The lattice for the sentence \"Mice like cheese\" and two possible tags, N and V*</div>\n",
    "\n",
    "\n",
    "In this case, we want to label “Mice like cheese\", and have an alphabet of only two tags, `N` and\n",
    "`V`. It is important to specify all the tags you want to use. We start from a designated start state and\n",
    "from there choose one of the tags with the respective probability $P(t)$. From each of those possible tag\n",
    "states, we can emit a word with the respective probability $P(w|t)$. Those are the horizontal lines in our\n",
    "lattice. Then, we choose the next tag with some probability $P(t_i|t_{i-1})$. Those are the crossing lines in\n",
    "the lattice.\n",
    "\n",
    "To visualize this, we list all tags as rows. If we have a ***dictionary*** that tells us that some words \n",
    "can only have certain tags, we simply set all other $P(w|t)$ for this word to $0$ (here, we could set\n",
    "($P($ `mice` | `verb` $) = 0$) and omit the arcs. If we don’t have a such a dictionary, we have to assume that all\n",
    "words can be emitted by all tags and let EM figure it out.\n",
    "\n",
    "Ultimately, we want to learn which tags follow one another, say `V` usually comes after `N`, and\n",
    "which words have which tags, e.g., like is most of the times a verb, never a noun. Only expressed as\n",
    "probabilities: what is $P($ `like`$|$ `V` $)$? So how do we assign those values? We reward good parameters, i.e.,\n",
    "transitions that increase $P(sentence)$, and we decrease bad ones. As a first step, instead of just taking\n",
    "whole counts of how often we see a transition, we \"weigh\" them by how likely the resulting sentence\n",
    "was ($P(sentence)$). This is called ***fractional counts***."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q9\"></a>**QUESTION**: What probabilistic event do you think is describe by $P(sentence)$?[$^9$](#answer9)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could try to just generate all possible taggings of a sentence (see Figure 9) and count how often\n",
    "we see `N` following `V` and “cheese\" was tagged as `N`, and then sum them all up to get $P(sentence)$. Here, we have two tags and three words, so we get $2^3=8$ possible paths through the lattice. But there is a problem...\n",
    "<img src=\"pics/naive.png\" width=\"600px\"/>\n",
    "<div align=\"center\">*Figure 9: Brute-force enumeration of all possible tag sequences*</div>\n",
    "\n",
    "Say each word has on avg. $2.5$ tags, and a sentence has $17$ words, like this one. That’s $2.5^{17}$, or\n",
    "$5,820,766$ possible paths through that lattice. For just one sentence! Imagine a sentence with 50\n",
    "words... Clearly, we cannot afford to do that. We will have to do something else. And that something\n",
    "is dynamic programming, see [Section 5.2](#secdp).\n",
    "\n",
    "An aside: since we multiply a lot of transitions to get through the lattice, the numbers can\n",
    "quickly become very small. To deal with this, you can use the logarithm of the probabilities.\n",
    "The smaller a number gets, the larger its negative log will be. Since the range of probabilities is\n",
    "between 0 and 1:0, this corresponds an interval between negative infinity and 0 in log world.\n",
    "If you do use logarithms, all multiplications shown here become additions (which is slightly faster),\n",
    "and all additions have to be log-additions, a special computation that unfortunately is relatively slow,\n",
    "and works like this (adapted from [Manning/Schütze 2000](#refMS)):\n",
    "\n",
    "```\n",
    "m = min(x,y)\n",
    "big = 10^30\n",
    "if y-x > log(big):\n",
    "    return (y)\n",
    "else if x- y > log(big):\n",
    "    return (x)\n",
    "else \n",
    "    return (m + log(exp(x-m)+exp(y-m)))\n",
    "```\n",
    "\n",
    "**NB:** In the remainder of the notebook, I use “normal” probabilities (not logs), to de-clutter the notation, but in the [implementation](#secrunnable), we will use logs, and the `scipy` function `logaddexp()`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5.2 Dynamic Programming\n",
    "--\n",
    "\n",
    "Back to our question: we wanted to know what $P($ `like` | `V` $)$ is. This now turns into the question of how\n",
    "likely it is that we end up at the node that has $P($ `like` | `V` $)$ as outgoing transition (in our example node\n",
    "(2,3)). And once we took that transition, what is the probability from the node we reach (i.e., (2,4))\n",
    "to the end?\n",
    "\n",
    "Since we have modeled the HMM as a lattice, we can use dynamic programming techniques. This\n",
    "allows us to compute how likely it is to arrive at each node (with the Forward algorithm), and to get to\n",
    "the end from there (with the Backward algorithm).\n",
    "\n",
    "So this is where we use the Forward-Backward algorithm! It consists of two parts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q10\"></a>**QUESTION**: Take a wild guess what the two parts are called.[$^{10}$](#answer10)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use Forward-Backward in order to efficiently compute for each sentence how often we see each\n",
    "transition and what the probability of that sentence is. We need both for the fractional counts. Forward-\n",
    "Backward is the E step in our EM implementation: we compute the expected counts given the current\n",
    "model parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q11\"></a>**QUESTION**: What was the E-step in clustering?[$^{11}$](#answer11)</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you model the lattice, it is a good idea to use a matrix or some such data structure in your\n",
    "implementation, so you can access the nodes directly. In the following, we’ll use a matrix with two\n",
    "columns for each word, to make everything easier to see. In practice, people use one column per word and later\n",
    "tease out the fractional counts for each type of parameter.\n",
    "\n",
    "**The Forward Algorithm **\n",
    "<img src=\"pics/alphalattice.png\" width=\"600px\"/>\n",
    "<div align=\"center\">*Figure 10: Computing the alphas in the Forward pass*</div>\n",
    "\n",
    "In the forward pass, we compute a new lattice with the same dimensions as the original one, which forward pass\n",
    "contains for each node the sum of all possible paths that lead up to there (see Figure 10). These values\n",
    "are also called alphas. $\\alpha[i, j]$ denotes the probability of all paths up to node $(i, j)$ (i.e., assigning tag $i$ to word $j$). \n",
    "\n",
    "$\\alpha[START]$ is always $1.0$. Each subsequent $\\alpha$ is just the sum of all transitions arriving there, each multiplied by the $\\alpha$ of the node where it (the transition) originated.\n",
    "\n",
    "$\\alpha[END]$ is the sum of all paths through the lattice, which is equal to $P(sentence)$. $P(data)$ is the\n",
    "sum of all $P(sentence)$ in the data. In each iteration, just add up all the $\\alpha[END]$ of the sentences.\n",
    "Remember, $P(data)$ has to increase with each iteration, or there is something wrong! EM guarantees\n",
    "that the likelihood of the data increases at each iteration over the data. Outputting $P(data)$ is thus a\n",
    "good way of debugging your code: if it does not increase, something went wrong..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def forward(alpha, tags, emissions, transitions, example):\n",
    "    \"\"\"\n",
    "    forward pass\n",
    "    @param alpha: initialized matrix. alpha[i,j] = the probability of ending up with tag i at word j\n",
    "    @param tags: tag set\n",
    "    @param emissions: emission parameters, nested dictionary\n",
    "    @param transitions: transition parameters, nested dictionary\n",
    "    @param example: list of words\n",
    "    @return: alpha\n",
    "    \"\"\"    \n",
    "    for i, tag in enumerate(tags):\n",
    "        if tag in start:\n",
    "            alpha[i][0] = start[ tag ]\n",
    "            if example[0] in emissions[tag]:\n",
    "                alpha[i][1] = alpha[i][0] + emissions[tag][ example[0] ]\n",
    "\n",
    "    for j in range(2, N, 2):\n",
    "        for i, tag1 in enumerate(tags):\n",
    "            for k, tag2 in enumerate(tags):\n",
    "                if tag1 in transitions[tag2]:\n",
    "                    # alpha[i,j] += P(t1|t2) * alpha[k, j-1]\n",
    "                    if alpha[i, j] == NINF:\n",
    "                        alpha[i,j] = transitions[tag2][tag1] + alpha[k,j-1]\n",
    "                    else:\n",
    "                        alpha[i,j] = logadd(alpha[i][j],\n",
    "                                            transitions[tag2][tag1]+alpha[k,j-1])\n",
    "            if example[int(j/2)] in emissions[tag1]:\n",
    "                # alpha[i, j+1] = alpha[i, j] * P(word_j|tag1)\n",
    "                alpha[i,j+1] = alpha[i,j] + emissions[tag1][ example[int(j/2)] ]\n",
    "\n",
    "    return alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The Backward Algorithm**\n",
    "<img src=\"pics/betalattice.png\" width=\"600px\"/>\n",
    "<div align=\"center\">*Figure 11: Computing the betas in the Backward pass*</div>\n",
    "\n",
    "The backward pass is almost the same as the forward pass, just backwards (note how the direction of backward\n",
    "the arrows is reversed in Figure 11). Again, we compute a new lattice, which contains for each node pass\n",
    "the sum of all possible paths that lead from that node to the end. These values are called ***betas***. $\\beta[i, j]$\n",
    "denotes the summed probability of all paths from node $(i, j)$ to the end. This time, however, we start at\n",
    "the end. $\\beta[END]$ is always $1.0$.\n",
    "\n",
    "A useful property for debugging is the fact that $\\beta[START] = \\alpha[END]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q12\"></a>**QUESTION**: Why are the two qunatities the same?[$^{12}$](#answer12)</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def backward(beta, tags, emissions, transitions, example):\n",
    "    \"\"\"\n",
    "    backward pass\n",
    "    @param alpha: initialized matrix\n",
    "    @param tags: tag set\n",
    "    @param emissions: emission parameters, nested dictionary\n",
    "    @param transitions: transition parameters, nested dictionary\n",
    "    @param example: list of words\n",
    "    @return: completed beta\n",
    "    \"\"\"    \n",
    "    # initialize beta starting from the back\n",
    "    for i, tag in enumerate(tags):\n",
    "        beta[i][-1] = sp.log(1.0)\n",
    "        if example[-1] in emissions[tag]:\n",
    "            beta[i][-2] = beta[i][-1] + emissions[tag][ example[-1] ]\n",
    "\n",
    "    for j in range(N-3,0,-2):\n",
    "        for i, tag1 in enumerate(tags):\n",
    "            for k, tag2 in enumerate(tags):\n",
    "                if tag2 in transitions[tag1]:\n",
    "                    if beta[i,j] == NINF:\n",
    "                        beta[i,j] = transitions[tag1][tag2] + beta[k,j+1]\n",
    "                    else:\n",
    "                        beta[i,j] = logadd(beta[i,j], transitions[tag1][tag2] + beta[k,j+1])\n",
    "            if example[int(j/2)] in emissions[tag1]:\n",
    "                beta[i, j-1] = beta[i, j] + emissions[tag1][example[int(j/2)]]\n",
    "\n",
    "    return beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Collecting fractional counts**\n",
    "<img src=\"pics/fraccount.png\" width=\"300px\"/>\n",
    "<div align=\"center\">*Figure 12: Collecting fractional counts for $P($ `like` | `V` $)$*</div>\n",
    "\n",
    "Once we have the alphas and betas, it is easy to compute for each transition how much it contributes\n",
    "to $P(sentence)$. So, once more, with conviction: how good is $P($ `like` | `V` $)$? Remember, we have to know\n",
    "the likelihood of all possible paths arriving at node $(2,3)$, and the probability – once we have taken the\n",
    "transition – from node $(2,4)$ to the end. See Figure 12.\n",
    "\n",
    "We used the forward algorithm to get the probability of arriving at node $(2,3)$, and the backward\n",
    "algorithm to compute how likely it is from node $(2,4)$ to the end. We divide that by the likelihood of\n",
    "the sentence ($= \\alpha[END]$), et voila!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def collect_counts(alpha, beta, tags, start, start_counts, emissions, emission_counts, transitions, transition_counts, example, example_likelihood):\n",
    "\n",
    "    for i, tag1 in enumerate(tags):\n",
    "        if tag1 in start:\n",
    "            # if no smoothing, start from scratch, otherwise just add up\n",
    "            if start_counts[tag1] == sp.log(0.0):\n",
    "                start_counts[tag1] = start[tag1] + beta[i, 0] - example_likelihood\n",
    "            else:\n",
    "                start_counts[tag1]=logadd(start_counts[tag1],\n",
    "                                          start[tag1]+beta[i,0]-example_likelihood)\n",
    "\n",
    "        for j, word in enumerate(example):\n",
    "            if word in emissions[tag1]:\n",
    "                if emission_counts[tag1][word] == sp.log(0.0):\n",
    "                    emission_counts[tag1][word] = alpha[i, j*2] + emissions[tag1][word] + beta[i, (j*2)+1] - example_likelihood\n",
    "                else:\n",
    "                    emission_counts[tag1][word]=logadd(emission_counts[tag1][word], alpha[i, j*2] + emissions[tag1][word] + beta[i, (j*2)+1] - example_likelihood)\n",
    "\n",
    "            for k, tag2 in enumerate(tags):\n",
    "                if j < len(example)-1 and tag2 in transitions[tag1]:\n",
    "                    if transition_counts[tag1][tag2] == sp.log(0.0):\n",
    "                        transition_counts[tag1][tag2] = alpha[i, j*2+1] + transitions[tag1][tag2] + beta[k, (j+1)*2] - example_likelihood\n",
    "                    else:\n",
    "                        transition_counts[tag1][tag2]=logadd(transition_counts[tag1][tag2], alpha[i, j*2+1] + transitions[tag1][tag2] + beta[k, (j+1)*2] - example_likelihood)\n",
    "   \n",
    "    return start_counts, transition_counts, emission_counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The M Step**\n",
    "\n",
    "Computing alphas and betas and collecting the fractional counts for all free parameter transitions over\n",
    "all examples is the ***E step***. This, as the name suggests, is one half of Forward-Backward EM, and in \n",
    "this case the bigger half. The ***M step*** is comparatively trivial: after having gone through all the data, we \n",
    "just normalize our fractional counts to get probabilities back (remember, probabilities are just normalized\n",
    "counts)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"blue\"><a name=\"q13\"></a>**QUESTION**: What do you need to compute conditional probabilities from counts?[$^{13}$](#answer13)</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normalize_counts(start_counts, start, emission_counts, emissions, transition_counts, transitions):\n",
    "    start_count_total = reduce(logadd, list(start_counts.values()))\n",
    "\n",
    "    print(\"\\nNew START probabilities:\")\n",
    "    print(('-' * 20))\n",
    "    for tag in start:\n",
    "        start[tag] = start_counts[tag] - start_count_total\n",
    "        print(('P(%s) = %.2f' % (tag, sp.exp(start[tag]))))\n",
    "\n",
    "    print(\"\\nNew EMISSION probabilities:\")\n",
    "    print(('-' * 20))\n",
    "    for tag, words in emissions.items():\n",
    "        emission_tag_total = reduce(logadd, list(emission_counts[tag].values()))\n",
    "        for word in words:\n",
    "            emissions[tag][word] = emission_counts[tag][word] - emission_tag_total\n",
    "            print(('P(%s|%s) = %.2f' % (word,tag, sp.exp(emissions[tag][word]))))\n",
    "\n",
    "    print(\"\\nNew TRANSITION probabilities:\")\n",
    "    print(('-' * 20))\n",
    "    for tag, tag_successors in transitions.items():\n",
    "        transition_tag_total = reduce(logadd, list(transition_counts[tag].values()))\n",
    "        for tag_successor in tag_successors:\n",
    "            transitions[tag][tag_successor] = transition_counts[tag][tag_successor] - transition_tag_total\n",
    "            print(('P(%s|%s) = %.2f' % (tag_successor,\n",
    "                                        tag,\n",
    "                                        sp.exp(transitions[tag][tag_successor]))))\n",
    "\n",
    "    print()\n",
    "    return start, emissions, transitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5.3 Putting It All Together\n",
    "--\n",
    "\n",
    "If we put the E step and the M step together, we end up with the Forward-Backward EM algorithm!\n",
    "\n",
    "The convergence criterion in this case is how much the data likelihood P(data) has improved since\n",
    "the last iteration. Once it starts to flatten out, we can assume that we reached a maximum on our curve\n",
    "and stop. It is also customary to set a maximum number of iterations (50) and stop even if EM has not\n",
    "converged, to avoid overfitting.\n",
    "\n",
    "And that’s it. No magic, no silver bullets, just counting and normalizing. EM will adjust all\n",
    "the free parameters to get the maximum data likelihood, and you can then use those probabilities to\n",
    "label data using the Viterbi algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secrunnable\"></a>5.4 Example Run\n",
    "--\n",
    "\n",
    "If you want to check whether your implementation is working, here is a little worked toy example. Let’s take\n",
    "the tags `N` and `V`, and the words *can* and *I* and initialize our model with the following numbers\n",
    "(this example is based on one used by Kevin Knight.)\n",
    "\n",
    "Transitions:\n",
    "\n",
    "$P($ `V` | `V` $) = 0.6$\n",
    "\n",
    "$P($ `N` | `V`$) = 0.4$\n",
    "\n",
    "$P($ `V` | `N`$) = 0.9$\n",
    "\n",
    "$P($ `N` | `N`$) = 0.1$\n",
    "\n",
    "$P($ `V` | `START` $) = P($ `V` $) = 0.6$\n",
    "\n",
    "$P($ `N` | `START` $) = P($ `N ` $) = 0.4$\n",
    "\n",
    "\n",
    "Emissions:\n",
    "\n",
    "$P($ `can` | `V` $) = 0.5$\n",
    "\n",
    "$P($ `I` | `V` $) = 0.5$\n",
    "\n",
    "$P($ `can` | `N` $) = 0.5$\n",
    "\n",
    "$P($ `I` | `N` $) = 0.5$\n",
    "\n",
    "\n",
    "As our only observed instance we use the sentences \"*can I can*\" and \"*I can can*\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "start = defaultdict(lambda: LINF)\n",
    "# P(x|y) = emissions[y][x] - this allows to sum over all values given y\n",
    "emissions = defaultdict(lambda: defaultdict(lambda: LINF))\n",
    "transitions = defaultdict(lambda: defaultdict(lambda: LINF))\n",
    "\n",
    "# initialize\n",
    "start['V'] = sp.log(0.6)\n",
    "start['N'] = sp.log(0.4)\n",
    "\n",
    "emissions['V']['I'] = sp.log(0.5)\n",
    "emissions['V']['can'] = sp.log(0.5)\n",
    "emissions['N']['I'] = sp.log(0.5)\n",
    "emissions['N']['can'] = sp.log(0.5)\n",
    "\n",
    "transitions['V']['V'] = sp.log(0.6)\n",
    "transitions['V']['N'] = sp.log(0.4)\n",
    "transitions['N']['N'] = sp.log(0.1)\n",
    "transitions['N']['V'] = sp.log(0.9)\n",
    "\n",
    "tags = ['V', 'N']\n",
    "\n",
    "examples = [\n",
    "    ['can', 'I', 'can'],\n",
    "    [ 'I', 'can', 'can']\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual code is two simple loops (one over the data, one until convergence) that call our functions defined above. You can run it and observe the output. Your $\\alpha$s and $\\beta$s for \"*can I can*\" after the first iteration should look like they do in Figure 13, where the $\\alpha$ values for each node are shown on a white\n",
    "background above the nodes, and $\\beta$ values with a grey background below the nodes. The arrows point\n",
    "forward, but you can just ignore them.\n",
    "\n",
    "<img src=\"pics/example_run.png\" width=\"600px\"/>\n",
    "<div align=\"center\">*Figure 13: $\\alpha$ (unshaded) and $\\beta$ (shaded boxes) values for example run*</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "====================\n",
      "iteration 1\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.6    0.3    0.36   0.18   0.171  0.085]\n",
      " [ 0.4    0.2    0.14   0.07   0.079  0.039]]\n",
      "BETA:\n",
      "[[ 0.125  0.25   0.25   0.5    0.5    1.   ]\n",
      " [ 0.125  0.25   0.25   0.5    0.5    1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.12499999999999997)\n",
      "0.00124sec\n",
      "ALPHA:\n",
      "[[ 0.6    0.3    0.36   0.18   0.171  0.085]\n",
      " [ 0.4    0.2    0.14   0.07   0.079  0.039]]\n",
      "BETA:\n",
      "[[ 0.125  0.25   0.25   0.5    0.5    1.   ]\n",
      " [ 0.125  0.25   0.25   0.5    0.5    1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.12499999999999997)\n",
      "0.00112sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.59\n",
      "P(N) = 0.41\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.34\n",
      "P(can|V) = 0.66\n",
      "P(I|N) = 0.36\n",
      "P(can|N) = 0.64\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.59\n",
      "P(N|V) = 0.41\n",
      "P(V|N) = 0.85\n",
      "P(N|N) = 0.15\n",
      "\n",
      "CHANGE:\n",
      "old: 0.00000 (-inf)\n",
      "new: 0.01562 (-4.15888308336)\n",
      "change:  0.015625\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 2\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.591  0.392  0.456  0.154  0.151  0.1  ]\n",
      " [ 0.409  0.264  0.199  0.071  0.073  0.047]]\n",
      "BETA:\n",
      "[[ 0.15   0.227  0.221  0.655  0.663  1.   ]\n",
      " [ 0.144  0.223  0.235  0.66   0.644  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14752815830794644)\n",
      "0.00113sec\n",
      "ALPHA:\n",
      "[[ 0.591  0.199  0.242  0.16   0.151  0.1  ]\n",
      " [ 0.409  0.146  0.103  0.066  0.075  0.048]]\n",
      "BETA:\n",
      "[[ 0.145  0.43   0.434  0.655  0.663  1.   ]\n",
      " [ 0.154  0.433  0.425  0.66   0.644  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14879225381021891)\n",
      "0.00114sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.58\n",
      "P(N) = 0.42\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.33\n",
      "P(can|V) = 0.67\n",
      "P(I|N) = 0.37\n",
      "P(can|N) = 0.63\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.59\n",
      "P(N|V) = 0.41\n",
      "P(V|N) = 0.81\n",
      "P(N|N) = 0.19\n",
      "\n",
      "CHANGE:\n",
      "old: 0.01562 (-4.15888308336)\n",
      "new: 0.02195 (-3.8189404333)\n",
      "change:  0.00632604717511\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 3\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.581  0.389  0.442  0.146  0.149  0.1  ]\n",
      " [ 0.419  0.265  0.212  0.078  0.075  0.048]]\n",
      "BETA:\n",
      "[[ 0.152  0.228  0.216  0.654  0.67   1.   ]\n",
      " [ 0.14   0.221  0.244  0.663  0.631  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14713268198536739)\n",
      "0.00113sec\n",
      "ALPHA:\n",
      "[[ 0.581  0.192  0.237  0.159  0.149  0.1  ]\n",
      " [ 0.419  0.154  0.109  0.069  0.079  0.05 ]]\n",
      "BETA:\n",
      "[[ 0.142  0.43   0.438  0.654  0.67   1.   ]\n",
      " [ 0.16   0.434  0.418  0.663  0.631  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14953405162679509)\n",
      "0.00120sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.57\n",
      "P(N) = 0.43\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.32\n",
      "P(can|V) = 0.68\n",
      "P(I|N) = 0.38\n",
      "P(can|N) = 0.62\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.58\n",
      "P(N|V) = 0.42\n",
      "P(V|N) = 0.77\n",
      "P(N|N) = 0.23\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02195 (-3.8189404333)\n",
      "new: 0.02200 (-3.81665164277)\n",
      "change:  5.02988888786e-05\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 4\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.57   0.386  0.431  0.139  0.146  0.099]\n",
      " [ 0.43   0.266  0.222  0.085  0.077  0.048]]\n",
      "BETA:\n",
      "[[ 0.155  0.228  0.21   0.653  0.678  1.   ]\n",
      " [ 0.136  0.22   0.254  0.665  0.618  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14678364688922718)\n",
      "0.00123sec\n",
      "ALPHA:\n",
      "[[ 0.57   0.183  0.234  0.158  0.147  0.1  ]\n",
      " [ 0.43   0.164  0.114  0.07   0.082  0.051]]\n",
      "BETA:\n",
      "[[ 0.138  0.43   0.443  0.653  0.678  1.   ]\n",
      " [ 0.166  0.436  0.411  0.665  0.618  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.15032583692873194)\n",
      "0.00126sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.56\n",
      "P(N) = 0.44\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.31\n",
      "P(can|V) = 0.69\n",
      "P(I|N) = 0.39\n",
      "P(can|N) = 0.61\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.58\n",
      "P(N|V) = 0.42\n",
      "P(V|N) = 0.74\n",
      "P(N|N) = 0.26\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02200 (-3.81665164277)\n",
      "new: 0.02207 (-3.81374566083)\n",
      "change:  6.40285020858e-05\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 5\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.557  0.383  0.422  0.132  0.144  0.099]\n",
      " [ 0.443  0.268  0.229  0.09   0.079  0.048]]\n",
      "BETA:\n",
      "[[ 0.157  0.229  0.204  0.653  0.687  1.   ]\n",
      " [ 0.133  0.219  0.263  0.666  0.606  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14639482670051945)\n",
      "0.00111sec\n",
      "ALPHA:\n",
      "[[ 0.557  0.174  0.231  0.159  0.145  0.1  ]\n",
      " [ 0.443  0.175  0.118  0.071  0.085  0.051]]\n",
      "BETA:\n",
      "[[ 0.134  0.43   0.449  0.653  0.687  1.   ]\n",
      " [ 0.172  0.437  0.403  0.666  0.606  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.15126587133971076)\n",
      "0.00110sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.54\n",
      "P(N) = 0.46\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.30\n",
      "P(can|V) = 0.70\n",
      "P(I|N) = 0.41\n",
      "P(can|N) = 0.59\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.58\n",
      "P(N|V) = 0.42\n",
      "P(V|N) = 0.72\n",
      "P(N|N) = 0.28\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02207 (-3.81374566083)\n",
      "new: 0.02214 (-3.81016426792)\n",
      "change:  7.91664544055e-05\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 6\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.542  0.378  0.415  0.126  0.142  0.099]\n",
      " [ 0.458  0.271  0.235  0.095  0.08   0.047]]\n",
      "BETA:\n",
      "[[ 0.16   0.229  0.198  0.653  0.697  1.   ]\n",
      " [ 0.13   0.219  0.272  0.668  0.593  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14588352130771209)\n",
      "0.00111sec\n",
      "ALPHA:\n",
      "[[ 0.542  0.164  0.23   0.16   0.145  0.101]\n",
      " [ 0.458  0.186  0.121  0.072  0.087  0.052]]\n",
      "BETA:\n",
      "[[ 0.13   0.43   0.455  0.653  0.697  1.   ]\n",
      " [ 0.179  0.439  0.396  0.668  0.593  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.15246971854455341)\n",
      "0.00167sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.53\n",
      "P(N) = 0.47\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.29\n",
      "P(can|V) = 0.71\n",
      "P(I|N) = 0.42\n",
      "P(can|N) = 0.58\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.58\n",
      "P(N|V) = 0.42\n",
      "P(V|N) = 0.70\n",
      "P(N|N) = 0.30\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02214 (-3.81016426792)\n",
      "new: 0.02224 (-3.80573604452)\n",
      "change:  9.82784135952e-05\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 7\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.526  0.372  0.409  0.12   0.14   0.099]\n",
      " [ 0.474  0.275  0.238  0.1    0.08   0.046]]\n",
      "BETA:\n",
      "[[ 0.162  0.229  0.192  0.654  0.707  1.   ]\n",
      " [ 0.127  0.218  0.281  0.669  0.58   1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14517934910133479)\n",
      "0.00197sec\n",
      "ALPHA:\n",
      "[[ 0.526  0.154  0.229  0.162  0.145  0.102]\n",
      " [ 0.474  0.199  0.124  0.072  0.089  0.052]]\n",
      "BETA:\n",
      "[[ 0.126  0.431  0.462  0.654  0.707  1.   ]\n",
      " [ 0.185  0.44   0.388  0.669  0.58   1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.15405545844788776)\n",
      "0.00115sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.51\n",
      "P(N) = 0.49\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.28\n",
      "P(can|V) = 0.72\n",
      "P(I|N) = 0.43\n",
      "P(can|N) = 0.57\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.58\n",
      "P(N|V) = 0.42\n",
      "P(V|N) = 0.69\n",
      "P(N|N) = 0.31\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02224 (-3.80573604452)\n",
      "new: 0.02237 (-3.80022803213)\n",
      "change:  0.000122851748897\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 8\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.508  0.364  0.405  0.114  0.138  0.099]\n",
      " [ 0.492  0.279  0.239  0.103  0.08   0.045]]\n",
      "BETA:\n",
      "[[ 0.164  0.229  0.185  0.655  0.717  1.   ]\n",
      " [ 0.124  0.218  0.29   0.671  0.568  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.1442267198147093)\n",
      "0.00129sec\n",
      "ALPHA:\n",
      "[[ 0.508  0.144  0.23   0.165  0.146  0.104]\n",
      " [ 0.492  0.213  0.126  0.072  0.091  0.052]]\n",
      "BETA:\n",
      "[[ 0.122  0.433  0.47   0.655  0.717  1.   ]\n",
      " [ 0.191  0.442  0.381  0.671  0.568  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.15613385857083142)\n",
      "0.00111sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.49\n",
      "P(N) = 0.51\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.27\n",
      "P(can|V) = 0.73\n",
      "P(I|N) = 0.44\n",
      "P(can|N) = 0.56\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.59\n",
      "P(N|V) = 0.41\n",
      "P(V|N) = 0.68\n",
      "P(N|N) = 0.32\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02237 (-3.80022803213)\n",
      "new: 0.02252 (-3.79341034629)\n",
      "change:  0.000153003090713\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 9\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.489  0.356  0.401  0.109  0.136  0.099]\n",
      " [ 0.511  0.284  0.238  0.106  0.079  0.044]]\n",
      "BETA:\n",
      "[[ 0.166  0.229  0.179  0.656  0.728  1.   ]\n",
      " [ 0.121  0.217  0.299  0.673  0.555  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.14298593453890276)\n",
      "0.00117sec\n",
      "ALPHA:\n",
      "[[ 0.489  0.133  0.232  0.169  0.147  0.107]\n",
      " [ 0.511  0.227  0.128  0.071  0.093  0.052]]\n",
      "BETA:\n",
      "[[ 0.118  0.435  0.478  0.656  0.728  1.   ]\n",
      " [ 0.198  0.444  0.373  0.673  0.555  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.15879933623596237)\n",
      "0.00110sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.47\n",
      "P(N) = 0.53\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.26\n",
      "P(can|V) = 0.74\n",
      "P(I|N) = 0.46\n",
      "P(can|N) = 0.54\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.59\n",
      "P(N|V) = 0.41\n",
      "P(V|N) = 0.67\n",
      "P(N|N) = 0.33\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02252 (-3.79341034629)\n",
      "new: 0.02271 (-3.78512292348)\n",
      "change:  0.000187397222172\n",
      "**************************************************\n",
      "\n",
      "====================\n",
      "iteration 10\n",
      "====================\n",
      "ALPHA:\n",
      "[[ 0.47   0.347  0.399  0.104  0.134  0.099]\n",
      " [ 0.53   0.288  0.236  0.108  0.078  0.042]]\n",
      "BETA:\n",
      "[[ 0.169  0.228  0.172  0.658  0.739  1.   ]\n",
      " [ 0.117  0.216  0.309  0.675  0.542  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.1414365590898522)\n",
      "0.00114sec\n",
      "ALPHA:\n",
      "[[ 0.47   0.123  0.236  0.174  0.15   0.111]\n",
      " [ 0.53   0.243  0.129  0.07   0.094  0.051]]\n",
      "BETA:\n",
      "[[ 0.114  0.437  0.487  0.658  0.739  1.   ]\n",
      " [ 0.205  0.447  0.366  0.675  0.542  1.   ]]\n",
      "\n",
      "('example likelihood: ', 0.16211732137101215)\n",
      "0.00110sec\n",
      "\n",
      "\n",
      "New START probabilities:\n",
      "--------------------\n",
      "P(V) = 0.45\n",
      "P(N) = 0.55\n",
      "\n",
      "New EMISSION probabilities:\n",
      "--------------------\n",
      "P(I|V) = 0.25\n",
      "P(can|V) = 0.75\n",
      "P(I|N) = 0.47\n",
      "P(can|N) = 0.53\n",
      "\n",
      "New TRANSITION probabilities:\n",
      "--------------------\n",
      "P(V|V) = 0.59\n",
      "P(N|V) = 0.41\n",
      "P(V|N) = 0.67\n",
      "P(N|N) = 0.33\n",
      "\n",
      "CHANGE:\n",
      "old: 0.02271 (-3.78512292348)\n",
      "new: 0.02293 (-3.77533900796)\n",
      "change:  0.000223244607723\n",
      "**************************************************\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dirkhovy/anaconda/envs/py3.4/lib/python3.4/site-packages/numpy/lib/scimath.py:262: RuntimeWarning: divide by zero encountered in log\n",
      "  return nx.log(x)\n"
     ]
    }
   ],
   "source": [
    "converged = False\n",
    "old_data_likelihood = sp.log(0.0)\n",
    "iteration = 1\n",
    "\n",
    "# repeat until convergence\n",
    "# (i.e., difference of data-log likelihood < THRESHOLD or MAX_NUM_ITERATIONS)\n",
    "while not converged:\n",
    "    print(('=' * 20))\n",
    "    print((\"iteration %s\" % (iteration)))\n",
    "    print(('=' * 20))\n",
    "    data_likelihood = sp.log(1.0)\n",
    "\n",
    "    # clear counts\n",
    "    start_counts = defaultdict(lambda: SMOOTHING)\n",
    "    emission_counts = defaultdict(lambda: defaultdict(lambda: SMOOTHING))\n",
    "    transition_counts = defaultdict(lambda: defaultdict(lambda: SMOOTHING))\n",
    "\n",
    "    ##########\n",
    "    # E-Step #\n",
    "    ##########\n",
    "    for example in examples:\n",
    "        starttime = time.time()\n",
    "        N = len(example) * 2\n",
    "        M = len(tags)\n",
    "\n",
    "        ################\n",
    "        # forward pass #\n",
    "        ################\n",
    "        alpha = forward(sp.ones((M,N)) * NINF, tags, emissions, transitions, example)\n",
    "        print(\"ALPHA:\")\n",
    "        print((sp.exp(alpha)))\n",
    "        \n",
    "        #################\n",
    "        # backward pass #\n",
    "        #################\n",
    "        beta = backward(sp.ones((M,N)) * NINF, tags, emissions, transitions, example)\n",
    "        print(\"BETA:\")\n",
    "        print((sp.exp(beta)))\n",
    "        print()\n",
    "\n",
    "        # sum of all paths through example\n",
    "        example_likelihood = reduce(logadd, alpha[:, -1])\n",
    "        print((\"example likelihood: \", sp.exp(example_likelihood)))\n",
    "\n",
    "        #############################\n",
    "        # collect fractional counts #\n",
    "        #############################\n",
    "        start_counts, transition_counts, emission_counts = collect_counts(alpha, beta, tags, start, start_counts, emissions, emission_counts, transitions, transition_counts, example, example_likelihood)\n",
    "\n",
    "        # add example likelihood to data likelihood\n",
    "        data_likelihood += example_likelihood\n",
    "\n",
    "        print((\"%.5fsec\" % (time.time()-starttime)))\n",
    "    print()\n",
    "\n",
    "    ############################\n",
    "    # M-Step: normalize counts #\n",
    "    ############################\n",
    "    start, emissions, transitions = normalize_counts(start_counts, start, emission_counts, emissions, transition_counts, transitions)\n",
    "\n",
    "    change = -(sp.exp(old_data_likelihood) - sp.exp(data_likelihood))\n",
    "    assert change > 0.0\n",
    "    print(\"CHANGE:\")\n",
    "    print(('old: %.5f (%s)' % (sp.exp(old_data_likelihood), old_data_likelihood)))\n",
    "    print(('new: %.5f (%s)' % (sp.exp(data_likelihood), data_likelihood)))\n",
    "    print(\"change: \", change)\n",
    "    print('*' * 50)\n",
    "    print()\n",
    "    if change <= THRESHOLD or iteration > MAX_NUM_ITERATIONS:\n",
    "        converged = True\n",
    "        print('Converged!')\n",
    "\n",
    "    old_data_likelihood = data_likelihood\n",
    "    if iteration == 10:\n",
    "        break\n",
    "        \n",
    "    iteration += 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secfinally\"></a>6 Finally...\n",
    "==\n",
    "\n",
    "CONGRATULATIONS! You have made it through. You now know EM and can go and try it out...\n",
    "You do not have to be constrained to HMMs with hidden inputs, like we used here. Maybe\n",
    "you do not have sequential data. In that case the HMM becomes a Bayesian network, which is like an\n",
    "HMM without the transitions ([Hovy et al. (2013)](#refMACE) is an example for training a Bayes net with EM to\n",
    "find out which annotators are more reliable than others, and what the most likely answer is).\n",
    "Maybe you know both the inputs and outputs, but not the hidden variable $X$ that connects them .\n",
    "In that case you just want to know the conditional probabilities $P(X|input)$ and $P(output|X)$ (in our\n",
    "diary example above, you might also have the weather reports from that time and want to compute the\n",
    "traffic probabilities). EM can help you in those cases, too.\n",
    "\n",
    "If you want to apply EM to other problems, always start with the inputs and outputs. What are\n",
    "the observed variables? Are the inputs hidden, or can they be observed, too? Is it a sequential model\n",
    "or a network? Write down the joint probability $P(x,y)$ and see how you can factor it. Drawing a\n",
    "graphical model helps.\n",
    "\n",
    "Think about how you can express the parameters as conditional probabilities. Do you need additional\n",
    "(hidden) variables? What is your E step, what is your M step, what is the data you iterate\n",
    "over? Once you have figured all this out, you can start to look into implementations. Do not focus on\n",
    "implementation details like Forward-Backward in the beginning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"sectrouble\"></a>7 Troubleshooting\n",
    "==\n",
    "\n",
    "There are some problems that can arise, and it is good to be aware of them and know how to deal with,\n",
    "or, better, avoid them.\n",
    "\n",
    "- **Problem**: What we maximize – $P(data)$ – is often not what we want to evaluate (say, the tagging\n",
    "accuracy as compared to gold data).\n",
    "\n",
    "- **Solution**: Design your model with this in mind and try to formulate the problem so that the two\n",
    "criteria are similar. Keep test data around and evaluate your models on it.\n",
    "\n",
    "\n",
    "- **Problem**: ***Local maxima***. EM improves P(data) at each iteration, but can get stuck in local \n",
    "maxima. Remember the graph from section 4: there are many suboptimal parameter configurations at\n",
    "which EM stops because it can not improve from there. The result is not the best model, though.\n",
    "\n",
    "- **Solution**: Restart EM 50 times or more with random initializations, and remember the model that\n",
    "got the best data likelihood. Each time you restart, you start at a different point along the curve, and\n",
    "hopefully eventually at one that leads to the global optimum.\n",
    "\n",
    "\n",
    "- **Problem**: EM does not care about semantics! Whether a label makes sense or not is irrelevant,\n",
    "as long as it explains the data! It uses this ***weird tag*** that only occurs once or twice, like `FW` (foreign \n",
    "word).\n",
    "\n",
    "- **Solution**: Use a dictionary that constrains EM’s $P(\\mathbf{w}|\\mathbf{t})$ choices to possible ones (all others will be\n",
    "0). If you do completely unsupervised training (no dictionary), the labels become meaningless, and\n",
    "the task is more like clustering. Each label becomes a cluster. In this case, you have to afterwards map\n",
    "each of the clusters to a tag in order to label data and get the accuracy. One mapping is many-to-1 \n",
    "(see [Johnson 2007](#refJohnson)). Also, keep your data clean. EM wants to use everything. Yes, also that weird tag\n",
    "that got in by accident... By assigning it to a frequent word, EM actually even thinks it has done a\n",
    "good job.\n",
    "\n",
    "\n",
    "- **Problem**: EM uses unlikely tags too frequently, i.e. we have a rather ***flat label distribution***. \n",
    "In language, however, most probability mass should go to one or two cases, the rest becomes less and \n",
    "less likely, i.e., we have a Zipfian distribution (see below).\n",
    "\n",
    "- **Solution**: Use smoothing and a dictionary. Additionally, you might want to try techniques like\n",
    "L$_0$ normalization ([Vaswani/Pauls/Chiang 2010](#refVaswani)) or Bayesian inference with Gibbs sampling and sparse\n",
    "priors.\n",
    "<img src=\"pics/distribution.png\" width=\"350px\"/>\n",
    "<div align=\"center\">*Figure 14: Difference between label distribution in EM and in reality*</div>\n",
    "\n",
    "\n",
    "\n",
    "- **Problem**: EM changes good initial parameters to make them worse. \n",
    "\n",
    "- **Solution**: Fix as many parameters as you can! That will guide EM to only optimize the right \n",
    "things. Also, add ***pseudo-counts*** before normalization to make changes less dramatic. The higher the\n",
    "pseudo-counts, the smaller the normalization effect. You can see in Table 1 how the normalized values\n",
    "of two parameters change when using different pseudo-counts. If you do not want to fix anything\n",
    "random restarts can help to find a good starting point.\n",
    "<img src=\"pics/smoothing.png\" width=\"400px\"/>\n",
    "<div align=\"center\">*Table 1: Influence of pseudo-counts on parameters*</div>\n",
    "\n",
    "\n",
    "- **Problem**: The resulting model is ***overfitting*** the training data (“Look, I can explain this data \n",
    "perfectly! And nothing else. . . \").\n",
    "\n",
    "- **Solution**: Again, use ***smoothing*** (add $n$ to the fractional counts before normalizing) and stop after a\n",
    "maximum amount of iterations (usually 50). This will prevent the data from being exactly modeled.\n",
    "\n",
    "\n",
    "- **Problem**: Ties. All transition options leaving from one node are equally good. EM doesn’t \n",
    "take a stance and just leaves all of them as is.\n",
    "\n",
    "- **Solution**: Start out randomly to avoid ties, and do restarts. This way, you can break the ties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secreading\"></a>8 Useful Reading\n",
    "==\n",
    "\n",
    "If you want to read the original, go for [Dempster/Laird/Rubin (1977)](#refDempster).\n",
    "[Rabiner/Juang (1986)](#refRabiner) is a general overview over ***parameter estimation***, but very math-heavy. [Manning/Schütze (2000)](#refMS) has a chapter on EM, based on clustering, but with an eye on other NLP applications.\n",
    "One of the most famous implementations for NLP is the tagging paper by [Merialdo (1994)](refMerialdo),\n",
    "which also gives you a good idea about the various parameters you can set.\n",
    "\n",
    "If you want to know more about the ***Forward-Backward*** procedure of calculating alphas and betas,\n",
    "check out Jason Eisner’s tutorial ([Eisner, 2002](#refEisner)), including a spreadsheet with the changing values.\n",
    "The second edition of the AI handbook ([Russell/Norvig, 2003, 724 – 733](refRN)) has a comprehensive\n",
    "section about EM.\n",
    "\n",
    "Kevin Knight has written a very compelling ***introduction*** ([Knight, 2009a](#refKnightA)). It is mainly about\n",
    "Bayesian Inference, but explains EM very nicely. You might also want to check out his tutorials\n",
    "for Carmel ([Knight, 2009b](#refKnightB)), a software that helps you implement graphical models as automata and\n",
    "train them with EM. Also, the workbook for MT ([Knight, 1999](#refKnightMT)) contains a useful section on EM.\n",
    "\n",
    "[Vaswani/Pauls/Chiang 2010](#refVaswani) show how using ***L$_0$ normalization*** can lead to smaller models and a L0\n",
    "sparser distribution, which improves language related tasks a lot, because it creates a more Zipfian normalization\n",
    "distribution (see [Hovy et al. 2011](#refHovy) for another application)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secadvanced\"></a>9 Advanced Topics\n",
    "==\n",
    "\n",
    "The goal of his section is mainly to give you an idea of some of the topics are out there to address\n",
    "problems with EM. It would lead too far to go into detail for each, so we will give a very high-level\n",
    "intuition and refer to related work for details.\n",
    "\n",
    "\n",
    "9.1 Scaling\n",
    "--\n",
    "\n",
    "In the examples above, we have used logarithms to prevent underflow of the probabilities. That turned\n",
    "all our multiplications into additions (which is good, because it’s fast), but it also turned every addition\n",
    "into a function call to a log addition (which unfortunately is very slow). There is another way to prevent\n",
    "underflow, called ***scaling***, but it requires some more bookkeeping. The idea is that we normalize each \n",
    "column in our lattice so that it sums to 1. That way, the probabilities won’t get too small, and we can\n",
    "still use multiplication and normal addition.\n",
    "\n",
    "Practically, we add a vector to our forward-backward procedure. It has the same length as out lattice,\n",
    "and for each position of the lattice contains the sum of the forward probabilities at that position, the\n",
    "***scaling factor***. Before we move on, we normalize the column by that factor, so that it sums to 1. The \n",
    "backward pass uses the same scaling parameters as the forward pass.\n",
    "We can then use the scaled forward and backward matrices as well as the scaling factor vector to\n",
    "compute the fractional counts. See [Shen (2008)](#refShen) for more details (abstract mathematical notation, but\n",
    "very good explanations).\n",
    "\n",
    "\n",
    "9.2 Variational Bayes\n",
    "--\n",
    "\n",
    "We have earlier seen the idea of adding pseudo-counts to the fractional counts, which helps with\n",
    "smoothing. Variational Bayes inference works similarly, only that now we define a prior for our prior\n",
    "parameters.\n",
    "\n",
    "A ***prior*** is essentially a curve that looks similar to what we would like to achieve (see the example\n",
    "above of how EM’s distribution differs from the real one: a prior would look like the real one). In\n",
    "practice, the prior here is a ***Dirichlet distribution***, which takes two parameters: a probability distribution and a vector of shape parameters. Both have the same number of elements (or dimensions). If \n",
    "we have only two dimensions in each, the prior is called a Beta function. Typically, the probability\n",
    "distribution are our transition parameters, and the shape parameters are something like pseudo-counts\n",
    "for each element.\n",
    "\n",
    "This sounds rather complicated, but is relatively easy to implement. The E-step stays the same\n",
    "as before. In the M-step, the elements of the shape parameter vector are added as pseudo-counts\n",
    "to the matching fractional counts, and the result is passed through a ***Digamma function*** and finally\n",
    "exponentiated. To normalize, we add the sum of all elements in the shape vector to our denominator\n",
    "and again run it though Digamma and exponentiation. This is like a softer version of smoothing with\n",
    "pseudo-counts, where we have separate counts for each parameter.\n",
    "\n",
    "\n",
    "9.3 Type and Token Constraints\n",
    "--\n",
    "\n",
    "So far, we have pretended that every word can be labeled with all tags. However, we know that that is\n",
    "not true: for example, “I” can never be a preposition or an adverb. One way to encode that knowledge\n",
    "is to add a ***dictionary***. For each word, it lists the allowed tags. Essentially, this sets all \n",
    "parameters for disallowed tags to 0. Using a dictionary encodes general type constraints for a word, type\n",
    "i.e., the ***constraints*** hold whenever we see the word anywhere in our data. Dictionaries are either built \n",
    "by hand or derived from some annotated data. In our example, we could restrict the legal tags for the\n",
    "word “I” to just nouns.\n",
    "\n",
    "There is another type of constraints, namely ***token constraints***. A token is a specific word in a \n",
    "specific context, so these constraints only hold for certain words in certain contexts. In our example,\n",
    "we just might happen to know that the first “can” is a verb. ***Type constraints*** can come from simple\n",
    "rules (e.g., something like “any word after the is either an adjective or a noun”), manual annotations\n",
    "(going through the data and tagging some sentences) or through transfer from another data set.\n",
    "\n",
    "Both types of constraints basically do the same thing, though: they reduce the number of nodes for a word\n",
    "in our lattice. Type constraints do that for all occurrences of the word, token constraints only for some.\n",
    "When we have an actual annotation for a word, the number of nodes is reduced to one. This greatly\n",
    "reduces the number of legal paths through a sentence, so not only does it make forward-backward\n",
    "more efficient, it also guides the algorithm because you can now only collect fractional counts for\n",
    "certain parameters. See [Täckström et al. (2013)](#refOscar) for examples of using type and token constraints in an\n",
    "application."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secfurther\"></a>10 Going Further...\n",
    "==\n",
    "\n",
    "If you want to get deeper into the matter, you could look into ***EM with features*** (Berg-Kirkpatrick et al., 2010). Instead of just using conditional probabilities, which often cannot capture \n",
    "useful properties (or only when encoded as additional states), you can add all the features you like in\n",
    "discriminative models (like word suffixes or capitalization) and still do unsupervised learning.\n",
    "\n",
    "Or explore ***structural EM***, which not only learns the parameter values, but also how many parameters there should be.\n",
    "\n",
    "If you do any of these things, or if you have questions, comments, or criticism, send me a mail—I’d\n",
    "be dead curious to know!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secacknowledgements\"></a>Acknowledgements\n",
    "==\n",
    "\n",
    "Thanks to Kevin Knight and David Chiang for the basics, Congxing Cai, Karl-Moritz Hermann, and\n",
    "Eduard Hovy for useful comments, and Ashish Vaswani, Victoria Fossum, and Taylor Berg-Kirkpatrick\n",
    "for many enlightening discussions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References<a name=\"secref\"></a>\n",
    "==\n",
    "\n",
    "<a name=\"refTaylor\"></a>BERG-KIRKPATRICK, TAYLOR ET AL. (2010): Painless Unsupervised Learning with Features. In: North American Chapter\n",
    "of the Association for Computational Linguistics.\n",
    "\n",
    "<a name=\"refBishop\"></a>BISHOP, C.M. (2006): Pattern recognition and machine learning. New York: Springer.\n",
    "\n",
    "<a name=\"refDempster\"></a>DEMPSTER, ARTHUR P./LAIRD, NAN M./RUBIN, DONALD B. (1977): Maximum likelihood from incomplete data via the\n",
    "EM algorithm. In: Journal of the Royal Statistical Society. Series B (Methodological), 39, Nr. 1, 1–38.\n",
    "\n",
    "<a name=\"refEisner\"></a>EISNER, JASON (2002): An interactive spreadsheet for teaching the forward-backward algorithm. In: Proceedings of the\n",
    "ACL-02Workshop on Effective tools and methodologies for teaching natural language processing and computational\n",
    "linguistics-Volume 1. Association for Computational Linguistics, 10–18.\n",
    "\n",
    "<a name=\"refMACE\"></a>HOVY, DIRK ET AL. (2013): Learning Whom to trust with MACE. In: Proceedings of NAACL HLT.\n",
    "\n",
    "<a name=\"refHovy\"></a>HOVY, DIRK ET AL. (2011): Unsupervised Discovery of Domain-Specific Knowledge from Text. In: Proceedings of the\n",
    "49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies. Portland,\n",
    "Oregon, USA: Association for Computational Linguistics, pages 1466–1475. [http://www.aclweb.org/anthology/P11-1147.pdf](http://www.aclweb.org/anthology/P11-1147.pdf)\n",
    "\n",
    "<a name=\"refJohnson\"></a>JOHNSON, MARK (2007): Why doesn’t EM find good HMM POS-taggers. In: Proceedings of the 2007 Joint Conference\n",
    "on Empirical Methods in Natural Language Processing and Computational Natural Language Learning (EMNLPCoNLL).,\n",
    "296–305.\n",
    "\n",
    "<a name=\"refKnightMT\"></a>KNIGHT, KEVIN (1999): A statistical MT tutorial workbook. JHU Summer Workshop, 1999.\n",
    "\n",
    "<a name=\"refKnightA\"></a>KNIGHT, KEVIN (2009a): Bayesian Inference with Tears. [http://www.isi.edu/natural-language/people/bayeswith-tears.pdf](http://www.isi.edu/natural-language/people/bayeswith-tears.pdf).\n",
    "\n",
    "<a name=\"refKnightB\"></a>KNIGHT, KEVIN (2009b): Training Finite-State Transducer Cascades with Carmel. 2009b.\n",
    "\n",
    "<a name=\"refMS\"></a>MANNING, C.D./SCHÜTZE, H. (2000): Foundations of statistical natural language processing. MIT Press.\n",
    "\n",
    "<a name=\"refMerialdo\"></a>MERIALDO, BERNARD (1994): Tagging English text with a probabilistic model. In: Computational linguistics, 20, Nr. 2,\n",
    "155–171.\n",
    "\n",
    "<a name=\"refRabiner\"></a>RABINER, L./JUANG, B. (1986): An introduction to hidden Markov models. In: IEEE ASSp Magazine, 3, Nr. 1 Part 1, 4–16.\n",
    "\n",
    "<a name=\"refRN\"></a>RUSSELL, S.J./NORVIG, P. (2003): Artificial intelligence: a modern approach. 2nd edition. Upper Saddle River, NJ: Prentice\n",
    "Hall.\n",
    "\n",
    "<a name=\"refShen\"></a>SHEN, DAWEI (2008): Some Mathematics for HMMs. [http://courses.media.mit.edu/2010fall/mas622j/ProblemSets/ps4/tutorial.pdf](http://courses.media.mit.edu/2010fall/mas622j/ProblemSets/ps4/tutorial.pdf) – visited on 3/21/2013.\n",
    "\n",
    "<a name=\"refOscar\"></a>TÄCKSTRÖM, OSCAR ET AL. (2013): Token and Type Constraints for Cross-Lingual Part-of-Speech Tagging. In: Transactions\n",
    "of the ACL.\n",
    "\n",
    "<a name=\"refVaswani\"></a>VASWANI, ASHISH/PAULS, ADAM/CHIANG, DAVID (2010): Efficient optimization of an MDL-inspired objective function\n",
    "for unsupervised part-of-speech tagging."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"secanswers\"></a>Answers\n",
    "===\n",
    "\n",
    "1. <a name=\"answer1\"></a>$P(x|y) = \\frac{P(x,y)}{P(y)}$ [back](#q1)\n",
    "\n",
    "2. <a name=\"answer2\"></a>$W$ is discrete, $T$ is discrete, but binary, and $L$ is boolean, and thus also binary. [back](#q2)\n",
    "\n",
    "3. <a name=\"answer3\"></a>We count the occurrences of each rule S!NP VP in a treebank and normalize by the number of times we have seen S. [back](#q3)\n",
    "\n",
    "4. <a name=\"answer4\"></a>Exactly, CKY parsing... You are getting good at this! [back](#q4)\n",
    "\n",
    "5. <a name=\"answer5\"></a>We compute the distance between the cluster centroids and each point and assign it to the closest centroid. [back](#q5)\n",
    "\n",
    "6. <a name=\"answer6\"></a>We divide $P(t) \\times P(w|t) = P(w) \\times P(t|w)$ by $P(w)$ [back](#q6)\n",
    "\n",
    "7. <a name=\"answer7\"></a>You already had your model and would not need EM... [back](#q7)\n",
    "\n",
    "8. <a name=\"answer8\"></a>Because words tend to be ambiguous. What is the most likely tag for \"can\"? It depends. It could be `NOUN`, `VERB`, or `AUX`. But what is the most likely tag if we see \"the can\"? Still `VERB` or `AUX`? Probably not... The context disambiguates it. That’s why we want the transition probabilities. [back](#q8)\n",
    "\n",
    "9. <a name=\"answer9\"></a>The probability that we end up with this sentence if we ran our model in generation mode. [back](#q9)\n",
    "\n",
    "10. <a name=\"answer10\"></a>Well, Forward and Backward... What did you think? Though you could of course call Forward \"Baum\" and Backward \"Welch\", but that seems a little unfair to Mr. Welch. [back](#q10)\n",
    "\n",
    "11. <a name=\"answer11\"></a>Assigning each data point to one of the centroids (= the current model parameters). [back](#q11)\n",
    "\n",
    "12. <a name=\"answer12\"></a>Because in both cases we walked through all possible paths of the lattice, summing them up. [back](#q12)\n",
    "\n",
    "13. <a name=\"answer13\"></a>The count of the two events together and the counts of the given part. [back](#q13)"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}