{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hidden Markov Models" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import daft\n", "from matplotlib import rc\n", "from IPython.display import Latex\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(TeX macro definitions)\n", "$ %TeX macros\n", "\\newcommand{\\vct}[1]{\\boldsymbol{#1}}\n", "\\newcommand{\\mtx}[1]{\\mathbf{#1}}\n", "\\newcommand{\\rvar}[1]{\\mathsf{#1}}\n", "\\newcommand{\\rvct}[1]{\\vct{\\rvar{#1}}}\n", "\\newcommand{\\set}[1]{\\mathcal{#1}}\n", "\\newcommand{\\prob}[1]{\\mathbb{P}\\left[#1\\right]}\n", "\\newcommand{\\fset}[1]{\\left\\lbrace #1 \\right\\rbrace}\n", "\\newcommand{\\lsb}{\\left[}\n", "\\newcommand{\\rsb}{\\right]}\n", "\\newcommand{\\lpa}{\\left(}\n", "\\newcommand{\\rpa}{\\right)}\n", "\\newcommand{\\lbr}{\\left\\lbrace}\n", "\\newcommand{\\rbr}{\\right\\rbrace}\n", "\\newcommand{\\tr}{^\\mathsf{T}}\n", "\\newcommand{\\gvn}{\\,|\\,}\n", "$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model definition\n", "\n", "A Hidden Markov model is a probablistic model for discrete sequence / time-dependent data. It is formed of hidden (or latent) variables $\\rvar{z}_t$ which make up a Markov chain (i.e. $\\rvar{z}_{t+1}$ is conditionally independent of $\\fset{\\rvar{z}_s}_{s=0}^{t-1}$ given $\\rvar{z}_t$) and take values in a finite discrete state space $\\rvar{z}_t \\in \\set{H} = \\fset{1,2,\\dots H}$ and observed variables $\\rvar{x}_t$ which only depend on the current hidden state $\\rvar{z}_t$ and take values in some finite discrete state space $\\rvar{x}_t \\in \\set{V} = \\fset{1,2,\\dots V}$.\n", "\n", "The joint distribution over hidden and observed variables can therefore be factored as\n", "\n", "$$\n", "\\prob{\\rvct{x}_{0:T} = \\vct{x}_{0:T}, \\rvct{z}_{0:T} = \\vct{z}_{0:T}} =\n", "\\prod_{t = 0}^{T} \\lbr \\prob{\\rvar{x}_t = x_t \\gvn \\rvar{z}_t = z_t} \\rbr\n", "\\prod_{t = 1}^{T} \\lbr \\prob{\\rvar{z}_t = z_t \\gvn \\rvar{z}_{t-1} = z_{t-1}} \\rbr\n", "\\prob{\\rvar{z}_0 = z_0}\n", "$$\n", "\n", "This can also be represented as a graphical model as follows" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": [ "iVBORw0KGgoAAAANSUhEUgAAAkUAAADxCAYAAADFusOLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n", "AAALEgAACxIB0t1+/AAAIABJREFUeJzs3XeYFFXWwOHfqR6GKCgmDBhAMaCOiIJpXdOCGRfDrq6r\n", "iBnDrgn5UEwgLphds4hhzWlVFgXMCTAgYA4EMZFEBAnDMH3P98etwXEkDFDdt7rnvM/TDwMzVJ2e\n", "6qp76txQoqoYY4wxxtR1UegAjDHGGGPSwJIiY4wxxhgsKTLGGGOMASwpMsYYY4wBLCkyxhhjjAEs\n", "KTLGGGOMASwpMsYYY4wBLCkyxhhjjAEsKTLGGGOMASwpMsYYY4wBLCkyxhhjjAEsKTLGGGOMASwp\n", "MsYYY4wBLCkyxhhjjAEsKTLGGGOMASwpMsYYY4wBLCkyxhhjjAEsKTLGGGOMASwpMsYYY4wBLCky\n", "xhhjjAEsKTLGGGOMASwpMsYYY4wBLCkyxhhjjAEsKTLGGGOMASwpMsYYY4wBLCkyxhhjjAEsKTLG\n", "GGOMASwpMsYYY4wBLCkyxhhjjAEsKTLGGGOMASwpMsYYY4wBLCkyxhhjjAEsKTLGGGOMAaAkdAD5\n", "JiICNAEaxC8FFgILVXVByNhM7YlIfaAx/hiWAOXxa56qupCxmdoRkQy/PRcX48/FBaq6KGRspvZE\n", "pDHQEH8MBX8eLgTmq6qGjM3UjojUw19PGwKlwCJ+PYaVIWPLNynmz2ycAG0GtI9fOwM74U/eBfgD\n", "L/HfGwMzgTHA+/GfY1R1Rt4DN78hIs2AdvjjV3UsN8Ufw4VAFqjPrxflccTHL359oarZ/EduqohI\n", "CbANvx6/9kAZ/tgtxJ+L9fDHsBEwiV+P3/vAWFWdl//ITXUisgG/vZ62B9bi13MRfj2G86lxPQW+\n", "tUQpLBFpiD/3ql9P2/DrjWUFv15PS4GP+e319GNVrch/5PlRlEmRiGwCnAp0x1eCaiY605byfwRo\n", "xe8TqB+AO4AHVHVOXt6AqaoEHQn0wJ/A4/ntifn50u5gRGRN/HGrfsKvAzwM3K6qH+XlDRgARGQn\n", "/DE8GpjK7xOduUv5P/WAbfltArU98C5wG/CMqi7OyxswiMjaQDfgNGBtfp/oTKmZ6MTX0w35fQJV\n", "AQwC7lbVqXl6C3VeXJU9ADgD2Bf4jN9eTz9aWnU2rgLWTKA2B57Dn4tvFVuSWzRJkYhEwP74C/Af\n", "gAeBO1T1s9XYpgB74T9InYHHgdtUdfzqR2yWRkQ2w198u+MToduAoavTCIrIRsDJ+ER5YrzNp4v5\n", "bickEWmAT4J6ABvgbyoGq+r01dhmKXB4vM02wN34hvW71Y/YLI2IdMD/vg/HN4K3A6NXtRGMr6c7\n", "AKcDfwVexJ+Lrxdbw5oWIrIu/lp6OvAj/vf92OoMFRGRtYDj8Z+NinibD6rqL6sfcXhFkRSJyD74\n", "E7YcuBV4WFXnJ7yPFviG9TTgS+A0VZ2Q5D7qsvjkvRGffD6AT2i/THgf9YDD8CfztsCFwEN2QU5G\n", "3OidCvTDVxJuA55PuutSRNrib1SOBZ4GLlDVn5PcR10mItsBd/LbhPbHhPfRFPg7/lzM4q+no5Lc\n", "R10mIo2AK4GTgGfwN/PvJbwPwVedegD7AAOBawt+DJKqFuwLP0jzVuA74FDiJC/H+ywB/onPuv8B\n", "RKF/D4X+Ao7Ad61cCzTO0z53Bj4EngU2CP07KPQXfozXi8B7wPZ52mdTfOL1LXBg6N9Bob/ia1tv\n", "/NjKU4FMHvYpwFHx+X8N0DD076HQX8Ae+Bv3R4D18rTPVsDLwDvAtqF/B6vzKtgp+XF16EP8AOnt\n", "VXWIxkcnl1S1UlVvBHbDj3l5XUS2yPV+i5GIrCMijwL9gSNU9QJNuMK3LKr6Pr8mRuNE5Lj4zses\n", "BPFOxVeGXgF20zyN21LVuaraAzgBuE1EBsdjysxKiqtDo4G9gfaqepfmYXKCek/gu9U2AT4QkV1z\n", "vd9iJCKNROQ64Emgl6oeo3maKKSqk/DDV+7Ft4kXxZMrCk7BJUUiEonIAOA/wNmq2k1VZ+c7DlX9\n", "Cvgj/gM4WkT+lu8YCpmI7Al8hK/y7aiqI/Mdg6pWqGof4ECgJ/B4PDPD1IKIrAEMBU4B9lbVqzVA\n", "6VxVX8E3quXARyLSPt8xFDIROQ14Fd9l1llVv8l3DKo6U1X/AvQB/isil9hNSu2JyJbAWPzg9u1V\n", "9el8xxAnuHcAuwCdgDdFZL18x7G6CmpMUZx53gVsDRyiqj8FDglYMsZhGPAvVb01dDxpJyIH4JPa\n", "v6nqiNDxwJLZbvcBLYDDtEgGDeZKPCPpefxg+B4hkqGlEZE/4xv3I1X1jdDxpFmcdPTCj5XspKoT\n", "A4cELBm/+QLwBnCu2rpjyyUiO+Dbn8tV9a7Q8cCSz9aV+K7RP6nqt4FDqrWCSYriKYUPAc2BP+er\n", "m6W2RGRz/JiK21X1utDxpJWIHIafOXS4pmxgZfwZuxW/JtKfdCnTxc2ShOg1fFLUKx/d1itDRPbH\n", "L8FwjKq+HDqetBKRK4Gu+IToh9DxVBd3g/4PP3X81LR9xtJCRNrhE6KzVfXx0PHUJCLnA2cBfwxR\n", "gVwVBZEUxVnn3fiFGA9R1fKwES2diGwMvAn0V9W7Q8eTNvE4sMeAgzXhmRBJiT9rt+GrkQep6sIV\n", "/Jc6Je4yexnf3ZK6hKiKiFR1bR+iqu+EjidtROQC/MykvVR1Zuh4lkZEmuBvNN8GLkzrZy0UEdkK\n", "f3NyZojustoSkX/iZ4v+IV9jnFZHoSRFl+MXntpfU76qbTzo+nX83c3Q0PGkhYhsiz+Bj1bV18JG\n", "s3xxxehB/Gyco+1i7MVrgb0AfA2cnvbfi4gcAtwD7Kqqk0PHkxYicgxwNb6RSnW3hog0x19P74kn\n", "uBiWVGs/AC5T1fsCh7NCInIFcAiwu6b8ET6pT4pEZBd8GbVMl7ISdRqJyB+AR/ED3lIx7imkeCzY\n", "SPyF7c7Q8dRGPMbofWCAqj4YOp40EJEe+LVl9szHzKQkiMhF+LWv9rexKSAiG+Ifg3OAqn4QOp7a\n", "iBd0fQ9f1VrlxXiLiYg8BMxQ1XNDx1IbcQX+GfzK2ZeEjmd5Up0UxQ3TGOAqVX0kdDwrQ0RuApqr\n", "6t9DxxKaiPQC9sOPXUjvB66GeBbTC/iEvE4/kiAeM/cePiH6PHQ8tRUn5G8D96nq7aHjCSlumJ7D\n", "P17l0tDxrAwROQP/qJE90jKoPxQRORy/plOZFtBDzMU/N28cfljCmNDxLEvak6L++IdIdi2kxhSW\n", "PDNmPHCeqj4XOp5Q4pl5rwE7q+qUwOGsNBHpi5/ufXihfQaTEnebvYxfnfqa0PGsLBHZBj/Wb5e6\n", "3I0mIscD5+N/DwX1iJv4M/giMFxVB4aOJ5S42+wjfLf+W6HjWVnx0jW98O1BKrvRUpsUFWK3WU0i\n", "shd+VdE62Y1WiN1mNVk3WmF2m9UkIj35dVxinetGq9Zt1llVx4aOZ1VYN1rhdZvVVAjdaGlOit7A\n", "N6b3h45ldYjIbcA8Ve0ZOpZ8E5ET8A8j3LuQqywisjMwBNi00O6wV1c8A2gKPiEq2IYoHjz/LtBX\n", "VZ8JHU++ichdwCxV/b/QsawOETkbv1zGYaFjyTcR6YifvbttIXWb1RR3o32CLxZ8HzqemlK5orWI\n", "bA+0xq81UuiuA06soysl9wCuKeSECJY8EuQz4M+hYwngWOCNQk6IAOIK1/X4z2SdEq/5cxRwU+hY\n", "EnAPsHtcNapregC3FHJCBBCPz3wUv2ho6qQyKcKvaXCXqi4OHcjqileJfR84OnQs+RRXV9bHD1Qu\n", "BrdRxxrUuNR9Jv69F4MngTIRaRM6kDw7HhhWqMMQqosTggeA00LHkk8isg5wGP7ZYsXgduBUEakX\n", "OpCaUpcUiUhT4K/4xRqLRZ1rUPHv945CHYOyFM8CW4h/cGZdsTvQED/IuuDFAzsHA6eHjiVf4sS2\n", "B8WT2ALcAXSPx/vVFScCz6rqrNCBJEH9Q6MnAl1Cx1JT6pIi4Djg5bQtO7+angdaxNWTohcvuPZn\n", "fKm7KMRVy7vxVcy6ogdwW5ENTL4TOF5EGoUOJE/2ARYDBTdTaVlU9Uv8zN4jQ8eSD/HMuzMorsQW\n", "/COVUlcsSGNSdDzFVSWqGs8wCD+Dpy44Aj91NpWPD1gNdwN/iwftFrV4DNxhQEFPdKhJVb/Gd2cf\n", "EjiUfPk7cHehj+tbirvxbUVdsAfwC37mXTH5L7C9iLQMHUh1qUqKRKQUvybM26FjyYE3gV1CB5En\n", "HfBPuE6MiOwkIm4Zr7wM2ItnSswC6sKYlDLgS1WdHTqQHHgD/xmtCxI/F1PiTWCXuHuw2HXAT3Yo\n", "qsQ2nsn7DpCqHpRUJUVAW2Cyqs5PaoNpaExjY/GDPEvyuM9QdsavRJ6Y+JEEO1V7tQfmALOBfD4d\n", "eky872KX+DFM0blYJ45hvJzCZvjpz0luN/hxjAeNLwA2z8f+AmtPgudiGo5fNWNIWVKUtgY60YMP\n", "vjEVkZ2q/ZMArwCOPDamqjpHRL7Hr9D9Ub72m28i0gDYCvgw6W2r6rhq+xkANMMvxjc36X0tR1WD\n", "WuwLObYHRie5wbSci/hjuJOIREU2XqqmHYFPkp7Fm7Lj2B6YlMd9htAe/wDfRKTo+IE/hqkaV5S2\n", "SlF7fH9/olR1XNULP7OtGXBUnhtT8O+t2O9Qtwe+UtWFudpBfEJfiF9l+pWlfP/FXO2bunEMIQc3\n", "KJCOc1FVfwR+xq+FVsxycgwhHceROlDxi2djb4RfJy0xKTl+EB/DNHWDprFSlLM78GU1piLSCj84\n", "uKqL5i5VnZODEKpO4vtysO202IkcXYireQIYU3N1XhHZD9/Q7ZfDfX8AtCvmKkM8yHoLcljRXM65\n", "2A5fTl8TPwbvohw9r6yqbP9VDradFjvhx97kzHKO407AWvjj+Kf4+0kfxzHAPxPeZtq0Az7M1UNw\n", "l3P8RuCvdRPx58kp+LWhqs7Ln1Q1iaUtfgCyQEvgmwS2t9rSlhS1AL7N4faX2pgCj6vqzgAi8j5+\n", "ZkMuFlv8FuiYg+3mVfyAzZnxHXdNOT2GInInfpzE7xIfVX0ZeFlE7sjV/lV1dnxT0xg/I6QgxYvB\n", "rY2v6tVM7tbGX/Ry+cDG352LItIM/6DIu+O/74d/COgWOdj/t/jFRQtanER+uYxxmLm+nsKyr6kv\n", "AZup6tx4iY4nSH7sSFEcwxXIe5sYf6Z6Vg1XiK93a1adl/G/nZLEzlVVRaTqOFpStBQNgPJcbHhZ\n", "jWmcKS95WGs89mf/XMSAf28NcrTtfPoXcIiIzMbfTbyO71Yag39/OSnBisiR+DuWU+Op1aFUHceC\n", "TYrwd6DDgMUi8jl+HZtR+GPoyNF5CMtNbFsDF/HrkhxjgFYi0jQHZf1iOReHAWuLyA/4Kdtv4M/F\n", "cfiFN0McR4gTovjr2UAuZk6V499jQRORDvgnxw/DX1OaVEtA8t4mAs3jG8wqR+OfuVZdksNc0nUu\n", "qmpqXvjGtFkOtnsk/kJ/8jK+93iNf5sA7JiDODrhLw7F9qqMf78KzAMuzMHvbk38xXUC/iTev9pr\n", "8xo/63L8Of0lBb/zXLwWx3+WA9/l6He3zHMx/v5m1b7eH/8Q01zEcXkKft+5eFXEfzr8uKmOIY5j\n", "jZ8dAeybgxhaAt/x63W1Zfzvhfb3u2ocw8X4pOP2+OtHAx8/R7XzMgexvAh0ytX2V/aVtkpRJZDo\n", "wnjxwxDvxs9QmFyjCjQR3++dLyX4O7j/5HGfudCPXysl9fB9wh/hK0bb4GczJO1U/GDApviTqLoB\n", "QD6f/l0OXEthV4o6AH8B5uMveg2Ar4GR+PPihKR3uKJzUVUn628rgKfiK4O5UA//XL6XcrT9fLkO\n", "f92cj6+azMEv//E6cAw56A2ozXGMf25zfOM7QpcyISIB9fBJw1j8Z7mqO7/Q/v4f4CTgC6A5fszp\n", "UHxbkQESfZxJbY9f/LM7AZM0t5X5EvxxTIfQWVmNjPEHYKOEt9kTf9HPxn9Wf12NH2A9osb/+Ync\n", "VIqOAJ4O/XtO4H1ciO9C+zOwCSDVvtcXuDRwfLmuFP0MrBX6OKzme9gJv8r6KfHXpdW+txkwJQf7\n", "XO65WONnTwG65vD9XwdcEPo4JPA+7gJ646sPa9f43svAfiGPY/zzv7vGJhTHtsCnoY9BQu+lSfzn\n", "+kBU7d//Bjwc6vjFP3t7jt/7aGC30Meg6pW2StE8fDXg+6Q2qKoDgYHL+n48qKz5Uv7fuKX8+Opq\n", "in+PBU1Vr1nOt+cB6+YrlnyLF99siF84rmCpXwxzWQu1zQOaiohofNVKaJ/LPRerxAOsJ2puqgtV\n", "mlEEM89U9dTlfLvqepr0Pld0TW0FHFHtOvEy8ISIbKbJVhya4StkBU9V58V/Tq/xrcSPYW3Pw9if\n", "8N14uZSqdjFt6xR9jH/MR96o6tjqf49P6Fytc1OGf4/FLO/HsIqItBORnoCKyL/ixjVp2wBfa25n\n", "ZoU2C1gEbJzvHVdNfKhKiOLB9bmwA3Yu5srm+BmMVVoBsxNOiMC/t6JdCDcW7Hoa248cdjHHy39s\n", "BnyZq32srLRViqrW8Xk0z/s9RUQuxPex7kLuxjG0B4bkaNtpMQbYOekqQ23ECe5Yan8XtCpytiBe\n", "WqiqikjVOj65ntK9RHxD8n78ddU/TwSeTHg/pcB2+DEbxWwM0D3fO1XVl0VkzXja9k/4akMublAS\n", "fxRNCk0CmojIeqo6I187jY/dn/ADv08VkZdrFhASUgZ8nqabzDQmRT3zvdNqjSnAU7nYh/gnq++I\n", "n8JetFR1mogsxGf/k1fw44Wo6JOiWNUNyn/ztUNVnUR+qtdt8dW+1JTsc2QMcGuIHatq9etoTq6p\n", "+M/nPTnadirENygf4N/rC3nc7938ujRGLqXuepq27rMlzyQKHUgOtAFmaHE+dbymYl5+PyePokkh\n", "O4aF7xugnohsGDqQpIlIfWBrYHzoWPKgmB8tZEnR8qjqTPyU0lysYBvaLqTs4OdQ6p58nAQRqYfv\n", "389FGTltUvdMogTVhW4X4u7rojwX8efhBM3hMxZTZAy+/ShGqTsXU5UUxYaRm0dshHY0/r3VBcOA\n", "o4qw4ncIME7z/9DEEL4HZgD7hA4kSfF4osOB4aFjyRO7nha+V4G9RCSfa+rlnIhsi5+pnKpqXxob\n", "rduB0+Kpz0UhXsRsV/I/gDyUd/Fr+XQOHUjCegC3hQ4iH+Iqw23491xMugKfqOrnoQPJk/uAg0Vk\n", "vdCBJCWesdQNuDNwKHkR96AMxb/nYnIGMEhVK0IHUl3qkqJ4faBv8HflxeI04AFVLei1bWqrGBtU\n", "EdkKX7LP1aDRNHoQ2FdENgodSILqTGILEI9hfJoAs9By6GjgPVWdGDqQPLoNOKNYqu8i0gS/MOVd\n", "oWOpKa2/4KJpUEWkAf6ClOsFsNLmEWA3EdkscBxJOR24J01TR3Mt7iZ8hNwtUZFXIrI9/qGzz4WO\n", "Jc9uA06PZ8AWgx4EmlUX0Cj8grG5WNoghL8Br6lq3pb8qK20JkVPAmXx3XmhOwoYq6oFv3ruyoir\n", "Yg/gS6QFTUQaA8dTR8r1NdyOX6ekNHQgCegB3KWq6XnOUh6o6hhgOnBw6FhWl4jsDKxH3RlPBPym\n", "+n5W6FhWVzx5I7UV21QmRfHd+I3ATYU8+0VEmgFX4R9YWhf9GzhJRLYMHchquhJ4QVWnhA4k31T1\n", "Y/zskF6hY1kdIrIj/uGkdTGxBf+swmvi8TgFKa503QwMVNVs6HgCeAhoV+MBroXoZKACyOVjfFaZ\n", "5HnR4VqLpz+PBm5T1YJcoEtE7sYn+ct7PlFRE5F/4KtlfyzEC5mI7IGvXG6vqj+u6OeLkYhsjF+G\n", "YH9VTdVMkdqIq1zvAjeq6n2BwwlGRB7HP+j3wtCxrAoROQ84DNhXVV3oeEIQkc74xH6HQpwFKyKb\n", "4G+y9olvuFIntUkRLBkD8AqwUxr7HpdHRDrhVwTdvhA/vEmJBwa+BjytqjcGDmelxHfV44Feqvp0\n", "6HhCEpETgXOADoXW/SQil+PXQzk034+eSRMRWRf4EOiqqqNCx7My4qEUbwEd45XP6ywRGQRUqurp\n", "oWNZGXGvz3D8WKL+oeNZllQnRQAicgmwJ3BgoVzQ4m6zD4GTVTVXD5ctGCKyBX6g4O6FNLZKRK4D\n", "NlTVY0LHElp8Qfsf8I6qXhk6ntqKu81GAO1U9fvQ8YQWP2C3H/73URALH8bdZm8CD6vqLaHjCa1a\n", "+3KSqubsYa1Ji5+ndiqwm6pWho5nWQohKarqRntAVW8KHc+KxJWRB4FfVPW00PGkhYicg59xsE8h\n", "LE0Ql6nvxZep62S3WU3x1PyxwJ9V9e3Q8axI3Hi8BVyrqveHjictROQxYDZwRiHcaIrIZfhFROts\n", "t1lN8fXpLmBXVZ0aOp4ViRdqfJ0Ud5tVSeVA6+riUv2RwAUiclzoeJYnvpu+EdgEOC9wOGlzC/A5\n", "8FTaZzKJyO7Af4CjLCH6VVxpOQ54Oq7ApJaINMJXtl7Dz4I0vzoV6ABcETqQFRGRM/AzP4+1hOhX\n", "qjocP7ZohIg0Dx3P8ohIK3y32T/TnhBBAVSKqsSZ5ivAWar6ZOh4aooTov7AAfhs+OfAIaVOvEr5\n", "E4ADjknbSqYAIrILvjE9Pr7wmBriLpibgU5pvMjFCdHT+MeUdLPG9PfiFa7fBO5V1X+FjmdpRKQb\n", "0BfYS1UnBw4ndeI2ZyDwR+AAVf0pcEi/IyKb4tvta1W1INbqS32lqIqqfopPOG4WkVStzhp3md0C\n", "dMI3FJYQLUXcj/xX/Ofu2bjxSg0R2Qe/nP7JlhAtW3xTcj7wkoh0CB1PdXGX2TBgJtDdEqKlU9UZ\n", "wL7A8SLSP21Ln8SzVq8E/mQJ0dLFXZ898cntayLSInBIvyEiW+Nju7FQEiIooKQIljwCZG/gsvhE\n", "rh84pKoZHU8C2+P7vGcGDinV4jWojsLfxb8qIm0Ch4R4pwOPAUer6pDQMaWdqj4CnAT8T0SOT0Oj\n", "Gs9WfQM/CPWENA/mTIO4O3Qv/M3cf0RkzcAhISKNROQm4EzgD3XoGXWrJE6MLgAeB96Ou/6DE5FD\n", "8Q+yvURV/x06npVRUEkRgKp+CXQEtgXej1c4DSLuRvgQmIivEM0JFUshiRurE/GLkY0UkXNDPYIg\n", "fgzJS/hHseytqq+FiKMQqepQfIN6PvCciGwYIg4RqSciF+PL9P8GzrYKUe3EY+b2BuYAH4nIQaFi\n", "EZE9gXHAOvgZSnVusdRVoV4/fNXoKRG5NtQinSLSXEQewI+t/auqFtx4voJLigBUdRrwZ+BqYKiI\n", "XJXPqpGIrBsvhNYPv+bHhapanq/9FwNVdap6M7Ar/li+ns+qUbXq0Hv4Kdu7x120ZiXE1dtdgA+A\n", "cSLy93xWjeLq0Gh8xaO9qg4qhBlVaaKq81T1TPyA5ltE5N58Vo3i6tAN+GpHT1X9m6rOytf+i4Wq\n", "PoXvsdgYfy7mtWoUV4c+An7Gz9p9PZ/7T0pBJkWwJDt+GCgD2gIfiMgJucyQRWQdEemJrw5Nwa/1\n", "UVCLoKWNqk7A36k+jq8aDRCRzXO1PxHJxCfv6/jq0B9VdYB1taw6Va1Q1cuAzvhS/nAR6Sw5fKK3\n", "iLQRkRvx1aHb8QNNv8nV/uoCVX0V2AH/4NGPReSfIrJWrvYnIo3jtWvG459ntr2qPpOr/dUFqvqj\n", "qv4V6I2vGg0SkbJc7S++ufyDiDyFrw4dq6rnqOr8XO0z1wo2KapSrWp0IfAX4BsRuUZEWiex/fig\n", "7xqXBL8CtgEOiqtDBbH4WdpVqxrtDNQD3hOR/4nIQUk1rCKynoj8HzAJuBgYhFWHEqWqY/FVoyfw\n", "z9r6QkTOS2rKsIiUiMjhIjICP4BzIbCjVYeSU61qdAR+2v4kEblbRNoltQ8R2ToeN/QN/iG1p1t1\n", "KFnVqkZT8L0pb4nIsUn1qIjIGnGlfTz+WvoGBVwdqq5gpuTXVpwMnQ50wz9j5aX4zw9qO+ZHRDYA\n", "2sevw4Bm+LvR++zEzb14Vtpf8YMt18I3su/jj+Pk2jSAccWwDH8MqwaTPgXcrv6p4SaH4i60XfFP\n", "wz4EP6vvLfwx/Kg23c1xQtwafwx3Bo4GvgNuBZ6MB+2bHBKR9fED6k8HvscfxzHAmHgGW2220RzY\n", "CX8cO+Er+4OAu6y6l3vil0I5FH8u7oCvyr+LP45faC2eSSl+EeW2+PNwV6ArfiD1bcArxXRTUnRJ\n", "UZW4UewC7I4/GcuAqfgPwgT8XWY5vlrWAGiMP+jtgfrxz43BH/iXbOBm/sUN6y74pRiqktRG+OPy\n", "ITAXfxwr8cewAdAy/rktgM/in30P34jOzvNbMCyZoXkE/li2B9oAX+CPzQ/4Y7gIXyVsAKyJP193\n", "wg8ArjoXn4+rUSbP4oa1M35NnPb4YzMPf1w+w3e5LQQUaBi/2sQ/uw5+APUY4G3guTSuUVYXiH+G\n", "XBd+vZ6ujz82Y4Gf8G1iBb4NbIDv1myPbxu/5tdz8SlV/S7P4edF0SZFNcWzm7bGH+BN8Qe8N76E\n", "+wD+hK66UE8ppsy3mMR3ru2B7YAmwG7A/vhB9+XAdH6tRlglIYXiG5Yd8A3revhzsRfwX/xAzXn4\n", "pPcDW+IineIbllb4c3FLfBJ0BtAcv4htOb6regzwpd1UplM8ZqwdsCPQFH8+ngEMwLeJP+ETpnGq\n", "Oi9UnPlUZ5KipRERBd5T1VQtQGdqT0TOA65T1eDr5JhVF5+LR2kKV6s3tSMiPwAb2LlYuOIlSibX\n", "5WNY8AOtjTHGGGOSYEmRMcYYYwyWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOM\n", "McYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wx\n", "xgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHG\n", "AJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYA\n", "lhQZY4xBD8gdAAAgAElEQVQxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOM\n", "McYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wx\n", "xgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgCWFBljjDHGAJYUGWOMMcYAlhQZY4wxxgB1NCkSkSYi\n", "8mX8111EZFjQgMwqEZGbgOviryeISFngkMxKEpE9RWRi/NcnRKRv0IDMKhGRt4AN4q8/EZGSwCGZ\n", "lSQiZwEfxF9PFpE/Bw4piDqZFKnqPGDNav9UGSoWs1p+ArLx15sBE5f9oyalJgEt468rgBkBYzGr\n", "rnpb0lBV7ZpaeKYB9eOv1wO+CRhLMHUyKYqNjf+sAN4IGYhZZe8B8+Ovp8XJrikgqvoDsCD+azkw\n", "JmA4ZtW9Drj463dCBmJWWfVzrx7wcahAQqrLSdHrwGL8hfi9wLGYVTOGX+9s7BgWrvHxn42AcSED\n", "MavsXWAesBC7ySxUXwMafz1FVRcFjCWYupwUvY8/gRsS96OawqKq0/GVosXYhbiQvYG/GP+gqgtW\n", "9MMmlcbgqwuLsWpfQVJVBT6N/zoqZCwh1eWkaAzQFJilqnNCB2NW2Xj8xdguxIXrXUCwal8h+xY/\n", "vq8x8GHgWMyqq7q5fCtoFAHV2aRIVWcBM7HGtNC9jq8yjF3RD5rUqjoHXw8ahVllcZXhI+AbVS0P\n", "HY9ZZaPjP+tsuyj+s1zcRCQD7AMcKiJ7RFG0EVBPVZuKSCW+C6Yim81OwF+Y/6uq1simjIi0Ao4A\n", "9slkMtsBDfAzXRqKyBxgsar+6Jx7BxgKDK2r/eJpJSINgS7AAVEUdRSR5kCJqq4JzBeRcmBBNpv9\n", "CHgFeFJVvw0YsqlBRATYBfgzsGcmk2mNr9Y2UdVIRH4BKpxz36rqW8BzwBtaFxqbAiIi6wNdgf0z\n", "mUw7fJWvvqquUe16Osc59z4wAt8uFn2vStEmRSKyr4j0E5HtnXNNRESbNm2a3WyzzTItW7aURo0a\n", "Ub9+fZxzlJeXM3/+fCZPnuy+/fZbt2DBghJAM5nM7Gw2+zZwnqpOCP2e6hoRWRu4NpPJHOCcW09V\n", "o/r162c33HBDWrdunWnatCkNGjQgk8lQUVFBeXk5U6dOZeLEiZWzZ8+OnHNRFEULgQnOuWtU9T+h\n", "31NdEzegp0RR9E9gc+dcg0wm45o3b+622GKLkvXWW48GDRpQv359KisrKS8vZ86cOUyYMCE7bdo0\n", "Fi1alBGRbBRF07PZ7BDgorpwYU4bEWkrItdFUdQxm802A6Rx48aVm2yySWbzzTeXhg0b0qhRIwAW\n", "LlzIwoUL+eabb/Trr7/Ozps3L6OqEkXRXFUdr6q9VHVk2HdU94hIfeDKTCbzV+fchqpaUq9evex6\n", "663HlltumVlzzTVp0KABpaWlLFq0iEWLFvHjjz8yYcKEyh9//DGqrKyMRKRCRL5xzg0CrlXV7Ir2\n", "W2iKKikSkXrAZVEUnemcW7NVq1Zut912i/bYYw9atmy5wv9fxTnHZ599xujRo3n77bezM2fOzERR\n", "9LVz7grgfrvjyS0R6RRF0UDnXNkaa6yR3X333TMdO3akrKyM0tLSWm/nxx9/5O233+add95xH3/8\n", "sYhIuXPuQXzDOjt378DEd6HXisjRIlKvrKyMjh07yu67785aa61V6+1UVFTwwQcf8M477zBq1Kjs\n", "/PnzoyiK3nfOna+qb+buHZg4oe0RRdH/Oec2atGiRXaPPfbI7Lbbbmy55ZZEUe1HX0yePJmRI0cy\n", "atSo7JQpUzJRFP3onLsBGFCMDWuaiEiZiNyoqn+sX7++23XXXTMdOnSgQ4cONGjQoNbbmTt3LqNH\n", "j+a9997T999/X7PZrFPV/wHnqurXOXsDeVYUSZGIrCEij6tqp9LSUt1nn30yxx9/PE2bNk1k+999\n", "9x2DBg3SsWPHoqoVqjoIONuSo2SJyAVRFF3inGvWpk0b171796ht27aJbLuiooLHH3+coUOHZufN\n", "mxdFUfSBc+5oVZ2UyA4M4CsKURQ94pzbvmnTppVdunQp6dq1KyUlySxwPHbsWO67777spEmTMlEU\n", "/eScu1hV70hk4wZYMtzgbhE5TkRKOnTowMknnyzrr79+Itv/+eefGTx4MG+++abLZrMKDFHVY1V1\n", "YSI7MACIyFFRFN3knNtgo402yh533HGZPffcM5FtO+d44YUXeOKJJ7KzZs3KRFE00Tl3YjHcqBR8\n", "UiQiJ4jI3U2bNo1OOumkzD777JOzfVVWVvLEE0/wxBNPuGw2O8c510lV38/ZDusIEdkoiqLXVLV1\n", "p06dpHv37ktK8bkwfvx47rjjjux3330XAVer6sU521kdEVcVbgbO3HzzzfWMM86Ittlmm5ztb86c\n", "OQwePJhXX31VReRj59x+qjozZzusI0TkjyIypLS0tPGxxx4bHX744StVEVoZzjlefPFF7r333uyC\n", "BQsqVfVvqvpUTnZWh4hI4yiKnnfO7bXrrrvqqaeeKuuuu27O9vf1119z++2366effirAo8BxhVz9\n", "K9ikSETWiKLoJedch0MOOYRTTjklZydvTeXl5Vx55ZXZjz76KAPcB3S3qtGqEZFewFUbbrihXnXV\n", "VZl11lknb/seMmQIgwYNUuA759zeVjVaNXF16GURWfess86K9t9//7zte+rUqfTu3Ts7a9YsVPUf\n", "qnpr3nZeREQkIyKPquqRu+yyi+vVq1e0Ml3Vq8M5x80336wvv/yyRFH0mnPuIKsarRoR+YuIPNCk\n", "SZPMlVdemdliiy3ytu93332XAQMGuMWLF89T1UMKtWpUkEmRiBwrIvc1a9Ys6tu3b2azzTYLEsdb\n", "b73F9ddfX1U12kdVx6/4fxkAEVk3iqKRqtr6uOOOk6OPPjpIHD///DOXXHJJdsqUKRFwpapeHiSQ\n", "AiUi1wHntmnTxvXt2zeTywrf8tx///089dRTVVWjPVT1lyCBFCAR2V1EhpWWljbu1atXtPPOOweJ\n", "44svvuDyyy/Pzp8/v1JVj1LVIUECKUAiUi8uEuzVqVMnPfPMMyVfRYLqKioq6N+/vxszZkwkIg+r\n", "6nGFVjAouKRIRC4EBh500EF62mmnBTnw1ZWXl3PppZe6zz//XFV1P1W1tVZWQEQ2FZFP1l577QbX\n", "XHNNXqtDyzJkyBDuvvtuBR5yzv09dDyFIIqiIcDB55xzjuSzOrQsU6dO5cILL8z+8ssvvzjn2lh3\n", "2oqJyKHAM+3ateOSSy7JW3VoWZxz3Hjjjfrqq68CnKaqdwcNqACISMMoir6qX7/+Bv3794/yWR1a\n", "lnfffZf+/furqo5xznUopMSooJIiEbkK6H3qqady6KGHhg7nNwYOHKhvvvmmAofbHc6yicjWIjJu\n", "k002KbnxxhszSQ3ATcK4ceO49NJLFRjqnEvXByxFRESiKHpTRHYfMGCAbLXVVqFDWqKiooIePXpk\n", "Z8yYsUhVt1XVKaFjSisR+RvwnwMOOIAzzzxTQsdT3aOPPspDDz0E0FNVrwkdT1qJSLMoiiY2adJk\n", "zTvvvDPTpEmT0CEt8f3333POOee4ysrKCc65bQtlnFHBJEUi0hu46txzz2XfffcNHc5S3XLLLTp8\n", "+HAF9rWK0e+JSEsR+apNmzb1Bg4cGIWu8i3Nl19+yYUXXqiq+qRzLkyfXsplMpkXRWS/f//737Iy\n", "S13kS2VlJeecc072u+++K1fVTePV6001InIY8MwRRxwh3bp1Cx3OUv3vf//jzjvvBDjDZhj+Xlwh\n", "mtKsWbPmd911V2Zlptfny6xZszjttNPc4sWLP3XO7VAIFaP0tUpLISKnAFedccYZqU2IAM466yz5\n", "wx/+ICLysoiUhY4nTURkLRH5tGXLliVpTYgA2rRpQ79+/QQ4UkRuCx1P2kRR9DCw3/XXX5/KhAig\n", "pKSEm2++ObPuuus2iKLoq3gVbRMTkT2B/x544IGkNSECOOSQQ/j73/8OcJuIHBU6njSJq7WfNWrU\n", "qPkdd9yRyoQIYO211+b222+PMpnMtlEUFcRDu9PZMlUjIhsBd/zlL3/hoIMOCh3OCvXs2VO22mor\n", "iaLo1XiasgGiKHqpWbNmDW+66aZMWhOiKttvvz29e/cW4AwRCT9YJiVE5ChVPaZfv37SqlWr0OEs\n", "V0lJCbfffnumYcOGTUXk+dDxpIWI1BORYe3bt6dHjx6pvz4dffTRHHjggYjIwyLSLHQ8KXJvJpNp\n", "eccddwSb3FBb6667LjfddFOkqnuIyAWh41mRdLdOQBRFr2ywwQZ63HHHhQ6l1vr27RtlMplmwKDQ\n", "saSBiJyiqjv1798/VWOIlmfXXXelY8eOTkSeiRezq9NEpKGI/Gf//ffX7bffPnQ4tVJaWsoVV1yR\n", "UdW9ReSI0PGkgYg8Wb9+/YaXXHJJ6q/9VXr06CFrrrmmRFH0UuhY0kBEOgAnXHDBBVGzZoWRJ7Zs\n", "2ZJjjjlGgAEiskHoeJYn1SeGiFykqlv279+/oBqlBg0acN5550VAdxEJM781JeJus9u6dOmyUo9a\n", "SYNevXpF9evXbygiT4aOJbQoil5o3Lhxydlnn5366kJ1W221Ffvuu6+KyEN1vRtNRDqp6mEXX3xx\n", "VCg3J1X69u2bcc7tLCInhY4lpLjbbNgOO+zgdt9999DhrJRjjjmGFi1aaBRFr4aOZXlSmxTF3Wb9\n", "jzvuOEnDlO2Vteeee7LddttloygaUZe70aIoemmttdaSk04qvGtZSUkJvXv3jlT18LrcjSYiRznn\n", "/njFFVekvutzaf7xj39I48aNS+pyN1rcbfb0brvtpjvuuGPocFbapptuSpcuXRCRO+p4N9q9mUym\n", "WZ8+fQrvRAT69euXUdU2ae5GS+0vNoqiF1q0aKGhFvVLwmWXXZaJu9H+HTqWEETkr865nfr161dQ\n", "lb7q2rVrR8eOHV0URXXy8QPiPbDffvtpmzZtQoezSqIoqt6N1il0PIE8UL9+/YY9e/Ys2Bu0k08+\n", "mbgbbWjoWEIQka2Ju83SOrB6RdZff/3q3WiNQ8ezNKlMikSkpXNu+969exdsYwq+G+2oo46K6mrJ\n", "N4qiq9u3b6+F1m1WU8+ePSNVXaOOzoDpISKl55xzTsE2puBnFW699dYuiqJrQ8eSb/EjPI7s3r17\n", "wXWb1XT++ednnHO718VqkYjc2KJFi2yhdZvVdMwxx9CwYUMF+oWOZWlSmRQBNzRv3rxy8803Dx3H\n", "ajviiCMQkfoi0i10LPkkIm2cc5udfPLJBd2Ygh+wu91222kURak8iXMpiqLeHTt2DL5yfBK6d+8e\n", "Oee2F5HcPR0znS7IZDJR586dQ8ex2srKymjSpIkDBoaOJZ9EpJ6q/unYY48t6EJBlU6dOmWiKDo5\n", "dBxLk7orXVyu73LkkUcW9i1NrLS0lB133JEoii4LHUs+icgN6623XuXGG28cOpREnHzyyZFzro2I\n", "FHbZayWIyPbOuQ1POumkgk9sAbbZZhuaNm1aCVwXOpZ8iqLovD333LMoEluAgw8+OBNFUV17FM8l\n", "paWlus8++4SOIxHHHXccqtpYRA4PHUtNaTxL/hlFUXTwwQeHjiMxp5xyijjnNhOR8A+lyQMRyahq\n", "52OOOaYoEluAVq1a0bx580rg+tCx5IuI3Ljhhhtm119//dChJKZLly4l8ZPEiyLRWxER6eCcW697\n", "9+5F836PPvpoVLVB/JiSOiGKorP33nvvoqgSgR9asu2222oURf1Dx1JT6pKiKIou3HXXXYvmrgZg\n", "4403Zt11181Sd+5QL6pXr56kefXxVdG1a9cSEekSOo58iBPbvf/2t78VzYUYoGvXrohIPeDE0LHk\n", "g4hcu8kmm2TXWmut0KEkprS0lLKyMq0r1XcR2dM5t9YJJ5wQOpRExdX3bdLWnZ2qzENExDnXomvX\n", "rkVzV1Nljz32yGQymT1Dx5EPItJl6623LqrEFuDQQw9FVevFs0CK3T4iInvttVfoOBJVUlLCpptu\n", "6oDCnda6EkRkp3333beoEluAgw8+OFLVdC+rnpzjmzdvXtm0adPQcSRqiy22oLS0NAv8JXQs1aWt\n", "1doFkC22KL5epg4dOpDNZovndm05RGTrtm3bJprY3nfffRx22GFMnDjxN/9+yimncOqppya5q2WK\n", "ooiGDRtWAnVhdeRDmzZtWhBPtV5Z22yzTSaTyRTeYj0rSUQyzrnGhT5baWl22mknVDVTF8b4RVG0\n", "W6tWrYousQU/RR9I1RpwaUuK/tykSZPKJCsMaWhMAdq2bQu+GFbUD4qNq31Nk74Qd+vWjbKyMi65\n", "5BLmz58PwC233MKMGTPo27dvovtano022igC9s7bDgMRkT0222yzRC/EaTkXd955Z5xzqSrZ50in\n", "KIp0gw2Se6pCWo5haWkp9evXzwJH5m2n4bTacccdE7vJTMsxBGjTpk0mk8nslNedrkDakqI/tGzZ\n", "MtELcVoa0yiKaNSoUSXQNW87DWNPEdFcLKfQq1cv1lhjDf71r3/x9ttvM2LECC666CLyORB4m222\n", "iTKZTFEntuCXVNhuu+0Srfal5Vxs164dqhrVgYkPhyRd7UvLMQRo0aIFQHENXKxBROo75xoleZOZ\n", "pmMY36Ck6lloqUqKMplM26S7XSAdjSnAxhtvHAF/zOtO869LkyZNctLt0rhxYy666CLGjx/PgAED\n", "OOCAA1jWxaJPnz65CKGqG3TtnGw8JeJq3xp77LFH4ttOw7lYUlJSVWUo6m7QKIp233zzzRPvdknD\n", "MQTYaqutUldlyIEDoyhy666bbGEzLcdw5513RlVL0vSQ2FRNmXbONd1hhx0S325VY3ruuecyfvz4\n", "pTam06ZNY+TIkbRu3ZqJEyfSuXNnGjdOdhXyVq1aRRMnTiz2u9MdNthgg5z1f7du3Zrdd9+dUaNG\n", "ccQRv2/Txo8fz9SpUxk/fnxO9h9/PiMRWUtVZ+dkJ+FtB+TkAb61ORcnTpzIhAkTmD9/Pl9++SXd\n", "unWrqgokZp111tHvv/++faIbTRkR2XTLLbdM/CazNsew6vjNmzePcePGccQRRyR+DNu2bctLL71U\n", "eA/GXDl7xItVJlrAqM0x7NOnD61bt2aDDTbgq6++YsSIEZx55plLzss11liDM888c7XiaNCgASUl\n", "Ja6ysvKPwKOrtbGEpKpSpKrSrFluVm+vakxFZKmN6YABA+jatStlZWV07tyZf/87+ceVxUlWaeIb\n", "TpdGDRs2zNnswbfffpuRI0dSVlbGv/71r999v6ysjAMOOCBXu6faeLdiHjTfXEQ0Vxtf3rk4f/58\n", "JkyYQOfOnenatSsHHnggl156aeIxNGjQIAIaJb7hFFHVeknf2FVZ0fW0T58+bLnlluyxxx5sscUW\n", "DBgwIPEYmjRpgqqmqg3LgcalpblpMpZ3DCdOnMiJJ55It27d6Ny585JjWXVe9urVi6QmREVRpEBq\n", "ptal7QMljRrl5jq1vMZ0woQJNGnSZMnfGzdunJNKQ/369QHqJb7hFBGRRvXq5eYtTps2jQEDBnDW\n", "WWdx0UUXMW3aNG699dac7Gt54oShmJOiZvGFKieWdy5OnTqVp5769dm7rVu3Ztq0aSxYsCDRGEpL\n", "SyOgYaIbTZ+SENdTgHvuuYeqfVe/tiapcePGdSEpapKr59Ut7xjOmzePVq1+XfHgrbfeYscdfzth\n", "M6mkKJPJWFK0PLlYaHZFjen06dN/11W2xhprMGnSpETjyGQyqGrRrcFUXfyYlpxsu6qc26lTJxo3\n", "bkyvXr0YPnw4w4cPz8n+VqAop8jGcnZdWNG5uMUWW/xmsGfVDUvSjXtc8SvmYwiEuZ4Cvzlew4YN\n", "o1u3bonHUWzroC1DTi6mKzqGZWW/nUsyfvz43yVFrVu3TjKk1JyLaftU6cKFCxPf6Ioa019++SXx\n", "fS7NwoULEZHFedlZIM65hYsXJ/8Wq2ZH9OrVa8m/lZWV0bVrV2677bbEE9jliRPbn/O2w/ybm6vk\n", "vTaJbfXBnsOGDeOss85KPI5FixY5INnyU/pky8vLE99obW9Opk2bxtNPP027du1+18gmYf78+eSy\n", "mzclFlRWVia+0ZW5wZwwYQItWrTI2SBs55wAc3Oy8VWQuqQo6QSlNo3pGmussWRqYpVcJErxBaqo\n", "kyJg4aJFixLf6FlnncWzzz77uxOzW7duPPvss78p9ebJnHzvMI/mxBeqRK1sYjt8+HD22muvZc4w\n", "XB0VFRUKJH8HliIisjjpbseVOYYtWrSga9eurL/++jmZDRq/t6JcYLSaxJOilT0Pl1YlSlI2mwWY\n", "l7MdrKRUJUVRFC388ssvE91mbRrTFi1aMG/e749J0g3tN998o8657xPdaPpMnD59evK3NilR7YIx\n", "I2QcOfYpwKxZsxLd6MoktuPHj6dFixY5SYgAZs2apUCyF5uUUdUZX3/9daLbrM0xrKoQVSkrK2P8\n", "+PFMnz490Vi++uoroijKT5k/nI/mzZuXaDu9sjeY48aNy1lS5JyjsrIyA7yXkx2sglQlRao64cMP\n", "P8x7ObRm3+i0adNy8iH4+uuvs6r6duIbTpfnf/7552Cfq4kTJ/LUU08hItx3332JD5gfOXIkURTN\n", "UdWiLdur6sIoispHjRoVZP9V44iqulzefjvZU8Y5x/z580uA/ya64ZRxzr07YcKEvN+gTJ8+nblz\n", "f+0NmTZtGk2aNEm8++XTTz912Wz240Q3mj7/raysjJZ2054vuawUffjhhwBOVT/PyQ5WQarWKVLV\n", "tydPntyWAHGdddZZPP3006y//vp89dVXnH322Ynv45dffskAzyW+4XR5wTkXzZw5k6QXHKuN1q1b\n", "07p166VOE07CJ598oqr6WU42ni5fjx8/futDDjkkrzudNm0a55133m/+rUWLFiS5kORXX30FoKTo\n", "7jRHhs2aNevYfO+0rKyMefPmMXz4cJo0acK4cePo169f4vv57rvvHPB64htOEVWdHUVRxejRo0v3\n", "3z+/jwgbPnw4Y8eORUQYNmwYZWVlSQ+u5t133yWTySRbkl5NqUqKgP/NmTPnjBA7rmpMgUQvwFWm\n", "TJlSNUD31cQ3niKquiiTySwcOXJkwy5duoQOJ3FTpkzJqupboePINefc6AkTJmxBnq8RLVq04Lnn\n", "cnvfMGrUKDKZzM+VlZVFW+2LPZfNZqM5c+aQq/XflqX6NTQX11PnHAsWLCgBnl7hDxc4Efnugw8+\n", "aJXvpKhz58507tw5p/v47LPPXDab/TCnO1lJqeo+A0Y452Tq1Kmh40hc3O3ySzF3u1QzKVcrSoc2\n", "b968DPBs6DjyYNjs2bPTdn1IxKeffqrZbPaT0HHkmqr+IiKLQnWD5tJnn30GoKpanBeaapxz7331\n", "1VdFOU7z+++/T121L1UXPVXNRlE0Z+jQoaFDSdzo0aOzqvpR6DjywTn36qefflp0s0Leeecd4py2\n", "2MeFQVxl+OKLL0LHkbiJEyc6YEToOPJBRCa88cYbLnQcSRs+fDhRFBXzZIfq/jt9+vRMLqbmhzRr\n", "1iwWLlxYAjwWOpbqUpUUATjn7h8xYkRRNahz585l0qRJGVW9JHQseXLp/Pnzo48+Kq4c8IEHHsiK\n", "yMi6UO2LB1t/cs899xRVg/rGG29QUVEhwMDQseSDc67vxx9/LBUVFaFDSYxzjrfeess55/K/nH0Y\n", "jwMVTz75ZOg4EjVo0CCNomi6qqZqFmjqkiKg98KFC6N33303dByJuf/++4miaLaqFvV4oirx4MAP\n", "Bg8eXDQN6uzZs/nmm28yqnrein+6ODjnLvr888+jpNe6Cemhhx7Kishrqpr8YloppKqPicjCRx55\n", "JHQoiXnllVeIx4P1Dx1LPqj3yJAhQ4qmWOCcY9SoUeqcS/6heKspdUmRqs4XkVEPPPBA0XwAXnvt\n", "taxzLvknzKaYc+7CCRMmRNWn5haywYMHaxRFM1S1eLL1FVDVoSLyy4MPPhg6lETMnDmTH374IaOq\n", "/wwdSz455+5/4YUXiuZ6+sgjj1QCw1S1aN5TLVwwd+7czCefFMdQuKFDh+Kcc8CNoWOpKXVJEYCq\n", "nj9lypTMzz8X/pMUXnnllapyffJzUlNMVV+Nomj2Aw88EDqU1RaX69U5d33oWPLNOXfXiy++WBSN\n", "T1yu/6GujO2r5v/mz58fjR07NnQcq+37779nxowZJXUtsVXVWVEUjb/33nuLovr+5JNPZlX12TQO\n", "RUhrUjQ6iqLp1113XUF/AJxzDB48OCsiw1W12B/v8TvOuRtfeuklLfTul4cffphsNuuAa0PHEkCf\n", "8vJyefbZwp5wN3v2bEaNGoVzrk7dnACo6hwReeeWW24p+OT2+uuvz0ZRNElVJ4SOJd+ccz2/+OKL\n", "KOlVyvNt9OjR/PTTTxng3NCxLE0qkyIA59yR48aNi955553Qoayyu+66S+fOnetU9S+hYwmkr6r+\n", "eNlllxVscjtt2jQef/xxVdUL61i5HvADroGrBg8erIVcub344ouzIjJFVW8PHUsIqnr4jBkzeOih\n", "h0KHsspeeuklvvzyy8g5d1joWEJQ1RFRFL3Xp0+frO95KjwVFRUMHDjQiciDqvpt6HiWJrVJkaq+\n", "JSKPDBw40BXizInJkyczdOhQUdVTVLXYn8+zVKqqzrn9Pv/8c3nllVdCh7NK4sb0M1VNXd93vqjq\n", "pcA3F198cUEmhc8++yzffvtt5JzbJ3QsoajqdKDnY489pkk/gywfFixYwC233OKAW1W1OAbWrALn\n", "3J/mzJnjBg0aFDqUVdK3b19XWVn5i6oeHzqWZUltUgSgqn9fvHjxvH79+hVUWuyco0+fPtkoit5R\n", "1ftDxxNSPH7j9ptvvtkVWjfagw8+yMyZM3HO7Rs6ltCcc/t+8803UaF1o82ePZvBgwcrvmr5deh4\n", "QlLV60Xki0JMbvv06ZNV1ZnAOaFjCUlV56jqqUOGDKHQutFGjx7NuHHjIlU9KI1jiaqkPSnKqurB\n", "Y2N9gsUAAArFSURBVMeOjUaPHh06nFq78847de7cuc4596fQsaTEWao669JLLy2Y5Hbq1KlV3WYX\n", "xHfZdZqqTgL6DR48WGfNStWjipYrTgCmqOploWNJA+fc3tOnT6eQZhSOGDGiqttsvzQ3pvmiqvcV\n", "WjfaggULuOaaa6q6zUaGjmd5Up0UwZJutDv79++vhbC67rPPPsvzzz8vqtqtrnab1VTVjfbFF19w\n", "ww03pP6iNnv2bM4++2wnIu/X5W6zmuJutC969OjhQj61u7b69evnvvvuO63L3WY1xQn+2Y899hiF\n", "0KU9duxYbrnlFgUG1OVus5ribrSKCy64IPVVv4qKCs4444xsZWXlj2nuNquS+qQIwDl3OvC/iy66\n", "SKdMmRI6nGV65ZVXiPt6z1fVh0PHkyZxN9qB8e8otYnRvHnzOP30093ixYsnO+c6ho4nbZxzO5SX\n", "l08//fTTs+Xl5aHDWaabbrpJ48ey/KGud5vVFA8273/DDTeQ5gr8F198wWWXXabAf1T1/0LHkyZx\n", "N9oOEyZMcH369EltuaiyspKzzz47O3v27PnOuTaFUOmTAohxiUwm85qI7DVgwADZaqutQofzG0OH\n", "DuWOO+4AuKoOPc5jpYnIX4BHDjjgAM4880wJHU91s2bNokePHq68vHyqc27zuriMQm2ISOMoiqY0\n", "adJkzTvvvDPTpEmT0CH9xsCBA/XNN99U4GBVHRY6nrQSkVtFpMe5557LPvukq5g2btw4Lr30UgWG\n", "OOe6hI4nrUSknYi8u80220RXX311FEXpqXNUVFTQo0eP7IwZM8pVdUtVLYgnvRdUUgQQRdEQ4ODe\n", "vXvLrrvuGjocwA/IfeyxxxT4P1VN3bLlaSMihwLPdOzYkd69e6fiRJ48eTIXXHCBq6ysnOica2sJ\n", "0fLFidEX9evX3+CGG26INtpoo9AhUVlZyaWXXuo+/vhjVdX9VDVVT99OIxH5F3DRiSeeSNeuXUOH\n", "A8Crr75a1c3+kHPu76HjSTsRaSsiH7Rs2TJz3XXXZRo0aBA6JGbPns0555yTnTt37i9xhWhm6Jhq\n", "q+CSIoAoiu5Q1dPatWvnLrnkkqi0tDRIHDNnzqR3797ZadOmCXCqqt4TJJACJCJ/EJEXGzVqVHLF\n", "FVdkQlX+nHPceeed+vzzz0sURW845/YuhBJvGohIvSiKRqnqTl27dpVu3boFi2X8+PH069fPLVq0\n", "aGHcZVb4yzfniYicB1y7ySabuKuuuiqz5pprBoljwYIFXHHFFe7TTz8V4BpVvShIIAVIRDaNomhM\n", "FEVrnX/++dGee+4ZLJZnnnmGe++9V4FJzrl2hTa2tiCTIgAR2VNEhtarV69Jz549o44d8zv84/HH\n", "H+fBBx9UEfkqbkgLojSYJiLSUESGqOp+++67r/7jH/+QfFaNJk+eTJ8+fbLxApsnqep/8rbzIiIi\n", "Z4jIv9dee2369euXyWfVqLKykoEDB7pRo0ZFIvKcqh5pVb6VJyKtoih6FWjZvXt36dIlvz1Wb7zx\n", "BjfccINzzv3knOtkSe3KExEBBgHd27Ztm7388svzWjWaPXs2F198cfbbb78V/BIYl+dt5wkq2KQI\n", "/IdARB5S1WN23HFHd95550VrrbVWTvc5adIkrr766qrqUC9VvSanO6wDRKSriDzcqFGjkgsuuCCz\n", "884753R/5eXl3H333TpixAiJomh0fBEuqLuZtBGRtaMoekVVtz/kkEOkW7du5LqCO3LkSG644Yaq\n", "6lBXVR2R0x3WASLSD+i9ySabuF69emVatmyZ0/3NnDmTa665xn322WcC3A2cbpXa1SMiO0dRNDyK\n", "ojVPO+206IADDsjp/pxzPPbYYzz66KMKTHHO7VPIkxsKOimqIiJ7RlH0pHNu/U022SR7wgknZDp0\n", "6JDY9p1zDBkyhKeffrryp59+Komi6CPnXGerDiVHRBpGUfSMc+5PjRo1cgceeGDm2GOPTbRhnTBh\n", "AoMGDXKffvqpiMh851wPqw4lS0TOFJFrgAY77LCDnnzyydFmm22W2PYXLFjAgw8+yP+3d3+vVddx\n", "HMdf7893G/PsyJAxhgarwfwxONJUEFY4AlEqZIGmF92U/0BYEFEQEVQIhdBdlJfebhM1iobUZBcD\n", "SR34C3ehaNOByUQ5J8/a9/3uYmckXaVMt47Pxx/wPd8vH75fnud9OJ/vyMhI/uDBg4Xp0F6mQ4un\n", "NjX60d3Xtbe35/v27ct27typxZzijo2N6ciRI/nU1FSWUvrd3QeYDi2e2tToe0nvNDY2Wn9/f9q/\n", "f79aW1sX7TOmp6d1+PDhOH36tCLir4g4WA/7gdVFFC0ws61mdigiXlqxYoXv2LEj27Vrl1avXv3I\n", "x3J3TU5Oanh4OMbHx8PdPSKOSnp/ub6zpR6Y2UpJB1NKb0dEoVQqxe7du1Nvb68aGhoe+XgzMzMa\n", "GRnR8ePH87t372YppUvu/nFtLfGEmNlbKaXP3L27ra0tHxgYyLZv3/5YD+XZ2VmdPXtWg4ODfvny\n", "5WRm9939O0mf1N7NhifAzLolHTKz17Mss23bttnAwIB1d3c/1vFu3LihEydO6OTJk3m1Wk1m9mtE\n", "vBcRE4t75lhgZo2SPkopvevubV1dXfmePXuyvr6+x/rCWalUNDo6qqGhoXx6enohaL+U9G29TPjq\n", "KooWmFmLpC9SSu+4e6uZRWtra97V1ZVt3LjRNmzYoGKxqEKhIHdXuVzW/fv3NTExoYsXL8b169fz\n", "crncIClSSrfc/StJ39TLov9fmNnelNLn7t4tKTU3N8+tWbMm9fT0pE2bNmnVqlVqaWlRU1OTyuWy\n", "KpWKrl69qomJCU1OTs7NzMykPM9TSulPdz+q+f2jmO49RWb2gqRDKaXX3L05yzJva2vztWvXNvT2\n", "9qqzs1MtLS0qFAqqVquqVCq6c+eOzp07p0uXLuU3b95UtVrNzCw3s8vu/mFE/LDEl/VMMbNM0gcp\n", "pQPu3mFmKhaLc52dnVmpVLJSqaRisahisSgzU7lcVrlc1oULF3T+/Pm4du1afu/evSwiLKU0Uwva\n", "TyOiutTX9iyp/bnla0mbI6Khqakp7+jo0Pr167MtW7aovb1dhUJBzc3NqlQqqlQqunXrls6cOaMr\n", "V67kt2/ftrm5uWRms5JORcSBetxQsy6j6GG1G3qHpF1m9rKZrY2IFRFhkhb2yQkz85TSvTzPL0g6\n", "JWlY0m+E0PJgZuskvSnplSzLXnT3tohIemgNNR+xs5Kuufu4pJ8kHWOasDzUpoADkl5NKW2V9Ly7\n", "N2p+DR++F/OU0h95np+V9IukwdprRrDEaj/L9El6w8z6zawnIlb++3mq+XuxHBFXImJM0jFJoxGx\n", "7HdgfhaY2XOS9kjanmXZZnfviIhM/2zoHJq/F+fMbMrdT0v6WdJQRMws0Wk/FXUfRQAAAP/F0u+a\n", "BwAAsAwQRQAAACKKAAAAJBFFAAAAkogiAAAASUQRAACAJKIIAABAElEEAAAgiSgCAACQRBQBAABI\n", "IooAAAAkEUUAAACSiCIAAABJRBEAAIAkoggAAEASUQQAACCJKAIAAJBEFAEAAEgiigAAACQRRQAA\n", "AJKIIgAAAElEEQAAgCSiCAAAQBJRBAAAIIkoAgAAkEQUAQAASCKKAAAAJBFFAAAAkogiAAAASUQR\n", "AACAJKIIAABAElEEAAAgiSgCAACQRBQBAABIkv4GINkG4JL0W1wAAAAASUVORK5CYII=\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# code for plotting HMM graphical model using daft\n", "rc(\"font\", family=\"serif\", size=20)\n", "gm = daft.PGM([5, 2], origin=[0.5, 0.5], grid_unit=4, node_unit=2)\n", "for t in range(4):\n", " gm.add_node(daft.Node('z_{0}'.format(t), r'$\\mathsf{{z}}_{0}$'.format(t), t+1, 2))\n", " gm.add_node(daft.Node('x_{0}'.format(t), r'$\\mathsf{{x}}_{0}$'.format(t), t+1, 1, observed=True))\n", " gm.add_edge('z_{0}'.format(t), 'x_{0}'.format(t))\n", " if t > 0:\n", " gm.add_edge('z_{0}'.format(t-1), 'z_{0}'.format(t))\n", "gm.add_node(daft.Node('z_T', r'$\\mathsf{z}_{T}$', t+2, 2))\n", "gm.add_node(daft.Node('x_T', r'$\\mathsf{x}_{T}$', t+2, 1, observed=True))\n", "gm.add_edge('z_T'.format(t), 'x_T')\n", "gm.add_edge('z_{0}'.format(t), 'z_T', plot_params={'linestyle' : 'dotted'})\n", "gm.render()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The joint distribution factors are parametrised in terms of matrices $\\mtx{A}$ and $\\mtx{B}$ and a vector $\\vct{\\pi}$ as\n", "\n", "\\begin{align*}\n", "\\prob{\\rvar{z}_t = i \\gvn \\rvar{z}_{t-1} = j} &= A_{ij} \\qquad &\\forall i \\in \\set{H}, ~j \\in \\set{H}\n", "\\\\\n", "\\prob{\\rvar{z}_0 = i} &= \\pi_i \\qquad &\\forall i \\in \\set{H}\n", "\\\\\n", "\\prob{\\rvar{x}_t = i \\gvn \\rvar{z}_t = j} &= B_{ij} \\qquad &\\forall i \\in \\set{V}, ~j \\in \\set{H}\n", "\\end{align*}\n", "\n", "To ensure the parameterd define valid normalised probability distributions it is required that\n", "\n", "\\begin{align*}\n", "\\sum_{i\\in\\set{H}} A_{ij} &= 1 \\qquad \\forall j \\in \\set{H} \n", "\\qquad &\n", "A_{ij} \\geq 0 \\qquad &\\forall i \\in \\set{H}, ~ j \\in \\set{H}\n", "\\\\\n", "\\sum_{i\\in\\set{V}} B_{ij} &= 1 \\qquad \\forall j \\in \\set{H} \n", "\\qquad &\n", "B_{ij} \\geq 0 \\qquad &\\forall i \\in \\set{V}, ~ j \\in \\set{H}\n", "\\\\\n", "\\sum_{i\\in\\set{V}} \\pi_i &= 1\n", "\\qquad &\n", "\\pi_i \\geq 0 \\qquad &\\forall i \\in \\set{H}\n", "\\end{align*}\n", "\n", "It is also required that observations lie in correct state space i.e.\n", "\n", "$$\n", "x_t \\in \\set{V} \\qquad \\forall t \\in \\set{T} = \\fset{0,1, \\dots T}\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def check_params(A, B, p, x):\n", " assert all(A.sum(0) == 1.) and all(A.flat >= 0.), \\\n", " 'Columns must sum to one and elements be non-negative'\n", " assert all(B.sum(0) == 1.) and all(B.flat >= 0.), \\\n", " 'Columns must sum to one and elements be non-negative'\n", " assert p.sum(0) == 1. and all(p >= 0.), \\\n", " 'Elements must sum to one and be non-negative'\n", " assert A.shape[0] == A.shape[1] and A.shape[0] == B.shape[1], \\\n", " 'Inconsistent A and/or B shapes'\n", " assert min(x) >= 0 and max(x) <= B.shape[0], \\\n", " 'Observations must be integers in range 0-{0}'.format(B.shape[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inference (forward and backward algorithms)\n", "\n", "Define \n", "$$\\alpha^{(t)}_i = \\prob{\\rvar{z}_t=i, \\rvct{x}_{0:t} = \\vct{x}_{0:t}} \\qquad \\forall i \\in \\set{H}$$\n", "\n", "then $\\alpha^{(t+1)}_i$ can be calculated recursively from $\\alpha^{(t)}_i$ using\n", "\n", "\\begin{align}\n", "\\alpha^{(t+1)}_i &= \n", "\\prob{\\rvar{z}_{t+1}=i, \\rvct{x}_{0:t+1} = \\vct{x}_{0:t+1}} \\\\\n", "&=\n", "\\prob{\\rvar{x}_{t+1} = x_{t+1} \\gvn \\rvar{z}_{t+1} = i} \n", "\\prob{\\rvar{z}_{t+1}=i, \\rvct{x}_{0:t} = \\vct{x}_{0:t}} \\\\\n", "&=\n", "\\prob{\\rvar{x}_{t+1} = x_{t+1} \\gvn \\rvar{z}_{t+1} = i} \n", "\\sum_{j \\in \\set{H}} \\lbr\n", " \\prob{\\rvar{z}_{t+1}=i \\gvn \\rvar{z}_{t}=j}\n", " \\prob{\\rvar{z}_{t}=j, \\rvct{x}_{0:t} = \\vct{x}_{0:t}}\n", "\\rbr \\\\\n", "&=\n", "B_{x_{t+1}i}\n", "\\sum_{j \\in \\set{H}} \\lbr A_{ij} \\alpha^{(t)}_j \\rbr\n", "\\end{align}\n", "\n", "or in matrix-vector notation\n", "$$ \\vct{\\alpha}^{(t+1)} = \\vct{B}_{x_{t+1}} \\circ \\mtx{A} \\vct{\\alpha}^{(t)} $$\n", "with $\\circ$ the elementwise product and $\\vct{B}_{x_{t+1}}$ the column vector formed by taking the transpose of the $x_{t+1}^\\text{th}$ row of $\\mtx{B}$. The recursion is initialised using\n", "$$ \\vct{\\alpha}^{(0)} = \\vct{B}_{x_{0}} \\circ \\vct{\\pi} $$\n", "\n", "Similarly define\n", "$$ \\beta^{(t)}_i = \\prob{\\rvct{x}_{t+1:T} = \\vct{x}_{t+1:T} \\gvn \\rvar{z}_t = z_t} \\qquad \\forall i \\in \\set{H} $$\n", "\n", "then $\\beta^{(t-1)}_i$ can be calculated recursively from $\\beta^{(t)}_i$ using\n", "\n", "\\begin{align}\n", "\\beta^{(t-1)}_i &= \n", "\\prob{\\rvct{x}_{t:T} = \\vct{x}_{t:T} \\gvn \\rvar{z}_{t-1} = z_{t-1}} \\\\\n", "&=\n", "\\sum_{j \\in \\set{H}} \\lbr\n", " \\prob{\\rvct{x}_{t:T} = \\vct{x}_{t:T}, \\rvar{z}_t = j \\gvn \\rvar{z}_{t-1} = z_{t-1}}\n", "\\rbr \\\\\n", "&=\n", "\\sum_{j \\in \\set{H}} \\lbr\n", " \\prob{\\rvar{x}_{t} = x_{t} \\gvn \\rvar{z}_{t} = j} \n", " \\prob{\\rvct{x}_{t+1:T} = \\vct{x}_{t+1:T}, \\rvar{z}_t = j \\gvn \\rvar{z}_{t-1} = z_{t-1}}\n", "\\rbr \\\\\n", "&=\n", "\\sum_{j \\in \\set{H}} \\lbr\n", " \\prob{\\rvar{x}_{t} = x_{t} \\gvn \\rvar{z}_{t} = j}\n", " \\prob{\\rvar{z}_t = j \\gvn \\rvar{z}_{t-1} = z_{t-1}}\n", " \\prob{\\rvct{x}_{t+1:T} = \\vct{x}_{t+1:T} \\gvn \\rvar{z}_t = j}\n", "\\rbr \\\\\n", "&=\n", "\\sum_{j \\in \\set{H}} \\lbr B_{x_t j} A_{ji} \\beta^{(t)}_j \\rbr\n", "\\end{align}\n", "\n", "or in matrix-vector notation\n", "$$ \\vct{\\beta}^{(t-1)} = \\mtx{A}\\tr \\lpa \\vct{B}_{x_{t+1}} \\circ \\vct{\\beta}^{(t)} \\rpa $$\n", "The recursion is initialised using\n", "$$ \\vct{\\beta}^{(T)} = \\vct{1} $$\n", "\n", "These two recursive update formulae for $\\vct{\\alpha}^{(t)}$ and $\\vct{\\beta}^{(t)}$ are called the forward and backward equations respectively (due to the direction of the recursion in each). Most conditional and marginal probabilities that might need to be inferred for a HMM can be expressed in terms of $\\vct{\\alpha}^{(t)}$ and $\\vct{\\beta}^{(t)}$ terms and the recursion formulae used to efficiently calculate them." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Functions for doing forward and backward passes given parameters and observations are" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def forward(A, B, p, x):\n", " check_params(A, B, p, x)\n", " alphas = np.empty((len(x), A.shape[0])) * np.nan\n", " alphas[0] = B[x[0]] * p\n", " for t in range(1, len(x)):\n", " alphas[t] = B[x[t]] * A.dot(alphas[t-1])\n", " return alphas\n", "\n", "def backward(A, B, p, x):\n", " check_params(A, B, p, x)\n", " betas = np.empty((len(x), A.shape[0])) * np.nan\n", " betas[-1] = np.ones(A.shape[0])\n", " for t in range(len(x)-2, -1, -1):\n", " betas[t] = A.T.dot(B[x[t]]*betas[t+1])\n", " return betas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem set up\n", "\n", "Consider problem on 3-element hidden state space $\\set{H} = \\fset{0,1,2}$, two element observed state space $\\set{V} = \\fset{0,1}$, with three observed time points $\\set{T} = \\fset{0,1,2}$. Model parameters (note 0 indexing for hidden/observed states and times) and observed states are as follows\n", "\n", "$$\n", "\\mtx{A} = \\lsb\n", "\\begin{array}{ccc}\n", "0.5 & 0.0 & 0.0 \\\\\n", "0.3 & 0.6 & 0.0 \\\\\n", "0.2 & 0.4 & 1.0\n", "\\end{array}\n", "\\rsb\n", "~~\n", "\\mtx{B} = \\lsb\n", "\\begin{array}{ccc}\n", "0.7 & 0.4 & 0.8 \\\\\n", "0.3 & 0.6 & 0.2\n", "\\end{array}\n", "\\rsb\n", "~~\n", "\\vct{\\pi} = \\lsb\n", "\\begin{array}{ccc}\n", "0.9 \\\\\n", "0.1 \\\\\n", "0.0 \n", "\\end{array}\n", "\\rsb\n", "~~\n", "\\vct{x}_{0:2} = \\lsb\n", "\\begin{array}{ccc}\n", "0 \\\\\n", "1 \\\\\n", "1 \n", "\\end{array}\n", "\\rsb\n", "$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "A = np.array([[0.5, 0.3, 0.2], [0.0, 0.6, 0.4], [0.0, 0.0, 1.0]]).T\n", "B = np.array([[0.7, 0.3], [0.4, 0.6], [0.8, 0.2]]).T\n", "p = np.array([0.9, 0.1, 0.0])\n", "x = [0,1,1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do forward pass to get $\\vct{\\alpha}^{(t)} ~\\forall t \\in \\set{T}$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{align*}\\vct{\\alpha}^{(0)} &= \\lsb 0.6300~~0.0400~~0.0000 \\rsb\\tr\\\\\\vct{\\alpha}^{(1)} &= \\lsb 0.0945~~0.1278~~0.0284 \\rsb\\tr\\\\\\vct{\\alpha}^{(2)} &= \\lsb 0.0142~~0.0630~~0.0197 \\rsb\\tr\\end{align*}" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "alphas = forward(A, B, p, x)\n", "# code for pretty printing\n", "Latex(r'\\begin{align*}' +\n", " r'\\\\'.join([r'\\vct{{\\alpha}}^{{({0})}} &= \\lsb {1} \\rsb\\tr'\n", " .format(t, '~~'.join(['{0:.4f}'.format(a) for a in alpha])) \n", " for t, alpha in enumerate(alphas)]) + r'\\end{align*}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do backward pass to get $\\vct{\\beta}^{(t)} ~\\forall t \\in \\set{T}$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{align*}\\vct{\\beta}^{(0)} &= \\lsb 0.2143~~0.1696~~0.1600 \\rsb\\tr\\\\\\vct{\\beta}^{(1)} &= \\lsb 0.3700~~0.4400~~0.2000 \\rsb\\tr\\\\\\vct{\\beta}^{(2)} &= \\lsb 1.0000~~1.0000~~1.0000 \\rsb\\tr\\end{align*}" ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "betas = backward(A, B, p, x)\n", "# code for pretty printing\n", "Latex(r'\\begin{align*}' +\n", " r'\\\\'.join([r'\\vct{{\\beta}}^{{({0})}} &= \\lsb {1} \\rsb\\tr'\n", " .format(t, '~~'.join(['{0:.4f}'.format(b) for b in beta])) \n", " for t, beta in enumerate(betas)]) + r'\\end{align*}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question: Find $\\prob{\\rvct{x}_{0:2} = \\vct{x}_{0:2}}$\n", "\n", "$$\n", "\\prob{\\rvct{x}_{0:2} = \\vct{x}_{0:2}}\n", "= \\sum_{z\\in\\set{H}} \\lbr \\prob{\\rvar{z}_2 = z, \\rvct{x}_{0:2} = \\vct{x}_{0:2}} \\rbr\n", "= \\sum_{z\\in\\set{H}} \\lbr \\alpha^{(2)}_{z} \\rbr\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$\\prob{\\rvct{x}_{0:2} = \\vct{x}_{0:2}} = 0.0969$" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P_x = alphas[2].sum()\n", "# code for pretty printing\n", "Latex(r'$\\prob{{\\rvct{{x}}_{{0:2}} = \\vct{{x}}_{{0:2}}}} = {0:.4f}$'.format(P_x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question: Find $\\prob{\\rvar{z}_1 = z_1 \\gvn \\rvct{x}_{0:2} = \\vct{x}_{0:2}}$\n", "\n", "\\begin{align*}\n", "\\prob{\\rvar{z}_1 = z_1 \\gvn \\rvct{x}_{0:2} = \\vct{x}_{0:2}}\n", "&= \n", "\\frac{\\prob{\\rvar{z}_1 = z_1, \\rvct{x}_{0:2} = \\vct{x}_{0:2}}}{\\prob{\\rvct{x}_{0:2} = \\vct{x}_{0:2}}}\n", "\\\\\n", "&= \n", "\\frac{\n", " \\prob{\\rvar{z}_1 = z_1, \\rvct{x}_{0:1} = \\vct{x}_{0:1}}\n", " \\prob{\\rvar{x}_2 = x_2 \\gvn \\rvar{z}_1 = z_1} \n", "}\n", "{\\prob{\\rvct{x}_{0:2} = \\vct{x}_{0:2}}}\n", "\\\\\n", "&= \\frac{\\alpha^{(1)}_{z_1}\\beta^{(1)}_{z_1}}{\\prob{\\rvct{x}_{0:2} = \\vct{x}_{0:2}}}\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$ \\prob{\\rvar{z}_1 = z_1 \\gvn \\rvct{x}_{0:2} = \\vct{x}_{0:2}} =\\begin{cases}0.3609 \\qquad : ~ z_1 = 0\\\\0.5804 \\qquad : ~ z_1 = 1\\\\0.0586 \\qquad : ~ z_1 = 2\\end{cases}$$" ], "text/plain": [ "" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P_z1_gvn_x = alphas[1]*betas[1] / P_x\n", "# code for pretty printing\n", "Latex(\n", " r'$$ \\prob{\\rvar{z}_1 = z_1 \\gvn \\rvct{x}_{0:2} = \\vct{x}_{0:2}} =' +\n", " r'\\begin{cases}' +\n", " r'\\\\'.join([r'{0:.4f} \\qquad : ~ z_1 = {1}'\n", " .format(p, z1) for z1, p in enumerate(P_z1_gvn_x)]) +\n", " r'\\end{cases}$$')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }