DIY/HMM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HMM三个基本问题\n",
    "\n",
    "----------\n",
    "#### 1.概率计算问题\n",
    "盒子和球模型λ=(A,B,π)，状态集合Q=｛1，2，3｝，观测集合V=｛红，白｝，\n",
    "\n",
    "A = [ 0.5 0.2 0.3           \n",
    "\n",
    "     0.3 0.5 0.2\n",
    "     \n",
    "     0.2 0.3 0.5]\n",
    "     \n",
    "B = [0.5 0.5\n",
    "\n",
    "     0.4 0.6\n",
    "     \n",
    "     0.7 0.3]\n",
    "     \n",
    "π=(0.2,0.4,0.4)^T\n",
    "\n",
    "设T=3，0=（红，白，红）。计算P(O|π)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.1       0.16      0.28    ]\n",
      " [ 0.077     0.1104    0.0606  ]\n",
      " [ 0.04187   0.035512  0.052836]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.130218"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "####### 1.1 前向算法 ########\n",
    "# 0.数值化\n",
    "A = np.array([[0.5,0.2,0.3],\n",
    "              [0.3,0.5,0.2],\n",
    "              [0.2,0.3,0.5]])\n",
    "B = np.array([[0.5,0.5],\n",
    "              [0.4,0.6],\n",
    "              [0.7,0.3]])\n",
    "pi = np.array([0.2,0.4,0.4]).reshape((3,1))\n",
    "lambd = [A,B,pi]\n",
    "Q = [0,1,2]\n",
    "V = ['r','w']\n",
    "V = [0,1]\n",
    "\n",
    "T = 3\n",
    "O = ('r','w','r')\n",
    "O = [0,1,0]\n",
    "\n",
    "def forwardCal(lambd,O):\n",
    "    T = len(O) #时间序列长度\n",
    "    Q = len(lambd[0][0,:])#隐含状态种类 \n",
    "    #V = len(lambd[1][0,:])#观测状态种类\n",
    "    alpha = np.zeros((T,Q))\n",
    "    # 1.计算初值\n",
    "    for i in range(T):\n",
    "        alpha[O[0]][i] = pi[i]*lambd[1][i][O[0]]\n",
    "    # 2.递推计算\n",
    "    for i in range(1,T):\n",
    "        for j in range(Q):\n",
    "            alpha[i][j] = np.sum(alpha[i-1,:]*lambd[0][:,j])*lambd[1][j,O[i]]\n",
    "    # 3.终止\n",
    "    print(alpha)\n",
    "    P = np.sum(alpha[-1,:])\n",
    "    return P\n",
    "\n",
    "forwardCal(lambd,O)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.      1.      1.    ]\n",
      " [ 0.54    0.49    0.57  ]\n",
      " [ 0.2451  0.2622  0.2277]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.130218"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "####### 1.2 后向算法 ########\n",
    "def backwardCal(lambd,O):\n",
    "    T = len(O) #时间序列长度\n",
    "    Q = len(lambd[0][0,:])#隐含状态种类 \n",
    "    #V = len(lambd[1][0,:])#观测状态种类\n",
    "    beta = np.zeros((T,Q))\n",
    "    # 1.赋初值\n",
    "    beta[0,:] = 1.0  \n",
    "    # 2.递推计算\n",
    "    for i in range(1,T):\n",
    "        for j in range(Q):\n",
    "            beta[i][j] = np.sum(beta[i-1,:]*lambd[0][j,:].T*lambd[1][:,O[-i]])\n",
    "    # 3.终止\n",
    "    print(beta)\n",
    "    P = np.sum(beta[-1,:]*pi.T*lambd[1][:,O[0]])\n",
    "    return P\n",
    "\n",
    "backwardCal(lambd,O)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.学习问题"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "####### 2.1 Baum-Welch算法 ########\n",
    "#同EM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.预测问题"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.1      0.16     0.28   ]\n",
      " [ 0.028    0.0504   0.042  ]\n",
      " [ 0.00756  0.01008  0.0147 ]]\n",
      "[0, array([2, 2, 2], dtype=int64), array([1, 1, 2], dtype=int64)]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[2, 2, 2]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "####### 3.1 维特比算法 ########\n",
    "#已知A,B,pi,O,求最优路径I*\n",
    "def viterbi(lambd,O):\n",
    "    A = lambd[0]\n",
    "    B = lambd[1]\n",
    "    pi = lambd[2]\n",
    "    T = len(O) #时间序列长度\n",
    "    Q = len(lambd[0][0,:])#隐含状态种类 \n",
    "    \n",
    "    delte = np.zeros((T,Q))\n",
    "    psi = [0]\n",
    "    \n",
    "    #初始化：\n",
    "    delte[0,:] = (pi*B[:,O[0]].reshape((3,1))).T\n",
    "    \n",
    "    #递推\n",
    "    for t in range(1,T):\n",
    "        delte[t,:] = np.max(delte[t-1,:]*A.T,axis=1)*B[:,O[t]]\n",
    "        psi_tmp = np.argmax(delte[t-1,:]*A.T,axis=1)\n",
    "        psi.append(psi_tmp)\n",
    "    #终止\n",
    "    print(delte)\n",
    "    print(psi)\n",
    "    P_end = np.max(delte[-1,:])\n",
    "    i_end = np.argmax(delte[-1,:])\n",
    "    I = []\n",
    "    for p in psi[1:]:\n",
    "        I.append(np.max(p))\n",
    "    I.append(i_end)\n",
    "    #回溯\n",
    "    return I\n",
    "    \n",
    "    \n",
    "viterbi(lambd,O)"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}