Example/Example 2 -  complete implementation of HMM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## complete implementation 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "python my HMM\n",
      "BaumWelch\n",
      "num iteration  2140\n",
      "==================================================\n",
      "HMM content: N = 2 ,M = 2\n",
      "hmm.A  [ 0.59077362  0.40922638]  hmm.B  [ 0.9636068  0.0363932]\n",
      "       [ 0.50971332  0.49028668]         [ 0.02315302  0.97684698]\n",
      "hmm.pi [ 0.33220246  0.66779754]\n",
      "==================================================\n"
     ]
    }
   ],
   "source": [
    "#-*-coding:UTF-8-*-\n",
    "'''\n",
    "Created on 2014年9月25日\n",
    "@author: Ayumi Phoenix\n",
    "'''\n",
    "import numpy as np\n",
    "\n",
    "class HMM:\n",
    "    def __init__(self, Ann, Bnm, pi1n):\n",
    "        self.A = np.array(Ann)\n",
    "        self.B = np.array(Bnm)\n",
    "        self.pi = np.array(pi1n)\n",
    "        self.N = self.A.shape[0]\n",
    "        self.M = self.B.shape[1]\n",
    "        \n",
    "    def printhmm(self):\n",
    "        print (\"==================================================\")\n",
    "        print (\"HMM content: N =\",self.N,\",M =\",self.M)\n",
    "        for i in range(self.N):\n",
    "            if i==0:\n",
    "                print (\"hmm.A \",self.A[i,:],\" hmm.B \",self.B[i,:])\n",
    "            else:\n",
    "                print (\"      \",self.A[i,:],\"       \",self.B[i,:])\n",
    "        print (\"hmm.pi\",self.pi)\n",
    "        print (\"==================================================\")\n",
    "\n",
    "    # 函数名称：Forward *功能：前向算法估计参数 *参数:phmm:指向HMM的指针\n",
    "    # T:观察值序列的长度 O:观察值序列    \n",
    "    # alpha:运算中用到的临时数组 pprob:返回值,所要求的概率\n",
    "    def Forward(self,T,O,alpha,pprob):\n",
    "    #     1. Initialization 初始化\n",
    "        for i in range(self.N):\n",
    "            alpha[0,i] = self.pi[i]*self.B[i,O[0]]\n",
    "    \n",
    "    #     2. Induction 递归\n",
    "        for t in range(T-1):\n",
    "            for j in range(self.N):\n",
    "                sum = 0.0\n",
    "                for i in range(self.N):\n",
    "                    sum += alpha[t,i]*self.A[i,j]\n",
    "                alpha[t+1,j] =sum*self.B[j,O[t+1]]\n",
    "    #     3. Termination 终止\n",
    "        sum = 0.0\n",
    "        for i in range(self.N):\n",
    "            sum += alpha[T-1,i]\n",
    "        pprob[0] *= sum\n",
    "    \n",
    "    # 带修正的前向算法\n",
    "    def ForwardWithScale(self,T,O,alpha,scale,pprob):\n",
    "        scale[0] = 0.0\n",
    "    #     1. Initialization\n",
    "        for i in range(self.N):\n",
    "            alpha[0,i] = self.pi[i]*self.B[i,O[0]]\n",
    "            scale[0] += alpha[0,i]\n",
    "        \n",
    "        for i in range(self.N):\n",
    "            alpha[0,i] /= scale[0]\n",
    "        \n",
    "    #     2. Induction\n",
    "        for t in range(T-1):\n",
    "            scale[t+1] = 0.0\n",
    "            for j in range(self.N):\n",
    "                sum = 0.0\n",
    "                for i in range(self.N):\n",
    "                    sum += alpha[t,i]*self.A[i,j]\n",
    "                \n",
    "                alpha[t+1,j] = sum * self.B[j,O[t+1]]\n",
    "                scale[t+1] += alpha[t+1,j]\n",
    "            for j in range(self.N):\n",
    "                alpha[t+1,j] /= scale[t+1]\n",
    "        \n",
    "    #     3. Termination\n",
    "        for t in range(T):\n",
    "            pprob[0] += np.log(scale[t])\n",
    "\n",
    "    # 函数名称：Backward * 功能:后向算法估计参数 * 参数:phmm:指向HMM的指针 \n",
    "    # T:观察值序列的长度 O:观察值序列 \n",
    "    # beta:运算中用到的临时数组 pprob:返回值，所要求的概率\n",
    "    def Backword(self,T,O,beta,pprob):\n",
    "    #     1. Intialization\n",
    "        for i in range(self.N):\n",
    "            beta[T-1,i] = 1.0\n",
    "    #     2. Induction\n",
    "        for t in range(T-2,-1,-1):\n",
    "            for i in range(self.N):\n",
    "                sum = 0.0\n",
    "                for j in range(self.N):\n",
    "                    sum += self.A[i,j]*self.B[j,O[t+1]]*beta[t+1,j]\n",
    "                beta[t,i] = sum\n",
    "                \n",
    "    #     3. Termination\n",
    "        pprob[0] = 0.0\n",
    "        for i in range(self.N):\n",
    "            pprob[0] += self.pi[i]*self.B[i,O[0]]*beta[0,i]\n",
    "    \n",
    "    # 带修正的后向算法\n",
    "    def BackwardWithScale(self,T,O,beta,scale):\n",
    "    #     1. Intialization\n",
    "        for i in range(self.N):\n",
    "            beta[T-1,i] = 1.0\n",
    "    \n",
    "    #     2. Induction\n",
    "        for t in range(T-2,-1,-1):\n",
    "            for i in range(self.N):\n",
    "                sum = 0.0\n",
    "                for j in range(self.N):\n",
    "                    sum += self.A[i,j]*self.B[j,O[t+1]]*beta[t+1,j]\n",
    "                beta[t,i] = sum / scale[t+1]\n",
    "    \n",
    "    # Viterbi算法\n",
    "    # 输入：A，B，pi,O 输出P(O|lambda)最大时Poptimal的路径I\n",
    "    def viterbi(self,O):\n",
    "        T = len(O)\n",
    "        # 初始化\n",
    "        delta = np.zeros((T,self.N),np.float)  \n",
    "        phi = np.zeros((T,self.N),np.float)  \n",
    "        I = np.zeros(T)\n",
    "        for i in range(self.N):  \n",
    "            delta[0,i] = self.pi[i]*self.B[i,O[0]]  \n",
    "            phi[0,i] = 0\n",
    "        # 递推\n",
    "        for t in range(1,T):  \n",
    "            for i in range(self.N):                                  \n",
    "                delta[t,i] = self.B[i,O[t]]*np.array([delta[t-1,j]*self.A[j,i]  for j in range(self.N)]).max()\n",
    "                phi[t,i] = np.array([delta[t-1,j]*self.A[j,i]  for j in range(self.N)]).argmax()\n",
    "        # 终结\n",
    "        prob = delta[T-1,:].max()  \n",
    "        I[T-1] = delta[T-1,:].argmax()\n",
    "        # 状态序列求取   \n",
    "        for t in range(T-2,-1,-1): \n",
    "            I[t] = phi[t+1,I[t+1]]\n",
    "        return I,prob\n",
    "    \n",
    "    # 计算gamma : 时刻t时马尔可夫链处于状态Si的概率    \n",
    "    def ComputeGamma(self, T, alpha, beta, gamma):\n",
    "        for t in range(T):\n",
    "            denominator = 0.0\n",
    "            for j in range(self.N):\n",
    "                gamma[t,j] = alpha[t,j]*beta[t,j]\n",
    "                denominator += gamma[t,j]\n",
    "            for i in range(self.N):\n",
    "                gamma[t,i] = gamma[t,i]/denominator\n",
    "    \n",
    "    # 计算sai(i,j) 为给定训练序列O和模型lambda时：\n",
    "    # 时刻t是马尔可夫链处于Si状态，二时刻t+1处于Sj状态的概率\n",
    "    def ComputeXi(self,T,O,alpha,beta,gamma,xi):\n",
    "        for t in range(T-1):\n",
    "            sum = 0.0\n",
    "            for i in range(self.N):\n",
    "                for j in range(self.N):\n",
    "                    xi[t,i,j] = alpha[t,i]*beta[t+1,j]*self.A[i,j]*self.B[j,O[t+1]]\n",
    "                    sum += xi[t,i,j]\n",
    "            for i in range(self.N):\n",
    "                for j in range(self.N):\n",
    "                    xi[t,i,j] /= sum\n",
    "                    \n",
    "    # Baum-Welch算法\n",
    "    # 输入 L个观察序列O，初始模型：HMM={A,B,pi,N,M}\n",
    "    def BaumWelch(self,L,T,O,alpha,beta,gamma):\n",
    "        print (\"run BaumWelch\")\n",
    "        DELTA = 0.01 ; round = 0 ; flag = 1 ; probf = [0.0]\n",
    "        delta = 0.0 ; deltaprev = 0.0 ; probprev = 0.0 ; ratio = 0.0 ; deltaprev = 10e-70\n",
    "        \n",
    "        xi = np.zeros((T,self.N,self.N))\n",
    "        pi = np.zeros((T),np.float)\n",
    "        denominatorA = np.zeros((self.N),np.float)\n",
    "        denominatorB = np.zeros((self.N),np.float)\n",
    "        numeratorA = np.zeros((self.N,self.N),np.float)\n",
    "        numeratorB = np.zeros((self.N,self.M),np.float)\n",
    "        scale = np.zeros((T),np.float)\n",
    "        \n",
    "        while True :\n",
    "            probf[0] = 0\n",
    "            # E - step\n",
    "            for l in range(L):\n",
    "                self.ForwardWithScale(T,O[l],alpha,scale,probf)\n",
    "                self.BackwardWithScale(T,O[l],beta,scale)\n",
    "                self.ComputeGamma(T,alpha,beta,gamma)\n",
    "                self.ComputeXi(T,O[l],alpha,beta,gamma,xi)\n",
    "                for i in range(self.N):\n",
    "                    pi[i] += gamma[0,i]\n",
    "                    for t in range(T-1): \n",
    "                        denominatorA[i] += gamma[t,i]\n",
    "                        denominatorB[i] += gamma[t,i]\n",
    "                    denominatorB[i] += gamma[T-1,i]\n",
    "                    \n",
    "                    for j in range(self.N):\n",
    "                        for t in range(T-1):\n",
    "                            numeratorA[i,j] += xi[t,i,j]\n",
    "                    for k in range(self.M):\n",
    "                        for t in range(T):\n",
    "                            if O[l][t] == k:\n",
    "                                numeratorB[i,k] += gamma[t,i]\n",
    "                            \n",
    "            # M - step\n",
    "            # 重估状态转移矩阵 和 观察概率矩阵\n",
    "            for i in range(self.N):\n",
    "                self.pi[i] = 0.001/self.N + 0.999*pi[i]/L\n",
    "                for j in range(self.N):\n",
    "                    self.A[i,j] = 0.001/self.N + 0.999*numeratorA[i,j]/denominatorA[i]\n",
    "                    numeratorA[i,j] = 0.0\n",
    "                \n",
    "                for k in range(self.M):\n",
    "                    self.B[i,k] = 0.001/self.M + 0.999*numeratorB[i,k]/denominatorB[i]\n",
    "                    numeratorB[i,k] = 0.0\n",
    "                \n",
    "                pi[i]=denominatorA[i]=denominatorB[i]=0.0;\n",
    "            \n",
    "            if flag == 1:\n",
    "                flag = 0\n",
    "                probprev = probf[0]\n",
    "                ratio = 1\n",
    "                continue\n",
    "            \n",
    "            delta = probf[0] - probprev\n",
    "            ratio = delta / deltaprev\n",
    "            probprev = probf[0]\n",
    "            deltaprev = delta\n",
    "            round += 1\n",
    "            \n",
    "            if ratio <= DELTA :\n",
    "                print (\"num iteration \",round)\n",
    "                break\n",
    "         \n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print (\"python my HMM\")\n",
    "  \n",
    "    A = [[0.8125,0.1875],[0.2,0.8]]\n",
    "    B = [[0.875,0.125],[0.25,0.75]]\n",
    "    pi = [0.5,0.5]\n",
    "    hmm = HMM(A,B,pi)\n",
    "    \n",
    "    O = [[1,0,0,1,1,0,0,0,0],\n",
    "     [1,1,0,1,0,0,1,1,0],\n",
    "     [0,0,1,1,0,0,1,1,1]]\n",
    "    L = len(O)\n",
    "    T = len(O[0])  # T等于最长序列的长度就好了\n",
    "    alpha = np.zeros((T,hmm.N),np.float)\n",
    "    beta = np.zeros((T,hmm.N),np.float)\n",
    "    gamma = np.zeros((T,hmm.N),np.float)\n",
    "    hmm.BaumWelch(L,T,O,alpha,beta,gamma)\n",
    "    \n",
    "    hmm.printhmm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## complete implementation 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "DELTA = 0.001\n",
    "\n",
    "class HMM:\n",
    "\n",
    "\n",
    "    def __init__(self, pi, A, B):\n",
    "        self.pi = pi\n",
    "        self.A = A\n",
    "        self.B = B\n",
    "        self.M = B.shape[1]\n",
    "        self.N = A.shape[0]\n",
    "        \n",
    "    def forward(self,obs):\n",
    "        T = len(obs)\n",
    "        N = self.N\n",
    "        \n",
    "        alpha = np.zeros([N,T])\n",
    "        alpha[:,0] = self.pi[:] * self.B[:,obs[0]-1]                                                                                                      \n",
    "    \n",
    "        for t in xrange(1,T):\n",
    "            for n in xrange(0,N):\n",
    "                alpha[n,t] = np.sum(alpha[:,t-1] * self.A[:,n]) * self.B[n,obs[t]-1]\n",
    "                     \n",
    "        prob = np.sum(alpha[:,T-1])\n",
    "        return prob, alpha\n",
    "        \n",
    "    def forward_with_scale(self, obs):\n",
    "        \"\"\"see scaling chapter in \"A tutorial on hidden Markov models and \n",
    "        selected applications in speech recognition.\" \n",
    "        \"\"\"\n",
    "        T = len(obs)\n",
    "        N = self.N\n",
    "        alpha = np.zeros([N,T])\n",
    "        scale = np.zeros(T)\n",
    "\n",
    "        alpha[:,0] = self.pi[:] * self.B[:,obs[0]-1]\n",
    "        scale[0] = np.sum(alpha[:,0])\n",
    "        alpha[:,0] /= scale[0]\n",
    "\n",
    "        for t in xrange(1,T):\n",
    "            for n in xrange(0,N):\n",
    "                alpha[n,t] = np.sum(alpha[:,t-1] * self.A[:,n]) * self.B[n,obs[t]-1]\n",
    "            scale[t] = np.sum(alpha[:,t])\n",
    "            alpha[:,t] /= scale[t]\n",
    "\n",
    "        logprob = np.sum(np.log(scale[:]))\n",
    "        return logprob, alpha, scale    \n",
    "        \n",
    "    def backward(self, obs):\n",
    "        T = len(obs)\n",
    "        N = self.N\n",
    "        beta = np.zeros([N,T])\n",
    "        \n",
    "        beta[:,T-1] = 1\n",
    "        for t in reversed(xrange(0,T-1)):\n",
    "            for n in xrange(0,N):\n",
    "                beta[n,t] = np.sum(self.B[:,obs[t+1]-1] * self.A[n,:] * beta[:,t+1])\n",
    "                \n",
    "        prob = np.sum(beta[:,0])\n",
    "        return prob, beta\n",
    "\n",
    "    def backward_with_scale(self, obs, scale):\n",
    "        T = len(obs)\n",
    "        N = self.N\n",
    "        beta = np.zeros([N,T])\n",
    "\n",
    "        beta[:,T-1] = 1 / scale[T-1]\n",
    "        for t in reversed(xrange(0,T-1)):\n",
    "            for n in xrange(0,N):\n",
    "                beta[n,t] = np.sum(self.B[:,obs[t+1]-1] * self.A[n,:] * beta[:,t+1])\n",
    "                beta[n,t] /= scale[t]\n",
    "        \n",
    "        return beta\n",
    "\n",
    "    def viterbi(self, obs):\n",
    "        T = len(obs)\n",
    "        N = self.N\n",
    "        psi = np.zeros([N,T]) # reverse pointer\n",
    "        delta = np.zeros([N,T])\n",
    "        q = np.zeros(T)\n",
    "        temp = np.zeros(N)        \n",
    "        \n",
    "        delta[:,0] = self.pi[:] * self.B[:,obs[0]-1]    \n",
    "        \n",
    "        for t in xrange(1,T):\n",
    "            for n in xrange(0,N):\n",
    "                temp = delta[:,t-1] * self.A[:,n]    \n",
    "                max_ind = argmax(temp)\n",
    "                psi[n,t] = max_ind\n",
    "                delta[n,t] = self.B[n,obs[t]-1] * temp[max_ind]\n",
    "\n",
    "        max_ind = argmax(delta[:,T-1])\n",
    "        q[T-1] = max_ind\n",
    "        prob = delta[:,T-1][max_ind]\n",
    "\n",
    "        for t in reversed(xrange(1,T-1)):\n",
    "            q[t] = psi[q[t+1],t+1]    \n",
    "            \n",
    "        return prob, q, delta    \n",
    "        \n",
    "    def viterbi_log(self, obs):\n",
    "        \n",
    "        T = len(obs)\n",
    "        N = self.N\n",
    "        psi = np.zeros([N,T])\n",
    "        delta = np.zeros([N,T])\n",
    "        pi = np.zeros(self.pi.shape)\n",
    "        A = np.zeros(self.A.shape)\n",
    "        biot = np.zeros([N,T])\n",
    "\n",
    "        pi = np.log(self.pi)        \n",
    "        A = np.log(self.A)\n",
    "        biot = np.log(self.B[:,obs[:]-1])\n",
    "\n",
    "        delta[:,0] = pi[:] + biot[:,0]\n",
    "\n",
    "        for t in xrange(1,T):\n",
    "            for n in xrange(0,N):\n",
    "                temp = delta[:,t-1] + self.A[:,n]    \n",
    "                max_ind = argmax(temp)\n",
    "                psi[n,t] = max_ind\n",
    "                delta[n,t] = temp[max_ind] + biot[n,t]   \n",
    "\n",
    "        max_ind = argmax(delta[:,T-1])\n",
    "        q[T-1] = max_ind            \n",
    "        logprob = delta[max_ind,T-1]      \n",
    "                \n",
    "        for t in reversed(xrange(1,T-1)):\n",
    "            q[t] = psi[q[t+1],t+1]    \n",
    "\n",
    "        return logprob, q, delta\n",
    "\n",
    "    def baum_welch(self, obs):\n",
    "        T = len(obs)\n",
    "        M = self.M\n",
    "        N = self.N        \n",
    "        alpha = np.zeros([N,T])\n",
    "        beta = np.zeros([N,T])\n",
    "        scale = np.zeros(T)\n",
    "        gamma = np.zeros([N,T])\n",
    "        xi = np.zeros([N,N,T-1])\n",
    "    \n",
    "        # caculate initial parameters\n",
    "        logprobprev, alpha, scale = self.forward_with_scale(obs)\n",
    "        beta = self.backward_with_scale(obs, scale)            \n",
    "        gamma = self.compute_gamma(alpha, beta)    \n",
    "        xi = self.compute_xi(obs, alpha, beta)    \n",
    "        logprobinit = logprobprev        \n",
    "        \n",
    "        # start interative \n",
    "        while True:\n",
    "            # E-step\n",
    "            self.pi = 0.001 + 0.999*gamma[:,0]\n",
    "            for i in xrange(N):\n",
    "                denominator = np.sum(gamma[i,0:T-1])\n",
    "                for j in xrange(N): \n",
    "                    numerator = np.sum(xi[i,j,0:T-1])\n",
    "                    self.A[i,j] = numerator / denominator\n",
    "                                   \n",
    "            self.A = 0.001 + 0.999*self.A\n",
    "            for j in xrange(0,N):\n",
    "                denominator = np.sum(gamma[j,:])\n",
    "                for k in xrange(0,M):\n",
    "                    numerator = 0.0\n",
    "                    for t in xrange(0,T):\n",
    "                        if obs[t]-1 == k:\n",
    "                            numerator += gamma[j,t]\n",
    "                    self.B[j,k] = numerator / denominator\n",
    "            self.B = 0.001 + 0.999*self.B\n",
    "\n",
    "            # M-step\n",
    "            logprobcur, alpha, scale = self.forward_with_scale(obs)\n",
    "            beta = self.backward_with_scale(obs, scale)            \n",
    "            gamma = self.compute_gamma(alpha, beta)    \n",
    "            xi = self.compute_xi(obs, alpha, beta)    \n",
    "\n",
    "            delta = logprobcur - logprobprev\n",
    "            logprobprev = logprobcur\n",
    "            # print \"delta is \", delta\n",
    "            if delta <= DELTA:\n",
    "                break     \n",
    "                \n",
    "        logprobfinal = logprobcur\n",
    "        return logprobinit, logprobfinal                \n",
    "            \n",
    "    def compute_gamma(self, alpha, beta):\n",
    "        gamma = np.zeros(alpha.shape)\n",
    "        gamma = alpha[:,:] * beta[:,:]\n",
    "        gamma = gamma / np.sum(gamma,0)\n",
    "        return gamma\n",
    "            \n",
    "    def compute_xi(self, obs, alpha, beta):\n",
    "        T = len(obs)\n",
    "        N = self.N\n",
    "        xi = np.zeros((N, N, T-1))\n",
    "            \n",
    "        for t in xrange(0,T-1):        \n",
    "            for i in xrange(0,N):\n",
    "                for j in xrange(0,N):\n",
    "                    xi[i,j,t] = alpha[i,t] * self.A[i,j] * \\\n",
    "                                self.B[j,obs[t+1]-1] * beta[j,t+1]\n",
    "            xi[:,:,t] /= np.sum(np.sum(xi[:,:,t],1),0)    \n",
    "        return xi\n",
    "\n",
    "def read_hmm(hmmfile):\n",
    "    fhmm = open(hmmfile,'r') \n",
    "\n",
    "    M = int(fhmm.readline().split(' ')[1])\n",
    "    N = int(fhmm.readline().split(' ')[1]) \n",
    "    \n",
    "    A = np.array([])\n",
    "    fhmm.readline()\n",
    "    for i in xrange(N):\n",
    "        line = fhmm.readline()\n",
    "        if i == 0:\n",
    "            A = np.array(map(float,line.split(',')))\n",
    "        else:\n",
    "            A = np.vstack((A,map(float,line.split(','))))\n",
    "        \n",
    "    B = np.array([])\n",
    "    fhmm.readline()\n",
    "    for i in xrange(N):\n",
    "        line = fhmm.readline()\n",
    "        if i == 0:\n",
    "            B = np.array(map(float,line.split(',')))\n",
    "        else:\n",
    "            B = np.vstack((B,map(float,line.split(','))))\n",
    "    \n",
    "    fhmm.readline()\n",
    "    line = fhmm.readline()\n",
    "    pi = np.array(map(float,line.split(',')))\n",
    "    \n",
    "    fhmm.close()\n",
    "    return M, N, pi, A, B \n",
    "    \n",
    "def read_sequence(seqfile):\n",
    "    fseq = open(seqfile,'r') \n",
    "    \n",
    "    T = int(fseq.readline().split(' ')[1])\n",
    "    line = fseq.readline()\n",
    "    obs = np.array(map(int,line.split(',')))\n",
    "    \n",
    "    fseq.close()\n",
    "    return T, obs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## complete implementation 3\n",
    "参考自：https://github.com/hankcs/hidden-markov-model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------viterbi test-----------------------------------\n",
      "            normal       cold      dizzy\n",
      "Healthy:   0.300000   0.084000   0.005880\n",
      "  Fever:   0.040000   0.027000   0.015120\n",
      "\n",
      "The most possible states and probability are:\n",
      "Healthy Healthy Fever 0.01512\n",
      "update\n",
      "update\n",
      "\n",
      "\n",
      "-----------------------baum_welch_train test-------------------------------\n",
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      "[0 1 1 1 1 1 1 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 1 1 0 0 0\n",
      " 0 1 1 0 0 0 0 1 0 0 0 1 1 1 1 1 0 0 1 1 0 1 1 1 1 0 0 0 0 1 0 0 0 0 0 0 1\n",
      " 0 0 1 1 1 0 0 0 0 0 1 0 0 1 0 0 1 1 1 0 0 0 0 0 0 0]\n",
      "Accuracy： 0.610000000\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "class HMM:\n",
    "    \"\"\"\n",
    "    Order 1 Hidden Markov Model\n",
    "    Attributes\n",
    "    ----------\n",
    "    A : numpy.ndarray\n",
    "        State transition probability matrix\n",
    "    B: numpy.ndarray\n",
    "        Output emission probability matrix with shape(N, number of output types)\n",
    "    pi: numpy.ndarray\n",
    "        Initial state probablity vector\n",
    "    Common Variables\n",
    "    ----------------\n",
    "    obs_seq : list of int\n",
    "        list of observations (represented as ints corresponding to output\n",
    "        indexes in B) in order of appearance\n",
    "    T : int\n",
    "        number of observations in an observation sequence\n",
    "    N : int\n",
    "        number of states\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, A, B, pi):\n",
    "        self.A = A\n",
    "        self.B = B\n",
    "        self.pi = pi\n",
    "\n",
    "    def _forward(self, obs_seq):\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "\n",
    "        F = np.zeros((N,T))\n",
    "        F[:,0] = self.pi * self.B[:, obs_seq[0]]\n",
    "\n",
    "        for t in range(1, T):\n",
    "            for n in range(N):\n",
    "                F[n,t] = np.dot(F[:,t-1], (self.A[:,n])) * self.B[n, obs_seq[t]]\n",
    "\n",
    "        return F\n",
    "\n",
    "    def _backward(self, obs_seq):\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "\n",
    "        X = np.zeros((N,T))\n",
    "        X[:,-1:] = 1\n",
    "\n",
    "        for t in reversed(range(T-1)):\n",
    "            for n in range(N):\n",
    "                X[n,t] = np.sum(X[:,t+1] * self.A[n,:] * self.B[:, obs_seq[t+1]])\n",
    "\n",
    "        return X\n",
    "\n",
    "    def observation_prob(self, obs_seq):\n",
    "        \"\"\" P( entire observation sequence | A, B, pi ) \"\"\"\n",
    "        return np.sum(self._forward(obs_seq)[:,-1])\n",
    "\n",
    "    def state_path(self, obs_seq):\n",
    "        \"\"\"\n",
    "        Returns\n",
    "        -------\n",
    "        V[last_state, -1] : float\n",
    "            Probability of the optimal state path\n",
    "        path : list(int)\n",
    "            Optimal state path for the observation sequence\n",
    "        \"\"\"\n",
    "        V, prev = self.viterbi(obs_seq)\n",
    "\n",
    "        # Build state path with greatest probability\n",
    "        last_state = np.argmax(V[:,-1])\n",
    "        path = list(self.build_viterbi_path(prev, last_state))\n",
    "\n",
    "        return V[last_state,-1], reversed(path)\n",
    "\n",
    "    def viterbi(self, obs_seq):\n",
    "        \"\"\"\n",
    "        Returns\n",
    "        -------\n",
    "        V : numpy.ndarray\n",
    "            V [s][t] = Maximum probability of an observation sequence ending\n",
    "                       at time 't' with final state 's'\n",
    "        prev : numpy.ndarray\n",
    "            Contains a pointer to the previous state at t-1 that maximizes\n",
    "            V[state][t]\n",
    "        \"\"\"\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "        prev = np.zeros((T - 1, N), dtype=int)\n",
    "\n",
    "        # DP matrix containing max likelihood of state at a given time\n",
    "        V = np.zeros((N, T))\n",
    "        V[:,0] = self.pi * self.B[:,obs_seq[0]]\n",
    "\n",
    "        for t in range(1, T):\n",
    "            for n in range(N):\n",
    "                seq_probs = V[:,t-1] * self.A[:,n] * self.B[n, obs_seq[t]]\n",
    "                prev[t-1,n] = np.argmax(seq_probs)\n",
    "                V[n,t] = np.max(seq_probs)\n",
    "\n",
    "        return V, prev\n",
    "\n",
    "    def build_viterbi_path(self, prev, last_state):\n",
    "        \"\"\"Returns a state path ending in last_state in reverse order.\"\"\"\n",
    "        T = len(prev)\n",
    "        yield(last_state)\n",
    "        for i in range(T-1, -1, -1):\n",
    "            yield(prev[i, last_state])\n",
    "            last_state = prev[i, last_state]\n",
    "\n",
    "    def simulate(self, T):\n",
    "\n",
    "        def draw_from(probs):\n",
    "            return np.where(np.random.multinomial(1,probs) == 1)[0][0]\n",
    "\n",
    "        observations = np.zeros(T, dtype=int)\n",
    "        states = np.zeros(T, dtype=int)\n",
    "        states[0] = draw_from(self.pi)\n",
    "        observations[0] = draw_from(self.B[states[0],:])\n",
    "        for t in range(1, T):\n",
    "            states[t] = draw_from(self.A[states[t-1],:])\n",
    "            observations[t] = draw_from(self.B[states[t],:])\n",
    "        return observations,states\n",
    "\n",
    "    def baum_welch_train(self, observations, criterion=0.05):\n",
    "        n_states = self.A.shape[0]\n",
    "        n_samples = len(observations)\n",
    "\n",
    "        done = False\n",
    "        while not done:\n",
    "            # alpha_t(i) = P(O_1 O_2 ... O_t, q_t = S_i | hmm)\n",
    "            # Initialize alpha\n",
    "            alpha = self._forward(observations)\n",
    "\n",
    "            # beta_t(i) = P(O_t+1 O_t+2 ... O_T | q_t = S_i , hmm)\n",
    "            # Initialize beta\n",
    "            beta = self._backward(observations)\n",
    "\n",
    "            xi = np.zeros((n_states,n_states,n_samples-1))\n",
    "            for t in range(n_samples-1):\n",
    "                denom = np.dot(np.dot(alpha[:,t].T, self.A) * self.B[:,observations[t+1]].T, beta[:,t+1])\n",
    "                for i in range(n_states):\n",
    "                    numer = alpha[i,t] * self.A[i,:] * self.B[:,observations[t+1]].T * beta[:,t+1].T\n",
    "                    xi[i,:,t] = numer / denom\n",
    "\n",
    "            # gamma_t(i) = P(q_t = S_i | O, hmm)\n",
    "            gamma = np.squeeze(np.sum(xi,axis=1))\n",
    "            # Need final gamma element for new B\n",
    "            prod =  (alpha[:,n_samples-1] * beta[:,n_samples-1]).reshape((-1,1))\n",
    "            gamma = np.hstack((gamma,  prod / np.sum(prod))) #append one more to gamma!!!\n",
    "\n",
    "            newpi = gamma[:,0]\n",
    "            newA = np.sum(xi,2) / np.sum(gamma[:,:-1],axis=1).reshape((-1,1))\n",
    "            newB = np.copy(self.B)\n",
    "\n",
    "            num_levels = self.B.shape[1]\n",
    "            sumgamma = np.sum(gamma,axis=1)\n",
    "            for lev in range(num_levels):\n",
    "                mask = observations == lev\n",
    "                newB[:,lev] = np.sum(gamma[:,mask],axis=1) / sumgamma\n",
    "\n",
    "            if np.max(abs(self.pi - newpi)) < criterion and \\\n",
    "                            np.max(abs(self.A - newA)) < criterion and \\\n",
    "                            np.max(abs(self.B - newB)) < criterion:\n",
    "                done = 1\n",
    "\n",
    "            self.A[:],self.B[:],self.pi[:] = newA,newB,newpi\n",
    "            #print(\"update\") \n",
    "            \n",
    "states = ('Healthy', 'Fever')\n",
    "\n",
    "observations = ('normal', 'cold', 'dizzy')\n",
    "\n",
    "start_probability = {'Healthy': 0.6, 'Fever': 0.4}\n",
    "\n",
    "transition_probability = {\n",
    "    'Healthy': {'Healthy': 0.7, 'Fever': 0.3},\n",
    "    'Fever': {'Healthy': 0.4, 'Fever': 0.6},\n",
    "}\n",
    "\n",
    "emission_probability = {\n",
    "    'Healthy': {'normal': 0.5, 'cold': 0.4, 'dizzy': 0.1},\n",
    "    'Fever': {'normal': 0.1, 'cold': 0.3, 'dizzy': 0.6},\n",
    "}\n",
    "\n",
    "##############转换数据类型###########\n",
    "def generate_index_map(lables):\n",
    "    index_label = {}\n",
    "    label_index = {}\n",
    "    i = 0\n",
    "    for l in lables:\n",
    "        index_label[i] = l\n",
    "        label_index[l] = i\n",
    "        i += 1\n",
    "    return label_index, index_label\n",
    "\n",
    "\n",
    "states_label_index, states_index_label = generate_index_map(states)\n",
    "observations_label_index, observations_index_label = generate_index_map(observations)\n",
    "\n",
    "\n",
    "def convert_observations_to_index(observations, label_index):\n",
    "    list = []\n",
    "    for o in observations:\n",
    "        list.append(label_index[o])\n",
    "    return list\n",
    "\n",
    "\n",
    "def convert_map_to_vector(map, label_index):\n",
    "    v = np.empty(len(map), dtype=float)\n",
    "    for e in map:\n",
    "        v[label_index[e]] = map[e]\n",
    "    return v\n",
    "\n",
    "\n",
    "def convert_map_to_matrix(map, label_index1, label_index2):\n",
    "    m = np.empty((len(label_index1), len(label_index2)), dtype=float)\n",
    "    for line in map:\n",
    "        for col in map[line]:\n",
    "            m[label_index1[line]][label_index2[col]] = map[line][col]\n",
    "    return m\n",
    "\n",
    "\n",
    "A = convert_map_to_matrix(transition_probability, states_label_index, states_label_index)\n",
    "# print A\n",
    "B = convert_map_to_matrix(emission_probability, states_label_index, observations_label_index)\n",
    "# print B\n",
    "observations_index = convert_observations_to_index(observations, observations_label_index)\n",
    "pi = convert_map_to_vector(start_probability, states_label_index)\n",
    "# print pi\n",
    "##############转换数据类型###########\n",
    "print('--------------------------viterbi test-----------------------------------')\n",
    "h = HMM(A, B, pi)\n",
    "V, p = h.viterbi(observations_index)\n",
    "print(\" \" * 7, \" \".join((\"%10s\" % observations_index_label[i]) for i in observations_index))\n",
    "for s in range(0, 2):\n",
    "    print (\"%7s: \" % states_index_label[s] + \" \".join(\"%10s\" % (\"%f\" % v) for v in V[s]))\n",
    "print('\\nThe most possible states and probability are:')\n",
    "p, ss = h.state_path(observations_index)\n",
    "for s in ss:\n",
    "    print(states_index_label[s],end=' ')\n",
    "print(p)\n",
    "\n",
    "simulength = 100\n",
    "# run a baum_welch_train\n",
    "observations_data, states_data = h.simulate(simulength)\n",
    "# print observations_data\n",
    "# print states_data\n",
    "guess = HMM(np.array([[0.5, 0.5],\n",
    "                          [0.5, 0.5]]),\n",
    "                np.array([[0.3, 0.3, 0.3],\n",
    "                          [0.3, 0.3, 0.3]]),\n",
    "                np.array([0.5, 0.5])\n",
    "                )\n",
    "guess.baum_welch_train(observations_data)\n",
    "states_out = np.array(list(guess.state_path(observations_data)[1]))\n",
    "print('\\n\\n-----------------------baum_welch_train test-------------------------------')\n",
    "print(states_out) #这个变量输出来竟然都是0，应该是程序有bug\n",
    "print(states_data)\n",
    "\n",
    "compare=states_out-states_data\n",
    "#输出准确率\n",
    "print('Accuracy： %.9f' % (len(compare[compare==0]) / simulength))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}