2018-05-13-viterbi-message-passing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inference in Hidden Markov Models continued: Viterbi Algorithm\n",
    "\n",
    "This is a short post that continues from [the more-detailed alpha recursion HMM post](2018-05-02-hmm-alpha-recursion.ipynb). In this post I'll implement the [Viterbi algorithm](https://en.wikipedia.org/wiki/Viterbi_algorithm) like Barber does in [\"Bayesian Reasoning and Machine Learning\"](http://www.cs.ucl.ac.uk/staff/d.barber/brml/). Like before, I'm porting the MatLab code from the textbook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import nb_code.hmm_alpha_recursion as prev_post"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# helper functions you can skip over :D\n",
    "SAVE = True\n",
    "def maybe_save_plot(filename):\n",
    "    if SAVE:\n",
    "        plt.tight_layout()\n",
    "        plt.savefig('images/' + filename, bbox_inches=\"tight\")\n",
    "\n",
    "def hide_ticks(plot):\n",
    "    plot.axes.get_xaxis().set_visible(False)\n",
    "    plot.axes.get_yaxis().set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viterbi algorithm in HMMs using message passing\n",
    "\n",
    "The Viterbi algorithm finds the most-likely path $h_{1:T}$ for the visibles $v_{1:T}$, where $T$ is the timestep of the last observed visible.\n",
    "The algorithm takes in the visibles $v_{1:T}$, the initial state $p(h_1)$, the transition probabilities $p(h_{t} \\mid h_{t - 1})$, and the emission probabilities $p(v_t \\mid h_t)$, and returns the list of most-likely hidden states $h_{1:T}$. \n",
    "\n",
    "I generated the required probabilities in the [alpha recursion HMM post](2018-05-02-hmm-alpha-recursion.ipynb).\n",
    "\n",
    "### Algorithm\n",
    "\n",
    "Barber frames the Viterbi algorithm as message passing using the max-product algorithm.\n",
    "\n",
    "This version of the algorithm begins at the end of the hidden states ($h_T$), and computes an incoming message from future states. The message is meant to represent the effect of maximizing over those states.\n",
    "Barber gives the messages as:\n",
    "\n",
    "$$\\mu(h_T) = 1$$\n",
    "\n",
    "$$\\mu(h_{t - 1}) = \\max_{h_t} p(v_t \\mid h_t)p(h_t \\mid h_{t - 1})\\mu(h_t).$$\n",
    "\n",
    "Once the messages are computed, the algorithm then computes the most-likely state for $h_1$, and uses that to compute the most-likely state for $h_2$ and so on. It basically maximizes the marginal of $p(h_t \\mid v_{1:T})$ and then uses the most-likely state for $h_t$ in the transition matrix for computing $p(h_{t + 1} \\mid v_{1:T})$ so it returns a valid path.\n",
    "\n",
    "$$h_1^* = \\max_{h_1} p(v_1 \\mid h_1)p(h_1)\\mu(h_1)$$\n",
    "\n",
    "$$h_t^* = \\max_{h_t} p(v_t \\mid h_t)p(h_t \\mid h_{t - 1}^*)\\mu(h_t).$$\n",
    "\n",
    "Now in Python!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[15  9  3  2  1  2  3  4 10  9]\n"
     ]
    }
   ],
   "source": [
    "def viterbi(visibles, p_hidden_start, p_transition, p_emission):\n",
    "    num_timestamps = visibles.shape[0]\n",
    "    num_hidden_states = p_transition.shape[0]\n",
    "    \n",
    "    # messages[t] corresponds to mu(h_t), which is the message coming into h_t\n",
    "    messages = np.zeros((num_timestamps, num_hidden_states))\n",
    "    \n",
    "    most_likely_states = np.zeros((num_timestamps,), dtype=int)\n",
    "\n",
    "    # The message coming into the last node is 1 for all states\n",
    "    messages[-1] = np.ones(num_hidden_states)\n",
    "\n",
    "    # normalize!\n",
    "    messages[-1] /= np.sum(messages[-1])\n",
    "    \n",
    "    # Compute the messages!\n",
    "    for t in range(num_timestamps - 1, 0, -1):\n",
    "        # use the data at time t to make mu[h_{t - 1}]\n",
    "        \n",
    "        # compute max p(v|h)p(h|h)mu(h)!\n",
    "        \n",
    "        # compute p(v|h)mu(h)\n",
    "        message_and_emission = messages[t] * p_emission[visibles[t]]\n",
    "        \n",
    "        # compute p(v|h)p(h|h)mu(h)\n",
    "        # message_and_emission.reshape(-1, 1): new_state x 1\n",
    "        # np.tile(...): new_state x old_state\n",
    "        # p_transition: new_state x old_state\n",
    "        # np.tile(...) * p_transition: new_state x old_state\n",
    "        all_h_ts = np.tile(\n",
    "            message_and_emission.reshape(-1, 1),\n",
    "            (1, num_hidden_states)\n",
    "        ) * p_transition\n",
    "        \n",
    "        # the message is the value from the highest h_t\n",
    "        messages[t - 1] = np.max(all_h_ts, axis=0)\n",
    "        \n",
    "        # and normalize\n",
    "        messages[t - 1] /= np.sum(messages[t - 1])\n",
    "    \n",
    "    # now from the beginning! compute h_t* using these messages\n",
    "    \n",
    "    # argmax will give us the state.\n",
    "    # argmax p(v_1|h_1)p(h_1)mu(h_1)\n",
    "    most_likely_states[0] = np.argmax(\n",
    "        p_hidden_start \n",
    "        * p_emission[visibles[0]] \n",
    "        * messages[0]\n",
    "    )\n",
    "    \n",
    "    for t in range(1, num_timestamps):\n",
    "        # argmax_h_t p(v_t|h_t)p(h_t|h_{t - 1})mu(h_t)\n",
    "        most_likely_states[t] = np.argmax(\n",
    "            p_emission[visibles[t], :]\n",
    "            * p_transition[:, most_likely_states[t - 1]] \n",
    "            * messages[t]\n",
    "        )\n",
    "    \n",
    "    return most_likely_states\n",
    "\n",
    "most_likely_states = viterbi(\n",
    "    prev_post.visibles, \n",
    "    prev_post.p_hidden_start,\n",
    "    prev_post.p_transition,\n",
    "    prev_post.p_emission,\n",
    ")\n",
    "\n",
    "print(most_likely_states)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing\n",
    "\n",
    "I can plot the most-likely states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIAAAAEYCAYAAAAkpY3TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHa9JREFUeJzt3XuUZWdZJ+DfW1Xp7lwgITdCEkgLhBAY1JlEbsIICwYS\nghJZLHQQIiAoumR0vIyIAxNmeZ1BHB1xQB25yEXQ0ahESIYlMkRASDuCMoFAkk4CCSFNEjrXvlR9\n80dVa9F2c77q7upTe9fzrFVrde16691f1Vt71zm/3vtUtdYCAAAAwHjNTHsBAAAAAKwuARAAAADA\nyAmAAAAAAEZOAAQAAAAwcgIgAAAAgJETAAEAAACMnAAIAAAAYOQEQAAAAAAjJwACAAAAGLm5aS9g\naKqqTXsN60VrrQ62x4knntg2b958CFZz+GzZsmXaSzhQ21prJx1MgyEeX+ecc860l3BAtmzZctDz\nSszscNm6dWu2bdt20OfEIc5rwJwTB8Q5cVicEwfpkBxjHtsfPp6LDU7XMSYAYtQ2b96cK6+8ctrL\nWJGqgz7XTsv1017ANAzt52uPqlqX80qGObNzzz132ktg5dblMTbE4ytxThwa58RBOiTHmMf2w2Je\nh1XXMeYWMAAAAICREwABAAAAjJwACAAAAGDkBEAAAAAAIycAAgAAABg5ARAAAADAyAmAAAAAAEZO\nAAQAAAAwcgIgAAAAgJETAAEAAACMnAAIAAAAYOQEQAAAAAAjJwACAAAAGDkBEAAAAMDICYAAAAAA\nRk4ABAAAADByAiAAAACAkRMAAQAAAIycAAgAAABg5A4qAKqqN1XVazrqPlNVT1n698VV9Y5vULu1\nqp5+MOsCAAAA4J/MfaMPVtUHknyitfbavbY/J8mbk5zeWts9aSettUcf1CoBAAAAOGCTrgB6W5IX\nVlXttf1FSd7ZE/4AAAAAMF2TAqBLkpyQ5Ml7NlTVA5I8O8nbq+qtVfXzS9tPrKr3VdUdVXVbVX2k\nqmaWPrb3bV2bquo9VXVnVf1tVX3LvnZeVTNV9aqquqaqvlpV762q45c+tqmq3rG0/Y6q+mRVPfDA\nvxUAAAAA4/QNA6DW2r1J3pvkomWbn5/ks621T+1V/pNJvpjkpCQPTPLqJG0/rZ+T5A+THJ/kXUku\nqaoj9lH3yiQXJvmOJKcmuT3JG5c+9v1Jjk3y4CyGVK9Icu83+noAAAAA1qOeF4F+W5LnVdWmpfcv\nWtq2t11JHpTkjNbartbaR1pr+wuAtrTW/qi1tivJG5JsSvL4fdS9IsnPtda+2FrbkeTipbXMLe3v\nhCQPb63Nt9a2tNa2d3w9AAAAAOvKxACotXZFkm1JLqyqhyV5bBav2tnbf03yhSSXV9W1VfWqb9D2\nxmX9F7J45dCp+6g7I8mfLN3idUeSq5LMZ/EKo99PclmSP6iqm6rqv+znKiIAAACAda33z8C/PYtX\n/rwwyWWttVv2Lmit3dla+8nW2kOTfFeSn6iqp+2n34P3/GPpdYJOT3LTPupuTHJ+a+24ZW+bWmtf\nWrrK6HWttUcleWIWX5foon30AAAAAFjXVhIAPT3Jy7Pv279SVc+uqocv/cWwr2XxSp2F/fQ7p6qe\nu3Qr148n2ZHk4/uoe1OSX6iqM5b2cdLSn6BPVT21qh5TVbNJtmfxlrD97Q8AAABg3eoKgFprW5N8\nNMnRSf5sP2VnJvlgkruSfCzJb7XWPrSf2j9N8j1ZfFHnFyV57tLrAe3t15f2d3lV3ZnFkOhxSx87\nJckfZTH8uSrJh7N4WxgAAAAAy8z1FrbWnrKPbS9e9u9fS/Jr+/nczcv+ffGE/SyvXcjii0S/YR91\n707y7gnLBgAAAFj3em8BAwAAAGCgBEAAAAAAIycAAgAAABg5ARAAAADAyAmAAAAAAEZOAAQAAAAw\ncgIgAAAAgJETAAEAAACMnAAIAAAAYOQEQAAAAAAjJwACAAAAGDkBEAAAAMDICYAAAAAARq5aa9Ne\nw6BU1a1Jrp/2OtaBM1prJx1sE/M6rA56ZuZ1WDnGhsW8hsc5cVgcY8NiXsNjZsNiXsPTNTMBEAAA\nAMDIuQUMAAAAYOQEQAAAAAAjJwACAAAAGDkBEAAAAMDICYAAAAAARk4ABAAAADByAiAAAACAkRMA\nAQAAAIycAAgAAABg5ARAAAAAACMnAAIAAAAYOQEQAAAAwMgJgAAAAABGTgAEAAAAMHICIAAAAICR\nEwABAAAAjJwACAAAAGDkBEAAAAAAIycAAgAAABg5ARAAAADAyM1NewHrQVW1aa9hira11k6a9iJW\nysyGNbPVmNc555zTVbdly5ZDveuVGty8kukeY72zTVZnvq21OuRNV9lQzomrdNwO7hgbyrx6rfCY\nHdy8kvHNbIUGNzPzMq8BGdy8EjNbzZkJgFht1097AayYmSW58soru+qqpv5c3rxWqHe2yZqYLyuw\nSsetY2zKVnjMmtfwmNmwmNewmNfwrOrM3AIGAAAAMHICIAAAAICREwABAAAAjJwACAAAAGDkBEAA\nAAAAIycAAgAAABg5ARAAAADAyAmAAAAAAEZubtoLgPWitdZVV1WrvBJ6mMN4me14me04mevw9D7m\nScwX4HByBRAAAADAyAmAAAAAAEZOAAQAAAAwcgIgAAAAgJETAAEAAACMnACoQ1U9paq+OO11AAAA\nAByIwQRAVfVXVXV7VW3sqN1cVa2q/Jl7AAAAYN0bRABUVZuTPDlJS/JdU10MAAAAwMAMIgBKclGS\njyd5a5Lv37Oxqo6sql+tquur6mtVdUVVHZnk/yyV3FFVd1XVE6rq4qp6x7LP/bqrhKrqJVV1VVXd\nWVXXVtUPHb4vDwAAAGD1DOUWqYuSvCHJ3yT5eFU9sLV2S5LXJ3l0kicm+XKSxyVZSPKvk1yX5LjW\n2u4kqapnTtjHV5I8O8m1S5///qr6ZGvtb1fh6zlorbWuuqpa5ZXQyywAgPXAY55h6X1ekZjtavL8\njsNhzV8BVFVPSnJGkve21rYkuSbJC6pqJslLk/xYa+1LrbX51tpHW2s7DmQ/rbVLW2vXtEUfTnJ5\nFm87AwAAABi0NR8AZfGWr8tba9uW3n/X0rYTk2zKYiB00Krq/Kr6eFXdVlV3JHnW0j4AAAAABm1N\n3wK29Ho+z08yW1VfXtq8MclxSR6U5L4kD0vyqb0+dV/Xz92d5Khl75+ybD8bk/yvLN5q9qettV1V\ndUkS19cBAAAAg7fWrwC6MMl8kkcl+dalt7OTfCSLYc3vJXlDVZ1aVbNLL/a8McmtWXwtoIcu6/V3\nSf51VT2kqo5N8rPLPrYhi8HSrUl2V9X5SZ6xul8aAAAAwOGx1gOg70/yltbaDa21L+95S/KbSb4v\nyauS/H2STya5LcmvJJlprd2T5BeS/HVV3VFVj2+t/e8k70ny6SRbkrxvz05aa3cm+XdJ3pvk9iQv\nSPJnh+uLBAAAAFhNtZJXfefAVNUh/yYP6FXit7TWzp32IlZqNWY2IIObmXkNa17J+p5Za23qJ+aV\nWs/zygCPMfMa1rwSMxvazMY2rxX+FTDzWiWr9PxucPNKhjOzVbKqM1vrVwABAAAAcJAEQAAAAAAj\nJwACAAAAGDkBEAAAAMDIzU17ARyYNfDizgAAwMB5XrE2mMPwDOgPM/0jVwABAAAAjJwACAAAAGDk\nBEAAAAAAIycAAgAAABg5ARAAAADAyAmAAAAAAEZOAAQAAAAwcgIgAAAAgJETAAEAAACM3Ny0FwAA\nAAAwJFU17SWsmCuAAAAAAEZOAAQAAAAwcgIgAAAAgJETAAEAAACMnAAIAAAAYOQEQAAAAAAjJwAC\nAAAAGDkBEAAAAMDICYAAAAAARk4ABAAAADByc9NewDqxLcn1017ElJwx7QUcIDMbFvManvU6M/Ma\nniHOzLyGx8yGxbyGxbyGx8xWSbXWVrM/AAAAAFPmFjAAAACAkRMAAQAAAIycAAgAAABg5ARAAAAA\nACMnAAIAAAAYOQEQAAAAwMgJgAAAAABGTgAEAAAAMHICIAAAAICREwABAAAAjJwACAAAAGDkBEAA\nAAAAIycAAgAAABg5ARAAAADAyAmAAAAAAEZOAAQAAAAwcnPTXsB6sGH2yHbkEcdOrDvzrDu6e37+\nc8cdzJIOev+9tnx6x7bW2kmHvPEq21Ab26YcPbmwVtC0HfBy9r/7I47oLOxf6PadtwxuZhtmNrUj\nZ46ZWNfmF7p77nxQx/yTtBWcRTfduruvcMfO7p7b222Dm1eydF6cm3xeTOs/cHacuKGvcAXH7cZb\n7uuqawt9P1v3tbuzs923kjPHmrChNrZNNfmY2HBW//8r7fpCb23/z8B9p27sqtt08wqOsV23Du4Y\nm73f0W3uxAdMrNu49Z7unjtP6zsnzvR/azN3691ddTXT/3O1feGrg5tXkmyoTV3H2I7TjuruefQx\nfeevh264s7vn5z/b9/iz7ez/Qbgztw9uZouPO+43sW7HyZu6e278auf3bAWPZe47rfOceNOO7p7b\nd28b3rw6H9fXbP+5pvsx5Qp+41dv8YbOx/9Jtu8Y3uP6ZAXPxaas5vqeCLTdnc8BsvrnRAHQYXDk\nEcfmiae/aGLdpZdd0t3zgiddeDBLOuj995p90BeuP+RND4NNOTqPm3n6xLqane3uuZIDv9fcKaf1\nFR7Rf6h/4NpfHdzMjpw5Jk845jkT6+a3b+/uecMPPrGrbscJ/Q/EHvlb27rq2vVf7O55+b3vGNy8\nkuTIuWPzxNO+b3Lhrv7j5roXb+6qm9/UHyg87A2f7apbuKfvifTHd7y/e99ryaY6Oo+fe+bEutPe\n1v9k5+bndD6RnZ/v7nnVzz2sq+5Rr7uhu+cHbvrNwR1jcyc+IA963Y9OrDvzJX/b3fO6H318V90x\nN3a3zElv/kRX3czR/aHH5dvfMrh5JUvH2BHnTay75sfO6e75bd/ed/561zd9qLvnBd8++Xdtkuze\n2n+MfXDhDwc3syNn7pcnHPvdE+u2vuzs7p7f9Pt9B0/bfld3z6tefWZX3dmvuaa752W3vnlw89qU\no/O42WdMrJs9pj9wmL+rN8DuT4B6w4SZU0/p7vmBa14/uHklSzOrpx3apr3/Ib6C/3ycPfHkrrr5\nW77S3fOD7Y9WdWZuAQMAAAAYOQEQAAAAwMgJgAAAAABGTgAEAAAAMHICIAAAAICREwABAAAAjJwA\nCAAAAGDkBEAAAAAAIzc37QWsB2eedUcuveySiXUXPOnCw7Cag9//pVdM/loGr5KanZ1Y1hZad8vr\nfvkJB7Oifbr6ov/RVfesf/M9h3zfa8rcEckpJ00su+a1j+pu+bCf/lhXXc/PyR4L1Ze5zzzim7p7\n5h/6S9eU3buz8NXbJ5bN/4uHdrd88C/2zWwl5lvfMV5zvb9O+88Za8n88Ufntu/8tsmFT97S3XP2\nhKO76lbyHXvEj3Tu/+QTV9B1eDZuvSdnvnjy9+Kym/6uu+czT12Fn92qrrKFu+859PteY3adfFRu\n/r5zJ9ad+fP9J/3b5+e76p557zndPbNwfVdZ/zkxyUJ/6ZrRWtrOnRPLNv/OF/pb7tzVV7iC7+0j\nfuSTXXV9PynDVbMzmT1m8u+chXv6zzXvu/ETXXXnv/Dl3T0/+M7f66o774zHdvccqpqdzez9j51Y\nN/+17d09Z775kV1173//u7t7nvddL+yqm13Bc8Z8pb/0QLgCCAAAAGDkBEAAAAAAIycAAgAAABg5\nARAAAADAyAmAAAAAAEZuTQdAVfWQqrqrqmaX3v+rqnrZIez/5Kr63KHqBwAAALAWrZkAqKq2VtW9\nS4HPXVV1V5LdrbVjWmv/7K8TVtWLq+qKg9lna+0jrbWzDqYHAAAAwFq3ZgKgJd+5FPjsebtptXZU\nVXOr1RsAAABgLVlrAdDXqarNVdX2Dmuq6uwkb0ryhKWrhe5Y2r6xql5fVTdU1S1V9aaqOnLpY0+p\nqi9W1c9U1ZeTvGXPtmV9t1bVT1XVp6vqa1X1nqratOzj/6Gqbq6qm6rqZUtre/hh+WYAAAAAHKBB\nXgXTWruqql6R5GWttSct+9AvJ3lYkm9NsivJu5K8NsnPLn38lCTHJzkji+HX4/bR/vlJzktyX5K/\nTvLiJG+qqvOS/ESSpyW5Lslv96738587Lhc86cKJdZdecUlvy1XRs8b14oRH7cwL/njrxLo/eMJj\nunue+Zs3HMSK9u3S522aXJSk7t1xyPe9lrQdOzJ/9TUT6x7xP1eQeZ9wfN++T39gd8t21eQ1Jkl2\n7uruOVRtYSELd98zsW7uM9f199y4sa9u9+7+np21tWFDX938mv5/l/2ave3uHP/OT06s+8ANV3b3\nfOpLX95Vd+THr+7umbbQVbaw/c7+nkNVNbHkvIec291u9oT7d9Ut3HV3d8+2o/d3U99ch+yIW+7O\nKf/tYxPrbv/efT103bfj3veZg1nSQVnJeXaQ5uYyc9IJE8tueN7p3S3vemjf9+yRP72CubbWVzcz\n299zgHZu3pAbXv/giXUP+Y/9j78uPPfErrojTr63u+cF/+qZXXVt1y3dPQdrbjbpOMZmawWPq778\n1a6y8854bHfLh3/08111n3/szu6eq22tPRK9pKruWHpbURpSVZXkB5P8+9baba21O5P8YpLvXVa2\nkOQ/tdZ2tNb2dzT+RmvtptbabUn+PIthUrIYDL2ltfaZ1to9SS5eyfoAAAAApmWtXQF0YWvtg3ve\nqarNK/jck5IclWRL/dP/elWS5ZH2ra21+yb0+fKyf9+T5NSlf5+aZPl/bd64grUBAAAATM1aC4BW\nYu9rGrcluTfJo1trX+r8nJW4Ocny6zgnX0cIAAAAsAastVvAVuKWJKdX1YYkaa0tJPmdJL9WVScn\nSVWdVlV9N1NO9t4kL6mqs6vqqCSvOUR9AQAAAFbVkAOgv0zymSRfrqptS9t+JskXkny8qrYn+WCS\nsw7Fzlpr70/yG0k+tGcfSx8a96vrAgAAAIO3Zm4Ba61t3se2rVl8HZ897z9l2b93Jrlgr/r7krx6\n6W3vXn+Vr7+F659t23sNrbWL93r/l5L8UvKPf4p+IYu3hgEAAACsWUO+Auiwq6rvrqqNVfWAJL+S\n5M9bayP/O5cAAADA0AmAVuaHknwlyTVJ5pP88HSXAwAAADDZmrkFbAhaa+dNew0AAAAAK+UKIAAA\nAICRcwXQGnLBky485D0vveKSQ95zPfjqZzbknY88fXJhbu9ventf7ewDHtDd8su7juuqa9vv7O45\nRI/45nty2WV/N7HugnNO6e65cNrJnYUL3T3brr6XDNv1wPt398zV/aVrSc3OZvbYyV/nX3zmQ909\nn/bCH+iqO+Ijf9/dM1WTa5Is3Nf3ByHbCn5e1pLasDEzm8+YWHf+WUd29zzqxK901e18zEO7e858\n5P921dXcuB/+1JGbMvPIsyfWffZHj+ru+fC39p2/Zv76U909u7V26HuuMQ95zF3575deMbHu2l3/\nr7vneb/ad146/6wnd/dcuHPcjye6zc+n3fG1iWWnXtH3OC1JZn/3+q66G175Ld09T/vlj/YVLsx3\n9xyiDdftzINftHVi3fw993T3rLkjuuraLbd29/zcr5/bVXfmK2/p7jlUbcfOzH/huo7CQ//7Ye5B\n/c8Xrn1a38/MzMa+x5NJknv7Sw+EK4AAAAAARk4ABAAAADByAiAAAACAkRMAAQAAAIycAAgAAABg\n5ARAAAAAACMnAAIAAAAYOQEQAAAAwMgJgAAAAABGbm7aC1gPzjzrjlx62SUT6y540oXdPXdfu7Wr\nbto9k9evoHaNqZpc01p/v5nZrrKbX3B2d8s/fMydXXXX/udHdPfMz/aXrhVXf/qoPPPUb51YN3v/\ne7p71vX3dtXN3/G17p67nnFuV90Rl1/Z3XOwZmdTx95vYtkFj3t2d8sjbv5UV13bvbu7Z7c2f+h7\nriFttrJw9KaJdTMPOrm75+6rr+mqm7nu+u6evea3bz/kPdeSdt99af9w9cS6R/7Gw7t7LvzD5zt3\nvoLfi/yjG686Lj/+2O+eWPfGT/xxd89n/cvndtUt3PmV7p4sOvPRd+YvLvvwxLpnPfqp3T3nt9/V\nVXfar3ysuyeLWmtpO3dNrJs56qj+njt2dNVd/cZzuns+4of/pruWlZu9//276ua33dbdc+b447rq\nFtbQ4w5XAAEAAACMnAAIAAAAYOQEQAAAAAAjJwACAAAAGDkBEAAAAMDICYAAAAAARk4ABAAAADBy\nAiAAAACAkRMAAQAAAIycAAgAAABg5Kq1Nu01jF5V3Zrk+mmvY0rOaK2dNO1FrJSZDWtm5jWseSXr\nembmNTyDm5l5DWteiZkNbWbmZV4DMrh5JWa2mjMTAAEAAACMnFvAAAAAAEZOAAQAAAAwcgIgAAAA\ngJETAAEAAACMnAAIAAAAYOQEQAAAAAAjJwACAAAAGDkBEAAAAMDICYAAAAAARk4ABAAAADByAiAA\nAACAkRMAAQAAAIycAAgAAABg5ARAAAAAACMnAAIAAAAYOQEQAAAAwMjNTXsB60FVtWmvYYq2tdZO\nmvYiVsrMhjWzac7rnHPO6a7dsmXLaixhcPNKVmdmvbNYpTl0a63VVBdwAHrntQaOh9UwuGNsPR9f\nGeC8Er/HhjazoTxOXKXj1ryGZXDzSsxsNWdWra3n7+3hsc5/gLe01s6d9iJWysyGNbNpzmsl59Cq\nVXneP7h5Jaszs95ZrNIcuo05AFoDx8NqGNwxtp6PrwxwXonfY0Ob2VAeJ67ScWtewzK4eSVmtpoz\ncwsYAAAAwMgJgAAAAABGTgAEAAAAMHICIAAAAICREwABAAAAjJwACAAAAGDkBEAAAAAAIycAAgAA\nABi5uWkvAPh6rbXu2qpaxZXQwwzWDrOYPjMYL7MdL7MdL7Mdlt7nAOa6usb+XMwVQAAAAAAjJwAC\nAAAAGDkBEAAAAMDICYAAAAAARk4ABAAAADBygw2AqupNVfWaQ9hvc1W1qtrnX0arqldX1e8eqv0B\nAAAAHC5r9s/AV9UHknyitfbavbY/J8mbk5zeWtu9tO0pSd7RWjt9tdbTWvvF1eoNAAAAsJrW8hVA\nb0vywqqqvba/KMk794Q/h8L+rvoBAAAAGIO1HABdkuSEJE/es6GqHpDk2UneXlVvraqfr6qjk7w/\nyalVddfS26lVNVNVr6qqa6rqq1X13qo6fqnPntu9fqCqbkjyl8v2+9Kquqmqbq6qn1q274ur6h2H\n4wsHAAAAOJTWbADUWrs3yXuTXLRs8/OTfLa19qlldXcnOT/JTa21Y5bebkryyiQXJvmOJKcmuT3J\nG/fazXckOTvJM5dte2qSM5M8I8nPVNXTD+kXdpi11rrfWBuqqvsNAAAYPo//14axPxdbswHQkrcl\neV5VbVp6/6KlbT1ekeTnWmtfbK3tSHLxUq/lt3td3Fq7eyls2uN1S9v+Pslbkvzbg/sSAAAAAKZr\nTb/2TWvtiqraluTCqvpkkscmeW7np5+R5E+qamHZtvkkD1z2/o37+Lzl265P8pgVLBkAAABgzVnr\nVwAlyduzeOXPC5Nc1lq7ZR81+7p/6cYk57fWjlv2tqm19qUJn/fgZf9+SJKbDnThAAAAAGvBUAKg\npyd5efZ/+9ctSU6oqmOXbXtTkl+oqjOSpKpOWvoT8pO8pqqOqqpHJ3lJkvcc+NIBAAAApm9N3wKW\nJK21rVX10STfkuTP9lPz2ap6d5Jrq2o2yaOS/HqSSnJ5VZ2a5CtZDHP+dMIuP5zkC1kMx17fWrv8\n0HwlAAAAANOx5gOgJGmtPWUf21681/sv3cenvmHpbe/P3ZrFcGh/2357H59zcd9qAQAAANaWIdwC\nBgAAAMBBEAABAAAAjJwACAAAAGDkBEAAAAAAIzeIF4HmwFXV5CIAAABg1FwBBAAAADByAiAAAACA\nkRMAAQAAAIycAAgAAABg5ARAAAAAACMnAAIAAAAYOQEQAAAAwMgJgAAAAABGTgAEAAAAMHICIAAA\nAICREwABAAAAjJwACAAAAGDkBEAAAAAAIycAAgAAABg5ARAAAADAyAmAAAAAAEZOAAQAAAAwcgIg\nAAAAgJETAAEAAACMnAAIAAAAYOTmpr2AdWJbkuunvYgpOWPaCzhAZjYs5jU863Vm5jU8Q5yZeQ2P\nmQ2LeQ2LeQ2Pma2Saq2tZn8AAAAApswtYAAAAAAjJwACAAAAGDkBEAAAAMDICYAAAAAARk4ABAAA\nADByAiAAAACAkRMAAQAAAIycAAgAAABg5ARAAAAAACP3/wGTk9cofOnEdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115a3d978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, all_axs = plt.subplots(4, prev_post.timesteps, figsize=(16, 4))\n",
    "all_axs = all_axs.T\n",
    "\n",
    "VISIBLES = 0\n",
    "TRUE_STATES = 1\n",
    "FILTERING = 2\n",
    "VITERBI = 3\n",
    "\n",
    "all_axs[0][VISIBLES].set_title('Visibles', x=-0.5, y=0.2)\n",
    "all_axs[0][TRUE_STATES].set_title('Actual', x=-0.5, y=0.4)\n",
    "all_axs[0][FILTERING].set_title('Filtering', x=-0.5, y=0.4)\n",
    "all_axs[0][VITERBI].set_title('Viterbi', x=-0.5, y=0.4)\n",
    "\n",
    "for i, (axs, hidden, visible, alpha, viterbi) in enumerate(zip(\n",
    "    all_axs, \n",
    "    prev_post.hiddens, \n",
    "    prev_post.visibles, \n",
    "    prev_post.alphas,\n",
    "    most_likely_states,\n",
    ")):\n",
    "    axs[VISIBLES].imshow([prev_post.map_visible_state_to_bump_creak[visible]], cmap='gray', vmin=0)\n",
    "    hide_ticks(axs[VISIBLES])    \n",
    "    \n",
    "    axs[TRUE_STATES].imshow(prev_post.plot_state_in_room(hidden), cmap='gray')\n",
    "    hide_ticks(axs[TRUE_STATES])\n",
    "    \n",
    "    axs[FILTERING].imshow(alpha.reshape(prev_post.height, prev_post.width))\n",
    "    hide_ticks(axs[FILTERING])      \n",
    "\n",
    "    axs[VITERBI].imshow(prev_post.plot_state_in_room(viterbi), cmap='gray')\n",
    "    hide_ticks(axs[VITERBI])   \n",
    "    \n",
    "maybe_save_plot('2018-05-13-viterbi')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## See Also\n",
    "\n",
    " - [This post](2018-05-13-viterbi-message-passing.ipynb) builds on this post to show the Viterbi algorithm.\n",
    " - [The more-detailed alpha recursion HMM post](2018-05-02-hmm-alpha-recursion.ipynb).\n",
    " - [This notebook](https://github.com/jessstringham/notebooks/blob/master/2018-05-13-hmm-check-results.ipynb) runs this code using the same example from Barber."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}