{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "对于HMM的参数学习，其实分两种情况，第一种是训练数据包括了隐状态序列和观测状态序列，比如$S$个序列对： \n", "\n", "$$\n", "\\{(O_1,I_1),(O_2,I_2),...,(O_S,I_S)\\}\n", "$$ \n", "\n", "另一种是训练数据仅包含了观测状态序列，而未包含隐状态序列，比如： \n", "\n", "$$\n", "\\{O_1,O_2,...,O_S\\}\n", "$$ \n", "\n", "显然，两者均可使用极大似然估计进行参数求解，只是后者需要用到EM算法，下面分别对其做介绍" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 一.有监督的方法\n", "这种情况就比较简单了，直接对数据做统计即可\n", "\n", "#### 1.1 转移概率$a_{ij}$的估计\n", "假设样本中任意时刻$t$处于状态$i$在下一时刻转移到状态$j$的频数为$A_{ij}$，那么： \n", "\n", "\n", "$$\n", "a_{ij}=\\frac{A_{ij}}{\\sum_{j=1}^NA_{ij}},i=1,2,...,N,j=1,2,...,N\n", "$$ \n", "\n", "#### 1.2 观测概率$b_j(k)$的估计\n", "假设样本中状态为$j$并且观测为$k$的频数为$B_{jk}$，那么： \n", "\n", "$$\n", "b_{jk}=\\frac{B_{jk}}{\\sum_{k=1}^MB_{jk}},j=1,2,...,N,k=1,2,...,M\n", "$$ \n", "\n", "#### 1.3 初始状态概率$\\pi$分布\n", "直接计算样本中初始隐状态$q_i,i=1,2,...,N$的频率即可" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 二.无监督的方法\n", "\n", "这种情况我们可以使用EM算法，隐变量我们自然可以选择隐状态$I$咯，假设目前所有的观测数据为$O=(o_1,o_2,...,o_T)$，所有的隐数据为$I=(i_1,i_2,...,i_T)$，令$\\lambda^-$是HMM参数的当前估计值，$\\lambda$是需要极大化的HMM参数，那么$Q$函数我们就可以写出来了： \n", "\n", "#### 2. 1 写出Q函数\n", "\n", "$$\n", "Q(\\lambda,\\lambda^-)=\\sum_IlogP(O,I\\mid\\lambda)P(O,I\\mid\\lambda^-)\n", "$$ \n", "\n", "我们首先看看$P(O,I\\mid\\lambda)$，显然可以直接由HMM的概率图模型直接写出： \n", "\n", "$$\n", "P(O,I\\mid\\lambda)=\\pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2)\\cdots a_{i_{T-1}i_T}b_{i_T}(o_T)\n", "$$ \n", "\n", "对于第二项，在Q函数中，它应该为$P(I\\mid O,\\lambda^-)$才对，但由于$P(I\\mid O,\\lambda^-)=\\frac{P(I,O\\mid\\lambda^-)}{P(O\\mid\\lambda^-)}$中$P(O\\mid\\lambda^-)$为常数，对极大化$Q$函数没有影响 \n", "\n", "接下里，对Q函数化简： \n", "\n", "$$\n", "Q(\\lambda\\mid\\lambda^-)=\\sum_Ilog\\pi_{i_1}P(O,I\\mid\\lambda^-)+\\sum_I(\\sum_{t=1}^{T-1}loga_{i_ti_{t+1}})P(O,I\\mid\\lambda^-)+\\sum_I(\\sum_{t=1}^Tlogb_{i_t}(o_t))P(O,I\\mid\\lambda^-)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2 极大化Q函数\n", "\n", "接下来便是极大化Q函数了，由于HMM的参数单独地出现在3个项中，所以我们需要对各项单独极大化即可 \n", "\n", "##### 初始状态$\\pi$的更新\n", "第一项公式可以写作： \n", "\n", "$$\n", "\\sum_Ilog\\pi_{i_1}P(O,I\\mid\\lambda^-)=\\sum_{i=1}^Nlog\\pi_iP(O,i_1=i\\mid\\lambda^-)\\\\\n", "s.t. \\sum_{i=1}^N\\pi_i=1\n", "$$ \n", "通过构造拉格朗日函数求KKT条件可以得到： \n", "\n", "$$\n", "\\pi_i=\\frac{P(O,i_1=i\\mid\\lambda^-)}{P(O\\mid\\lambda^-)},i=1,2,...,N\n", "$$\n", "##### 隐状态概率矩阵$A$的更新\n", "\n", "对于第二项： \n", "\n", "$$\n", "\\sum_I(\\sum_{t=1}^{T-1}loga_{i_ti_{t+1}})P(O,I\\mid\\lambda^-)=\\sum_{i=1}^N\\sum_{j=1}^N\\sum_{t=1}^{T-1}loga_{ij}P(O,i_t=i,i_{t+1}=j\\mid\\lambda^-)\\\\\n", "s.t. \\sum_{j=1}^Na_{ij}=1,i=1,2,...,N\n", "$$ \n", "同理，通过构造拉格朗日函数并求KKT条件可得： \n", "\n", "$$\n", "a_{ij}=\\frac{P(O,i_t=i,i_{t+1}=j\\mid\\lambda^-)}{\\sum_{t=1}^{T-1}P(O,i_t=i\\mid\\lambda^-)}\n", "$$ \n", "\n", "##### 观测状态概率矩阵$B$的更新\n", "第三项： \n", "\n", "$$\n", "\\sum_I(\\sum_{t=1}^Tlogb_{i_t}(o_t))P(O,I\\mid\\lambda^-)=\\sum_{j=1}^N\\sum_{k=1}^M\\sum_{t=1}^Tlogb_j(k)P(O,i_t=j\\mid\\lambda^-)I(o_t=v_k)\\\\\n", "s.t.\\sum_{k=1}^Mb_j(k)=1,j=1,2,..,M\n", "$$ \n", "\n", "同理，求得： \n", "\n", "$$\n", "b_j(k)=\\frac{\\sum_{t=1}^TP(O,i_t=j\\mid\\lambda^-)I(o_t=v_k)}{\\sum_{t=1}^TP(O,i_t=j\\mid\\lambda^-)}\n", "$$ \n", "\n", "好的，到这里参数的求解就推导完了，接下里需要对上面的计算进行优化，我们发现上式的求解主要包含两种类型，一种是$P(O,i_t=i\\mid\\lambda^-)$的形式，另一种是$P(O,i_t=i,i_{t+1}=j\\mid\\lambda^-)$的形式，下面对其求解做介绍 \n", "\n", "#### 2.3 $P(O,i_t=i\\mid\\lambda^-)$和$P(O,i_t=i,i_{t+1}=j\\mid\\lambda^-)$的求解\n", "\n", "这里其实可以使用上一节介绍的前向概率$\\alpha_t(i)$和后向概率$\\beta_t(i)$求解（符号定义与上一节相同）\n", "\n", "##### $P(O,i_t=i\\mid\\lambda^-)$的求解\n", "对于$P(O,i_t=i\\mid\\lambda^-)$，其实上一节已经能给出答案了，它的求解如下图： \n", "![avatar](./source/12_HMM前向后向.png) \n", "\n", "所以只需要将$t$时刻的前向概率和后向概率乘起来即可，令： \n", "\n", "$$\n", "\\gamma_t(i)=P(O,i_t=i\\mid\\lambda^-)=\\alpha_t(i)\\beta_t(i)\n", "$$\n", "##### $P(O,i_t=i,i_{t+1}=j\\mid\\lambda^-)$的求解\n", "对于$P(O,i_t=i,i_{t+1}=j\\mid\\lambda^-)$的求解，可以让后向概率往后退一步，如下图： \n", "![avatar](./source/12_HMM前向后向2.png) \n", "\n", "这样，我们可以方便将结果写出来，令： \n", "\n", "$$\n", "\\xi_t(i,j)=P(O,i_t=i,i_{t+1}=j\\mid\\lambda^-)=\\alpha_t(i)a_{ij}b_j(o_j)\\beta_{t+1}(j)\n", "$$ \n", "\n", "这便是Baum-Welch算法，接下来对上面的无监督学习过程做总结，并梳理出算法流程\n", "\n", "#### 2.4 Baum-Welch算法\n", "已知观测数据$O=(o_1,o_2,...,o_T)$： \n", "\n", ">（1）初始化：对$n=0$，随机初始化 $a_{ij}^{(0)},b_j(k)^{(0)},\\pi_i^{(0)}$（注意满足概率的约束条件） \n", "\n", ">（2）递推，对$n=1,2,...$ \n", "\n", "$$\n", "a_{ij}^n=\\frac{\\sum_{t=1}^{T-1}\\xi_t(i,j)}{\\sum_{t=1}^{T-1}\\gamma_t(i)}\\\\\n", "b_j(k)^n=\\frac{\\sum_{t=1}^T\\gamma_t(j)I(o_t=v_k)}{\\sum_{t=1}^T\\gamma_t(j)}\\\\\n", "\\pi_i^n=\\gamma_1(i)\n", "$$ \n", "\n", ">（3）根据收敛条件进行终止" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 三.代码实现\n", "让我们接着上一节的代码继续写..." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "\"\"\"\n", "隐马尔科夫模型实现，封装到ml_models.pgm\n", "\"\"\"\n", "\n", "class HMM(object):\n", " def __init__(self, hidden_status_num=None, visible_status_num=None):\n", " \"\"\"\n", " :param hidden_status_num: 隐状态数\n", " :param visible_status_num: 观测状态数\n", " \"\"\"\n", " self.hidden_status_num = hidden_status_num\n", " self.visible_status_num = visible_status_num\n", " # 定义HMM的参数\n", " self.pi = None # 初始隐状态概率分布 shape:[hidden_status_num,1]\n", " self.A = None # 隐状态转移概率矩阵 shape:[hidden_status_num,hidden_status_num]\n", " self.B = None # 观测状态概率矩阵 shape:[hidden_status_num,visible_status_num]\n", "\n", " def predict_joint_visible_prob(self, visible_list=None, forward_type=\"forward\"):\n", " \"\"\"\n", " 前向/后向算法计算观测序列出现的概率值\n", " :param visible_list:\n", " :param forward_type:forward前向，backward后向\n", " :return:\n", " \"\"\"\n", " if forward_type == \"forward\":\n", " # 计算初始值\n", " alpha = self.pi * self.B[:, [visible_list[0]]]\n", " # 递推\n", " for step in range(1, len(visible_list)):\n", " alpha = self.A.T.dot(alpha) * self.B[:, [visible_list[step]]]\n", " # 求和\n", " return np.sum(alpha)\n", " else:\n", " # 计算初始值\n", " beta = np.ones_like(self.pi)\n", " # 递推\n", " for step in range(len(visible_list) - 2, -1, -1):\n", " beta = self.A.dot(self.B[:, [visible_list[step + 1]]] * beta)\n", " # 最后一步\n", " return np.sum(self.pi * self.B[:, [visible_list[0]]] * beta)\n", "\n", " def fit_with_hidden_status(self, visible_list, hidden_list):\n", " \"\"\"\n", " 包含隐状态的参数估计\n", " :param visible_list: [[],[],...,[]]\n", " :param hidden_list: [[],[],...,[]]\n", " :return:\n", " \"\"\"\n", " self.pi = np.zeros(shape=(self.hidden_status_num, 1))\n", " self.A = np.zeros(shape=(self.hidden_status_num, self.hidden_status_num))\n", " self.B = np.zeros(shape=(self.hidden_status_num, self.visible_status_num))\n", " for i in range(0, len(visible_list)):\n", " visible_status = visible_list[i]\n", " hidden_status = hidden_list[i]\n", " self.pi[hidden_status[0]] += 1\n", " for j in range(0, len(hidden_status) - 1):\n", " self.A[hidden_status[j], hidden_status[j + 1]] += 1\n", " self.B[hidden_status[j], visible_status[j]] += 1\n", " self.B[hidden_status[j + 1], visible_status[j + 1]] += 1\n", " # 归一化\n", " self.pi = self.pi / np.sum(self.pi)\n", " self.A = self.A / np.sum(self.A, axis=0)\n", " self.B = self.B / np.sum(self.B, axis=0)\n", "\n", " def fit_without_hidden_status(self, visible_list=None, tol=1e-5, n_iter=10):\n", " \"\"\"\n", " 不包含隐状态的参数估计:Baum-Welch算法\n", " :param visible_list: [...]\n", " :param tol:当pi,A,B的增益值变化小于tol时终止\n", " :param n_iter:迭代次数\n", " :return:\n", " \"\"\"\n", " # 初始化参数\n", " if self.pi is None:\n", " self.pi = np.ones(shape=(self.hidden_status_num, 1)) + np.random.random(size=(self.hidden_status_num, 1))\n", " self.pi = self.pi / np.sum(self.pi)\n", " if self.A is None:\n", " self.A = np.ones(shape=(self.hidden_status_num, self.hidden_status_num)) + np.random.random(\n", " size=(self.hidden_status_num, self.hidden_status_num))\n", " self.A = self.A / np.sum(self.A, axis=0)\n", " if self.B is None:\n", " self.B = np.ones(shape=(self.hidden_status_num, self.visible_status_num)) + np.random.random(\n", " size=(self.hidden_status_num, self.visible_status_num))\n", " self.B = self.B / np.sum(self.B, axis=0)\n", " for _ in range(0, n_iter):\n", " # 计算前向概率\n", " alphas = []\n", " alpha = self.pi * self.B[:, [visible_list[0]]]\n", " alphas.append(alpha)\n", " for step in range(1, len(visible_list)):\n", " alpha = self.A.T.dot(alpha) * self.B[:, [visible_list[step]]]\n", " alphas.append(alpha)\n", " # 计算后向概率\n", " betas = []\n", " beta = np.ones_like(self.pi)\n", " betas.append(beta)\n", " for step in range(len(visible_list) - 2, -1, -1):\n", " beta = self.A.dot(self.B[:, [visible_list[step + 1]]] * beta)\n", " betas.append(beta)\n", " betas.reverse()\n", " # 计算gamma值\n", " gammas = []\n", " for i in range(0, len(alphas)):\n", " gammas.append((alphas[i] * betas[i])[:, 0])\n", " gammas = np.asarray(gammas)\n", " # 计算xi值\n", " xi = np.zeros_like(self.A)\n", " for i in range(0, self.hidden_status_num):\n", " for j in range(0, self.hidden_status_num):\n", " xi_i_j = 0.0\n", " for t in range(0, len(visible_list) - 1):\n", " xi_i_j += alphas[t][i][0] * self.A[i, j] * self.B[j, visible_list[t + 1]] * \\\n", " betas[t + 1][j][0]\n", " xi[i, j] = xi_i_j\n", " loss = 0.0 # 统计累计误差\n", " # 更新参数\n", " # 初始概率\n", " for i in range(0, self.hidden_status_num):\n", " new_pi_i = gammas[0][i]\n", " loss += np.abs(new_pi_i - self.pi[i][0])\n", " self.pi[i] = new_pi_i\n", " # 隐状态转移概率\n", " for i in range(0, self.hidden_status_num):\n", " for j in range(0, self.hidden_status_num):\n", " new_a_i_j = xi[i, j] / np.sum(gammas[:, i][:-1])\n", " loss += np.abs(new_a_i_j - self.A[i, j])\n", " self.A[i, j] = new_a_i_j\n", " # 观测概率矩阵\n", " for j in range(0, self.hidden_status_num):\n", " for k in range(0, self.visible_status_num):\n", " new_b_j_k = np.sum(gammas[:, j] * (np.asarray(visible_list) == k)) / np.sum(gammas[:, j])\n", " loss += np.abs(new_b_j_k - self.B[j, k])\n", " self.B[j, k] = new_b_j_k\n", " # 归一化\n", " self.pi = self.pi / np.sum(self.pi)\n", " self.A = self.A / np.sum(self.A, axis=0)\n", " self.B = self.B / np.sum(self.B, axis=0)\n", " if loss < tol:\n", " break" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "#随便造一组数据\n", "O = [\n", " [1, 2, 3, 0, 1, 3, 4],\n", " [1, 2, 3],\n", " [0, 2, 4, 2],\n", " [4, 3, 2, 1],\n", " [3, 1, 1, 1, 1],\n", " [2, 1, 3, 2, 1, 3, 4]\n", "]\n", "I = O" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.16666667]\n", " [0.33333333]\n", " [0.16666667]\n", " [0.16666667]\n", " [0.16666667]]\n", "[[0. 0.125 0.16666667 0. 0. ]\n", " [0. 0.375 0.33333333 0.5 0. ]\n", " [0. 0.375 0. 0.33333333 0.33333333]\n", " [1. 0.125 0.33333333 0. 0.66666667]\n", " [0. 0. 0.16666667 0.16666667 0. ]]\n", "[[1. 0. 0. 0. 0.]\n", " [0. 1. 0. 0. 0.]\n", " [0. 0. 1. 0. 0.]\n", " [0. 0. 0. 1. 0.]\n", " [0. 0. 0. 0. 1.]]\n" ] } ], "source": [ "# 有监督学习\n", "hmm = HMM(hidden_status_num=5, visible_status_num=5)\n", "hmm.fit_with_hidden_status(visible_list=O, hidden_list=I)\n", "print(hmm.pi)\n", "print(hmm.A)\n", "print(hmm.B)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.50829342]\n", " [0.32216634]\n", " [0.00678242]\n", " [0.09996799]\n", " [0.06278983]]\n", "[[0.17883858 0.18538523 0.24457604 0.18546919 0.2047107 ]\n", " [0.18031489 0.22859986 0.24448691 0.16848722 0.17755413]\n", " [0.16895248 0.21542431 0.16450599 0.23683223 0.21459622]\n", " [0.19926632 0.17417557 0.19592024 0.213706 0.21671724]\n", " [0.27262773 0.19641504 0.15051082 0.19550536 0.18642172]]\n", "[[0.14183119 0.28448812 0.11769356 0.18773283 0.18558097]\n", " [0.15028341 0.25009173 0.18529553 0.16628994 0.18546453]\n", " [0.26203855 0.10616088 0.24203416 0.28671766 0.17614 ]\n", " [0.18980262 0.18713587 0.2094639 0.1681282 0.27618863]\n", " [0.25604423 0.1721234 0.24551285 0.19113137 0.17662586]]\n" ] } ], "source": [ "# 无监督学习\n", "hmm = HMM(hidden_status_num=5, visible_status_num=5)\n", "hmm.fit_without_hidden_status(O[0] + O[1] + O[2] + O[3] + O[4] + O[5])\n", "print(hmm.pi)\n", "print(hmm.A)\n", "print(hmm.B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ps：想想无监督学习的意义...，也许可以拿来作半监督学习，先让监督学习学一会儿，然后再无监督学习" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }