2018-05-02-hmm-alpha-recursion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inference in discrete state Hidden Markov Models using numpy\n",
    "\n",
    "This demo shows exact inference on a Hidden Markov Model with known, discrete transition and emission distributions that are fixed over time. It does alpha recursion, which is a sum-product algorithm on HMMs.\n",
    "\n",
    "The problem is taken from Example 23.3 from [\"Bayesian Reasoning and Machine Learning\"](http://www.cs.ucl.ac.uk/staff/d.barber/brml/), and this is more or less a port of the provided Matlab implementation with my commentary.\n",
    "\n",
    "I'm trying to understand how my probabilistic graphical modelling class can be represented in code, so I try to go through this painfully slow and show how the `numpy` arrays map to the discrete probability distributions.\n",
    "\n",
    "## Problem\n",
    "\n",
    "The idea is that there's someone moving around a room. You can't see them, but you can hear bumps and creaks from them moving around. You know the probability of someone bumping or creaking at each location in the room. You also know how one moves around the room.\n",
    "Using the sequence of bumps and/or creaks, the goal is to figure out the locations in the room.\n",
    "\n",
    "Below is what I will end up with:\n",
    "\n",
    "![Three sets of 10 images, second row showing a dot moving around a grid, third row showing highlighted area where model thinks the dot is](images/2018-05-02-filtering.png)\n",
    "\n",
    "The first row represents the bumps/creaks you hear. The second row represents where the person actually is. The third row represents where alpha recursion guesses the person is.\n",
    "\n",
    "### HMMs\n",
    "\n",
    "![HMM diagram](images/2018-05-02-hmm.png)\n",
    "\n",
    "This problem can be modeled using a Hidden Markov Model. At each timestep $t$ in the sequence, there's the visible state $v_t$ (the bump/creak combination) and the hidden state $h_t$ (the location in the room.)\n",
    "\n",
    "The goal of *filtering* is to show the probability distribution over the locations in the room at a given timestep $t$ given the bumps/creaks from the current and previous timesteps, $p(h_t\\mid v_{1:t})$. Plotting this distribution at each timestep gives a heatmap of where in the room we think the person is. \n",
    "\n",
    "In this model, there is a transition distribution $p(h_t \\mid h_{t - 1})$ to show that each hidden state depends on the previous timestep's hidden state. There also is the emission distribution $p(v_t \\mid h_t)$ to show each visible state depends on the corresponding hidden state. There's also the probability of where in the room the person starts, $p(h_1)$.\n",
    "\n",
    "The rest of the notebook shows:\n",
    " - Setting up the known distributions $p(h_1)$, $p(h_t \\mid h_{t - 1})$, and $p(v_t \\mid h_t)$. This is a lot of the notebook!\n",
    " - Generating some data for the locations ($h_{1:t}$) and sounds ($v_{1:t}$) from that distribution.\n",
    " - Finally I get to alpha recursion. I try to predict $h_{1:t}$ based on $v_{1:t}$, $p(h_1)$, $p(h_t \\mid h_{t - 1})$, and $p(v_t \\mid h_t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# helper functions you can skip over :D\n",
    "SAVE = False\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": [
    "## Distributions\n",
    "\n",
    "Before anything, I'll set up the room. I want a $6 \\times 5$* room with a total of 30 possible locations. \n",
    "\n",
    "In HMM, this means there are 30 hidden states, and at each timestep $t$, $h_t \\in \\{0,...,29 \\}$. Heads up, I'm using index-at-1 notation for the variables to match the book's notation and index-at-0 for states and code.\n",
    "\n",
    "In `numpy`, I'll set a `width` and `height` for the room and `num_hidden_states` as the total possible locations in the room. While I'm here, I'll set up `map_x_y_to_hidden_state_id` as a convenient way to map from an `x` and `y` coordinate to the index of the state. \n",
    "\n",
    "*The original example used a $5\\times5$ room, but I didn't want to accidentally transpose my math, so I made it a rectangular room instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "width = 6\n",
    "height = 5\n",
    "num_hidden_states = width * height"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial distribution (starting location) $p(h_1)$\n",
    "\n",
    "The initial distribution $p(h_1)$ gives the probability over the hidden states at the first timestep.\n",
    "\n",
    "In this case, it's the probability over locations in the room where the person is at the first timestep. The Barber book uses a uniform distribution for this problem. This means at the first timestep, they could be anywhere in the room.\n",
    "\n",
    "To represent $p(h_1)$ as a `numpy` array `p_hidden_start`:\n",
    " - The distribution is uniform, so all values in the array should be the same.\n",
    " - The distribution's domain is the possible hidden states, so the array should have the same size `(30,)`.*\n",
    " - The array represents a probability distribution so it should sum to 1 and each value should be greater than or equal to 0.\n",
    " \n",
    "*Technically this distribution can be specified with only 29 values because it has to sum to 1. But I'm not going to worry about that here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prob of starting starting locations\n",
    "p_hidden_start = np.ones(num_hidden_states) / num_hidden_states\n",
    "\n",
    "# verify it's a valid probability distribution\n",
    "assert np.all(np.isclose(np.sum(p_hidden_start), 1))\n",
    "assert np.all(p_hidden_start >= 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transition (moving around the room) $p(h_t \\mid h_{t-1})$\n",
    "\n",
    "Given the previous location, the transition distribution $p(h_t \\mid h_{t-1})$ gives the probability of moving to a new location. \n",
    "\n",
    "For this problem, the person can only move up, down, left, or right and they can't stand still. This can be encoded in a probability distribution, which is cool! I'll set the value to 0 when it's an illegal move and otherwise evenly distribute it among the legal moves.\n",
    "\n",
    "Using this I can define $p(h_t, h_{t - 1})$, which I can use to find $p(h_t \\mid h_{t-1})$. Translating to `numpy`:\n",
    " - This distribution defines probabilities for every combination of two hidden states, so it will have $30 \\times 30$ values*. To make things easier when I compute $p(h_t \\mid h_{t - 1})$, I'll use a `(30, 30)` array. I'll have `p_transition_joint[new_state][old_state]` as the probability $p(h_t=\\texttt{new_state}, h_{t - 1}=\\texttt{old_state})$.\n",
    " - It's a joint probability, so summing over the entire 2D array should give 1 and every value should be greater than or equal to 0.\n",
    "\n",
    "*Like with $p(h_1)$, technically this can be specified with $30 \\times 30 - 1$ values, because the entire table sums to 1. And like before, I'm not going to worry about that in code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_transition_joint(width, height):\n",
    "    num_hidden_states = width * height\n",
    "    \n",
    "    # begin by building an unnormalized matrix with 1s for all legal moves.\n",
    "    unnormalized_transition_joint = np.zeros((num_hidden_states, num_hidden_states))\n",
    "    \n",
    "    # This will help me map from height and width to the state\n",
    "    map_x_y_to_hidden_state_id = np.arange(num_hidden_states).reshape(height, width).T\n",
    "    \n",
    "    for x in range(width):\n",
    "        for y in range(height):\n",
    "            h_t = map_x_y_to_hidden_state_id[x, y]\n",
    "    \n",
    "            # hax to go through each possible direction\n",
    "            for d in range(4):\n",
    "                new_x = x\n",
    "                new_y = y\n",
    "                if d // 2 == 0:\n",
    "                    # move left or right!\n",
    "                    new_x = x + ((d % 2) * 2 - 1)\n",
    "                else:\n",
    "                    # move up or down!\n",
    "                    new_y = y + ((d % 2) * 2 - 1)\n",
    "    \n",
    "                # make sure they don't walk through walls\n",
    "                if any((\n",
    "                    new_x > width - 1,\n",
    "                    new_x < 0,\n",
    "                    new_y > height - 1,\n",
    "                    new_y < 0\n",
    "                )):\n",
    "                    continue\n",
    "    \n",
    "                h_t_minus_1 = map_x_y_to_hidden_state_id[new_x, new_y]\n",
    "                unnormalized_transition_joint[h_t_minus_1][h_t] = 1\n",
    "    \n",
    "    # normalize!\n",
    "    p_transition_joint = unnormalized_transition_joint / np.sum(unnormalized_transition_joint)\n",
    "    \n",
    "    # make sure this is a joint probability\n",
    "    assert np.isclose(np.sum(p_transition_joint), 1)\n",
    "    # not super necessary, but eh\n",
    "    assert np.all(p_transition_joint >= 0)\n",
    "\n",
    "    return p_transition_joint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing $p(h_t \\mid h_{t - 1})$ in math\n",
    "\n",
    "Now I have $p(h_t, h_{t - 1})$. To get the transition distribution $p(h_t \\mid h_{t - 1})$, from this, I need\n",
    "\n",
    "\\begin{align}\n",
    "p(h_t \\mid h_{t - 1}) &= \\frac{p(h_t, h_{t - 1})}{p(h_{t - 1})} \\\\\n",
    " &= \\frac{p(h_t, h_{t - 1})}{\\sum_{h_t} p(h_t, h_{t - 1})}.\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Translating to `numpy`:\n",
    " - The distribution maps from a hidden state to another hidden state, so it should be of size `(30, 30)` to represent each possible transition. The way Barber does it is by having a 2D array, where the first axis is the `new_state` $h_t$ and the second is the `old_state` $h_{t-1}$.*\n",
    " - The distribution is a conditional probability distribution, so given a value for $h_{t-1}$, the values for $h_t$ should sum to 1.\n",
    " \n",
    "*One more time because I need to know this for my exam next week! Technically this can be specified with $30 \\times (30 - 1)$ values because each of the 30 columns sums to 1. And like before I'm not going to worry about that in code. \n",
    "\n",
    "#### Computing $p(h_t \\mid h_{t - 1})$ in `numpy`\n",
    "\n",
    "The element `p_transition_joint[new_state][old_state]` is the probability \n",
    "\n",
    "$$p(h_t=\\texttt{new_state}, h_{t - 1}=\\texttt{old_state}).$$\n",
    "\n",
    "The column `p_transition_joint[:, old_state]` s the 30-element vector representing the unnormalized distribution $\\tilde{p}(h_t , h_{t - 1} = \\texttt{old_state})$. If I marginalize over all possible values of $h_t$, I end up with \n",
    "\n",
    "$$\\sum_{h_t} p(h_t, h_{t - 1}=\\texttt{old_state}) = p(h_{t - 1} = \\texttt{old_state}).$$\n",
    "\n",
    "This is like summing over the rows or columns of a probability table to get the marginal.\n",
    "\n",
    "| `p(x, y)` | `x = 0` | `x = 1` | `p(y)` |\n",
    "|--|--|--|--|\n",
    "| `y = 0` | 0.4 | 0.2 | **0.6** |\n",
    "| `y = 1` | 0.3 | 0.1 | **0.4** |\n",
    "| `p(x)` | **0.7** | **0.3** | |\n",
    "\n",
    "I'll build $p(h_t \\mid h_{t - 1})$ column-wise. Each column will be the normalized distribution ${p}(h_t \\mid h_{t - 1} = \\texttt{old_state})$.\n",
    " For each `old_state` I'll compute $p(h_t \\mid h_{t - 1}=\\texttt{old_state})$ as follows:\n",
    " - Sum over the column `old_state` in `p_transition_joint`. By doing this, I'm summing over $h_t$. This gives $p(h_{t - 1}=\\texttt{old_state}) = \\sum_{h_t} p(h_t, h_{t - 1}=\\texttt{old_state})$.\n",
    " - Divide the column `old_state` in `p_transition_joint` which represents $p(h_t, h_{t - 1}=\\texttt{old_state})$ by that sum $p(h_{t - 1}=\\texttt{old_state})$. That gives $p(h_t \\mid h_{t - 1}=\\texttt{old_state})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_transition(width, height):\n",
    "    p_transition_joint = create_transition_joint(width, height)\n",
    "    \n",
    "    num_hidden_states = width * height\n",
    "\n",
    "    p_transition = np.zeros((num_hidden_states, num_hidden_states))\n",
    "    \n",
    "    for old_state in range(num_hidden_states):\n",
    "        p_transition[:, old_state] = p_transition_joint[:, old_state] / np.sum(p_transition_joint[:, old_state])\n",
    "    \n",
    "    # verify it's a conditional distribution\n",
    "    assert np.all(np.sum(p_transition, axis=0)) == 1\n",
    "    \n",
    "    return p_transition\n",
    "    \n",
    "p_transition = create_transition(width, height)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualizing\n",
    "\n",
    "This one is weird to visualize but it makes a fun pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO4AAAD7CAYAAABt9agKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADGZJREFUeJzt3XusZWdZx/HvU4ZOy/SWYtVOKW1lLJ1pGQfRlsZgG4lO\nqBo0BG8VbY1KVf4QE6NEJUSrxqhpVEyARltDhOINI6IexToDBWnDZQq0SFq0ODItnXYo04tTaPv4\nx3pP2Z6evfecfeacvZ51vp9kJ/vsd6193vXO+u11eeY9OzITSbUcN+8OSFo5gysVZHClggyuVJDB\nlQoyuFJBgw9uRGyOiDsj4sx596VPImJPRFw+47oZEdvGtF0VEbesqnNrKCJ2RsSH5t2P1Rp8cIGf\nBt6fmffOuyPjrHRnj4hzW3g2rWW/5iEiboyIa1ew/IrGLjM/ATwUEd87Uwd7YiME9xrg7fPuhHrl\nz4HXzrsTq5KZpR/APcAbgDuBLwI3ACe0tucD/wtsGln+RuCPgfcCDwO3Ai8Yab8A+BfgEPAZ4Afa\n6+cBDwHHtZ+vB+4fWe/twM9P6etVwH+23/tfwJXAduAI8CTwCPBQW/a7gY8Dh4H9wJtG3ue/gWzL\nPwJc2l7/CeDTbRwWgHMm9GUPcPmE9p8C7m7j8HfA1pG2BLa1589t7YeB24DfAG6ZMg4BXAfc39b7\nJHAR3dnRV4Avt+16T1v+l4HPtnG7E/j+9vq4sdsM/F4bpy8AbwFOHPn9Z7X9YvO899+Z9/t5d2DV\nG9AF91PA2cDpwAeBa0d2/juWLH8j8CBwMbCJ7tP3pta2pYXk6tb2YuABYEd+NTAvac8/00K4faTt\nxRP6uaXtpC9sP58JXNieX7V0ZwcuB15Ed1a0s+2A39fazm3hGf1AemUL2vbW918FPjShP2ODC3xH\n2+5vbiH4I7rLjeWCexPwF237LgI+fxTB3Q18FDithXg7cObIv8+1S5Z/NbC1jcUPAo+OLL/c2F1H\n92FyOnAy8B7gt5cscxjYOe/9d9bHUE6V35yZ+zPzEPCbwA+310+j+5Re6t2ZeVtmPkEX3F3t9e8B\n7snMGzLzicz8OPDXdDsOwF7gsoj4+vbzX7WfzwNOAW6f0s+ngIsi4sTMvDcz7xi3YGbuycxPZuZT\n2V2XvRO4bMJ7X0O3c366bddvAbsi4pwpfVrOlcCfZubHMvNxujOaSyPi3NGFIuJZwKuAN2bmo5n5\nKeDPjuL9v0IXqAuAaH0eew8iM/8yMw+0sXgXcBfdB+8zRETQHblfn5mHMvNhurH4oSWLPky3f5Q0\nlODuH3n+ObpPZ+hOGU9eZvn7Rp4/BpzUnp8DXBIRDy0+6HbixaDupTsSfjvwfrqj1mXt8YHMfGpc\nBzPzUbqjxTXAvRHx3oi4YNzyEXFJRPxbRByMiC+19b5m3PKt738w0u9DdEezsyasM85WunFc7Psj\ndGcpS9/rDLqj+9LxnygzbwbeTHfJcn9EvC0iThm3fET8WETsG9m2ixg/FmcAzwE+OrL8P7XXR51M\nd+lT0lCCe/bI8+cDB9rzTwDnreDu635gb2aeNvI4KTN/prXvBV5GF969wC3At9EFd++0N8/Mhcz8\nTrrT5P+gu06G7tRzqXfQne6dnZmn0l2nxYTl9wOvXdL3EzNzltLHAboPAgAiYgvdteznlyx3EHiC\nZ47/VJn5h5n5EmAHcD7wi4tNo8u1M4brgdcBz83M0+gujcaNxQN0168XjozDqZl50sh7ngUcT3e5\nU9JQgvtzEfG8iDgd+BXgXQCZ+T90133LnlYt4++B8yPiNRHx7Pb41ojY3t7vLrqd4kfpAn6Y7trz\nVUwJbkR8XUS8soXgcbqbKYtH6C8Az4uI40dWORk4lJlHIuJi4EdG2g62db9h5LW3AG+IiAvb7zs1\nIl7NbN4JXB0RuyJiM92p5q2Zec/oQpn5JPA3wJsi4jkRsQP48Wlv3sb0koh4Nt316hH+/1iMbtcW\nunAebOteTXfEZWT5p8eunfVcD1wXEV/b1jkrInaPrHMZcHO7DChpKMF9B/DPdDeLPguM1gHfCrzm\naN6kXQ99F9310AG6U+rfobtBs2gv8GBm7h/5OYCPTXn744BfaO97iG7nWTyS3wzcAdwXEQ+0134W\n+PWIeBh4I90NoMV+PkZ3Lf/Bdjr40sx8d+vrTRFxmO6o9Iqj2e6lMvN9wK/RXd/fC7yAZ14jLnod\n3aXGfXQ3lm44il9xCl24vkh3av0g8Lut7U+AHW27/jYz7wR+H/h3upC+iO4G5KLlxu6X6D6wP9zG\n4n3AC0fWuZLug66saHfYyoqIe4CfbDvbcu2b6coqL590A2SjiYg9dCWmPXPuyrqKiJ3AWzPz0nn3\nZTUG9z9vlmqnQzvm3Q/1Q7tDXzq0sAGCu94i4pExTa/IzA+sa2cmu5GuBr4mIuJlwD8u1zZ6o0iz\nKX+qLG1EQ7k5JW0oBlcqaEXXuMfH5jyBLcu2PfmNm5d9fdGz7lr/ktmkPs2jP9I0R3iUL+fjMW25\nFV3jnhKn5yXx8mXbvvQPy86rftqpV9x91L/nWJnUp3n0R5rm1vxXDuehqcH1VFkqyOBKBRlcqSCD\nKxW0orvK5+98jIWFfcu27d667MtPm8eNoknve/d1Lx3btu31H16L7kjHjEdcqSCDKxVkcKWCDK5U\nkMGVCjK4UkEGVyqoF5MM+jYZwBqv5sVJBtKAGVypIIMrFWRwpYIMrlSQwZUKOmbloGlmLflUKhWB\n5SKtjuUgacAMrlSQwZUKMrhSQQZXKsjgSgWtWzlokllLPgsHlv+Lk4t2b901c59m5cwirYblIGnA\nDK5UkMGVCjK4UkEGVyrI4EoFrehLv9bKrLODpn3R2DxMKvn0baaT6vKIKxVkcKWCDK5UkMGVCjK4\nUkEGVyrI4EoFrWha37d80wl528LZy7bNYwrdavRt+p01XoHT+qRBM7hSQQZXKsjgSgUZXKkggysV\ndMz+yuOkcgbMp6Qxa4mlb6Uiv2hs47AcJA2YwZUKMrhSQQZXKsjgSgUZXKmgdfvSr77Nfpm15NPH\n0kzfyleaneUgacAMrlSQwZUKMrhSQQZXKsjgSgWtWzlokqGUila77lroW380meUgacAMrlSQwZUK\nMrhSQQZXKsjgSgUZXKmgXtRxJ6lU44XJtdEhbYvWhnVcacAMrlSQwZUKMrhSQQZXKsjgSgX1vhw0\nycKBfRPbd2/dtU49+aohTaPrW/lqI7AcJA2YwZUKMrhSQQZXKsjgSgUZXKmg0uWgaYZSzujjLJ6h\njG3fWA6SBszgSgUZXKkggysVZHClggyuVNCmeXdgLU0qS/RtFs+k8sq2K/r3RWOVxnaIPOJKBRlc\nqSCDKxVkcKWCDK5UkMGVCjK4UkGDntY3q77VIVczrW9I27IROK1PGjCDKxVkcKWCDK5UkMGVCjK4\nUkGWg1aoj+WMWUs+Q9qWobAcJA2YwZUKMrhSQQZXKsjgSgUZXKkgy0HHWN/KGavpT9++2KtvY7sW\nLAdJA2ZwpYIMrlSQwZUKMrhSQQZXKsjgSgVZx11HleqiUKs22rexnZV1XGnADK5UkMGVCjK4UkEG\nVyrI4EoFWQ7qiaGUM6B/0+8WDuyb2L5766516sl0loOkATO4UkEGVyrI4EoFGVypIIMrFbRp3h1Q\nZ1LJp4+zeCaVr7Zd0a8vGptW7qlYivOIKxVkcKWCDK5UkMGVCjK4UkEGVyrI2UED0LfZOKvpT9+2\nZb1LRc4OkgbM4EoFGVypIIMrFWRwpYIMrlSQwZUKso47cH2ri65mWl/ftmUtarzWcaUBM7hSQQZX\nKsjgSgUZXKkggysVZDloA+vjX4+cteTTt7/UOKk/ML5PloOkATO4UkEGVyrI4EoFGVypIIMrFWQ5\nSGP1rcTSt9lBq7FwYN+yr1+8ez8fuf2I5SBpiAyuVJDBlQoyuFJBBlcqyOBKBVkO0kz6VipajT5t\ni7ODpAEzuFJBBlcqyOBKBRlcqSCDKxVkcKWCNs27A6ppUn2zb9Pvxk2hW7R76/i2PtV4R3nElQoy\nuFJBBlcqyOBKBRlcqSCDKxXktD6tqz5+0disJZ+1KBU5rU8aMIMrFWRwpYIMrlSQwZUKMrhSQZaD\n1Ct9m1m0FqWiSetaDpIGzOBKBRlcqSCDKxVkcKWCDK5UkOUglTGUUhGM/wN2F+/ez0duP2I5SBoi\ngysVZHClggyuVJDBlQoyuFJBBlcqyDquBqFvX87ltD5Jz2BwpYIMrlSQwZUKMrhSQQZXKmhF5aCI\nOAh8bu26I21452TmGdMWWlFwJfWDp8pSQQZXKsjgSgUZXKkggysVZHClggyuVJDBlQoyuFJB/weQ\nWtpLkm61gwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f92b550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(p_transition)\n",
    "ax.set_title('p(new_state | old_state)')\n",
    "hide_ticks(ax)\n",
    "maybe_save_plot('2018-05-02-transition')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sample paths\n",
    "\n",
    "I can also generate a few path examples to see if they look right. Later, after I have the emission distribution, I can use this code to help generate the visibles.\n",
    "\n",
    "The `np.random.choice` comes in handy a few times. The `p` argument allows me to specify a list of probabilities for each state. I use it to sample from the `p_hidden_start` distribution and the `p_transition` distribution for a given hidden state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_state_in_room(state_id, width=width, height=height):\n",
    "    h = np.zeros(width * height)\n",
    "    h[state_id] = 1\n",
    "    return h.reshape(height, width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAGOCAYAAACUpwdsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC3FJREFUeJzt3E1u28oWhdFTFx5C3HbmIM1/BNYc4r7mULfhFwSvEZK6\n+iGr9lqAW2LsOi4y+CDB1XrvBQCQ6p+9FwAAsCcxBABEE0MAQDQxBABEE0MAQDQxBABEE0MAQDQx\nBABEE0MAQLS3Wy5urQ19XHXvvS29Pvp8VXXtvb8vXTD6jGt7WGXGEXgWx5/Rffpt9hlHn682PItV\n3hmazdfeCwCqyrMIR7HpWRRDAEA0MQQARBNDAEA0MQQARBNDAEA0MQQARBNDAEA0MQQARBNDAEA0\nMQQARBNDAEA0MQQARBNDAEA0MQQARHvbewEAHE/v/a+vnc/nF66Ee4y+j0vrr6pqrT3k53hnCACI\nJoYAgGhiCACIJoYAgGhiCACIJoYAgGhiCACI5pwh4GavOvvjqBLmn2GGNaOfwbPF6Pv4qvV7ZwgA\niCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiDbVOUOjnxmRcHbJmtH3sCpjH2eY4R7p88/CPvKb\nd4YAgGhiCACIJoYAgGhiCACIJoYAgGhiCACIJoYAgGi3njN0raqvZyzkEVbOjPjY8C12ne8BZ14c\nfsY1D9jDqrH3cYgZ7zT8fbrB7DO6T/+YfcaR56vauI9t7YA4AICZ+ZgMAIgmhgCAaGIIAIgmhgCA\naGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGII\nAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgm\nhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgmhgCA\naGIIAIj2dsvFrbX+rIVscTqdFl+/XC6Lr/fe29Lrz57v3vVvcO29vy9dcO+ML5hh0doeVq3PuPcM\na7bM+OPHj/7z58+/vr73DGv2fhZf4OnP4t4e8SwenRnHn682PItVVa337XPu/UtZW2try/ft3pt+\n7/o3uPTezys/464ZXzDD2s+/+z+nvWdYs2XG8/ncPz8///r63jOs2ftZfIGnP4t7EwrfZp9x9Plq\nw7NY5WMyACCcGAIAookhACCaGAIAookhACDaUDHUWlv8OrrR119lhqO4XC7Dz7DkdDpV7/2vXwCP\nNFQMAQA8mhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCDmf2c5SqnKUERyKGAIBoYggAiCaG\nAIBoYggAiCaGAIBoYggAiCaGAIBoYghgB85SyjhLaWn+EX4HKXsohgCAaGIIAIgmhgCAaGIIAIgm\nhgCAaGIIAIgmhgCAaG97LwCAOf0+SynZ6POn7KF3hgCAaGIIAIgmhgCAaGIIAIgmhgCAaGIIAIgm\nhgCAaLeeM3Stqq9nLOQFPjZcM/J8VfPPuGW+KjMe3ez3adX8M7pP/5h9xpHnq9q4j633/uyFAAAc\nlo/JAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaG\nAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBo\nYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggAiCaGAIBoYggA\niCaGAIBoYggAiCaGAIBoYggAiCaGAIBob7dc3Frrz1rIK/Te29Lro89XVdfe+/vSBaPPuLaHVWYc\ngWdx/Bndp99mn3H0+WrDs1jlnaHZfO29AKCqPItwFJueRTEEAEQTQwBANDEEAEQTQwBANDEEAEQT\nQwBANDEEAEQTQwBANDEEAEQTQwBANDEEAEQTQwBANDEEAEQTQwBANDEEAEQTQwBANDEEAEQTQwBA\nNDEEAEQTQwBANDEEAEQTQwBAtLe9F8B2vffF11trL1oJ6dLvxfT5Z7G0j+fz+YUr+e9GvxePsn7v\nDAEA0cQQABBNDAEA0cQQABBNDAEA0cQQABBNDAEA0aY6Z2iGMyOWHP28iEeYYQ+Pcm7GM80wwz0S\n5p/hWVwzwz6OPsNR1u+dIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAg2q3nDF2r6usZC3mE\nlfMKPjZ8i0PPt8HwMz5gD6t2nvHOczOGmPFOw9+nGww/4wzP4p3M+G3k+ao27mNbOyAOAGBmPiYD\nAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJ\nIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAg\nmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgCAKKJIQAgmhgC\nAKKJIQAgmhgCAKKJIQAgmhgCAKK93XJxa60vvX46nRb//eVyueXHPVzvvS29vjbfAK699/elC0af\ncW0Pq8w4As/i+DO6T7/NPuPo89WGZ7Hqxhha8/n5ufh6a6v3Fff52nsBQFV5FuEoNj2LPiYDAKKJ\nIQAgmhgCAKKJIQAgmhgCAKI99K/JZv9rsd6X/8Jw9vkBYEbeGQIAookhACCaGAIAookhACCaGAIA\nookhACCaGAIAoomhG7TWFr9m0Htf/BrB6XQafgYAXkcMAQDRxBAAEE0MAQDRxBAAEE0MAQDRxBAA\nEE0MAQDRxBD/Z4azlC6Xy/AzrHGW0vjs4Rznmt3L7+AYxBAAEE0MAQDRxBAAEE0MAQDRxBAAEE0M\nAQDRxBAAEE0MwYBmP0sp4Qye2fdwixnONbvX7L+DUc5REkMAQDQxBABEE0MAQDQxBABEE0MAQDQx\nBABEE0MAQLS3G6+/VtXXMxbyAh8brhl5vqr5Z9wyX5UZj251xsvlcm2tjTpflWfxNzMe29Pv0wOc\nlbRpH9uRDj0CAHg1H5MBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQ\nTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwB\nANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHE\nEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHEEAAQTQwBANHebrm4tdaftZAtTqfT4uuXy2Xx9d57\nW3p97/ke4Np7f1+6YO8Zn72HVfvPeC8zjj9fDfAs3muG+3Tp/6Nfv37V9XodfsY1nsVvrfftc+79\nS1lba2vL923Apl967+elC/ae8dl7+L/vMfQ+mnH8+WqAZ/FeM9ynS/8fnc/n+vz8HH7GNZ7Fbz4m\nAwCiiSEAIJoYAgCiiSEAIJoYAgCi3fSn9Xtb+0uj0d37l1YjmGGGNTPs4wwz3CN9/qr1v7SagX08\n/j6+6ln0zhAAEE0MAQDRxBAAEE0MAQDRxBAAEE0MAQDRxBAAEG2oc4Zml3DmRYIZ9nGGGe6RPn+V\n38EsRt/HV63fO0MAQDQxBABEE0MAQDQxBABEE0MAQDQxBABEE0MAQDQxBABEE0MAQDQxBABEE0MA\nQDQxBABEE0MAQDQxBABEE0MAQDQxBABEE0MAQDQxBABEE0MAQDQxBABEE0MAQDQxBABEE0MAQLS3\nG6+/VtXXMxbyAh8brhl5vqr5Z9wyX5UZj272+7Rq/hndp3/MPuPI81Vt3MfWe3/2QgAADsvHZABA\nNDEEAEQTQwBANDEEAEQTQwBANDEEAEQTQwBANDEEAEQTQwBAtH8Bsvp7Pvrp7TwAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fa76588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# code for plots\n",
    "timesteps = 10\n",
    "example_count = 5\n",
    "\n",
    "figs, all_axs = plt.subplots(example_count, timesteps, figsize=(timesteps, 1.5 * example_count))\n",
    "\n",
    "for ex_i in range(example_count):\n",
    "    # choose the starting state using `p_hidden_start`\n",
    "    curr_hidden_state = np.random.choice(num_hidden_states, p=p_hidden_start)\n",
    "    \n",
    "    # plot it\n",
    "    all_axs[ex_i][0].imshow(plot_state_in_room(curr_hidden_state), cmap='gray')\n",
    "    hide_ticks(all_axs[ex_i][0])\n",
    "\n",
    "    for t in range(1, timesteps):\n",
    "        # use the transition matrix to get the probability distribution of the\n",
    "        # next state given the current one\n",
    "        curr_hidden_state = np.random.choice(\n",
    "            num_hidden_states,\n",
    "            p=p_transition[:, curr_hidden_state]\n",
    "        )\n",
    "        \n",
    "        all_axs[ex_i][t].imshow(plot_state_in_room(curr_hidden_state), cmap='gray')\n",
    "        hide_ticks(all_axs[ex_i][t])\n",
    "maybe_save_plot('2018-05-02-sample-paths')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Emission distribution (bumps and creaks) $p(v_t \\mid h_t)$\n",
    "\n",
    "Now I need the emission distribution $p(v_t \\mid h_t)$, which is the probability of a visible state given the hidden state.\n",
    "\n",
    "In this problem, this distribution is the combination of the probability of a bump given the location, and the probability of a creak given the location. \n",
    "\n",
    "$$p(v_t \\mid h_t) = p(v_t^{bump} \\mid h_t)p(v_t^{creak} \\mid h_t)$$\n",
    "\n",
    "Both of these have two states, so combining the two gives 4 visible states for $v_t$:\n",
    " - bump and creak\n",
    " - bump and no creak\n",
    " - no bump and creak\n",
    " - no bump and no creak"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bumps and Creaks map, $p(v^{bump}=True \\mid h_t)$\n",
    "\n",
    "First I'll set up the bumpy and creaky aspects of the room. For example, `prob_bump_true_given_location` gives the probability distribution\n",
    "\n",
    "$$p(v^{bump}=True \\mid h_t)$$\n",
    "\n",
    " - There are `num_hidden_states = 30` values for $h_t$, so I need 2 sets of 30 probabilities. \n",
    " - This isn't a probability distribution. They just need to be between 0 and 1.\n",
    "\n",
    "I try to implement this in a similar way to Barber. Ten random spots in the room have a higher probability, and the rest has a low probability. This introduces some noise to the problem.\n",
    "Here I use `np.random.choice` in a new way. In this case, I'm using it to select 10 spots to have a high probability of making a noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_sound_map():\n",
    "    NUM_SOUNDS = 10\n",
    "    LOW_PROB = 0.1\n",
    "    HIGH_PROB = 0.9\n",
    "\n",
    "    # everything has at least LOW_PROB of triggering the sound\n",
    "    grid = LOW_PROB * np.ones(num_hidden_states)\n",
    "    # select NUM_BUMP_CREAKS to make HIGH_PROB\n",
    "    locs = np.random.choice(\n",
    "        num_hidden_states, \n",
    "        size=NUM_SOUNDS, \n",
    "        replace=False\n",
    "    )\n",
    "    grid[locs] = HIGH_PROB\n",
    "    \n",
    "    return grid\n",
    "\n",
    "prob_bump_true_given_location = make_sound_map()\n",
    "prob_creak_true_given_location = make_sound_map()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing bumps and creaks\n",
    "\n",
    "Let's plot it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAACgCAYAAAA7FRY2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACxRJREFUeJzt3HuMXGUdxvHngUJLqaYpRaVQWrzUAFoqpEVuAcNlA4rE\nIF5ajIAV8RL8B/CGCoJBE0SNiCgJRSgoRLyDLq2koAIKSilRtEql3VKQ0lJS2nLtzz/O2Xo63d2Z\n2d2Z+W37/SST7My8Z973nfnNM+e8Z2YdEQIAdN5OnR4AAKBAIANAEgQyACRBIANAEgQyACRBIANA\nEiMykG2Ptv1323uV16+zfWmnx5WV7eNsL6xcn2/7on7a7mb7n7b3aNsAsUVtbbex32Nsr2xnn0Nl\n+1LbF1aur7R9TD9tD7b9+7YNbpBGZCBLOlvS3RHxRKcH0izb+9p+rnIJ2xsq14/q5PgiYpOkH0q6\noJPj2IGNyNq2PadSw5tsb67WeafHFxF/lbTJ9omdHstARmognyPphk4PYjAiYkVEjOu9lDcfVLlt\nm09x2zu3eZg3SjrT9i5t7heDrG3bo1owloZFxI2Vmj5R0qo+6nyLDo33Rkkf60C/DUsbyLYfs/25\n8vDtGdvzbI+xva+k10v6U80mE20vsL3e9l22p5SPM7XcCx1VeexFtueWf59h+4+2v2l7ne1ltg8v\nb++x/ZTtD1e2vc721X31NYxzn2/7u7Z/a3uDpKNs/8H2GZU2c20vqlw/wPZC22tt/8P2qXW6mWD7\nN+Uc7rW9X+8dEbFc0gZJs4ZzXig0U9vlEtI3bC+3/WxZB7tV6vojtldIurNs/3bb95S1/FD1EN72\nmbYfKV/zZbb7DSfb55bj22cY573S9vm2H5a0wfaocg5TK222Wk6z/e5yHuvKub+lTjcH2364fK5+\nZHt05b5Fko7PvKORNpBLcyR1SXqDpGmSLpT0VknLIuLlPtpeImmipMUqPg0bdaikJZL2kHSTpB9L\nminpjZJOl3Sl7eqnfL992V5SFk9fl6uaGNNsSRdLepWkewdqWI5tgaTrJb2mHN8PbL+5zuN/UdIE\nSSvK+VQ9IumgJsaL5jRa25dLOkTS4Speqwskba7cf7Sk/SV12d5b0m2SLi3bnifpVtt7lm2fkvQu\nSa+WdKakb9o+uHZgtr8k6QxJR0fEShfLbP3V9Drbs5uY9wdU7EGPr9fQ9kxJ10iaq+K9ea2kX9je\ndYDN3ifpeBUfbIdI+lDvHeWOhiW9qYnxtlX2QL4yInoiYq2kr0r6oIoXcn0fbW+LiLsj4gVJX5B0\nmO3JDfbzn4iYFxGvSLpZ0mRJX4mIFyLiDkkvqgjnun1FxPSIGN/P5RNNzP1nEXFvRGwu+xnIKZKW\nRsT1EfFyRPxF0s8lvXeAbX4SEQ9ExEsqPlBm1Ny/Xg28aTBodWvb9k6SzpL06Yh4PCJeiYh7aurh\noojYUK79ny7p9oi4vaybBZIekHSSJEXEbRHxaBTuknSHpOo5C9u+QtIJkt4REavL7VYMUNPjI+Km\nJub97YhYWY63nrMlXRUR95dzv7a8feYA23wrIp6MiDWSfq0RVtfZA7mn8vdySZMkPaNir7HfthHx\nnKS1ZftG/Lfy96byMWpvq+4hD6WvRvXUb7LFFElHVPdaJL1f0kBn6p+s/L1RW89PKp7jdU2MAc1p\npLYnShoj6dEGH2eKpNNq6uBIlXVg+0Tb95XLWutUBPXEyvbjVYTgZRHx7CDnVU+zdf2ZmvnsJWnv\nAbYZ0XWdPZCre7j7SlqlYmlhP297UmBL2/IQfkLZfkN589hK29cN17hq+pLtv3nrb1FUL1c30Uft\nv+HboP7n0CPpdzV7LeMi4lPNTWsr+0t6aAjbY2CN1PbTkp5XsazRn2qd9Ei6oaYOdo+Ir5Vrqbeq\nWAJ5bUSMl3S7ikP4Xs+oWNKYZ/uI3hu97TeDai9zmpj3lvGWSzMvaOC6vrhmPmMj4pYm+tvC/z/X\n86/BbN8O2QP5k7b3sT1BxdLAzRGxUtK/te0Jp5NsH1muL10i6b7ykHC1pMclnW57Z9tnaeACb0Sf\nfUlSRBxYPbtcczlnCH0ulnRqeUJnmopD2V6/lHSg7dm2dykvs+qsIffLxcmlcZLuH8J4MbC6tR0R\nm1Wsm15he1JZv4fVnKiqmi/pZNtdZdsxLr5fvI+kXSWNlrRa0ssuvv51Qu0DRMQiFevbP7XdO44V\nA9T0uIho5nxNrYckzSnH+04Ve/S9rimfp5kujLN9su3dB9nX0ZIWlst0KWUP5JtUrHMtU3HY1vvj\nj++rslhfaftlFcsHh6hYT+v1UUnnS1oj6UBJ9wzDuPrrq1UuV7F38ZSKN+n83jvKw8uuchxPqDhs\nu0zFG3Aw5kiaFxEvDmXAGFCjtX2epIdVfDiulfR19fO+LXcKTpH0eRXB26Oi7neKiPWSzpV0i4o9\n4dkqPsj7epwFKj7wf9XXSb9hdq6k96hYRjitOqaIuE/SxyV9rxzzUg3tvTZHUjNHqW3nSPoP6m0/\nJmluRCzs477Rkh6UdGy0+Qv0tq+TtDIiLqzXNgvbx0n6bEQc10Db3VTsjR8REU+3fHA7oKy1PdK4\n+HXu8xFR91e6tt8m6TsRcWS9tp3U0S+TD1Z5lvmATo9je1Se/R7UUgeGjtpujYh4UFsvh6Q0IgMZ\nTVum4jvKwPbkTkm1v0cY0dIuWQDAjib7ST0A2GE0tWSxq0fHGA32GycYyaZN39jyPh7reUlPr33F\n9VsOr4kTdo6pk1v77w2WLhlbv9EwaMfr1K65tEO2um4qkMdodx3qYwc3Koxo3d2LW97HrK5mfsQ1\nfKZO3kV/7m70V/aD0zWp9he8rdGO16ldc2mHbHXNkgUAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0AS\nBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIA\nJEEgA0ASBDIAJEEgA0ASBDIAJDGq0wPolO5Vi1veR9ekGS3vA0O3dMnY7ea12l7m0S7teL6WxpqG\n27KHDABJEMgAkASBDABJEMgAkASBDABJEMgAkASBDABJEMgAkASBDABJEMgAkASBDABJEMgAkASB\nDABJEMgAkASBDABJEMgAkASBDABJEMgAkASBDABJEMgAkASBDABJEMgAkASBDABJEMgAkASBDABJ\njOr0AIBOmzZ9o7q7F7e0j65JM1r6+O3Uvaq1z1U7ZXtd2EMGgCQIZABIgkAGgCQIZABIgkAGgCQI\nZABIgkAGgCQIZABIgkAGgCQIZABIgkAGgCQIZABIgkAGgCQIZABIgkAGgCQIZABIgkAGgCQIZABI\ngkAGgCQIZABIgkAGgCQIZABIgkAGgCQIZABIYlQzjadN36ju7sWtGgsS65o0o+V9LI01Le+jz36X\njG3L/LYX7XquulfteFnDHjIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEg\nA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0ASBDIAJEEgA0AS\nBDIAJEEgA0ASBDIAJDGq0wOo1TVpRlv66V61uC39AO20PdV1O7KgHc/XrK6NDbdlDxkAkiCQASAJ\nAhkAkiCQASAJAhkAkiCQASAJAhkAkiCQASAJAhkAkiCQASAJAhkAkiCQASAJAhkAkiCQASAJAhkA\nkiCQASAJAhkAkiCQASAJAhkAkiCQASAJAhkAkiCQASAJAhkAkiCQASAJR0Tjje3Vkpa3bjjYwU2J\niD3b3Sl1jRZruK6bCmQAQOuwZAEASRDIAJAEgQwASRDIAJAEgQwASRDIAJAEgQwASRDIAJAEgQwA\nSfwPusS1xeigWYEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111496438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (bump, creak) = plt.subplots(1, 2)\n",
    "\n",
    "bump.imshow(prob_bump_true_given_location.reshape(height, width))\n",
    "bump.set_title('p(bump=True|h)')\n",
    "hide_ticks(bump)\n",
    "\n",
    "creak.imshow(prob_creak_true_given_location.reshape(height, width))\n",
    "creak.set_title('p(creak=True|h)')\n",
    "hide_ticks(creak)\n",
    "\n",
    "maybe_save_plot('2018-05-02-bumps-creaks')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Emission distribution $p(v \\mid h)$\n",
    "\n",
    "The goal is to define the distribution $p(v \\mid h)$. In `numpy`\n",
    " - The distribution gives probability of the 4 visible states given each of the 30 hidden states, so the array should be of size `(4, 30)`.\n",
    " - The distribution is conditional over the 4 states, so each column should add to 1.\n",
    "\n",
    "To get the emission distribution in math, I'd multiply \n",
    "\n",
    "$$p(v \\mid h) = p(v^{bump} \\mid h)p(v^{creak} \\mid h).$$\n",
    "\n",
    "It's a little more complicated in code. First to get $p(v^{bump} \\mid h_t)$,  I need to add a second state for $p(v^{bump}=False \\mid h_t) = 1 - p(v^{bump}=True \\mid h_t)$.\n",
    "But I also need to flatten out the visible states from 2 2-state visible variables to one 4-state visible variable. So like Barber, I'll define the four states' probabilities by hand. \n",
    "\n",
    "I also define a `map_visible_state_to_bump_creak` to map from the visible state to if there was a bump and/or creak."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_visible_states = 4\n",
    "\n",
    "def get_emission_matrix(prob_bump_true_given_location, prob_creak_true_given_location):\n",
    "    # prob_bump_given_state[v][state] = p(v | state)\n",
    "    p_emission = np.vstack((\n",
    "        prob_bump_true_given_location * prob_creak_true_given_location,\n",
    "        prob_bump_true_given_location * (1 - prob_creak_true_given_location),\n",
    "        (1 - prob_bump_true_given_location) * prob_creak_true_given_location,\n",
    "        (1 - prob_bump_true_given_location) * (1 - prob_creak_true_given_location),\n",
    "    ))\n",
    "    \n",
    "    assert np.all(np.sum(p_emission, axis=0)) == 1\n",
    "    \n",
    "    return p_emission\n",
    "\n",
    "p_emission = get_emission_matrix(prob_bump_true_given_location, prob_creak_true_given_location)\n",
    "    \n",
    "# 1 means True. ex: [1, 0] means bump=True, creak=False\n",
    "map_visible_state_to_bump_creak = np.vstack((\n",
    "    [1, 1],\n",
    "    [1, 0],\n",
    "    [0, 1],\n",
    "    [0, 0],\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing $p(\\textbf{v} \\mid h_t)$\n",
    "\n",
    "This one's a little tricker to visualize. I'll show it as four plots of the probability each visible state.\n",
    "\n",
    "When I do alpha recursion the data below, the result of filtering for the first timestep will be exactly the same as one of these maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAD7CAYAAADn/2ZiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFJ9JREFUeJzt3Xu0ZGV95vHvQ99osLk0InJpaC5ClKg9UUCR0Tg4C4TM\nGicuNYqXHgPGWRkva7yF6CgrahINE3RCoolRCQbEUSLLQbAzxKBBAhgSZBwxINBNd0BCN93cwy3v\n/LH3wU111XvqnN6Hc7rP97PWWV1V+633fXedXz21967dZ6eUgiRpuJ1mewKSNJcZkpJUYUhKUoUh\nKUkVhqQkVRiSklTxlIdkkiVJfpRk3/b+OUk+9lTPY3uX5Iwkfz5i2TuSfOKpntOOypqde5JcnuTU\nEcsuTPLKvsaajS3JtwHfLaXcMQtjT0n75jg3yeYk/5zk07M9pzF9DjglyTNmeyI7CGt2+/IJoLcP\nsdkIybcDX5qFcadjNfALwCHAwcBFfXWcZGFffQ0qpfwLcCnw5pkaY56xZnuWZMFM9V1KuQbYLckL\n++iv95BMsjbJ6e3uyeYkX0yyc7vsQJpf3tUDT3t6kv+T5L4k30lyUNt+ZZLSDZTuZnaS1Um+l+Ss\nJFuS3JLk2Pbx9e0n6Vs6zz0nyWeHjTXCo8A9pZTNpZQHSil/Pcb6r0jyF0nuSrIpydlD5roJOKN9\n/K1JbmhfqzXd+ST5dLse9ya5Nsm/HTHmoiRfbnczFrcPXw6cPNl8Zc226//eJNcnuSfJVybWv11+\nWpKfJLk7yTeS7Ffp67gkV7brtj7J6s56fCbJJUkeAF7ebvWemeS2JHe267m0bb9nkovb99Hm9vYB\nI8bct537+zoPX05P9T9TW5KnACcAhwKHAx9qH38ucEsp5bEh7T8KPB24DjhvCmMdA1wP7AWcD1wA\nHAUcBrwRODvJ06Y51rXAi5L81jgTSfPpeDGwDlgJ7N/OpzvXW4B9gI8n+Y/AbwK/DOwN/A3w5U77\n7wOrgOXtun21W7ztmEtpthYeBl5bSnmkXXQD8Pxx5i1gntZsx2uBE2m2Pp9Hs0VKkn8H/E67fF+a\n2r5gWAdteF8K/AFNPa9q5zvhDcDHgWXAFcDv0rzWq2jWfX/gw23bnYAvAgcBBwIPAWcPGfNg4DvA\n2aWU3+ss6q/+Sym9/gBrgbd37p8E3NzePgW4aqD9OcAFnftPAx4HVtAETQEWdpZfDpza3l4N3NRZ\n9ty2/T6dxzYBqyYba8h6LAduoymcq4EzOss2AM8d8pwXA3d159tZthq4beCxS4Ff7dzfCXgQOGjE\na7sZeH57+wzgGzQF8j+BDLR9FvB437/fHfFnPtdsZ/3f2Ln/SeCz7e3PA58cGP9RYOWQfk4Hvj5i\njHOAczv3AzwAHNp57MXArSOevwrYPPCa/n4799cPaX8a8O0+6mOmjout79xeB0xsnm+m+RQZ2b6U\ncn+Su9vn3DnGWN02D7V9DD7W/VQeNVZ3zgCvAW4opXwryfeB7yaB5pe9EPjhkLmsANaVrbc6thq7\ndRDw6ST/o/NYaD5R1yV5L/Cr7fwKsBvN1sSEFwGLaIpk8C+VLAPuGTEPbW2+1uyEn3ZuP8jP1n8/\n4O8Hxt9EU6NrB/pYAdxcGaM7372BXYBr2zlCU/sLAJLsApxFE/h7tsuXJVlQSnm8vX8K8BPga0PG\nWgZsqcxlbDO1u72ic/tA4Pb29vXAwdn6S4sn2re7Gcvb5zzQPrxLp+0z+5rbwFiDFtIEEKWUTcC/\nB94CrAHOHBJK0BTBgUPWb8Lgc9YDv1ZK2aPzs7SUcmV7/PH9NLs5e5ZS9qAJvXSe/5c0u0J/lWSf\ngb6fDfxgxDy0tflas5O5nebDfGL8XWkOE/zTkLbraQ5XjNIdfyPNh8GRndrfvZQy8eHwHuAI4JhS\nym7ASyem0OnjjLaf87P1F0G91f9MheSvJzkgyXLgg8BXAEopG2iS/+iB9ie1B3wX0xx7uaqUsr6U\nchfNL+ONSRYkeSv1X8I4ho41pN0lwFFJfi3JIppdjCtpjqE8OKLva4A7gN9NsmuSnZO8pDKXzwKn\nJzkSIMnuSV7TLlsGPEa7+57kwzRbkk9SSvkkzXGtv0rS3cp8Gc3uvMYzX2t2Ml8G/nOSVUmWAL8N\nXF1KWTuk7XnAK5K8NsnCJHslWTWs01LKv9KcqnZW2lPVkuyf5IS2yTKaEN3S/k4+MqSbR2m2nncF\nzk3SzbPe6n+mQvJ8mq2cW2g2v7vnLP0x8KYh7T8C3A28gObg9YTTgPfRHKc5kuaXvq1zGzXWE0op\ntwKvpDmNZhPNp9KdwMuBTyQ5cchzHgf+A81B6NtojgO9btRESilfpzmn64Ik99LsDk2cBLsG+BZw\nI83u37+w9e7VRD8fpfny5rIky9svd04C/mzU2NrKvKzZyZRSLgP+O3AhzQbAocCvjGh7G03dvaed\n63XUvzz5AM0H0FVt/V9Gs/UI8ClgKc2W4lU074VhYz5C88XnPsAXkuyU5Cjg/tKcCrTNMr0t8EqH\nyVqag9SXjVi+BPgH4PjyFJ+cm+QcYEMp5UOTtd2eJXkHzYH998/2XLYH1uyOJcmFwOdLKZf00d+M\nndA8SinlYeA5T/W480kp5Q9mew47Emt2+1JKeXWf/fkHLiSpovfdbUnakbglKUkVhqQkVUzpi5vF\nC5aWpYt2723w8vAjkzeagixZPHmjMfU9t74d/rzpnva2tbXrH2Xj3Y9n8pY7psVZUnZm197667MO\nod9a7Htuc929D9+5sZSy97b0MaWQXLpod449YPB0sel77Ja1vfUFsPCAlb311ffc+rZmzXWTNxrT\n0ScMPf1y3tiZXTkmx/fWX591CP3WYt9zm+u+dfOZ67a1D3e3JanCkJSkCkNSkioMSUmqMCQlqcKQ\nlKQKQ1KSKgxJSaowJCWpwpCUpApDUpIqDElJqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQlKSKKV2+\noTz8yJy/rMF8cfJxr+qtr5s2fKm3vrZHWbJ4Tl/WYM3t/V2q4+TjVvbW10z45hUX9drfgn23vQ+3\nJCWpwpCUpApDUpIqDElJqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQlKQKQ1KSKgxJSaowJCWpwpCU\npApDUpIqDElJqjAkJanCkJSkCkNSkiqmdI2bvq8FMpevl7PwkJWzPQU9Rfq+dtNcrp25/J6bq9yS\nlKQKQ1KSKgxJSaowJCWpwpCUpApDUpIqDElJqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQlKQKQ1KS\nKgxJSaowJCWpwpCUpApDUpIqpnT5ht4Hn0d/5r7vde1zfqU80ltf26PDn/cga9Zc12OPffbVr7lc\nhwAnH/eqXvuDM7e5B7ckJanCkJSkCkNSkioMSUmqMCQlqcKQlKQKQ1KSKgxJSaowJCWpwpCUpApD\nUpIqDElJqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQlKQKQ1KSKmb1Gjd9Xx9jze39XVuk72ttzOVr\n5mTD4t760kxcp6U/c7kOof/59cEtSUmqMCQlqcKQlKQKQ1KSKgxJSaowJCWpwpCUpApDUpIqDElJ\nqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQlKQKQ1KSKgxJSaowJCWpwpCUpIqUUsZvnNwFrJu56WiW\nHFRK2Xu2JzFbrOsd2jbX9pRCUpLmG3e3JanCkJSkil5DMsmSJD9Ksm+f/bZ9r03yir773dEluTzJ\nqSOWXZjklU/1nOYSa3bHlqQkOWzI40uS/DjJpMcr+96SfBvw3VLKHT33O6clOTjJd5Lcl2RdkjfP\n9pzG9AngY7M9iVlmzW5fNduLUsrDwBeA35isbd8h+XbgSz33uT34bWAtsBw4BvhRXx0nWdBXX4NK\nKdcAuyV54UyNsR2wZnuu2b4kWTjDQ5wPvCXJklqjKYVku/twert7sjnJF5Ps3C47EDgEuLrT/pwk\nf5jkm+0n1tVJDu0sPzbJ95Pc0/577CRTOGrE2KuTXDEw1yc2s9t5/FGSS5Pcn+R7SZ6Z5FNtXz9O\n8m/GWc8RHgU2lFIeLaX8tJTyd2O8lscluTLJliTrk6zuzPUzSS5J8gDw8nbX4MwktyW5M8lnkyxt\n2++Z5OIkd7VzvTjJASPG3DfJ9Une13n4cuDkyea7vbJmR5pSzSZZ2c7vLW0dbkzywc7yJe3cbm9/\nPlULnySnJbmhfY1/lOQXOuvxgSTXAw8kWZhkvzSHhu5KcmuSd3b6OTrJ37bvozuSnJ1k8Ygxj2vf\na78IUErZAGwGXlRbd0opY//QfPL8EFhB8wn0PeBj7bKTgf830P4cYBNwNLAQOA+4oF22vJ3gm9pl\nr2/v7zWNsVcDVwy0L8BhnXlsBF4A7Ax8G7gVeDOwgGaX86/HGWvE3N4BPAycOObreBBwX7vOi4C9\ngFWdud4DvITmQ2xn4CzgG+1clgH/G/idtv1ewKuBXdplXwUu6ox1OXAqcDBwI/C2gbn8N+AvplIH\n29OPNdtbza5s5/c5YCnw/Pb5z26X/xZwFfAMYG/gSuCjI/p6DfBPwFFAgMNozmecWI/r2vVY2r4H\nrgU+DCym+VC7BTihbf8CmpBb2M7xBuDdg68pcCKwHjh6YC7fAN5ZXfdpFNzbO/dPAm5ub58CXDWk\n4P50oP2P29tvAq4ZaP+3wOppjD1OwX1uoEBu6Nx/LrBlnLGGzOslbfG+DNgwUXTtL2Yj7bmoA885\nHfj6iP7OAc7t3A/wAHBo57EXA7eOeP4qYHPn/uXA77fr9Poh7U8Dvj2VOtiefqzZ3mp2ZTu/AzqP\nXQP8Snv7ZuCkzrITgLUjxl8DvKvymr21c/8Y4LaBNqcDXxzx/Hd331vtnE+n+c8CPz+k/XnAh2s1\nNJ19/vWd2+uA/drbm2m2ZAb9tHP7QeBp7e392Pp/OawD9p/G2OO4s3P7oSH3n/bk5mOP9V+BPyml\nfCfJfwIuSfImYB+aT/phZ+uvoCmqUbpj702zlXhtkonHQrM1QZJdaLY0TwT2bJcvS7KglPJ4e/8U\n4CfA14aMtQzYUpnLjsCafbLp1OyEcV+b2vhTqf+DgP2SdGt0AfA3AEkOp9kIeCHN+2QhzZZn17tp\nNjx+OGSsSet/Ol/crOjcPhC4vb19PXBwxj/YejvNC9B1IM1m+FTHfoDmBQIgyTPHnEPNqLEGLaTZ\nZaaU8n3gdcAFwBnA7414znrg0BHLoPn0m7CR5g1xZCllj/Zn91LKRHG+BzgCOKaUshvw0vbxdPo4\no+3n/Gz9RdCzgR9U5rIjsGafbDo1O5nB16Y2/lTqfz3NXtMenZ9lpZST2uWfAX4MPKut/9/kybUP\nze79q5K8a8hYk9b/dELy15MckGQ58EHgK/DEQdCf0BzLGcclwOFJ3tAenH0d8Bzg4qmOTbOSRyZZ\n1R6sPmPKazX+WIO+CrwzyUuT7ATcQbPL8EzgsRHPOQ94RZLXtuu+V5JVwxqWUv6V5jjQWUmeAZBk\n/yQntE2W0YTolnauHxnSzaM0hbIrcG47zwkvAy4dMc8dhTX7ZNOp2cl8GfhQkr2TPJ3mGOKfj2j7\np8B7k7wgjcOSDH74TLgGuK/9MmdpkgVJfj7JUe3yZcC9wP1Jfg74L0P6uB04HnhXkieWJ9mf5vjt\nVbUVm05Ing/8Jc3B05t58nl2f0xz3GZSpZRNwC/RbAltAt4P/FIpZeNUxy6l3Ehz4Pgy4CbgilEd\nTEFtPbvr8b9ozrX6E5ovY75O8zq8D7g4zTeog8+5jeaY0XuAu2kOVD+/MpcP0LyZr0pyL816HtEu\n+xTNAe6NNL/sb42Y5yPAL9PsUn0hyU5tod1fmlOBdmTW7JPXY8o1O4aPAX9Hs3X+f4G/r4z/VeDj\n7XzvAy6iCathbR+nec1X0RxH3UgTsru3Td4LvKHt53OM+GBo33PHA7+Rn/3nijcAf1aacyZHmupf\nAVoLnFpKuWzE8iXAPwDHl+345NzJ1nNHkeRC4POllEtmey4zxZrVMO3v/QfAS0sp/1xr2+vJmm0i\nP6fPPjVzSimvnu05zDZrdn5qf+8/N05b/8CFJFX49yQlqcItSUmqmNIxyacvX1BWrljU2+A3/eMe\nvfUF8Kwj+jsnuu+5zWUPPXoPjzz+0OC5ZfPG4gVLy9JFu0/ecEx91iH0W4t9z61vN16/y+SNpuA+\nNm8s23j5himF5MoVi7hmzYrJG47p5ONe1VtfAN9cc1FvffU9t7nsyg3z8Y/g/MzSRbtz7AFjnQU0\nlj7rEPqtxb7n1rcT9ht6uvC0XVa+ts3XLnJ3W5IqDElJqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQ\nlKQKQ1KSKgxJSaowJCWpwpCUpApDUpIqDElJqjAkJanCkJSkCkNSkioMSUmqmNLlG276xz3m9GUN\n+vzT7wsP6a2rGfHYLWt766uUR3rra3v0rCO2zPnLGswXCw9Z2W+HN297F25JSlKFISlJFYakJFUY\nkpJUYUhKUoUhKUkVhqQkVRiSklRhSEpShSEpSRWGpCRVGJKSVGFISlKFISlJFYakJFUYkpJUYUhK\nUoUhKUkVhqQkVUzpGjd9XwtkLl8v55tXzO1rnvR5PZ/5ru9rN83l2pnL77m5yi1JSaowJCWpwpCU\npApDUpIqDElJqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQlKQKQ1KSKgxJSaowJCWpwpCUpApDUpIq\nDElJqjAkJaliSpdvuPH6XXq+bMDaHvvqV99/5r7vP+m/8JCVvfWVDYt760tz+xIJfdfhXH+fLNh3\n2/twS1KSKgxJSaowJCWpwpCUpApDUpIqDElJqjAkJanCkJSkCkNSkioMSUmqMCQlqcKQlKQKQ1KS\nKgxJSaowJCWpwpCUpApDUpIqDElJqjAkJaliSte46X3wHq/T0rf5dC2Qo0/Y0ltf6r92+ryu1Fyu\nQ5iJ6wOduc09uCUpSRWGpCRVGJKSVGFISlKFISlJFYakJFUYkpJUYUhKUoUhKUkVhqQkVRiSklRh\nSEpShSEpSRWGpCRVGJKSVGFISlKFISlJFYakJFUYkpJUkVLK+I2Tu4B1MzcdzZKDSil7z/YkZot1\nvUPb5tqeUkhK0nzj7rYkVRiSklRhSEpShSEpSRWGpCRVGJKSVGFISlKFISlJFYakJFX8f61hNSZ2\n7WqSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110b68c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(2, 2)\n",
    "for ax, probs, (bump, creak) in zip(axs.flatten(), p_emission, map_visible_state_to_bump_creak):\n",
    "    ax.imshow(probs.reshape(height, width))\n",
    "    ax.set_title(\n",
    "        'p('\n",
    "        + ('' if bump else 'no ') + 'bump & '\n",
    "        + ('' if creak else 'no ') + 'creak)'\n",
    "    )\n",
    "    hide_ticks(ax)\n",
    "maybe_save_plot('2018-05-02-v-given-h')    \n",
    "plt.show()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating data: $\\textbf{v}_{1:t}$\n",
    "\n",
    "Now I have all of the probabilities I need. On to part 2: generating some data!\n",
    "\n",
    "I'll do it for 10 timesteps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "timesteps = 10\n",
    "\n",
    "hiddens = np.zeros(timesteps, dtype=int)\n",
    "visibles = np.zeros(timesteps, dtype=int)\n",
    "\n",
    "hiddens[0] = np.random.choice(num_hidden_states, p=p_hidden_start)\n",
    "visibles[0] = np.random.choice(\n",
    "    num_visible_states,\n",
    "    p=p_emission[:, hiddens[0]]\n",
    ")\n",
    "\n",
    "for t in range(1, timesteps):\n",
    "    hiddens[t] = np.random.choice(\n",
    "        num_hidden_states,\n",
    "        p=p_transition[:, hiddens[t - 1]]\n",
    "    )\n",
    "\n",
    "    visibles[t] = np.random.choice(\n",
    "        num_visible_states,\n",
    "        p=p_emission[:, hiddens[t]]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing\n",
    "\n",
    "Let's take a peak at where they were moving and what sounds they made.\n",
    "\n",
    "Visibles are plotted as a rectangle with two colors. The left is black if there was a bump and white if not. The right is black if there was a creak and white if not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8sAAAC+CAYAAAAC9p1AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC3NJREFUeJzt3V+MpXddx/HPt13tohdUaFNowN0AXggxmqwRJII1Ek1N\nAw0xmCC0amLkxj+Rm6rR1EQu/FdiogZvFLCCVIyFYJRqoljEaneIaEy9sJVKIWBru1AUm/75enFO\nw/hlZnfO2fbM2enrlUwye+Y55/zmO0+e2Xee58yp7g4AAADwZRcd9gIAAABg24hlAAAAGMQyAAAA\nDGIZAAAABrEMAAAAg1gGAACAQSwDAADAIJYBAABgEMsAAAAwHDvsBWybyy67rE+ePHnYy/gKOzs7\nh72E/TzQ3Zevc8eq6qd6MUfc2rO2X6+mu2ud+9mnV3bk9ulttbOzc+SO1adOnTrsJezpKM56i5n1\nhqz7e3Fbj9Xb+v+P+L24MQc9Vovl4eTJkzl9+vRhL+MrVK11jNqEew97Ac8ga8/afs2WOnL79Laq\nqiN3rN7Wn/9RnPUWM+stt63H6i3+/4ffixty0GO1y7ABAABgEMsAAAAwiGUAAAAYxDIAAAAMYhkA\nAAAGsQwAAACDWAYAAIBBLAMAAMAglgEAAGAQywAAADCIZQAAABjEMgAAAAxiGQAAAIbziuWqekdV\n/fwBtvuXqrpq+fmNVXXzWbb9ZFW95nzWBQAAAOfj2Nm+WFV/nuQfuvsXxu2vS/I7SV7Q3Y+d60m6\n+2XntUoAAADYoHOdWX5XkjdVVY3b35zkDw4SygAAAHChOVcs35rkuUle9eQNVfV1Sa5J8u6qemdV\n/dLy9suq6kNVdaaqHqyq26vqouXX5qXVx6vqfVX1cFV9vKq+ea8nr6qLquqGqrq7qv6rqm6pqucs\nv3a8qm5e3n6mqu6sqivWHwUAAAAsnDWWu/tLSW5Jct2um9+Q5F+7+xNj87cmuS/J5UmuSPKzSXqf\nh35dkj9K8pwk70lya1V91R7b/XiSa5N8Z5IrkzyU5LeWX7s+ybOTvDCLoH9Lki+d7fsBAACAgzjI\nH/h6V5Lvr6rjy39ft7xtejTJ85Oc6O5Hu/v27t4vlne6+/3d/WiSm5IcT/KKPbZ7S5Kf6+77uvuR\nJDcu13Js+XzPTfKS7n68u3e6+wsH+H4AAADgrM4Zy9390SQPJLm2ql6c5NuyOBs8/WqSf0tyW1Xd\nU1U3nOVhP7Xr8Z/I4oz0lXtsdyLJnywvsz6T5K4kj2dx5vr3k3w4yR9W1Weq6lf2OTsNAAAAKzno\nW0e9O4szym9K8uHu/tzcoLsf7u63dveLkrw2yU9X1Xfv83gvfPKT5euaX5DkM3ts96kkV3f3pbs+\njnf3p5dnr3+xu1+a5JVZvI76uj0eAwAAAFaySiy/JsmPZu9LsFNV11TVS5Z/OfvzWZwBfmKfxztV\nVa9fXk79U0keSXLHHtu9I8nbqurE8jkuX75tVarqu6rqm6rq4iRfyOKy7P2eDwAAAA7sQLHc3Z9M\n8rEkX5vkg/ts9g1J/jLJF5P8XZLf7u6/2mfbDyT5gSz+YNebk7x++frl6TeWz3dbVT2cRVC/fPm1\n5yV5fxahfFeSj2RxaTYAAACcl2MH3bC7r9rjth/a9fnbk7x9n/ue3PX5jed4nt3bPpHFHwC7aY/t\n3pvkvedYNgAAAKzsoJdhAwAAwDOGWAYAAIBBLAMAAMAglgEAAGAQywAAADCIZQAAABjEMgAAAAxi\nGQAAAAaxDAAAAINYBgAAgEEsAwAAwCCWAQAAYKjuPuw1bJWquj/JvYe9jgvIie6+fJ07mvXKzHoz\nzHlzzHpzzHpzzHpzzHozzHlzzHpzDjRrsQwAAACDy7ABAABgEMsAAAAwiGUAAAAYxDIAAAAMYhkA\nAAAGsQwAAACDWAYAAIBBLAMAAMAglgEAAGAQywAAADCIZQAAABjEMgAAAAxiGQAAAAaxDAAAAINY\nBgAAgEEsAwAAwCCWAQAAYBDLAAAAMIhlAAAAGMQyAAAADGIZAAAABrEMAAAAg1gGAACAQSwDAADA\nIJYBAABgEMsAAAAwiGUAAAAYxDIAAAAMYhkAAAAGsQwAAACDWAYAAIBBLAMAAMAglgEAAGAQywAA\nADCIZQAAABjEMgAAAAxiGQAAAAaxDAAAAINYBgAAgEEsAwAAwCCWAQAAYBDLAAAAMIhlAAAAGMQy\nAAAADGIZAAAABrEMAAAAg1gGAACAQSwDAADAIJYBAABgEMsAAAAwiGUAAAAYxDIAAAAMYhkAAAAG\nsQwAAACDWAYAAIBBLAMAAMAglgEAAGAQywAAADCIZQAAABjEMgAAAAxiGQAAAAaxDAAAAINYBgAA\ngEEsAwAAwCCWAQAAYBDLAAAAMBw77AVsm6rqw17DBeaB7r58nTua9crMekO6u9a5nzmvzD69OWa9\nOWa9OWa9IX4vbox9enMONGtnljlf9x72Ap5BzJqjxj69OWa9OWa9OWbNUWOf3pwDzVosAwAAwCCW\nAQAAYBDLAAAAMIhlAAAAGMQyAAAADGIZAAAABrEMAAAAg1gGAACAQSwDAADAIJYBAABgEMsAAAAw\niGUAAAAYxDIAAAAMYhkAAAAGsQwAAACDWAYAAIBBLAMAAMAglgEAAGAQywAAADCIZQAAABjEMgAA\nAAwXZCxX1VVVdd9hrwMAAICjae1Yrqq/rqqHquqSA2x7sqq6qo6t+3wAAACwKWvFclWdTPKqJJ3k\ntU/hegAAnhbdvdIHAF+26jH0KBxH1z2zfF2SO5K8M8n1T95YVc+qql+vqnur6vNV9dGqelaSv1lu\ncqaqvlhV315VN1bVzbvu+//OPlfVD1fVXVX1cFXdU1U/tuZaAQAAYCXrXhZ9XZKbkvx9kjuq6oru\n/lySX0vysiSvTPLZJC9P8kSSVyf59ySXdvdjSVJV33uO5/jPJNckuWd5/z+rqju7++NrrhkAAAAO\nZOUzy1X1HUlOJLmlu3eS3J3kjVV1UZIfSfKT3f3p7n68uz/W3Y+ss7Du/tPuvrsXPpLktiwu/QYA\nAICn1TqXYV+f5LbufmD57/csb7ssyfEs4vm8VdXVVXVHVT1YVWeSfN/yOQAAAOBptdJl2MvXH78h\nycVV9dnlzZckuTTJ85P8b5IXJ/nEuOter+7+7yRfs+vfz9v1PJck+eMsLvf+QHc/WlW3JqlV1gsA\nAADrWPXM8rVJHk/y0iTfsvz4xiS3ZxG2v5vkpqq6sqouXv4hr0uS3J/Fa5dftOux/jHJq6vq66vq\n2Ul+ZtfXvjqLCL8/yWNVdXWS71n5uwMAAIA1rBrL1yf5ve7+j+7+7JMfSX4zyQ8muSHJPye5M8mD\nSX45yUXd/T9J3pbkb6vqTFW9orv/Isn7kvxTkp0kH3rySbr74SQ/keSWJA8leWOSD57H9wkAAAAH\nVkfh/a+eSlVlIKvZ6e5vXeeOZr0ys96Q7l7rJR/mvDL79OaYdbLye35WrXUoMOvNMesN8XtxY7Z6\nn16nG9c8jm7CgWa97vssAwAAwJEllgEAAGAQywAAADCs9NZRzwSnTp3K6dOnV7rPFl+Lf+Rs6PVm\nrOGIvY7lSPGz2W6Oa5tjdptjv16P4/XqzGxznolzc2YZAAAABrEMAAAAg1gGAACAQSwDAADAIJYB\nAABgEMsAAAAwiGUAAAAYxDIAAAAMYhkAAAAGsQwAAACDWAYAAIBBLAMAAMBw7LAXsG12dnZSVYe9\nDPbhZ7O9/Gy2l5/NdvPz4SiyX6/H3FZnZjydnFkGAACAQSwDAADAIJYBAABgEMsAAAAwiGUAAAAY\nxDIAAAAMYhkAAAAGsQwAAACDWAYAAIBBLAMAAMAglgEAAGAQywAAADAcO+wFbKEHktx72Iu4gJw4\nj/ua9WrMejPMeXPMenPMenPMenPMejPMeXPMenMONOvq7qd7IQAAAHBBcRk2AAAADGIZAAAABrEM\nAAAAg1gGAACAQSwDAADAIJYBAABgEMsAAAAwiGUAAAAYxDIAAAAM/wcnrRUJUmnrFAAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110b68860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, all_axs = plt.subplots(2, timesteps, figsize=(16, 4))\n",
    "all_axs = all_axs.T\n",
    "\n",
    "VISIBLES = 0\n",
    "TRUE_STATES = 1\n",
    "\n",
    "all_axs[0][0].set_title('Visibles', x=-0.5, y=0.2)\n",
    "all_axs[0][1].set_title('Actual', x=-0.5, y=0.4)\n",
    "\n",
    "for i, (axs, hidden, visible) in enumerate(zip(all_axs, hiddens, visibles)):\n",
    "    axs[VISIBLES].imshow([map_visible_state_to_bump_creak[visible]], cmap='gray', vmin=0, vmax=1)\n",
    "    hide_ticks(axs[VISIBLES])    \n",
    "    \n",
    "    axs[TRUE_STATES].imshow(plot_state_in_room(hidden), cmap='gray')\n",
    "    hide_ticks(axs[TRUE_STATES])\n",
    "    \n",
    "maybe_save_plot('2018-05-02-v-over-time')\n",
    "plt.show()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering: $p(h_t \\mid v_{1:t})$\n",
    "\n",
    "Now for part 3: filtering using alpha recursion. Filtering finds the probability of the current hidden state given the visibles from all previous timesteps, as well as the transition, emission, and start probabilities.\n",
    "\n",
    "Taking from Barber, alpha recursion gives\n",
    "\n",
    "$$\\alpha(h_t) = p(h_t, v_{1:t}) \\propto p(h_t \\mid v_{1:t})$$\n",
    "\n",
    "The values of $\\alpha$ are computed with\n",
    "\n",
    "$$\\alpha(h_t) = p(v_t \\mid h_t) \\sum_{h_{t - 1}} p(h_t \\mid h_{t - 1})\\alpha(h_{t - 1})$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\alpha(h_1) = p(h_1, v_1) = p(v_1 \\mid h_1)p(h_1).$$\n",
    "\n",
    "In code I'll compute an $\\alpha$ for each timestep. I'll use these to give the distribution over all states for that timestep given the visibles up until that timestep. Because I want the filtered posterior, I can rescale, so like Barber's code, I'm normalizing each $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def alpha_recursion(visibles, p_hidden_start, p_transition, p_emission):\n",
    "    num_timestamps = visibles.shape[0]\n",
    "    num_hidden_states = p_transition.shape[0]\n",
    "    \n",
    "    # There will be one alpha for each timestamp \n",
    "    alphas = np.zeros((num_timestamps, num_hidden_states))\n",
    "\n",
    "    # alpha(h_1) = p(h_1) * p(v_1 | h_1)\n",
    "    alphas[0] = p_hidden_start * p_emission[visibles[0]]\n",
    "\n",
    "    # normalize to avoid overflow\n",
    "    alphas[0] /= np.sum(alphas[0])\n",
    "    for t in range(1, num_timestamps):\n",
    "        # p(v_s | h_s)\n",
    "        # size: new_states\n",
    "        corrector = p_emission[visibles[t]]\n",
    "        \n",
    "        # sum over all hidden states for the previous timestep and multiply the \n",
    "        # transition prob by the previous alpha\n",
    "        # transition_matrix size: new_state x old_state\n",
    "        # alphas[t_minus_1].T size: old_state x 1\n",
    "        # predictor size: new_state x 1,\n",
    "        predictor = p_transition @ alphas[t - 1, None].T\n",
    "\n",
    "        # alpha(h_s)\n",
    "        alphas[t, :] = corrector * predictor[:, 0]\n",
    "\n",
    "        # normalize\n",
    "        alphas[t] /= np.sum(alphas[t])\n",
    "\n",
    "    return alphas\n",
    "\n",
    "alphas = alpha_recursion(\n",
    "    visibles, \n",
    "    p_hidden_start,\n",
    "    p_transition,\n",
    "    p_emission,\n",
    ")\n",
    "\n",
    "assert np.all(np.isclose(np.sum(alphas, axis=1), 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing\n",
    "\n",
    "Again, but with the predictions!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA80AAADhCAYAAADyMFclAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFy5JREFUeJzt3Xt0ZWd5H+DfK8mj8YzxfbgMxjP4QgIsEhMTwFmQOoGF\nwUAgSSFtSAyB0tA/2gTSlZBr3aY0TUNhJU26aGkKTmlCDGkgF2KDQ6B2sak9WVxqXFJsZrAxtmew\nB4+Zi2ekr39IA/J4tq3vaHQky8+zltaSjt53f5++s88+56e9dVSttQAAAAAPNrHSEwAAAIDVSmgG\nAACAAUIzAAAADBCaAQAAYIDQDAAAAAOEZgAAABggNAMAAMAAoRkAAAAGCM0AAAAwYGqlJ7BanX76\n6W3r1q0rPY0H2bZt20pP4ahaazVKX1W1Yz2XNW5Xa23TKI326W4jr7X9utua269Xo+3bt2fXrl0j\nHauT1btfn3/++Ss9hQdZylqv1nVexRyrx2fNHavX4muQ1brWq1HPsVpoHrB169bccMMNKz2NB6ka\n+fUOa8OOURvt091GXmu6rbn9ejV61rOetdJTWBar8f5fq2u9SjlWj8+aO1avxdcgq3WtV6OeY7XL\nswEAAGCA0AwAAAADhGYAAAAYIDQDAADAAKEZAAAABgjNAAAAMEBoBgAAgAFCMwAAAAwQmgEAAGCA\n0AwAAAADhGYAAAAYIDQDAADAAKEZAAAABhyT0FxV76qqX11E3Y1VdeH855dW1fseonZ7Vb3wWMwP\nAAAARjG1mKKquiLJ/26t/doRt78iyX9KckZr7dDDbae19vSRZgkAAAArYLFnmi9L8hNVVUfc/pNJ\n/vtiAjMAAAA80iw2NH8oyWlJnn/4hqo6JcnLkvxBVb23qv71/O2nV9VfVNXuqrq7qq6uqon57x15\nyfX6qvrjqtpTVX9bVd99tMGraqKq3lpVN1fV16vq8qo6df5766vqffO3766q66vqcf1LAQAAAA+0\nqNDcWtuX5PIklyy4+dVJ/m9r7bNHlP9cktuSbEryuCS/lKQNbPoVST6Q5NQkf5jkQ1V13FHq/mmS\nVyb5e0k2J7knye/Nf++1SU5K8qTMBfs3Jdm3mJ8LAAAAHkrPG4FdluTvV9X6+a8vmb/tSAeTPCHJ\nltbawdba1a21odC8rbX2wdbawSTvSLI+yXOPUvemJL/cWruttXYgyaXzc5maH++0JOe01mZaa9ta\na/d2/FwAAABwVIsOza21a5LsSvLKqjo7ybMzd3b4SL+V5EtJPlpVt1TVWx9is7cu2P5s5s5Qbz5K\n3ZYkfzp/+fXuJDclmcncmez/luTKJO+vqtur6t8NnK0GAACALr3/cuoPMneG+SeSXNlau/PIgtba\nntbaz7XWzkryQ0neUlUvGNjekw5/Mv93z2ckuf0odbcmeUlr7eQFH+tba1+dP5v9L1trT0vyfZn7\nO+tLjrINAAAA6DJKaH5hkjfm6Jdmp6peVlXnzL/T9jcyd0Z4dmB751fVj8xfZv2zSQ4kue4ode9K\n8raq2jI/xqb5f3eVqvqBqnpGVU0muTdzl2sPjQcAAACL1hWaW2vbk3wqycYkfzZQdm6Sq5Lcl+Ta\nJP+xtfY3A7UfTvJjmXtjr59M8iPzf998pN+eH++jVbUnc8H6OfPfe3ySD2YuMN+U5JOZu2QbAAAA\nlmSqt6G1duFRbnvdgs/fmeSdA71bF3x+6cOMs7B2NnNvFPaOo9T9UZI/ephpAwAAQLfey7MBAADg\nUUNoBgAAgAFCMwAAAAwQmgEAAGCA0AwAAAADhGYAAAAYIDQDAADAAKEZAAAABgjNAAAAMEBoBgAA\ngAFCMwAAAAwQmgEAAGBAtdZWeg6rUlXtTLJjpefxCLGltbZplEbr3M1aj4+1Hh9rPR4jr3NirTvZ\np8fHWo+PtR4faz0ei15noRkAAAAGuDwbAAAABgjNAAAAMEBoBgAAgAFCMwAAAAwQmgEAAGCA0AwA\nAAADhGYAAAAYIDQDAADAAKEZAAAABgjNAAAAMEBoBgAAgAFCMwAAAAwQmgEAAGCA0AwAAAADhGYA\nAAAYIDQDAADAAKEZAAAABgjNAAAAMEBoBgAAgAFCMwAAAAwQmgEAAGCA0AwAAAADhGYAAAAYIDQD\nAADAAKEZAAAABkyt9ARWq6pqKz2HR5LWWo3SZ5277WqtbRql0Vp3s9bjY63HZNRjdWKte3leHBvH\nj/Gx1uNjrcdkscdqZ5rhkWXHSk/gUcRaj4+1Bkbl+DE+1np8rPUqIzQDAADAAKEZAAAABgjNAAAA\nMEBoBgAAgAFCMwAAAAwQmgEAAGCA0AwAAAADhGYAAAAYIDQDAADAAKEZAAAABgjNAAAAMEBoBgAA\ngAFCMwAAAAwQmgEAAGCA0AwAAAADhGYAAAAYIDQDAADAAKEZAAAABgjNAAAAMEBoBgAAgAGP6NBc\nVRdW1W0rPQ8AAADWpiWH5qr6RFXdU1XTi6jdWlWtqqaWOi4AAAAstyWF5qramuT5SVqSHzoG8wEA\nAIBVY6lnmi9Jcl2S9yZ57eEbq+r4qvr3VbWjqr5RVddU1fFJ/ud8ye6quq+qLqiqS6vqfQt6H3A2\nuqp+qqpuqqo9VXVLVf30EucMADyKtda6PgD4tt5j6Fo4ji71MulLkrwjyaeTXFdVj2ut3Znk7Ume\nnuT7ktyR5DlJZpN8f5IvJzm5tXYoSarqoocZ464kL0tyy3z/X1XV9a21v13i3AEAAOAhjXymuaqe\nl2RLkstba9uS3Jzkx6tqIsnrk/xMa+2rrbWZ1tqnWmsHRhmntfaXrbWb25xPJvlo5i4JBwAAgGW1\nlMuzX5vko621XfNf/+H8bacnWZ+5EL1kVfWSqrququ6uqt1JLp4fAwAAAJbVSJdnz/998quTTFbV\nHfM3Tyc5OckTkuxPcnaSzx7RerQL2r+ZZMOCrx+/YJzpJH+SucvAP9xaO1hVH0pSo8wbAAAAeox6\npvmVSWaSPC3JefMfT01ydeYC7n9N8o6q2lxVk/Nv+DWdZGfm/rb5rAXb+kyS76+qM6vqpCS/uOB7\n6zIXxncmOVRVL0nyohHnDAAAAF1GDc2vTfKe1tpXWmt3HP5I8rtJXpPkrUk+n+T6JHcn+c0kE621\nvUneluR/VdXuqnpua+1jSf44yeeSbEvyF4cHaa3tSfLPklye5J4kP57kz0acMwAAAHSptfAW4Muh\nqixMh9baSJfMW+du21przxql0Vp3s9bjY63HZNRjdbK21rr3tU9V/7J5Xhwbx4/xsdbjs6rXepT8\nOMpxdBwWe6xe6v9pBgAAgDVLaAYAAIABQjMAAAAMEJoBAABgwEj/p/nR4Pzzz88NN9zQ1bNa/8B9\nrRnHG7gwmrX0xhBrkftndXNsGx9rNx726dE4Vo/Guo3Po3HdnGkGAACAAUIzAAAADBCaAQAAYIDQ\nDAAAAAOEZgAAABggNAMAAMAAoRkAAAAGCM0AAAAwQGgGAACAAUIzAAAADBCaAQAAYIDQDAAAAAOm\nVnoCq9W2bdtSVSs9DY7C/bJ6uW9WN/fP6ub+Ya2xT4/Guo3GurGcnGkGAACAAUIzAAAADBCaAQAA\nYIDQDAAAAAOEZgAAABggNAMAAMAAoRkAAAAGCM0AAAAwQGgGAACAAUIzAAAADBCaAQAAYIDQDAAA\nAAOmVnoCq9iuJDtWehKPEFuW0Gud+1jr8bHW42Otx2Mp65xY6x726fGx1uNjrcfHWo/Hote5WmvL\nOREAAAB4xHJ5NgAAAAwQmgEAAGCA0AwAAAADhGYAAAAYIDQDAADAAKEZAAAABgjNAAAAMEBoBgAA\ngAFCMwAAAAwQmgEAAGCA0AwAAAADplZ6AqvVuppu62tjV0+tm+6qb/cf6KofZYxx2HfwG7l/Zm+N\n0ruuptv6dK7z5OQoQ/WZ7P99UpvqnNfe/d1j7Mk9u1prm7obM+JaV9/d2lrrqk+Smlj+3921Ntvd\ns6eNd63TudYjPeA6xxilp83MdA8x7v16raipvqfwfTN7cv/svpF2nWQ8az3K8aDN9j++l9v+fDP3\ntwNjfF7sW7c2079mozz3jnI86LW048f6dvzECZ1Nx/XVHzrUV59k9vj+13r3n9y3u03f+s3uMcb+\nGqRznxtlf6v1I7yuPth5n070HwruPbhz9LWe6N+ve48Jhzb1PxdM7ezf53pf7DzlGXu76rffejC7\n7p5Z1ChC84D1tTHPnbqoq2di65au+tlbdnTVjzLGOFy7/bKRe9dnY54z+aKunsmTTuwbZITAVL1j\nJJnZdFJXfdt2Y/cYV81+oH+nmbc+G/OcekFXz8T69V31s/v7fxEwsWH5A0870P8Lqo8dfP9Y17qO\nW9dXP8IvdjLCi946ru9pYmb3N7rHuKp9cGlrPfHCvqbeX+6M8suGEX6B1Gvy9Md21V+76wNLGm8c\naz1xwmP6tp9kds+e7p5MdD4WZvtejH+6/XXf9hcY5fgxeULfc9bMvfd21SfJ5EmndPfM7N7d1zDC\n42Ypx4/jJ07Icze8rKunztzcN8hdX++rT7L/mU/u7tn+8r4wf+7PXNc9xpKP1b37dec+172/JZk8\n+9zunnxtZ1d5Tfc9vyfJFV/7vSXt1xec8Iqunt5jws5XXdBVnySb3tW/z/X+4uSKK2/oqn/2Rbcu\nutbl2QAAADBAaAYAAIABQjMAAAAMEJoBAABggNAMAAAAA4RmAAAAGCA0AwAAwAChGQAAAAYIzQAA\nADBAaAYAAIABUys9gdWq1k1nYuuWZR3jiq/c0N1z8YXLO6ck+cgn/qSr/tkX3TPyWDU5mcmTTuzq\nueXN39lVf9Y7vtBVnyQ7/uGTunv2PmG2q/47d5zePUbu6m85rKoysX59X8+6dV31O37xe7rqk+TJ\nb/98d0+v+155fn/TB94/+oBVqenprpb9L/iurvr1f/25rvokOfC8p3b3TH3zUFf95PU3dY+R/f0t\nD9BaX31VV/nkYx7Tt/0kM/fe293Ta/bu3V31bWZm6YNW7+/a+46Ls3v2dG4/3ffn3EDHYC2WS1Xq\nuL5jbzvU9zgdac1a332ZJJOP3dRVP3PnEp7kRjG9Ljm77/n+qz94alf95nff3lWfJF/+sf5zWmd+\nuO/+mdi4sXuM3NffclhNTmbyxJP6mib71uHVX/ha3/aT/I+Lz+zuaU/o26/3be5/Dkn/j/ItbWa2\n+zloovN5bma6/xgyynPp7L6+Fwgvfvlruur/7kv/edG1zjQDAADAAKEZAAAABgjNAAAAMEBoBgAA\ngAFCMwAAAAwQmgEAAGCA0AwAAAADhGYAAAAYIDQDAADAgJFCc1WdWVX3VdXk/NefqKp/dKwmVVXP\nr6ovHqvtAQAAwCgeNjRX1faq2jcfku+rqvuSHGqtndBamzlK/euq6pqlTKq1dnVr7TuWsg0AAABY\nqsWeaX75fEg+/HH7ck2oqqaWa9sAAADQY6SAWlVbk3w5yXGttUMLbn9qknclOW7BGemTq2o6yduS\nvDrJdJI/TfLm1tq+qrowyfuS/Ickb07ysar6/STva62dMb/d7Ul+N8klSbYkuSLJa1tr++e///Pz\nvS3JryV5d5JzW2tfGuXnS5J2/4HM3rKjq2firC2jDrdovXNa9SrJZN9fCZz+7Dv7hti4sas+Sc58\n143dPZmc7O8Zo9ZaZg8c6OqZ6PyZpnd3lSdJZu+7r7+pta7yAyeO+e0bWks7eOjh6xa47TV99ed+\nsv/wPf3xz3X3ZKL66mf77pulqqmpTJ66qatnZufOvkHWHddXn2Riw4bunt59pvu+WQET09Nd9bP7\n93ePMXnaqd09NdX3+Dl0R9/zzpK0lnbw/r6WmQdd+HfMzXzj3mUfIzXCPr2EQ07bvz/tC30vFU95\n0jO76mf39e/Tx+3sP75v+OT/6aofZV5L0WZmuvehiRNO6Kp/w0l3dNUnyeXbv9LdU1N9zwnrNjyl\ne4ylqMmJTJ5wYl/Pxr7nrOnd/Q+8m3/h6d09Z7+379h78PjO5+uO59Fj+kqytXZTkjcluXb+jPTJ\n89/6t0mekuS8JOckeWLmwu1hj09yauYC8T8e2Pyrk7w4yZOTfFeS1yVJVb04yVuSvHB+2xcesx8I\nAACAR7XFhuYPVdXu+Y8P9QxQVZW5IPzm1trdrbU9Sf5Nkn+woGw2yb9orR1ore0b2NTvtNZub63d\nneTPMxfAk7kw/Z7W2o2ttb1JLu2ZHwAAAAxZ7PUfr2ytXXX4i/nLsxdrU5INSbbVty+7qSQLr/vc\nefhS64ew8JqLvUk2z3++OckNC753a8fcAAAAYNByvOnWkRe570qyL8nTW2tfXWRPj68lOWPB109a\nwrYAAADgW5bj3XHuTHJGVa1LktbabObemOudVfXYJKmqJ1bVRcdovMuT/FRVPbWqNiT51WO0XQAA\nAB7lliM0fzzJjUnuqKpd87f9QpIvJbmuqu5NclWSY/J/mFtrf5Xkd5L8zeEx5r/V9zbBAAAAcISH\nvTy7tbb1KLdtz9zfJR/++sIFn9+f5KVH1O9P8kvzH0du6xN54OXVD7rtyDm01i494uvfSPIbybf+\n7dVs5i7bBgAAgJGN+Z+XLo+q+uGqmq6qU5L8ZpI/X/j/owEAAGAUayI0J/npJHcluTnJTJJ/srLT\nAQAAYC1YjnfPHrvW2otXeg4AAACsPWvlTDMAAAAcc2viTPNyeMoz9uaKK2/o7Oqt7zdx1pbuntlb\ndnTVX3zhj3bV/7/tl3XVP8DkZHLKSV0tJ72p88/Vqx6+5gh3vepp3T2nv+f6rvodv/zs7jFyaX/L\nYbV+OpNnn9vV85GrLu+qv/AN/euWGuV3d7Nd1ae+57qHLzqGanpdJrf0PVbPef0XuuonTjm5qz5J\nsuH47pZDt2zvqp993nndY+Tq/pbD2qFDmdm16+ELF/Zc8N1d9bPX39hVn8zNq9fEhg1d9bN79/YN\n0Fpf/VG30ffYm92/v6u+pvpflsx8/e7unpqc7GuY6Kyf6Stfqt6fpx3sn+CVt3+mu+eiJz6zr2Gk\n54PRVap7n5v+SN9rve2//tyu+iQ557dv7u75yBf7DqQv/d6Lu8fIbf0th1VVJqanu3p+5bOf7Kp/\n6fe8pKs+Se655KzunlMuu7ar/tCJfT/3Uu3ffHy++M/7Xo99x6U3ddW/6Gev6apPkuvP6zyOJrnr\njRd01Z/2Xzpf783uW3SpM80AAAAwQGgGAACAAUIzAAAADBCaAQAAYIDQDAAAAAOEZgAAABggNAMA\nAMAAoRkAAAAGCM0AAAAwQGgGAACAAUIzAAAADBCaAQAAYMDUSk9gLbn4wh9d9jFmb9nR3TNx1pZl\nHaMdOtBV/wAzM8ndu7ta/vLzH++qv+iJz+yqT5LT33Nnd0+q73dQZ/76p7uH+Lvujm87dEay67da\nV89FZ5zfVb/xcbd21SfJodmZ7p7Vrh24PzM39z2OrrxtW1f9RZvP66pPklQte8/ENZ/pH2MJav10\nJree3dUzc+1nu+qvuL3/Z+p97CTJ7N69XfU1Pd03wIER7v8jtb5jyMR5T+uqn/3MF7rqk+TKUe6f\n3ueFzp97SSqpqb6XZ+3g/X1DdG4/Ge2Y8+6vXN1V/8Yzn9c9xlK01jK7f39Xz32vek5X/dZfubar\nPkk+MsI+ffHTf6Cr/qa3b+4eI6/vbzns4GkbsvOH+x53/+qsvvvm0A+e0VWfJKdc1n//fP0NF3TV\nn/b7/WMsxfSt38w5b7m+r2nTaV3ln/r5vsdBktz7huO6e057d9/aHXjJ93bVt2s+tehaZ5oBAABg\ngNAMAAAAA4RmAAAAGCA0AwAAwAChGQAAAAYIzQAAADBAaAYAAIABQjMAAAAMEJoBAABggNAMAAAA\nA4RmAAAAGCA0AwAAwIBqra30HFalqtqZZMdKz+MRYktrbdMojda5m7UeH2s9PtZ6PEZe58Rad7JP\nj4+1Hh9rPT7WejwWvc5CMwAAAAxweTYAAAAMEJoBAABggNAMAAAAA4RmAAAAGCA0AwAAwAChGQAA\nAAYIzQAAADBAaAYAAIABQjMAAAAM+P8VXmKfOpxPcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1120b41d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, all_axs = plt.subplots(3, timesteps, figsize=(16, 4))\n",
    "all_axs = all_axs.T\n",
    "\n",
    "VISIBLES = 0\n",
    "TRUE_STATES = 1\n",
    "FILTERING = 2\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",
    "\n",
    "for i, (axs, hidden, visible, alpha) in enumerate(zip(all_axs, hiddens, visibles, alphas)):\n",
    "    axs[VISIBLES].imshow([map_visible_state_to_bump_creak[visible]], cmap='gray', vmin=0, vmax=1)\n",
    "    hide_ticks(axs[VISIBLES])    \n",
    "    \n",
    "    axs[TRUE_STATES].imshow(plot_state_in_room(hidden), cmap='gray')\n",
    "    hide_ticks(axs[TRUE_STATES])\n",
    "    \n",
    "    axs[FILTERING].imshow(alpha.reshape(height, width), vmin=0)\n",
    "    hide_ticks(axs[FILTERING])      \n",
    "    \n",
    "maybe_save_plot('2018-05-02-filtering')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## See Also\n",
    "\n",
    " - [This post](2018-05-13-viterbi-message-passing.ipynb) builds on this post to show the Viterbi algorithm.\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": {},
   "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
}