200010036_Lab_10_Hidden_Markov_Model.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zZDMbDP0uJN0"
      },
      "source": [
        "# **LAB 11 : Hidden Markov Model**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "MKFudLjzt_K5"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "64j_rFZeF57i"
      },
      "source": [
        "Please refer to the following [article](http://www.adeveloperdiary.com/data-science/machine-learning/introduction-to-hidden-markov-model/) to understand Hidden Markov Model\n",
        "\n",
        "Here we will be dealing with 3 major problems : \n",
        "  \n",
        "  1. Evaluation Problem\n",
        "  2. Learning Problem\n",
        "  3. Decoding Problem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wEasU5TNGNbE"
      },
      "source": [
        "1. Evaluation Problem : Implementation of Forward and Backward Algorithm"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vjUe9joAGXvC",
        "outputId": "0a1ef95b-9653-4518-adbc-d288b6343da5"
      },
      "outputs": [],
      "source": [
        "data = pd.read_csv('data_python.csv') ## Read the data, change the path accordingly\n",
        "\n",
        "V = data['Visible'].values\n",
        "\n",
        "# Transition Probabilities\n",
        "a = np.array(((0.54, 0.46), (0.49, 0.51)))\n",
        "\n",
        "# Emission Probabilities\n",
        "b = np.array(((0.16, 0.26, 0.58), (0.25, 0.28, 0.47)))\n",
        "\n",
        "# Equal Probabilities for the initial distribution\n",
        "initial_distribution = np.array((0.5, 0.5))\n",
        "\n",
        "\n",
        "def forward(V, a, b, initial_distribution):\n",
        "\t# Probability that hidden state at time t is i given the observations up to time t\n",
        "\talpha = np.zeros((V.shape[0], a.shape[0]))\n",
        "\t# Define alpha for t=0\n",
        "\talpha[0, :] = initial_distribution * b[:, V[0]]\n",
        "\n",
        "\t# Define other alpha values\n",
        "\tfor t in range(1, V.shape[0]):\n",
        "\t\tfor j in range(a.shape[0]):\n",
        "\t\t\talpha[t, j] = np.dot(alpha[t-1, :], a[:, j]) * b[j, V[t]]\n",
        "\n",
        "\treturn alpha\n",
        "\n",
        "\n",
        "alpha = forward(V, a, b, initial_distribution)\n",
        "\n",
        "\n",
        "def backward(V, a, b):\n",
        "\t# Probability that the hidden state at time t is i given the observations after time t\n",
        "\tbeta = np.zeros((V.shape[0], a.shape[0]))\n",
        "\n",
        "\t# Set beta to one for last time step\n",
        "\tbeta[V.shape[0]-1, :] = np.ones(a.shape[0])\n",
        "\n",
        "\t# Define other beta values\n",
        "\tfor t in range(V.shape[0]-2, -1, -1):\n",
        "\t\tfor j in range(a.shape[0]):\n",
        "\t\t\tbeta[t, j] = np.dot(beta[t+1, :] * b[:, V[t+1]], a[j, :])\n",
        "\n",
        "\treturn beta\n",
        "\n",
        "\n",
        "beta = backward(V, a, b)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5Iyc4mIlGzPk"
      },
      "source": [
        "2. Learning Problem : Implementation of Baum Welch Algorithm"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "H6eaCtNEG3FQ",
        "outputId": "143855ca-81f4-4bfb-97f9-5a40596c15eb"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[0.53810394 0.48523844]\n",
            " [0.46245188 0.51417766]]\n",
            "[[0.16352973 0.26235368 0.57411658]\n",
            " [0.25062945 0.27803505 0.4713355 ]]\n"
          ]
        }
      ],
      "source": [
        "def baum_welch(V, a, b, initial_distribution, n_iter=100):\n",
        "\t# Number of hidden states\n",
        "\tM = a.shape[0]\n",
        "\t# Number of emission states\n",
        "\tK = b.shape[1]\n",
        "\t# Number of observations\n",
        "\tT = len(V)\n",
        "\n",
        "\tfor _ in range(n_iter):\n",
        "\t\t# Calculate alpha and beta\n",
        "\t\talpha = forward(V, a, b, initial_distribution)\n",
        "\t\tbeta = backward(V, a, b)\n",
        "\n",
        "\t\t# Calculate xi\n",
        "\t\t# ie. probability of being in state i at time t and state j at time t+1 given the observations and the model\n",
        "\t\txi = np.zeros((M, M, T-1))\n",
        "\n",
        "\t\tfor t in range(T-1):\n",
        "\t\t\tfor i in range(M):\n",
        "\t\t\t\tfor j in range(M):\n",
        "\t\t\t\t\txi[i, j, t] = alpha[t, i] * a[i, j] * b[j, V[t+1]] * beta[t+1, j]\n",
        "\t\t\t\t\n",
        "\t\t# Normalize xi\n",
        "\t\tfor t in range(T-1):\n",
        "\t\t\txi[:, :, t] /= np.sum(xi[:, :, t])\n",
        "\n",
        "\t\t# Calculate gamma\n",
        "\t\t# ie. probability of being in state i at time t given the observations and the model\n",
        "\t\tgamma = np.sum(xi, axis=1)\n",
        "\n",
        "\t\t# Recalculate transition probabilities\n",
        "\t\ta = np.sum(xi, axis=2) / np.sum(gamma, axis=1)\n",
        "\n",
        "\t\t# Add the last gamma ie. at time t=T\n",
        "\t\tgamma = np.hstack((gamma, np.sum(xi[:, :, T-2], axis=0).reshape((-1, 1))))\n",
        "\n",
        "\t\t# Recalculate emission probabilities\n",
        "\t\tdenom = np.sum(gamma, axis=1)\n",
        "\n",
        "\t\tfor k in range(K):\n",
        "\t\t\tb[:, k] = np.sum(gamma[:, V == k], axis=1)\n",
        "\n",
        "\t\tb /= denom.reshape((-1, 1))\n",
        "\n",
        "\treturn (a,b)\n",
        "\n",
        "data = pd.read_csv('data_python.csv')\n",
        "\n",
        "V = data['Visible'].values\n",
        "\n",
        "# Transition Probabilities\n",
        "a = np.ones((2, 2))\n",
        "a = a / np.sum(a, axis=1)\n",
        "\n",
        "# Emission Probabilities\n",
        "b = np.array(((1, 3, 5), (2, 4, 6)))\n",
        "b = b / np.sum(b, axis=1).reshape((-1, 1))\n",
        "\n",
        "# Equal Probabilities for the initial distributionhidden markov model\n",
        "initial_distribution = np.array((0.5, 0.5))\n",
        "\n",
        "a,b = baum_welch(V, a, b, initial_distribution, n_iter=100)\n",
        "\n",
        "print(a)\n",
        "print(b)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A40-TPXZHTHE"
      },
      "source": [
        "3. Decoding Problem : Implementation of Viterbi Algorithm"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nBG74dknHYgM",
        "outputId": "1ac04c53-fea0-4047-aae0-768c6c6e0780"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "['A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A']\n"
          ]
        }
      ],
      "source": [
        "def viterbi(V, a, b, initial_distribution):\n",
        "\t# Number of hidden states\n",
        "\tM = a.shape[0]\n",
        "\t# Number of observations\n",
        "\tT = V.shape[0]\n",
        "\n",
        "\t# Calculate omega\n",
        "\t# ie. at time t probability of being in state i given the previous observations and the model\n",
        "\t# Note we use the log scale to avoid underflow\n",
        "\tomega = np.zeros((T, M))\n",
        "\tomega[0, :] = np.log(initial_distribution * b[:, V[0]])\n",
        "\n",
        "\t# prev matrix store the most likely previous state\n",
        "\tprev = np.zeros((T-1, M))\n",
        "\n",
        "\tfor t in range(1, T):\n",
        "\t\tfor j in range(M):\n",
        "\t\t\t# Calculate probability of being in state j coming from i at time t given the previous observations and the model\n",
        "\t\t\tprobabilities = omega[t-1, :] + np.log(a[:, j]) + np.log(b[j, V[t]])\n",
        "\n",
        "\t\t\t# Store the most likely previous state\n",
        "\t\t\tprev[t-1, j] = np.argmax(probabilities)\n",
        "\n",
        "\t\t\t# Store the probability of being in state j at time t given the previous observations and the model\n",
        "\t\t\tomega[t, j] = np.max(probabilities)\n",
        "\n",
        "\t# Calculate the most likely path\n",
        "\tpath = np.zeros(T, dtype=int)\n",
        "\n",
        "\t# The last state is the one with the highest probability\n",
        "\tlast = np.argmax(omega[T-1, :])\n",
        "\n",
        "\t# Calculate the rest of the path\n",
        "\tfor t in range(T-2, -1, -1):\n",
        "\t\tpath[t] = prev[t, path[t+1]]\n",
        "\n",
        "\t# Convert state indices to names\n",
        "\tstates = {0: 'A', 1: 'B'}\n",
        "\tresult = [states[int(i)] for i in path]\n",
        "\n",
        "\treturn result\n",
        "\n",
        "\n",
        "data = pd.read_csv('data_python.csv')\n",
        "\n",
        "V = data['Visible'].values\n",
        "\n",
        "# Transition Probabilities\n",
        "a = np.ones((2, 2))\n",
        "a = a / np.sum(a, axis=1)\n",
        "\n",
        "# Emission Probabilities\n",
        "b = np.array(((1, 3, 5), (2, 4, 6)))\n",
        "b = b / np.sum(b, axis=1).reshape((-1, 1))\n",
        "\n",
        "# Equal Probabilities for the initial distribution\n",
        "initial_distribution = np.array((0.6, 0.4))\n",
        "\n",
        "a, b = baum_welch(V, a, b, initial_distribution, n_iter=100)\n",
        "\n",
        "result1 = viterbi(V, a, b, initial_distribution)\n",
        "\n",
        "print(result1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6wTG6SkvKOHl"
      },
      "source": [
        "4. Use the built-in **hmmlearn** package to fit the data and generate the result using the decoder"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "id": "R1iy8tj4KrwF"
      },
      "outputs": [],
      "source": [
        "#%pip install hmmlearn\n",
        "from hmmlearn import hmm"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "id": "Pw-zo9opLLlh"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "['B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'A', 'B', 'A', 'B', 'B', 'A', 'B', 'A', 'B', 'A', 'B', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'B', 'A', 'B', 'B', 'B', 'A', 'B', 'B', 'B', 'A', 'B', 'B', 'B', 'A', 'B', 'A', 'B', 'B', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'B', 'A', 'B', 'B', 'B', 'A', 'A', 'B', 'A', 'A', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'B', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'B', 'A', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'B', 'A', 'B', 'A', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'B', 'A', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'B', 'A', 'B', 'A', 'B', 'A', 'B', 'B', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'B', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'A', 'B', 'A', 'B', 'B', 'B', 'A', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'B', 'A', 'B', 'B', 'A', 'B', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'B', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'B', 'B', 'A', 'A', 'B', 'B', 'A', 'B', 'B', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'B', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'A', 'B', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'B', 'B', 'A', 'B', 'A', 'B', 'A', 'B', 'B', 'A', 'A', 'B', 'A', 'A', 'B', 'A', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'B', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'B', 'A', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'A', 'B', 'A', 'A', 'A', 'A', 'B', 'A', 'B', 'A', 'B', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'B', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'A', 'B', 'A', 'A', 'A', 'B', 'A', 'B', 'A', 'A', 'B', 'B', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'A', 'B', 'B', 'B', 'A', 'B', 'A', 'A', 'B', 'B', 'B', 'A', 'B', 'B', 'A', 'A', 'B', 'A']\n"
          ]
        }
      ],
      "source": [
        "# Load the data\n",
        "data = pd.read_csv('data_python.csv')\n",
        "\n",
        "V = data['Visible'].values\n",
        "\n",
        "a = np.ones((2, 2))\n",
        "a = a / np.sum(a, axis=1)\n",
        "\n",
        "b = np.array(((1, 3, 5), (2, 4, 6)))\n",
        "b = b / np.sum(b, axis=1).reshape((-1, 1))\n",
        "\n",
        "initial_distribution = np.array((0.5, 0.5))\n",
        "\n",
        "# Use the hmmlearn library to train the model\n",
        "model = hmm.CategoricalHMM(n_components=2)\n",
        "model.startprob_ = initial_distribution\n",
        "model.transmat_ = a\n",
        "model.emissionprob_ = b\n",
        "\n",
        "# Predict the most likely hidden states\n",
        "log_prob, result2 = model.decode([V])\n",
        "states = {0: 'A', 1: 'B'}\n",
        "result2 = [states[i] for i in result2]\n",
        "\n",
        "print(result2)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "name": "Lab_11_Hidden_Markov_Model.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3.10.5 64-bit",
      "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.10.8"
    },
    "vscode": {
      "interpreter": {
        "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}