{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Dishonest Casino Example\n\nWe'll use the ubiquitous dishonest casino example to demonstrate\nhow to train a Hidden Markov Model (HMM) on somewhat realistic\ntest data (e.g. http://www.mcb111.org/w06/durbin_book.pdf\nChapter 3).\n\nIn this example, we suspect that a casino is trading out a fair\ndie (singular or dice) for a loaded die. We want to figure out\n1) when the loaded die was used (i.e. the most likely path) 2) how\noften the loaded die is used (i.e. the transition probabilities) and\n3) the probabilities for each outcome of a roll for the loaded die\n(i.e. the emission probabilities).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, import necessary modules and functions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom hmmlearn import hmm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, let's act as the casino and exchange a fair die for a loaded one\nand generate a series of rolls that someone at the casino would\nobserve.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# make our generative model with two components, a fair die and a\n# loaded die\ngen_model = hmm.CategoricalHMM(n_components=2, random_state=99)\n\n# the first state is the fair die so let's start there so no one\n# catches on right away\ngen_model.startprob_ = np.array([1.0, 0.0])\n\n# now let's say that we sneak the loaded die in:\n# here, we have a 95% chance to continue using the fair die and a 5%\n# chance to switch to the loaded die\n# when we enter the loaded die state, we have a 90% chance of staying\n# in that state and a 10% chance of leaving\ngen_model.transmat_ = np.array([[0.95, 0.05],\n                                [0.1, 0.9]])\n\n# now let's set the emission means:\n# the first state is a fair die with equal probabilities and the\n# second is loaded by being biased toward rolling a six\ngen_model.emissionprob_ = \\\n    np.array([[1 / 6, 1 / 6, 1 / 6, 1 / 6, 1 / 6, 1 / 6],\n              [1 / 10, 1 / 10, 1 / 10, 1 / 10, 1 / 10, 1 / 2]])\n\n# simulate the loaded dice rolls\nrolls, gen_states = gen_model.sample(30000)\n\n# plot states over time, let's just look at the first rolls for clarity\nfig, ax = plt.subplots()\nax.plot(gen_states[:500])\nax.set_title('States over time')\nax.set_xlabel('Time (# of rolls)')\nax.set_ylabel('State')\nfig.show()\n\n# plot rolls for the fair and loaded states\nfig, ax = plt.subplots()\nax.hist(rolls[gen_states == 0], label='fair', alpha=0.5,\n        bins=np.arange(7) - 0.5, density=True)\nax.hist(rolls[gen_states == 1], label='loaded', alpha=0.5,\n        bins=np.arange(7) - 0.5, density=True)\nax.set_title('Roll probabilities by state')\nax.set_xlabel('Count')\nax.set_ylabel('Roll')\nax.legend()\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, let's see if we can recover our hidden states, transmission matrix\nand emission probabilities.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# split our data into training and validation sets (50/50 split)\nX_train = rolls[:rolls.shape[0] // 2]\nX_validate = rolls[rolls.shape[0] // 2:]\n\n# check optimal score\ngen_score = gen_model.score(X_validate)\n\nbest_score = best_model = None\nn_fits = 50\nnp.random.seed(13)\nfor idx in range(n_fits):\n    model = hmm.CategoricalHMM(\n        n_components=2, random_state=idx,\n        init_params='se')  # don't init transition, set it below\n    # we need to initialize with random transition matrix probabilities\n    # because the default is an even likelihood transition\n    # we know transitions are rare (otherwise the casino would get caught!)\n    # so let's have an Dirichlet random prior with an alpha value of\n    # (0.1, 0.9) to enforce our assumption transitions happen roughly 10%\n    # of the time\n    model.transmat_ = np.array([np.random.dirichlet([0.9, 0.1]),\n                                np.random.dirichlet([0.1, 0.9])])\n    model.fit(X_train)\n    score = model.score(X_validate)\n    print(f'Model #{idx}\\tScore: {score}')\n    if best_score is None or score > best_score:\n        best_model = model\n        best_score = score\n\nprint(f'Generated score: {gen_score}\\nBest score:      {best_score}')\n\n# use the Viterbi algorithm to predict the most likely sequence of states\n# given the model\nstates = best_model.predict(rolls)\n\n# plot our recovered states compared to generated (aim 1)\nfig, ax = plt.subplots()\nax.plot(gen_states[:500], label='generated')\nax.plot(states[:500] + 1.5, label='recovered')\nax.set_yticks([])\nax.set_title('States compared to generated')\nax.set_xlabel('Time (# rolls)')\nax.set_xlabel('State')\nax.legend()\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's check our learned transition probabilities and see if they match.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f'Transmission Matrix Generated:\\n{gen_model.transmat_.round(3)}\\n\\n'\n      f'Transmission Matrix Recovered:\\n{best_model.transmat_.round(3)}\\n\\n')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, let's see if we can tell how the die is loaded.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f'Emission Matrix Generated:\\n{gen_model.emissionprob_.round(3)}\\n\\n'\n      f'Emission Matrix Recovered:\\n{best_model.emissionprob_.round(3)}\\n\\n')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this case, we were able to get very good estimates of the transition and\nemission matrices, but decoding the states was imperfect. That's because\nthe decoding algorithm is greedy and picks the most likely series of states\nwhich isn't necessarily what happens in real life. Even so, our model could\ntell us when to watch for the loaded die and we'd have a better chance at\ncatching them red-handed.\n\n"
      ]
    }
  ],
  "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.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}