{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we implement the Viterbi algorithm as described in these notes:\n",
    "[nvit.pdf](https://courses.cit.cornell.edu/info2950_2017sp/nvit.pdf)<br>\n",
    "`T[i,j]` is the i->j transition probability, `e[i][a]` is the emission probability for observable `a` from state `i`, `w` are the initial state weights, `snames` is a list of the names to call the states (instead of the integers `i,j`), and `O` is a string of observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def viterbi(T,e,w,snames,O):\n",
    "  states=range(len(snames))  #label states as integers\n",
    "  #initialize probability for first observation:\n",
    "  v = [w[j]*e[j][O[0]] for j in states]\n",
    "  path=snames   # initialize paths to state names\n",
    "\n",
    "  for o in O[1:]:\n",
    "    #  nv[j]=v[i]*T[i,j] array with probabilities for arriving from i to j:\n",
    "    nv = [v*T[:,j] for j in states]\n",
    "    #  update maximum probability to state i for next step:\n",
    "    v = [max(nv[j])*e[j][o] for j in states]\n",
    "    #  path[j] keeps track of max prob path to get to state j:\n",
    "    path = [path[np.argmax(nv[j])] + snames[j] for j in states]\n",
    "\n",
    "  return path[np.argmax(v)],max(v)  #return path to max probability final state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First consider the simple case of three urns 1,2,3 and three observations RBR as in the above notes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('121', 0.016666666666666663)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T=np.array([[.3,.6,.1],[.5,.2,.3],[.4,.1,.5]])\n",
    "e=[{'B':1./2,'R':1./2},{'B':2./3,'R':1./3},{'B':1./4,'R':3./4}]\n",
    "w=np.ones(3)/3.\n",
    "state_names=['1','2','3']\n",
    "viterbi(T,e,w,state_names,'RBR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now return to the \"slightly dishonest casino\", defined by this HMM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQEAYABgAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0a\nHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIy\nMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCACxASsDASIA\nAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQA\nAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3\nODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWm\np6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEA\nAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSEx\nBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElK\nU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3\nuLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+gnA\nyelFR3EXn20sQbaXQrn0yMUAYOmeMrDVL+C3jtryGG73fYruaNRDeYBJ8shiegJG4DIBIyKt6z4g\nh0ia3tls7u/vbhWeK0s1UyMi43N8zKoA3LySOoxmuP0W1v7mDwdocmmXltN4edGvZ5YWSE+XbyQj\ny3IxJuLg/L0Gc4PFXteubia80zX9PtdWtLhYp7UyHTWuCqGRMo8AIf5igZXHAAyeCMgGrP4ztBYW\nVzZafqGovdyyQrbW0aCVGjzvDh2UDaQQeetbWm3j39hFdSWVzZO+c290FEiYJHO1mHOM8E8EV5ra\n+HvsGk2k+v2OvziW7vLtm0+WQSxGaQECRIMMSQATtyFOQfWu38HxalB4djj1NpzIJZfJFy26YQeY\n3lCQ93Cbck8+vOaAN2iiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooA\nKKKKACs7Wdat9EtY5Zo5p5ZpBDBbwKGkmkIJ2qCQM4BPJAABJNaNcz4sguI7vQtYhtprmLS71pp4\nYELyGNoZI9yqOWILg4HOM4z0oA1dG1q31u1klhjmglhkMM9vOoWSGQAHawBIzgg8Egggg1nWnjKw\nvL6GFbW+jtrmUwW19JCBBPIN2QpzkfcbBYAHsTkVS8PzTW19qmrT6ffxx63qcf2eI27b40WBEEkq\n9YwTG3Xplc8nFc1baJqtzf6NokDanb6Zpeox3QhlshEsMcbMyp54O2VT8qhVAIB+Y0AeqUUUUAFF\nFFABSMwVSzEBQMknoKWvGPE2qap8U/GVx4L0G6e00CwIGrXsZ5lOcFB6jIIA7kEngUAbmvfGrQtP\nv20zQrS78Q6kMjyrBcpkdt+Dnr1UNVFfHnxPuAs0Hw32QnnbLdYfGffBH5V33hrwlonhLT1s9HsY\n4FwA8mMySn1ZuprboA8li+NUukXKW/jPwlqeibmCrcBTJEeM55Az/wAB3V6bpOsadrunx3+l3kN3\nayfdkibI+h9D7Hmp7uztr+1ktby3iuLeVdskUqBlYehB4NeLeIvDuo/CDWG8WeE1eXw9IwGpaWWJ\nEYJ+8vt6Hqp9QTgA9voqjo+rWeu6PaapYSeZa3UYkjbvg9j6EdCPUVeoAKKKKACiiigAooooAKKK\nKACiiigAooooAKKKKACiiigAooooAKKKKACiimTSx28Mk0zqkUal3djgKAMkmgBl3d21haS3V5PH\nb28SlpJZWCqoHck9K8w1D43WVxfPYeEdB1HxFcoQC0CFIxzjOcFse5AHvWFDDqPxy8STTzzTWngf\nT5dkcaZRrxx6/oSewIA5JI9l0rR9O0OwSx0uyhtLZOkcKBR9T6n3NAHmn/CdfFAHzT8OF8nGdouh\nv6fX+lSaZ8b9Pjvo9P8AFmiaj4cu3OAblC0XXHXAYfXbgeteqVQ1jRNM8Qae9jq1lDd2z9UlXOPc\nHqD7jmgC1b3EF3bx3FtNHNDIu5JI2DKw9QR1qWvDh9v+B3ie3ja4lu/BGpzbAJDuazkPP+J46gHu\nK9vR1kRXRgyMAVYHII9aAHUUUUAc54+1x/DfgTWdVicJNDbERMe0jfKh/wC+mFYnwc0CPQvhxp0h\nANzqK/bZ5McsX5X8l2/r60z42wyTfCXWfLUsUMLkAZOBMmT+HWul8GXEV34G0CeGNY430+ArGvRP\n3a/KPp0/CgDcooooARmCqWYgKBkk9BWRear4ev7KezutS02W3njaOWNrhCGUjBB59Kvan/yCrz/r\ng/8A6Ca+UfhT8NbT4h/2v9q1Gez+w+Tt8pA27fvznPps/WgD1L4Ma1a6GPEfhS91KAQ6XfsbSWWd\nAJI2JHynOD93dx/fr1X/AISDRf8AoL2H/gSn+NeP/wDDNWk/9DDe/wDfhP8AGj/hmrSf+hhvf+/C\nf40Aewf8JBov/QXsP/AlP8aP+Eg0X/oL2H/gSn+NeP8A/DNWk/8AQw3v/fhP8aP+GatJ/wChhvf+\n/Cf40Aewf8JBov8A0F7D/wACU/xo/wCEg0X/AKC9h/4Ep/jXjh/Zw0QTrAfE10JmUuI/KTcVBAJx\nnoCRz7in/wDDNWk/9DDe/wDfhP8AGgD2D/hINF/6C9h/4Ep/jVu1vLW9iMlpcw3EYO0tE4cA+mR3\n5FfPvir4B6b4e8K6nrEWuXc0lnbtMsbQqAxA6E5rqP2cf+Seah/2FZP/AEVFQB7BRRRQAUUUUAFF\nFFABRRRQAEgAknAFZ3/CQaL/ANBew/8AAlP8auXP/HrN/uH+VfJPwq+G9p8Q5tVS61Caz+xLEVMS\nBt2/d1z/ALtAH1V/wkGi/wDQXsP/AAJT/Gj/AISDRf8AoL2H/gSn+NeP/wDDNWk/9DDe/wDfhP8A\nGj/hmrSf+hhvf+/Cf40Aewf8JBov/QXsP/AlP8aP+Eg0X/oL2H/gSn+NeP8A/DNWk/8AQw3v/fhP\n8aY37OGipKkTeJboSPkohiQFsdcDPOKAPY/+Eg0X/oL2H/gSn+Ned/GrxbawfD2ax0y/t57rUpkt\nMQTKxCHLNnnoQu3/AIFWH/wzVpP/AEMN7/34T/Gj/hmrSf8AoYb3/vwn+NAHpPhZvDvhjwxp+jW2\nq6eEtYQrH7QnzP1Zjz1LEn8a6G1vbW+jMlpcw3CKcFopA4B9MivAfEvwA0zQvDGqatHrl3K9layT\nrG0KgMVUnB59q6H9m/8A5EbU/wDsJN/6KjoA9kooooAwvGfh6HxT4Q1LR5VBM8LeUxGdkg5RvwYC\nuX+COvNrnwzsklbdNp7tZNxjhQCn/jjKPwr0RmVELMQFUZJPYV5P8AXNx4U1q9WNkhudYlkj3DGQ\nUT/9X4UAes0UUUAUda0qDXNEvtKuR+5u4Hhf2DDGfw615r8GNelsrW88B6ywi1jRpXWNGbPmRE5+\nXPUAn/vkrXrFeffEL4dy+I7m21/QLv8As7xNYjMFwp2iUDorn8wDz1IOQaAO+lmigTfLIka5A3Ow\nAyelPr578b+OvEd74B1Pw54u8H6hZ3rxp/ptuuYCyOrBicEAZXnDH8K8l0Tx34p8OhV0rXbyCNRh\nYi++MD2Rsr+lAH2pqf8AyCrz/rg//oJrwn9mX/maf+3T/wBrVteBfGnxE8UaTcS6noVkdLMDf6dJ\nutmYbTkqPm3+2AB71i/sy/8AM0/9un/tagD07xdb3SeKPCd2upXItn1RYWsgEERPkTtvJxuJ4HBO\nPbPNdRCL7+0LoztbmyKp9nVFYSA8795JwR93GAO+c0Xmm2l/LZyXUXmPZz/aIDuI2SbWTPB5+V2G\nDxzRDp9tBqF1fRowuLpUWVjIxBCZ24UnA+8egGe+cCgC1XG+Lre6TxR4Tu11K5Fs+qLC1kAgiJ8i\ndt5ONxPA4Jx7Z5rsqq3mm2l/LZyXUXmPZz/aIDuI2SbWTPB5+V2GDxzQBwviy7ums9Yv7XzD5t7a\naPHsl8smMyqJSrD7pZpGTd1G0elaHg6NdL1zU9Gk0yHTLhIIbn7PZXJmtSjtIA67kUq5KsGGOcKa\n6X+xdNOlS6YbONrKUu0kLDIYuxZic9yxJz6mmaRoOnaFFIlhA6GQgvJLM8sj46ZdyWIHYE8UAY/x\nL/5Jn4j/AOvGT+VcT+zj/wAk81D/ALCsn/oqKu2+Jf8AyTPxH/14yfyrif2cf+Seah/2FZP/AEVF\nQB7BRRRQAUUUUAFFFFABRRRQBFc/8es3+4f5V4B+zP8A8fXiX/ctv5yV7/c/8es3+4f5V4B+zP8A\n8fXiX/ctv5yUAen/ABBtLiT+wrpdRuY4ItZsFe0QJ5cpNynLHbu47AEDjpXVgX39qMS1v/Z/kjao\nVvN83ccknONu3HGM5zzRf6baanFFHeReYkM8dwg3FcSRsHQ8HswBx09aBp9sNUbUgjfamhEBbzGx\nsDFgNuduck84z78UAcv4ut7pPFHhO7XUrkWz6osLWQCCInyJ23k43E8DgnHtnmtHyn1bxHqbw3Lw\nfYrYWUU0agtFLIBJIw3AqSF8nqD0Na95ptpfy2cl1F5j2c/2iA7iNkm1kzwefldhg8c0+2tILMSi\nBNvmyNK+WJLMxyTz/kdBQB5nbyXVt8P7q2kvriUyeJmsbm6dsSNE1+I3JIxjKkjjGM8Yro/D1nba\nJ431bRdLiW30tLC1uRbR8JFK7zK20dtwRSR7Z7mt9NB0xNOu9PFojWl3LLNPE5LB3kYu55PGWJPH\nTtijSNC07QoJItPgZBI26R5JXldz0G53JY4HAyeKAMv4hf8AJOPEn/YNn/8AQDXAfs3/APIjan/2\nEm/9FR13/wAQv+SceJP+wbP/AOgGvAvh14i8a+HvhzqM/hfRIL63F+xnnbMjxHy06RggkYxzz34o\nA+oqZHNFKXEciOUO1grA7T6H0r4u1v4l+MfEG5b7XrvymJ/dQN5KYPYhMZH1zXefC3xxqfhvwbLp\nWg+FNR1fUbq8aUTLGfIUlVUZIB4G3nOOp5FAHrHxc8XDw54TksLNi+s6tm1s4UPz/NwzAewOB7kV\nt/D/AMNf8Ij4I0zR32/aIo99wVOQZWO5ue4BOB7AVyngj4e6q/iA+M/HFwt1rzj9xbDBjsx2Axxk\nZ4xwOTkk5r1CgAooooAKKKKAMzxHpTa54Z1XSVdY3vbSW3V2GQpZSAT9Cc1xHg34KeGfCxiuruP+\n1tRTB865UeWjeqR9Bz0JyR616VRQBU1IY0m8A/54P/6Ca8K/Zl/5mn/t0/8Aa1e7an/yCrz/AK4P\n/wCgmvCf2Zf+Zp/7dP8A2tQB6P4tBhvbl7rXdRiuJYwNJ0/S5GWUuBguUUfP8zDO7KAAZxzXTaRN\nqzxiLVLSKNkt4SZ45Q3mylf3g24+UKQMcnOe2KzL7wd9r1651i31/V7C4uY0idbUwbdqZwBviY4y\nScZ6k1sWGnSWU0kj6je3ZeKKPbcMpAKAguAqjDNnLdsgYAoAvVyXii51m08TeGfIv4o9LudQEEtu\nkREjt5Mzcvuxt+QfLt69+1dbVDUtIt9Un02ad5VbT7oXcQQgBn2OmGyDxiQ9Mc4oA4/xVq11anVd\nStGlDpNa6PbvEgLxmSRTK6Z4LYkQDPG6OtDwfdOmo6hpk51e3mijinFlqsqzyIrFx5iyq77lYqRt\nz8pU9Aa2H8N6fLok2lSiV4ZZmuGcyHzBK0hl3huxDnIx0wPSjR/D8OkTTXLXl3f3s6rHJd3jKZCi\nklV+VVUAFm6AdTnNAGZ8S/8AkmfiP/rxk/lXE/s4/wDJPNQ/7Csn/oqKu2+Jf/JM/Ef/AF4yfyri\nf2cf+Seah/2FZP8A0VFQB7BRSM6ohd2CqBkknAFchq/xU8E6IWW68QWryKxUx2xM7A+h2A4/HFAF\nrx3NLF4dijSR4o7m/s7ad0baRFJOiPz2ypIz71V8PWdtonjfVtF0uJbfS0sLW5FtHwkUrvMrbR23\nBFJHtnua5DUPjn4S1VJ9Mt9E1fV4pEIdI7ZdrL9C2705x6Vm2/jmyh02fT7T4aeKXtZpEkldmmaW\nRlIZd0hy5wQON2McdOKAN/XraC9tPHmt3Zxqmiu402cn5rQJaxyJs9NzsSfXODkcVDeab9l1fVvF\nms+G9J1K0W7jlN0bnddW6JFEmEjCEZVlY43A89M9ce/+I2h3mopqWs/D3xHHcRkBm8twjbeV3rlV\nfHbcDjtU9p8TfhTq2p/2jeQ3VldSSLI63MUnlvICMO6RsyMw2qQxBIwOeKAPaaKyNH8VaBr+f7J1\nmxvGABKQzKXGfVeo/EVmt4/0KLx23hCW4EeoCJHVmI2M7ZPl57Nt2nHfNAHSXP8Ax6zf7h/lXgH7\nM/8Ax9eJf9y2/nJXv9z/AMes3+4f5V4B+zP/AMfXiX/ctv5yUAel+LQYb25e613UYriWMDSdP0uR\nllLgYLlFHz/MwzuygAGcc102mzasfs0OoWkQ/wBCjea5SUc3HR0CY4A6hs98Y4rMvvB32vXrnWLf\nX9XsLi5jSJ1tTBt2pnAG+JjjJJxnqTWxa6dJb3Uc76je3G21S3KSsuxipJMpCqPnOeSMDjgCgDA8\nUXOs2nibwz5F/FHpdzqAglt0iIkdvJmbl92NvyD5dvXv2q5MLjUPE12bSWJH0608mJ5VLos8uGO5\nQQThVjPUcOa0tS0i31SfTZp3lVtPuhdxBCAGfY6YbIPGJD0xzipbLT4rA3LRs7NcztPIz4yWOB2A\n4ACgewFAHA22s60nge4+0aiZNRl8QNpj3aLtEateeSWQEnbhc7Rk4461u6B52leLdS8Pi8u7uzjs\nre8ha7naaSNneVGUuxJI/dggE8ZPbpfXwppw0a90t/Okt7u7lvGLOAySvKZcqQBja5yvcYHWptF8\nPw6M9xObu7vry52ia7vGVpHVc7V+VVAAycAAdSepoAofEL/knHiT/sGz/wDoBrgP2b/+RG1P/sJN\n/wCio67/AOIX/JOPEn/YNn/9ANcB+zf/AMiNqf8A2Em/9FR0AdX4y+E/hjxlvnntfsWoNz9stQFd\nj/tjo/Tvz7itXwH4Xfwb4Os9CkuVuXt2lJmVSofdIzDjtwRXSUUAFFFFABRRRQAUUUUAFFFFAFXU\n/wDkFXn/AFwf/wBBNeE/sy/8zT/26f8AtavdtT/5BV5/1wf/ANBNeE/sy/8AM0/9un/tagD0bxzq\n+taDbX2qpqhtLeCL/QbaGz88XEgQsfPbafLTOBkFQAMlucDrbLUGu7iWE2syCOONxcEAwy7wTiNg\ncnGOeB1HWsrVfD+p3d5dSWGvPZW97GI7iFrcTYwpXdESQEJBGchhxnHXOhpWmSaUv2WK4U6bDBDB\naWwiwYAilTl8/NkbewxjvmgDSrl/EWsazp3iPw/a28FsumXt6IJpzITKx8qV9gTbgD5Ad27PbHeu\norK1jRv7WutIm+0eV/Z18LzGzd5mIpI9vUY/1mc89OnNAHPeJvElzpUmqX0EwWO1Nvp0O4M6LcTO\nu52QH5tqvGcdfvDvV3whq01/NfW02qTXpgCMFvbJrS6j3bvvRlFBQ4G1gOzA5IqzP4WiutBn06e7\nkaaW7N6LoKAyTCXzUIHQhSFAB6hcGpNG0K5stRudU1PUEvtSuIUt2lig8iNY0LMqhNzd3YkknPtQ\nBR+Jf/JM/Ef/AF4yfyrwn4U+MvE2m+Gbjw74T8PNqGoT3j3DXUh/cwKyRqM9Bn5SeSO3WvdviX/y\nTPxH/wBeMn8q4n9nH/knuof9hWT/ANFRUAEPwk8Q+KGFz4/8W3d0GIb7BYtshXByOox+Sg+9dno/\nwx8GaEqfY/D9m0inIluE858+uXzj8K62igBkcUcKbIo1Rf7qjAp9FFABVDUNE0nVo2j1LTLO8Ruq\n3ECuP1FX6KAPLPGHwk8Bw6Te6y1tNpJs4mnaaxlK42jPCnK54xwBXy1NeXNxeG7muZpLksG86SQs\n+R0O4854Ffek9vDdQtDcRRyxN95JFDKfqDWK/gfwlI7O/hbRGdjlmbT4iSfU/LQB5l8Kvi8mv2P/\nAAj/AIhnVNVWMrb3LnAuQB0P+3/P69cb9mf/AI+vEv8AuW385K9fufAnhD7PIy+FtFVlUsrJYxKQ\nceoWvIP2Z/8Aj68S/wC5bfzkoA9T8ZajqmlCW/8A7VfTtOiiAh+zWf2l5psOSJRsOyMALzle+WHG\nd3S9SuboWsN3YyRzvZR3Es0ZDW4duGjV85JBGemMEHNUdV8P6nd3t1JYa89nb3sYjuIXtxNjCld0\nRJGwkEZyGHGcdc6Gn6ZJpj29tbXCrpVtZx20Np5XzKU43b88jbtGMds55oAyfEWsazp3iPw/a28F\nsumXt6IJpzITKx8qV9gTbgD5Ad27PbHerdzc3s3iKeOxHmCxsi3ktKY45ppD8gZgDjAQ9jjzAasa\nxo39rXWkTfaPK/s6+F5jZu8zEUke3qMf6zOeenTmp9O04WMt9M0nmS3lyZ3bbjHyhVXqeiqo9+Tx\nmgDj7bxVrTeDZ7u5FumrSa02lxhPmjhJuvIBHA3BRzzjOO2a2NCvNRtfEd/4f1K/bUWgtYbyG7ki\nSN2WRpFKsEAXgx8EAcH2yRfB8R0G+0yS8cm41GXUYp0Ta0MjTmZMDJztbH1x0GataJoVxYXt3qWp\nX63+p3SRxPNHB5KLGhYqqpubHLsScnJP0FAFb4hf8k48Sf8AYNn/APQDXAfs3/8AIjan/wBhJv8A\n0VHXf/EL/knHiT/sGz/+gGuA/Zv/AORG1P8A7CTf+io6APZKKKKACiiigAooooAKKKKACiiigCrq\nf/IKvP8Arg//AKCa8J/Zl/5mn/t0/wDa1e7an/yCrz/rg/8A6Ca8J/Zl/wCZp/7dP/a1AHrHjDxW\n/h2O2is7cXF3NNCHDA7YonmSIuxHfLgAZ5OTyFNb0OoW0+oXVjG7G4tVRpVMbAAPnbhiMH7p6E47\n4yK4/wAU/D641l765sNf1G3uLy5tpZIHePyQInQ8ZiZhgKWUZxuPPBNdZYR6hBNJBcyRzWkcUSwT\ns2Z5WAIcyAKFBzgjb6ngcUAXq5zXPEd3pfiHRdMi0uV4NQuhC967L5a/JI+0ANuLfJ3GMd88V0dY\nuvaRcapeaFNA8Srp+oi7lDkgsnkyphcA85kHXHGaAM3X/FEuj3GoTrsaC0jgt0SRwiSXUzgKGbBK\nhQUJPo5ODirPhPXp9cjui97pOoRRMAt3pcmYyTnKFSzFWAxznBznjpVe+8Jzaj4du7S4mg+3Tah/\naAk2Exl0lDxqw6kBERD6gGrOh6PqMet3uuauLOK8ureK28iydnjVY2chizKpZjvPYYAA5oAr/Ev/\nAJJn4j/68ZP5V4b8K/ilY+AfB81rqOk6hPFcahJIlxAq7M+XGCvJHzDAOPRhXuXxL/5Jn4j/AOvG\nT+VeY/Bnwxpvi74P6vpOpxboZdUk2uAN8T+TFh1J6EUAa8f7R3hJiBJputJk9RDEQP8AyJXofhPx\ndpvjPS31LSlufsqSmLfPFs3MACceo5618leJfh5rfhrxhD4dmi86a6kVLOZBhLgMcAj0Oeo7fTBP\nseg+JvGPww0i20TW/BDXOl2uVW80sluOpZgMgkknJO3rQB7hRXnWk/HDwLqnlq+pS2Mr/wAF3Ay4\n+rDKj866mz8ZeGNQz9k8RaVMQMkJeRkge4zkUAblFVf7TsP+f62/7+r/AI1n3njDwzp+PtniHS4C\negku4wT9BnJoA2qK871f43+BdKDqmpyX8qj/AFdnCz5+jHC/rWG/j74heLv3XhDwi+nWzHA1DVOM\nDHUKcDr6b/pQB6Z4g1vTNB0ie81W9htIAhG6VsbjjoB1J9hXh/7M/wDx9eJf9y2/nJXVWHwaF3cv\nrPjnWbjX9RClliLFYE7gY6kZHQbR7Vyv7M//AB9eJf8Actv5yUAetePfFc/hXw9c3On28d1qSwvN\nHDJnYqJyzvjnaMgdskgcZyOgGoWx1RtNDt9qWETlfLbGwsVB3Y25yDxnPtzXJeMvAM/iGDWJ7DW7\n62vL+0Ft5DPH9nYKDhWzGzhckk7T3rp7KDUbW4S3kmS5sY7ZQLiZs3EkuTksFUJjGOR3zwKAMzXP\nEd3pfiHRdMi0uV4NQuhC967L5a/JI+0ANuLfJ3GMd88VZvNSuk117e2iknitLJrieCILvldmxGql\niAD8knUjtnijXtIuNUvNCmgeJV0/URdyhyQWTyZUwuAecyDrjjNWdM0+W1utRurh1aa7ufMG0khY\n1UIi8gdlyR6setAHOW3jW7m8JT6rJpohvTqbabBaO/3ZDceSgcgkcE5bGehxmtPQ9W1J9YvdE1oW\njX9tBFcrNaIyRyxyF1HyszEEGNgeTng+woL4Qum8OX1g9zCl0+sS6payLllVvtJnjDcA+gbHvg96\nv6HpOpJrF7retG0W/uYIrZYbR2eOKOMuw+ZlUkkyMTwMcD3IBH8Qv+SceJP+wbP/AOgGuA/Zv/5E\nbU/+wk3/AKKjrv8A4hf8k48Sf9g2f/0A1wH7N/8AyI2p/wDYSb/0VHQB7JRRRQAUUUUAFFFFABRR\nRQAUUUUAVdT/AOQVef8AXB//AEE14T+zL/zNP/bp/wC1q921P/kFXn/XB/8A0E14T+zL/wAzT/26\nf+1qAPWNU8SatD4jl0fR9Gtb14LaKeWS41D7OAZGkVVA8tsn92T26iuggvIJ5pLYTQm7hVGngSQM\n0W7ONw6gHBwT1wa43xT4ebUNVvJT4OsNVluLdY7a/aVFeBgG++X+ZQCQQ0eTz04ro9EgurTNpd2x\neSG2t0bUmZCb1wpDEgfMCCO/97jvQBr1g6t4ptdL17S9Ha3uZLnUJhGriFxEg2u2TJjaT8h+UHPf\npW9WB4k027v7/wAOSWsXmJZ6qLic7gNkfkTJnk8/M6jA55oAbq3iYaVeXgaEywWkEWVQfPLPK+yK\nJSSBk45z/fU5FW9F1S+v2ni1DTPsU0WCGjm86KRTn7r7R8wxyMcZHXrWBq/he91nwzfpPDH9uuNT\nW+8hpPlkSKVdkZYdN0Uaj0BY1N4Q0OXT9V1C/j0WPQbG5hijTTIzHjzEL7pSIyUBIZV4OSFGaAJf\niX/yTPxH/wBeMn8q4n9nH/knmof9hWT/ANFRV23xL/5Jn4j/AOvGT+VcT+zj/wAk81D/ALCsn/oq\nKgD1e5sLS8ntp7i2illtXMkDuoJjYgqSD24JFWaKKAMjVPCvh/W2Danomn3bjOHmt1Zh+OM1xGv/\nAAt+GOnRfa9Q0r7KJW2IsE05Z2xnCRqSWOAThR2Jr06uV8T3MWl+JNA1i+l8rTLdbmKaVlykTuq7\nHY/wjCuuTxlx60Acrp/wY+G+q2SXljbTT275AdbuTqDggjOQQQQQeQRirVl8JPhrDqzWKaWs99DE\nszxSXUrbUYkAkbsckHr6Vt+FNRtopbxpJ/LTWdVuJdNjdChlRUXcVBHQlHcE9Q2aj0nTbHTfirrC\n2NpDbCfSbaeURIF8yRp7jc5x1JwOfagDb0vwn4e0QltM0TT7Rz1eK3VWP1bGa2KKKAIrn/j1m/3D\n/KvAP2Z/+PrxL/uW385K9/uf+PWb/cP8q8A/Zn/4+vEv+5bfzkoA9c1TxJq0PiOXR9H0a1vXgtop\n5ZLjUPs4BkaRVUDy2yf3ZPbqK6CO8ga5Fo00K3oiEz24kBdVJI3Y67cgjOMcGuN8U+Hm1DVbyU+D\nrDVZbi3WO2v2lRXgYBvvl/mUAkENHk89OK6PS4Lq0mtrS6tjcPDYRJJqrMm6aQcMpH3u27OMfN65\noAg1bxTa6Xr2l6O1vcyXOoTCNXELiJBtdsmTG0n5D8oOe/SrF7rP2TV/s2zdBBZvdXRSNpHUZAQK\nq5JJxJ0BPy8VW8Sabd39/wCHJLWLzEs9VFxOdwGyPyJkzyefmdRgc81a0qzuI9Q1S+uk2yXM4WIZ\nBxCihV6ep3t/wPtQBlW3ji0uvC8+uJZ3Cot41lDbuu2SWTzfKQEEDaWYjg9O/Sr2ia7cX97d6bqV\ngthqdqkcrwxz+cjRuWCsr7VzyjAjAwR9DWCvhjVH8L31uIkjvU1+XVLZHcbZFW8MyAkZxuUY56Z5\nHFamhWeo3XiO/wDEGpWDac09rDZw2kkqSOqxtIxZihK8mTgAnge+AAL8Qv8AknHiT/sGz/8AoBrg\nP2b/APkRtT/7CTf+io67/wCIX/JOPEn/AGDZ/wD0A1wH7N//ACI2p/8AYSb/ANFR0AeyUUUUAFFF\nFABRRRQAUUUUAFFFFAFXU/8AkFXn/XB//QTXhP7Mv/M0/wDbp/7Wr3bU/wDkFXn/AFwf/wBBNeE/\nsy/8zT/26f8AtagD2PVvFdppN7JafZL69lgiE9yLOISfZ4juw78g87WwFyxx0railSaJJY2DI6hl\nYdweleb+INImt/FGt3s8HiK4a9jifTTpcsixpKse0q3lkAHKg5l+TDcfxV12gRCO6nN3Zuusm1tR\nf3ghKxXDhW+43Q7Tu4HTcM9aAN2sjUPE+k6brNho893GdRvnCQ2yMDJjDHeVzkL8p59eK165rxVZ\nzXGqeFpYLZ5fJ1cPK6IW8tPs84yxHQZIGTxkigC5qPiO10u7uo7gEQ2lqtxNIuWYF2KxoqgZZmKt\ngDnIAwc1LpGuR6s88LWd5Y3UG1nt7tAr7WztYbSQQdp6HjHOK5XXNG1DWPDWo3S290l1Lq0V0IVO\n2UwQTIAE9CUjLqD3arvg6GdNU1CS0TVodCeGLyItVaYy+eC/mFRMS4Xb5fXjIJFAE/xL/wCSZ+I/\n+vGT+VcT+zj/AMk81D/sKyf+ioq7b4l/8kz8R/8AXjJ/KuJ/Zx/5J5qH/YVk/wDRUVAHsFFFFABR\nRRQAUUUUAFFFFAEVz/x6zf7h/lXgH7M//H14l/3Lb+cle/3P/HrN/uH+VeAfsz/8fXiX/ctv5yUA\nez6t4rtNJvZLT7JfXssEQnuRZxCT7PEd2HfkHna2AuWOOlbUUqTRJLGwZHUMrDuD0rzfxBpE1v4o\n1u9ng8RXDXscT6adLlkWNJVj2lW8sgA5UHMvyYbj+Kuu0qIJqVu1/Zu2uDS4Vur5ISIXwTlA3TO/\nc23qAR60AbtZGoeJ9J03WbDR57uM6jfOEhtkYGTGGO8rnIX5Tz68Vr1zXiqzmuNU8LSwWzy+Tq4e\nV0Qt5afZ5xliOgyQMnjJFAFvWfEtrol1FFco7I0RkZkGSCXREQKOrOz4A9jT9I8QQardXFm1pd2N\n7bqsj212iq+xiQrjazAqSrdD25xWHqt1d6fDrWtjSJr25a5htbOJbV5WWNCAJCigsQHeV/lHIAx1\nzUng4w3N9e38/wDa0+rTRRpcXN9p0tomxSxVIldQAoLMcZJ55J4oAtfEL/knHiT/ALBs/wD6Aa4D\n9m//AJEbU/8AsJN/6Kjrv/iF/wAk48Sf9g2f/wBANcB+zf8A8iNqf/YSb/0VHQB7JRRRQAUUUUAF\nFFFABRRRQAUUUUAVtRVn0y7RFLM0LgADJJwa+VvBF18SfAP27+yPB+oSfbfL837TpU7Y2bsYxj++\nf0r6yooA+e/+FnfGL/oSpP8AwT3P/wAVR/ws74xf9CVJ/wCCe5/+Kr6ErzDXL7xbq3xZm8MaL4n/\nALGtItJW9z9giuMt5gQj5sHncO/bpzQBxX/CzvjF/wBCVJ/4J7n/AOKo/wCFnfGL/oSpP/BPc/8A\nxVeswag/gbw7Lc+NfFkV9mclLuS0W342jEaomdx+Vjxyc+1RSeOdC1/whq2o6B4khhFrAzPdiBna\n164domG49DgEc4oA8r/4Wd8Yv+hKk/8ABPc//FUf8LO+MX/QlSf+Ce5/+Kr1OP4g+HNH0vT49Y8R\nRyXLaVFftcNbPGLiJsL5oULgFm/gHzc9K7BWDoGU5VhkGgD5t1vxx8Wte0S80m78GXC293E0Uhj0\ni4DBT1wSTzXf/ALSdS0bwLfW+qafd2M7anI6x3ULRMV8qIZAYA4yCM+xr1SigAooooAKKKKACiii\ngAooooAjuAWtpQASShAA+lfKHgmf4keAZL19I8H6jIbwIJPtOlztjbuxjGP7xr6zriPiDr+p6Hf+\nEo9OufIS/wBbgtLkeWreZEx+ZfmBxn1GDQB5l/ws74xf9CVJ/wCCe5/+Ko/4Wd8Yv+hKk/8ABPc/\n/FV9AyyxwRPLK6xxopZnY4CgdSTXP6L498LeItRk0/Sdatrq7jzmJcgtjqVyBuHuM0AePf8ACzvj\nF/0JUn/gnuf/AIqj/hZ3xi/6EqT/AME9z/8AFV6vefE/wVp+tNpF14htY7xWKOpDFEYZBDOBsUgg\n5BIrQ1rxn4d8O3trZ6tq1va3F0QIo3JJOe5wDtHucCgDxj/hZ3xi/wChKk/8E9z/APFUf8LO+MX/\nAEJUn/gnuf8A4qu98AePBJ8LrXxF4v1eFHeaVHuJVVN2GIACoBk4HQDPFdZpvjDw9q+iTazY6vbS\n6fACZpy20RcZ+YNgrx2IoA8J1jx38W9b0a90u68F3C295C8Ehj0i5DBWGDgknnmu5+AOkano3g3U\nYNU067sZn1BnWO6gaJmXy0GQGA4yDz7VZvPiVbav448H2XhjW4rjT7ye6iv40jGW2IpTO5dy87sE\nYz716fQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAV5PrfgWz8WfGu4bXtGuLrR10VfLmPmxxecJR\ngb1Iy20t8ufwr1iigDzPxbosnhbRvD+neGdEm/sm2vGeaW0the3Vnk5DwrKG5JZstyQOg9OWGha4\n0/xKlk0vWn/tHS4haPeQq0tyQhU/6obS2R90cgEZAr3WigDyDwvoWoR/EPwhd3WlXSQWvg+3geaW\n3YLDOCcoSRhXAJ46816/RRQAUUUUAFFFFABRRRQAUUUUAFFFFABXm3xchv8AHhS+stLv9RGn61Dd\nzRWUDSvsTk8D6Y54zXpNFAHmOra3qPxG8L634dsvDfiHRrqayLx3GpWnkRSEOv7vcSeWBI+mfSqc\nUOq+Kda8I28PhjUtDXQJ1lurq6jRE2qmDFEQTvVuhI7V61RQB8+eIrDxTrHgvxFp8uha1BftdeYd\nPsLCGOxK+YreYrhQ8zEDnBJ3HkYrrxb6n4b+IU+sSeHtR1e01fTbe2ia2iVjbOq4ZHDEbEOMknuf\nwr1SigD580rwl4gX4d+Erv8As7V4JNJ1C5mubK2jVLlUcna6JKCCy44BHRj9a0L3wNqPiHwX4rvN\nPstfi1HVDbt5etyQpLciFgf9XGqhDjgbjyfTrXudFAHj91caj4k+IPgW/h8H6vplrYSXEc8t1a7d\nnyAAfKThM9C2M9u9ewUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFAB\nRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFAH/9k=\n",
      "text/plain": [
       "<IPython.core.display.Image at 0x109f12410>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image('casino.jpg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def roll(die='F'):\n",
    "  if die=='F': return np.random.randint(1,7)\n",
    "  elif die== 'L': return np.random.choice([1,2,3,4,5,6,6,6,6,6])\n",
    "  else: return None\n",
    "    \n",
    "def trial(state):\n",
    "  r=np.random.rand()\n",
    "  if state=='F' and r < .05: state = 'L'\n",
    "  elif state=='L' and r < .1: state = 'F'\n",
    "  return state,str(roll(state))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And generate a series of 300 hidden state transitions and observations starting from the 'F' state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "q='F'\n",
    "states=q\n",
    "rolls=str(roll(q))\n",
    "for i in xrange(299): #range in python3\n",
    "    q,o=trial(q)  #update state q, and get observed o\n",
    "    states += q\n",
    "    rolls += o"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4242515462351161555223451263544366336116515561154145424646262651632521164555362346662536566466643511\n",
      "FFFFFFFFFFFFFLLLFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLFFFLLLLLLLFLLLLLLLLLLFFFFFF \n",
      "\n",
      "4556213455631362425354166562666661626416623151424364234166226656653231266666621665546336665546622165\n",
      "FFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLFFFFFFFFFLLFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLFLLLLLLLLLLLLLLLLL \n",
      "\n",
      "2641313326461151232644165264466616655623346562552425232663226433134456133666541325651625451255323413\n",
      "LLFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLFFFFFLLLFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF \n",
      "\n"
     ]
    }
   ],
   "source": [
    "for i in 0,1,2:\n",
    "  print rolls[100*i:100*(i+1)]   #add parentheses to print statements for python3\n",
    "  print states[100*i:100*(i+1)],'\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above the hidden states are known, but suppose we try to infer the state transitions knowing only the obervations in the middle line above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "observed: 4556213455631362425354166562666661626416623151424364234166226656653231266666621665546336665546622165\n",
      "  hidden: FFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLFFFFFFFFFLLFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLFLLLLLLLLLLLLLLLLL \n",
      "\n",
      "inferred: FFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL\n"
     ]
    }
   ],
   "source": [
    "O='4556213455631362425354166562666661626416623151424364234166226656653231266666621665546336665546622165'\n",
    "Q='FFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLFFFFFFFFFLLFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLFLLLLLLLLLLLLLLLLL'\n",
    "\n",
    "T=np.array([[.95,.05],[.1,.9]])\n",
    "snames=['F','L']\n",
    "\n",
    "e = [{str(o): 1/6. for o in range(1,7)},\n",
    "     {str(o): 1/10. if o<6 else 1/2. for o in range(1,7) }]\n",
    "\n",
    "w=[.5,.5]\n",
    "\n",
    "print 'observed:',O   #add parentheses to print statements for python3\n",
    "print '  hidden:',Q,'\\n'\n",
    "print 'inferred:',viterbi(T,e,w,snames,O)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that it gets the essential large-scale features of the state transitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}