viterbi_algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probabilistic Reasoning over Time\n",
    "---\n",
    "# Finding the Most Likely Sequence with Viterbi Algorithm\n",
    "\n",
    "## Introduction\n",
    "An ***Hidden Markov Model*** (HMM) network is parameterized by two distributions:\n",
    "\n",
    "- the *emission or sensor probabilties* giving the conditional probability of observing evidence values for each hidden state;\n",
    "- the *transition probabilities* giving the conditional probability of moving between states during the sequence. \n",
    "\n",
    "Additionally, an *initial distribution* describes the probability of a sequence starting in each state.\n",
    "\n",
    "At each time $t$, $X_t$ represents the *hidden state* and $E_t$ represents an *observation* at that time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from probability import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;32mclass\u001b[0m \u001b[0mHiddenMarkovModel\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m\"\"\"A Hidden markov model which takes Transition model and Sensor model as inputs\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mdef\u001b[0m \u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtransition_model\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msensor_model\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprior\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransition_model\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtransition_model\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msensor_model\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msensor_model\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprior\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mdef\u001b[0m \u001b[0msensor_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mev\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32mif\u001b[0m \u001b[0mev\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msensor_model\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msensor_model\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%psource HiddenMarkovModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding the Most Likely Sequence\n",
    "\n",
    "There is a linear-time algorithm for finding the most likely sequence: the easiest way to think about the problem is to view each sequence as a path through a graph whose nodes are the possible states at each time step. Now consider the task of finding the most likely path through this graph, where the likelihood of any path is the product of the transition probabilities along the path and the probabilities of the given observations at each state. There is a recursive relationship between most likely paths to each state $x_{t+1}$ and most likely paths to each state $x_t$ . We can write this relationship as an equation connecting the probabilities of the paths:\n",
    "\n",
    "$$ \n",
    "\\begin{align*}\n",
    "m_{1:t+1} &= \\max_{x_{1:t}} \\textbf{P}(\\textbf{x}_{1:t}, \\textbf{X}_{t+1} | \\textbf{e}_{1:t+1}) \\\\\n",
    "&= \\alpha \\textbf{P}(\\textbf{e}_{t+1} | \\textbf{X}_{t+1}) \\max_{x_t} \\Big(\\textbf{P}\n",
    "(\\textbf{X}_{t+1} | \\textbf{x}_t) \\max_{x_{1:t-1}} P(\\textbf{x}_{1:t-1}, \\textbf{x}_{t} | \\textbf{e}_{1:t})\\Big)\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The *Viterbi algorithm* is a dynamic programming algorithm for *finding the most likely sequence of hidden states*, called the Viterbi path, that results in a sequence of observed events in the context of HMMs.\n",
    "This algorithms is useful in many applications, including *speech recognition*, where the aim is to find the most likely sequence of words, given a series of sounds and the *reconstruction of bit strings transmitted over a noisy channel*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;32mdef\u001b[0m \u001b[0mviterbi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHMM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mev\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m\"\"\"\u001b[0m\n",
       "\u001b[0;34m    [Equation 15.11]\u001b[0m\n",
       "\u001b[0;34m    Viterbi algorithm to find the most likely sequence. Computes the best path and the\u001b[0m\n",
       "\u001b[0;34m    corresponding probabilities, given an HMM model and a sequence of observations.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mev\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mev\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mev\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mev\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minsert\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mev\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m# the recursion is initialized with m1 = forward(P(X0), e1)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mforward\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHMM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mHMM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprior\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mev\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m# keep track of maximizing predecessors\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mbacktracking_graph\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0melement_wise_product\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHMM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msensor_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mev\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m                                    \u001b[0;34m[\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melement_wise_product\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHMM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransition_model\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m                                     \u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melement_wise_product\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHMM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransition_model\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mbacktracking_graph\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melement_wise_product\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHMM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransition_model\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m                                   \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melement_wise_product\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHMM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransition_model\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m# computed probabilities\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mml_probabilities\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0.0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mev\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m# most likely sequence\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mml_path\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mev\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m# the construction of the most likely sequence starts in the final state with the largest probability, and\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m# runs backwards; the algorithm needs to store for each xt its predecessor xt-1 maximizing its probability\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mi_max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mml_probabilities\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi_max\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mml_path\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mi_max\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0mi_max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbacktracking_graph\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi_max\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0mml_path\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mml_probabilities\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%psource viterbi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Umbrella World\n",
    "---\n",
    "\n",
    "> You are the security guard stationed at a secret under-ground installation. Each day, you try to guess whether it’s raining today, but your only access to the outside world occurs each morning when you see the director coming in with, or without, an umbrella.\n",
    "\n",
    "In this problem $t$ corresponds to each day of the week, the hidden state $X_t$ represent the *weather* outside at day $t$ (whether it is rainy or sunny) and observations record $E_t$ whether at day $t$ the security guard sees the director carrying an *umbrella* or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Observation Emission or Sensor Probabilities $P(E_t := Umbrella_t | X_t := Weather_t)$\n",
    "We need to assume that we have some prior knowledge about the director's behavior to estimate the emission probabilities for each hidden state:\n",
    "\n",
    "| |  $yes$  | $no$ |\n",
    "| --- | --- | --- |\n",
    "| $Sunny$ |   0.10  | 0.90 |\n",
    "| $Rainy$ | 0.80 | 0.20 |\n",
    "\n",
    "#### Initial Probability $P(X_0 := Weather_0)$\n",
    "We will assume that we don't know anything useful about the likelihood of a sequence starting in either state. If the sequences start each week on Monday and end each week on Friday (so each week is a new sequence), then this assumption means that it's equally likely that the weather on a Monday may be Rainy or Sunny. We can assign equal probability to each starting state:\n",
    "\n",
    "| $Sunny$ | $Rainy$ |\n",
    "| --- | ---\n",
    "| 0.5 | 0.5 |\n",
    "\n",
    "#### State Transition Probabilities $P(X_{t} := Weather_t | X_{t-1} := Weather_{t-1})$\n",
    "Finally, we will assume that we can estimate transition probabilities from something like historical weather data for the area. Under this assumption, we get the conditional probability:\n",
    "\n",
    "| |     $Sunny$ | $Rainy$ |\n",
    "| --- | --- | --- |\n",
    "|$Sunny$| 0.70 | 0.30 |\n",
    "|$Rainy$| 0.30 | 0.70 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "umbrella_transition = [[0.7, 0.3], [0.3, 0.7]]\n",
    "umbrella_sensor = [[0.9, 0.2], [0.1, 0.8]]\n",
    "umbrellaHMM = HiddenMarkovModel(umbrella_transition, umbrella_sensor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from graphviz import Digraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"218pt\" height=\"332pt\"\n",
       " viewBox=\"0.00 0.00 218.30 331.60\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 327.5952)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-327.5952 214.2976,-327.5952 214.2976,4 -4,4\"/>\n",
       "<!-- I -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>I</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"31.2976\" cy=\"-292.2976\" rx=\"27.1222\" ry=\"27.1222\"/>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"31.2976\" cy=\"-292.2976\" rx=\"31.0965\" ry=\"31.0965\"/>\n",
       "<text text-anchor=\"middle\" x=\"31.2976\" y=\"-288.5976\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Start</text>\n",
       "</g>\n",
       "<!-- R -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>R</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"83.2976\" cy=\"-192\" rx=\"31.6951\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"83.2976\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Rainy</text>\n",
       "</g>\n",
       "<!-- I&#45;&gt;R -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>I&#45;&gt;R</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M45.7825,-264.3591C53.2052,-250.0422 62.1971,-232.6986 69.5166,-218.5807\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"72.7814,-219.8877 74.277,-209.3989 66.567,-216.6657 72.7814,-219.8877\"/>\n",
       "<text text-anchor=\"middle\" x=\"72.2976\" y=\"-231.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.5</text>\n",
       "</g>\n",
       "<!-- S -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>S</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"109.2976\" cy=\"-105\" rx=\"33.2948\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"109.2976\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Sunny</text>\n",
       "</g>\n",
       "<!-- I&#45;&gt;S -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>I&#45;&gt;S</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M21.7433,-262.0794C15.6507,-237.2103 11.1169,-201.6765 24.2976,-174 34.9509,-151.6305 56.7941,-134.1515 75.6692,-122.3825\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"77.5458,-125.3388 84.3574,-117.2238 73.972,-119.3198 77.5458,-125.3388\"/>\n",
       "<text text-anchor=\"middle\" x=\"33.2976\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.5</text>\n",
       "</g>\n",
       "<!-- R&#45;&gt;R -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>R&#45;&gt;R</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M111.9831,-199.8269C123.4881,-200.2691 133.1448,-197.6602 133.1448,-192 133.1448,-188.1086 128.5805,-185.6594 121.9892,-184.6524\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"122.1391,-181.1557 111.9831,-184.1731 121.8041,-188.1476 122.1391,-181.1557\"/>\n",
       "<text text-anchor=\"middle\" x=\"142.1448\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.6</text>\n",
       "</g>\n",
       "<!-- R&#45;&gt;S -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>R&#45;&gt;S</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M90.7759,-174.2292C93.0448,-168.4731 95.4288,-162.0202 97.2976,-156 99.5965,-148.5943 101.695,-140.4653 103.4717,-132.916\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"106.8989,-133.628 105.6835,-123.1032 100.0702,-132.0888 106.8989,-133.628\"/>\n",
       "<text text-anchor=\"middle\" x=\"109.2976\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.2</text>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"156.2976\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"156.2976\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Yes</text>\n",
       "</g>\n",
       "<!-- R&#45;&gt;Y -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>R&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M110.7113,-182.4727C135.7407,-172.1362 171.2122,-152.9757 187.2976,-123 201.0417,-97.3873 187.3021,-64.9049 174.0104,-42.9007\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"176.9457,-40.994 168.6151,-34.4478 171.0451,-44.7602 176.9457,-40.994\"/>\n",
       "<text text-anchor=\"middle\" x=\"201.2976\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.8</text>\n",
       "</g>\n",
       "<!-- N -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>N</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"83.2976\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"83.2976\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">No</text>\n",
       "</g>\n",
       "<!-- R&#45;&gt;N -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>R&#45;&gt;N</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M72.1427,-174.9693C64.0159,-161.4905 53.6917,-141.8589 49.2976,-123 45.6668,-107.4174 45.6668,-102.5826 49.2976,-87 52.8233,-71.8682 60.1668,-56.2389 67.1297,-43.6896\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"70.3231,-45.1595 72.3106,-34.7528 64.2672,-41.6487 70.3231,-45.1595\"/>\n",
       "<text text-anchor=\"middle\" x=\"58.2976\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.2</text>\n",
       "</g>\n",
       "<!-- S&#45;&gt;R -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>S&#45;&gt;R</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M90.3466,-119.9126C84.3564,-125.7707 78.4737,-133.0114 75.2976,-141 72.4293,-148.2144 72.4181,-156.386 73.6589,-164.0408\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"70.305,-165.0714 75.9762,-174.0206 77.1236,-163.488 70.305,-165.0714\"/>\n",
       "<text text-anchor=\"middle\" x=\"84.2976\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.4</text>\n",
       "</g>\n",
       "<!-- S&#45;&gt;S -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>S&#45;&gt;S</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M139.2431,-112.8442C150.8213,-113.1929 160.4446,-110.5781 160.4446,-105 160.4446,-101.165 155.8961,-98.7307 149.2909,-97.6971\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"149.417,-94.1989 139.2431,-97.1558 149.0403,-101.1888 149.417,-94.1989\"/>\n",
       "<text text-anchor=\"middle\" x=\"169.4446\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.8</text>\n",
       "</g>\n",
       "<!-- S&#45;&gt;Y -->\n",
       "<g id=\"edge9\" class=\"edge\">\n",
       "<title>S&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M118.8087,-87.3943C125.5213,-74.9689 134.628,-58.1119 142.1848,-44.1236\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"145.4467,-45.4493 147.1204,-34.9875 139.288,-42.1221 145.4467,-45.4493\"/>\n",
       "<text text-anchor=\"middle\" x=\"144.2976\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.1</text>\n",
       "</g>\n",
       "<!-- S&#45;&gt;N -->\n",
       "<g id=\"edge10\" class=\"edge\">\n",
       "<title>S&#45;&gt;N</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M103.9104,-86.9735C100.3308,-74.9958 95.5623,-59.0396 91.5133,-45.4912\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"94.8058,-44.2848 88.5889,-35.7057 88.0989,-46.2892 94.8058,-44.2848\"/>\n",
       "<text text-anchor=\"middle\" x=\"106.2976\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">0.9</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7fb3800bcb50>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dot = Digraph()\n",
    "\n",
    "dot.node('I', 'Start', shape='doublecircle')\n",
    "dot.node('R', 'Rainy')\n",
    "dot.node('S','Sunny')\n",
    "\n",
    "dot.edge('I', 'R', label='0.5')\n",
    "dot.edge('I', 'S', label='0.5')\n",
    "\n",
    "dot.edge('R', 'S', label='0.2')\n",
    "dot.edge('S', 'R', label='0.4')\n",
    "\n",
    "dot.node('Y', 'Yes')\n",
    "dot.node('N', 'No')\n",
    "\n",
    "dot.edge('R', 'R', label='0.6')\n",
    "dot.edge('R', 'Y', label='0.8')\n",
    "dot.edge('R', 'N', label='0.2')\n",
    "\n",
    "dot.edge('S', 'S', label='0.8')\n",
    "dot.edge('S', 'Y', label='0.1')\n",
    "dot.edge('S', 'N', label='0.9')\n",
    "\n",
    "dot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $[true, true, false, true, true]$ is the umbrella sequence for the security guard’s first five days on the job. What is the weather sequence most likely to explain this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import rounder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1, 1, 0, 1, 1], [0.8182, 0.5155, 0.1237, 0.0334, 0.021])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "umbrella_evidence = [True, True, False, True, True]\n",
    "\n",
    "rounder(viterbi(umbrellaHMM, umbrella_evidence))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}