HMMBasicVer1.2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Theory Behind HMM (Jurafsky Notes):\n",
    "\"\"\" A Markov chain is useful when we need to compute a probability for a sequence of\n",
    "events that we can observe in the world. In many cases, however, the events we are\n",
    "interested in may not be directly observable in the world. For example, in part-of-\n",
    "speech tagging, we don’t observe part-of-speech tags in the world; we see words and have to infer the correct \n",
    "tags from the word sequence. Hence, we call the part-of-speech tags hidden because they are not observed.\n",
    "\"\"\"\n",
    "\"\"\"Imagine that you are a climatologist in the year 2799 studying the history of global\n",
    "warming. You cannot find any records of the weather in Baltimore, Maryland, for\n",
    "the summer of 2007, but you do find Jason Eisner’s diary, which lists how many ice\n",
    "creams Jason ate every day that summer. Our goal is to use these observations to\n",
    "estimate the temperature every day. We’ll simplify this weather task by assuming\n",
    "there are only two kinds of days: cold (C) and hot (H). So the Eisner task is as\n",
    "follows:\n",
    "Given a sequence of observations O, each observation an integer cor-\n",
    "responding to the number of ice creams eaten on a given day, figure\n",
    "out the correct ‘hidden’ sequence Q of weather states (H or C) which\n",
    "caused Jason to eat the ice cream.\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Terminology\n",
    "\"\"\" \n",
    "Q = q1 q2 . . . qN             : Set of N states\n",
    "A = a[i][j] s.t. 1<=i,j<=N     : Transition probability matrix A,s.t a[i][j] is prob of moving from state i to j\n",
    "a[i][j]=P(qj|qi)\n",
    "pi = pi[1],pi[2],...,pi[N]     : Initial prob distr over states. pi[i] is prob of Markov chain starting in state i.\n",
    "O = o[1]o[2]...O[T]            : A Sequence of T Observations, drawn from a Vocab V=v[1]v[2]...,v[v]\n",
    "\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Assumptions:\n",
    "\"\"\"\n",
    "I have Assumped the Observations,States to be 0,1,2.... i.e Natural Numbers:\n",
    "So, please make sure to encode textual/categorical/arbitrary numeric data to Continuous Integers starting from 0 \n",
    "\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"A first-order hidden Markov model instantiates two simplifying assumptions.\n",
    "First, as with a first-order Markov chain, the probability of a particular state depends\n",
    "only on the previous state:\n",
    "    Markov Assumption: P(qi |q1 ...qi−1) = P(qi|qi−1)\n",
    "Second, the probability of an output observation o[i] depends only on the state that\n",
    "produced the observation qi and not on any other states or any other observations:\n",
    "    Output Independence: P(o[i]|q1 . . qi , . . , qT , o[1] , . . , o[i] , . . , o[T])= P(o[i]|qi)\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"There is a (non-zero) probability of transitioning between any two states. Such an HMM is called \n",
    "a fully connected or ergodic HMM. Sometimes, however, we have HMMs in which many of the transitions \n",
    "between states have zero probability.\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"hidden Markov models should be characterized by three fundamental problems:\n",
    "    Problem 1 (Likelihood): Given an HMM λ = (A, B) and an observation sequence O, determine the \n",
    "    likelihood P(O|λ ).\n",
    "    Problem 2 (Decoding): Given an observation sequence O and an HMM λ =(A, B), discover the best hidden state \n",
    "    sequence Q.\n",
    "    Problem 3 (Learning): Given an observation sequence O and the set of states in the HMM, learn the HMM \n",
    "    parameters A and B.\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from math import log2\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Likelihood Computation: The Forward Algorithm  O(N*N*T)\n",
    "#Our first problem is to compute the likelihood of a particular observation sequence say 3 1 3?\n",
    "\"\"\"For a Markov chain, where the surface observations are the same as the hidden\n",
    "events, we could compute the probability of 3 1 3 just by following the states labeled\n",
    "3 1 3 and multiplying the probabilities along the arcs. For a hidden Markov model,\n",
    "things are not so simple. We want to determine the probability of an ice-cream\n",
    "observation sequence like 3 1 3, but we don’t know what the hidden state sequence\n",
    "is!\"\"\"\n",
    "#O: Observations, a:State Transition Prob, b:Emission Prob, pi: Initial Prob. \n",
    "#a[i][j]: Prob of transition from state i to state j\n",
    "#b[i][j]: prob of emitting observation j at state i\n",
    "#pi[i]: Initial Prob.of state i\n",
    "#P(O|a,b,pi)=P(O|X=x1,a,b,pi)+P(O|X=x2,a,b,pi)+......+P(O|X=xn,a,b,pi) s.t. x1,x2,....,xn are all possible sequences\n",
    "#alpha[t][i]:Partial Observation sequence upto time t is generated and we are state i\n",
    "#alpha[t][i]= alpha[t-1][0]*a[0][i]*b[i][o[t]]+\n",
    "#    alpha[t-1][1]*a[1][i]*b[i][o[t]]+....+alpha[t-1][n-1]*a[n-1][i]*b[i][o[t]]#Assuming we have n states\n",
    "def computeAlpha(observations,a,b,pi,alphaDP):\n",
    "    statesC=a.shape[0]\n",
    "    timePts=observations.shape[0]\n",
    "    if timePts<1:\n",
    "        return\n",
    "    for i in np.arange(statesC):\n",
    "        alphaDP[0][i]=pow(2,log2(pi[i])+log2(b[i][observations[0]]))#pi[i]*b[i][observations[0]]\n",
    "    for t in np.arange(1,timePts):\n",
    "        for i in np.arange(statesC):\n",
    "            for j in np.arange(statesC):\n",
    "                alphaDP[t][i]+=pow(2,log2(alphaDP[t-1][j])+log2(a[j][i])+log2(b[i][observations[t]]))\n",
    "                #alphaDP[t][i]=alphaDP[t][i]-(log10(alphaDP[t-1][j])+log10(a[j][i])+log10(b[i][observations[t]]))\n",
    "                #alphaDP[t][i]+=alphaDP[t-1][j]*a[j][i]*b[i][observations[t]]\n",
    "def observationProb(observations,a,b,pi,alphaDP):\n",
    "    computeAlpha(observations,a,b,pi,alphaDP)\n",
    "    print(alphaDP)\n",
    "    timePts=observations.shape[0]\n",
    "    stateC=a.shape[0]\n",
    "    ans=0.0\n",
    "    for i in np.arange(stateC):\n",
    "        ans+=alphaDP[timePts-1][i]\n",
    "    return ans\n",
    "def observationsLikelihood(alphaDP):\n",
    "    timePts=alphaDP.shape[0]\n",
    "    stateC=alphaDP.shape[1]\n",
    "    ans=0.0\n",
    "    for i in np.arange(stateC):\n",
    "        ans+=alphaDP[timePts-1][i]\n",
    "    return ans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  6.00000000e-02   2.80000000e-01]\n",
      " [  6.16000000e-02   3.72000000e-02]\n",
      " [  5.80000000e-03   2.85600000e-02]\n",
      " [  7.74200000e-03   1.88760000e-03]\n",
      " [  6.17444000e-04   2.41861200e-03]\n",
      " [  5.59862240e-04   3.27280080e-04]\n",
      " [  2.61407800e-04   3.64326720e-05]]\n",
      "0.000297840472\n",
      "0.000297840472\n"
     ]
    }
   ],
   "source": [
    "observations=np.array([0,1,0,2,0,1,2])\n",
    "A=np.array([[0.7,0.3],[0.4,0.6]])\n",
    "B=np.array([[0.1,0.4,0.5],[0.7,0.2,0.1]])\n",
    "pi=np.array([0.6,0.4])\n",
    "alphaDP=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "#Can Easily Optimize Space Requirement of above alphaDP; Only need 2*count_of_Hidden_States Space\n",
    "print(observationProb(observations,A,B,pi,alphaDP))\n",
    "print(observationsLikelihood(alphaDP))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"Problem 2 (Decoding): The Viterbi Algorithm -\n",
    "Decoding: Given as input an HMM λ = (A, B, pi) and a sequence of observations O = o1 , o2 , ..., oT\n",
    "find the most probable sequence of states Q = q1 q2 q3 . . . qT \n",
    "\"\"\"\n",
    "\"\"\"Note that the Viterbi algorithm is identical to the forward algorithm except that it takes the max over the\n",
    "previous path probabilities whereas the forward algorithm takes the sum. Note also\n",
    "that the Viterbi algorithm has one component that the forward algorithm doesn’t have: backpointers. The reason \n",
    "is that while the forward algorithm needs to produce an observation likelihood, the Viterbi algorithm must \n",
    "produce a probability and also the most likely state sequence. We can compute this best state sequence by keeping\n",
    "track of the path of hidden states that led to each state, and then at the end backtracing \n",
    "the best path to the beginning (the Viterbi backtrace).\n",
    "\"\"\"\n",
    "def ViterbiAlgo(observations,a,b,pi,viterbiDP,viterbiBackPtr):\n",
    "    statesC=a.shape[0]\n",
    "    timePts=observations.shape[0]\n",
    "    if statesC<1:\n",
    "        return\n",
    "    for i in np.arange(statesC):\n",
    "        viterbiDP[0][i]=pow(2,log2(pi[i])+log2(b[i][observations[0]]))#pi[i]*b[i][observations[0]]\n",
    "    #Doesn't Handle the Case when there is single observation(1 Time Point)\n",
    "    for t in np.arange(1,timePts):\n",
    "        for i in np.arange(statesC):\n",
    "            maxSoFar=0\n",
    "            maxJ=-1\n",
    "            for j in np.arange(statesC):\n",
    "                tempVal=pow(2,log2(viterbiDP[t-1][j])+log2(a[j][i]))#viterbiDP[t-1]*a[j][i]\n",
    "                if tempVal>maxSoFar:\n",
    "                    maxSoFar=tempVal\n",
    "                    maxJ=j\n",
    "            viterbiDP[t][i]=tempVal*b[i][observations[t]]\n",
    "            #At State i,Produce observation b[i][observations[t]]\n",
    "            viterbiBackPtr[t][i]=maxJ\n",
    "def getOptimalHiddenStates(observations,a,b,pi,viterbiDP,viterbiBackPtr):\n",
    "    ViterbiAlgo(observations,a,b,pi,viterbiDP,viterbiBackPtr)\n",
    "    print(viterbiDP)\n",
    "    print(viterbiBackPtr)\n",
    "    timePts=observations.shape[0]\n",
    "    stateC=a.shape[0]\n",
    "    ans=0.0\n",
    "    endState=-1\n",
    "    for i in np.arange(stateC):\n",
    "        if ans<alphaDP[timePts-1][i]:\n",
    "            ans=alphaDP[timePts-1][i]\n",
    "            endState=i\n",
    "    stateSeq=np.zeros(shape=(timePts))\n",
    "    stateSeq[timePts-1]=endState#This is where we are adding the Final Optimal State for the Given Observation\n",
    "    i=timePts-2\n",
    "    prevState=endState\n",
    "    while i>=0:\n",
    "        prevState=viterbiBackPtr[i+1][prevState]\n",
    "        stateSeq[i]=prevState\n",
    "        i=i-1\n",
    "    return {\"Probability\":ans,\"Sequence\":stateSeq}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.06        0.28      ]\n",
      " [ 0.0448      0.0336    ]\n",
      " [ 0.001344    0.014112  ]\n",
      " [ 0.0028224   0.00084672]]\n",
      "[[ 0.  0.]\n",
      " [ 1.  1.]\n",
      " [ 0.  1.]\n",
      " [ 1.  1.]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alok/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:50: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'Probability': 0.0077420000000000067, 'Sequence': array([ 1.,  1.,  1.,  0.])}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observations=np.array([0,1,0,2])\n",
    "A=np.array([[0.7,0.3],[0.4,0.6]])\n",
    "B=np.array([[0.1,0.4,0.5],[0.7,0.2,0.1]])\n",
    "pi=np.array([0.6,0.4])\n",
    "viterbiDP=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "viterbiBackPtr=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "getOptimalHiddenStates(observations,A,B,pi,viterbiDP,viterbiBackPtr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Problem 3 (HMM Parameter Training): The Forward-Backward/Baum-Welch Algorithm(Special Case of EM Algo) -\n",
    "Given an observation sequence O and the set of possible states in the HMM, learn the HMM \n",
    "parameters λ = (A, B) \n",
    "\"\"\"\n",
    "\"\"\"The input to such a learning algorithm would be an unlabeled sequence of ob-\n",
    "servations O and a vocabulary of potential hidden states Q. Thus, for the ice cream\n",
    "task, we would start with a sequence of observations O = {1, 3, 2, ..., } and the set\n",
    "of hidden states H and C. For the part-of-speech tagging task, we would start with\n",
    "a sequence of observations O = {w1 , w2 , w3 . . .} and a set of hidden states NN, NNS,\n",
    "VBD, IN,... and so on.\n",
    "\"\"\"\n",
    "\"\"\"\n",
    "Let us begin by considering the much simpler case of training a Markov chain rather than a hidden Markov model. Since \n",
    "the states in a Markov chain are observed, we can run the model on the observation sequence and directly see which path \n",
    "we took through the model and which state generated each observation symbol. A Markov chain of course has no emission \n",
    "probabilities B (alternatively, we could view a Markov chain as a degenerate hidden Markov model where all the b proba-\n",
    "bilities are 1.0 for the observed symbol and 0 for all other symbols). Thus, the only probabilities we need to train are \n",
    "the transition probability matrix A.\n",
    "\"\"\"\n",
    "\"\"\"\n",
    "We get the maximum likelihood estimate of the probability a[i][j] of a particular transition between states i and j by \n",
    "counting the number of times the transition was taken, lets say C(i → j), and then normalizing by the total count of all\n",
    "times we took any transition from state i:\n",
    "    a[i][j]=C(i → j)/(Sum(C(i → q))forAll q∈Q)\n",
    "\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "For an HMM, we cannot compute these counts directly from an observation sequence since we don’t know which path of states \n",
    "was taken through the machine for a given input. The Baum-Welch algorithm uses two neat intuitions to solve this problem.\n",
    "The first idea is to iteratively estimate the counts. We will start with an estimate for the transition and observation \n",
    "probabilities and then use these estimated probabilities to derive better and better probabilities. The second idea is \n",
    "that we get our estimated probabilities by computing the forward probability for an observation and then dividing that \n",
    "probability mass among all the different paths that contributed to this forward probability.\n",
    "\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Lets define a useful probability related to the forward probability, called as Backward Probability(β).\n",
    "β[t][i]: Probability of seeing the observations from t + 1 to end(T), given \n",
    "            that we are in state i at time t (and given λ )\n",
    "β[t][i] = P(o[t+1] , o[t+2] . . . o[T] |q[t] = i, λ=(A,B) )\n",
    "\"\"\"\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Similar to alpha, beta also follows the Recursive SubProperty, containing lot of recomputations\n",
    "beta[t][i]= a[i][0]*b[0][o[t+1]]*beta[t+1][0] + a[i][1]*b[1][o[t+1]]*beta[t+1][1] + ......... +\n",
    "                a[i][N-2]*b[N-2][o[t+1]]*beta[t+1][N-2] + a[i][N-1]*b[N-1][o[t+1]]*beta[t+1][N-1]\n",
    "#Assuming we have n states, from 0 to n-1\n",
    "\"\"\"\n",
    "\"\"\"\"\"\"\n",
    "def computeBeta(observations,a,b,pi,betaDP):\n",
    "    statesC=a.shape[0]\n",
    "    timePts=observations.shape[0]\n",
    "    if timePts<1:\n",
    "        return\n",
    "    for state in np.arange(statesC):\n",
    "            betaDP[timePts-1][state]=1\n",
    "    for t in np.arange(timePts-2,-1,-1):\n",
    "        for i in np.arange(statesC):\n",
    "            for j in np.arange(statesC):\n",
    "                #betaDP[t][i]=betaDP[t][i]+a[i][j]*b[j][observations[t+1]]*betaDP[t+1][j]\n",
    "                betaDP[t][i]+=pow(2,log2(a[i][j])+log2(b[j][observations[t+1]])+log2(betaDP[t+1][j]))\n",
    "    return betaDP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Compute Transition Probability: a[i][j]\n",
    "\"\"\"\n",
    "a[i][j]=Expected No of Transitions from State i to State j/Expected No of Transitions from State i=aNum(i,j)/aDenom(i)\n",
    "diGamma(t,i,j) : Prob. of being in state i at time t and state j at time t+1\n",
    "diGamma(t,i,j)=P(q[t]=i,q[t+1]=j|O,λ)=p(q[t]=i,q[t+1]=j,O|λ)/P(O|λ)=diGammaNum(t,i,j)/diGammaDenom()\n",
    "diGammaNum(t,i,j)=p(q[t]=i,q[t+1]=j,O|λ)=alphaDP[t][i]*a[i][j]*b[j][observations[t+1]]*betaDP[t+1][j]\n",
    "diGammaDenom=P(O|λ)=observationsLikelihood(alphaDP) i.e. Sum of alphaDP[t-1][i] over all i s.t. i is a state\n",
    "aNum(i,j)=Sum of DiGamma(t,i,j) over all t's s.t. 0<=t<T-1\n",
    "aDenom(i)=Sum of DiGamma(t,i,k) over all t's and all k's s.t. 0<=t<T-1, 0<=k<=N-1\n",
    "\"\"\"\n",
    "def computeDiGammaNum(t,i,j,alphaDP,detaDP,a,b,observations):\n",
    "    return alphaDP[t][i]*a[i][j]*b[j][observations[t+1]]*betaDP[t+1][j]\n",
    "def computeDiGammaDP(alphaDP,betaDP,a,b,observations):\n",
    "    observationsC=alphaDP.shape[0]\n",
    "    statesC=alphaDP.shape[1]\n",
    "    diGammaDP=np.ndarray(shape=(statesC,statesC),dtype=float)\n",
    "    diGammaDenom=observationsLikelihood(alphaDP)\n",
    "    for i in np.arange(statesC):\n",
    "        for j in np.arange(statesC):\n",
    "                for t in np.arange(observationsC-1):\n",
    "                    diGammaDP[i][j]+=(computeDiGammaNum(t,i,j,alphaDP,betaDP,a,b,observations)/diGammaDenom)\n",
    "    return diGammaDP\n",
    "def computeTransitionProbabilityA(alphaDP,betaDP,a,b,observations):\n",
    "    statesC=alphaDP.shape[1]\n",
    "    newlyComputedTransitionProbA=np.ndarray(shape=(statesC,statesC),dtype=float)\n",
    "    diGammaDP=computeDiGammaDP(alphaDP,betaDP,a,b,observations)\n",
    "    diGammaDPSumGrpByJ=np.ndarray(shape=(statesC),dtype=float)\n",
    "    for i in np.arange(statesC):\n",
    "        for j in np.arange(statesC):\n",
    "            diGammaDPSumGrpByJ[i]+=diGammaDP[i][j]\n",
    "    for i in np.arange(statesC):    \n",
    "        for j in np.arange(statesC):\n",
    "            newlyComputedTransitionProbA[i][j]=diGammaDP[i][j]/diGammaDPSumGrpByJ[i]\n",
    "    return newlyComputedTransitionProbA   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Compute Observation Probability: b[i][O[k]]~ b[i][vk] i.e. the Probability of observing the symbol vk at state i\n",
    "\"\"\"\n",
    "b[i][vk]=Expected no of times in state i and observing vk/Expected no of times in State i\n",
    "            =ObsrProbNum(i,vk)/ObsrProbDenom(i)\n",
    "gamma[t][j]=Probability of being in state j at time t\n",
    "gamma[t][j]=p(q[t]=j|O,λ)=p(q[t]=j,O|λ)/p(O|λ)=gammaNum(t,j)/gammaDenom()\n",
    "gammaNum(t,j)=p(q[t]=j,O|λ)=alphaDP[t][j]*betaDP[t][j]\n",
    "gammaDenom()=observationsLikelihood(alphaDP) i.e. Sum of alphaDP[t-1][i] over all i s.t. i is a state\n",
    "obsrProbDenom(i)= Sum of gamma[t][i] over all t s.t. 0<=t<T\n",
    "obsrProbNum(i,vk)=Sum of gamma[t][i] over all t s.t. 0<=t<T and O[t]=vk\n",
    "\"\"\"\n",
    "def computeGammaNum(t,j,alphaDP,detaDP):\n",
    "    return alphaDP[t][j]*betaDP[t][j]\n",
    "def computeGammaDP(alphaDP,betaDP):\n",
    "    observationsC=alphaDP.shape[0]\n",
    "    statesC=alphaDP.shape[1]\n",
    "    gammaDP=np.ndarray(shape=(statesC,observationsC),dtype=float)\n",
    "    gammaDenom=observationsLikelihood(alphaDP)\n",
    "    for i in np.arange(statesC):\n",
    "        for t in np.arange(observationsC):\n",
    "            gammaDP[i][t]=computeGammaNum(t,i,alphaDP,betaDP)/gammaDenom  \n",
    "    return gammaDP\n",
    "def computeObsrProbNum(gammaDP,i,vk,observations):\n",
    "    observationC=len(observations)\n",
    "    obsrProbNum=0.0\n",
    "    for t in np.arange(observationC):\n",
    "        if observations[t]==vk:\n",
    "            obsrProbNum+=gammaDP[i][t]\n",
    "    return obsrProbNum\n",
    "def computeTransitionProbabilityB(alphaDP,betaDP,a,b,observations):\n",
    "    statesC=a.shape[0]\n",
    "    observationsC=b.shape[1]\n",
    "    newlyComputedObsrProbB=np.ndarray(shape=(statesC,observationsC),dtype=float)\n",
    "    gammaDP=computeGammaDP(alphaDP,betaDP) \n",
    "    for i in np.arange(statesC):\n",
    "        obsrProbDenom =np.sum(gammaDP[i])\n",
    "        for vk in observations:\n",
    "            newlyComputedObsrProbB[i][vk]=computeObsrProbNum(gammaDP,i,vk,observations)/obsrProbDenom\n",
    "    return newlyComputedObsrProbB   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.42772562  0.37468667]\n",
      " [ 0.38746073  0.47755799]]\n",
      "[[ 0.23425023  0.45724387  0.3085059 ]\n",
      " [ 0.70813797  0.23628532  0.05557671]]\n"
     ]
    }
   ],
   "source": [
    "observations=np.array([0,1,0,2,1,0])\n",
    "A=np.array([[0.7,0.3],[0.4,0.6]])\n",
    "B=np.array([[0.1,0.4,0.5],[0.7,0.2,0.1]])\n",
    "pi=np.array([0.6,0.4])\n",
    "alphaDP=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "betaDP=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "computeAlpha(observations,A,B,pi,alphaDP)\n",
    "computeBeta(observations,A,B,pi,betaDP)\n",
    "print(computeTransitionProbabilityA(alphaDP,betaDP,A,B,observations))\n",
    "print(computeTransitionProbabilityB(alphaDP,betaDP,A,B,observations))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Change Convergence Criteria to be more reasonable/Useful\n",
    "def isConverged(count):\n",
    "    if count>=5:\n",
    "        return True\n",
    "    return False\n",
    "def Forward_Backward_EM_Algo(observations,A,B,pi):\n",
    "    count=0\n",
    "    while isConverged(count)==False:\n",
    "        print(\"State Transition Matrix (A) ===>\")\n",
    "        print(A)\n",
    "        print(\"Observation Probability Matrix (B) ===>\")\n",
    "        print(B)\n",
    "        #Expectation(E)-Step\n",
    "        alphaDP=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "        betaDP=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "        computeAlpha(observations,A,B,pi,alphaDP)\n",
    "        computeBeta(observations,A,B,pi,betaDP)\n",
    "        #Maximization(M)-Step\n",
    "        newA=computeTransitionProbabilityA(alphaDP,betaDP,A,B,observations)\n",
    "        newB=computeTransitionProbabilityB(alphaDP,betaDP,A,B,observations)\n",
    "        A=newA\n",
    "        B=newB\n",
    "        count=count+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "State Transition Matrix (A) ===>\n",
      "[[ 0.7  0.3]\n",
      " [ 0.4  0.6]]\n",
      "Observation Probability Matrix (B) ===>\n",
      "[[ 0.1  0.4  0.5]\n",
      " [ 0.7  0.2  0.1]]\n",
      "State Transition Matrix (A) ===>\n",
      "[[ 0.53304969  0.46695031]\n",
      " [ 0.44792178  0.55207822]]\n",
      "Observation Probability Matrix (B) ===>\n",
      "[[ 0.23425023  0.45724387  0.3085059 ]\n",
      " [ 0.70813797  0.23628532  0.05557671]]\n",
      "State Transition Matrix (A) ===>\n",
      "[[ 0.45285207  0.54714793]\n",
      " [ 0.4825146   0.5174854 ]]\n",
      "Observation Probability Matrix (B) ===>\n",
      "[[ 0.38237174  0.37125943  0.24636883]\n",
      " [ 0.70945176  0.25178477  0.03876347]]\n",
      "State Transition Matrix (A) ===>\n",
      "[[ 0.4256403   0.5743597 ]\n",
      " [ 0.51099754  0.48900246]]\n",
      "Observation Probability Matrix (B) ===>\n",
      "[[ 0.50922725  0.27045219  0.22032056]\n",
      " [ 0.6781742   0.29022522  0.03160057]]\n",
      "State Transition Matrix (A) ===>\n",
      "[[ 0.40804421  0.59195579]\n",
      " [ 0.53939866  0.46060134]]\n",
      "Observation Probability Matrix (B) ===>\n",
      "[[ 0.614906    0.18607394  0.19902006]\n",
      " [ 0.61632299  0.35584658  0.02783043]]\n"
     ]
    }
   ],
   "source": [
    "observations=np.array([0,1,0,2,1,0])\n",
    "A=np.array([[0.7,0.3],[0.4,0.6]])\n",
    "B=np.array([[0.1,0.4,0.5],[0.7,0.2,0.1]])\n",
    "pi=np.array([0.6,0.4])\n",
    "Forward_Backward_EM_Algo(observations,A,B,pi)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}