白板推导系列Pytorch-隐马尔可夫模型-学习问题

HMM 博客汇总

• 问题三合一（为了防止部分读者厌烦，故而单独发了三个问题的博客）
• 概率计算问题
• 学习问题
• 解码问题
• 词性标注

隐马尔可夫模型的学习问题分为监督学习和非监督学习问题，监督学习采用极大似然估计，非监督学习采用Baum-Welch算法，我们直接先讲Baum-Welch算法。

Baum-Welch算法（EM算法）

事实上，Baum-Welch算法是EM算法在HMM学习问题中的应用。

我们要估计的$\hat{\lambda }$应该为
$$\hat{\lambda }=\underset{\lambda }{\mathrm{argmax}}P(O|\lambda )$$
我们还是先把EM算法的公式写出来先
$${\theta }^{t+1}=\underset{\theta }{\mathrm{argmax}}{\int }_{z}logP(X,Z|\theta )\cdot P(Z|X,\theta )dz$$
其中Z表示隐变量，X表示观测变量，那在HMM中，隐变量是I，观测变量是O，且是离散的
$$\begin{array}{rl}{\lambda }^{t+1}& =\underset{\theta }{\mathrm{argmax}}\sum _{I}logP(O,I|\lambda )\cdot P(I|O,{\lambda }^t)\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{I}logP(O,I|\lambda )\cdot \frac{P(O,I|{\lambda }^t)}{P(O|{\lambda }^t)}\\ & =\underset{\theta }{\mathrm{argmax}}\frac{1}{P(O|{\lambda }^t)}\sum _{I}logP(O,I|\lambda )\cdot P(O,I|{\lambda }^t)\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{I}logP(O,I|\lambda )\cdot P(O,I|{\lambda }^t)\end{array}$$
其中
$$\begin{array}{rl}F(T)=P(O,I\mid \lambda )& =P(o_1,o_2,...,o_T,i_1,i_2,...,i_T|\lambda )\\ & =P(o_T\mid o_1,o_2,...,o_{T-1},i_1,i_2,...,i_T,\lambda )\cdot P(o_1,o_2,...,o_{T-1},i_1,i_2,...,i_T|\lambda )\\ & =P(o_T\mid i_T,\lambda )\cdot P(o_1,o_2,...,o_{T-1},i_1,i_2,...,i_{T-1}|\lambda )\cdot P(i_T|o_1,o_2,...,o_{T-1},i_1,i_2,...i_{T-1},\lambda )\\ & =b_{i_T}(o_T)\cdot F(T-1)\cdot {\alpha }_{i_{T-1}{i}_T}\end{array}$$
故（简单的等比数列）
$$\begin{array}{rl}F(T)& =\prod _{t=2}^T{b}_{i_t}(o_t)\cdot {\alpha }_{i_{t-1}{i}_t}\cdot F(1)\\ & =[\prod _{t=2}^T{b}_{i_t}(o_t)\cdot {\alpha }_{i_{t-1}{i}_t}]\cdot P(o_1,i_1\mid \lambda )\\ & =[\prod _{t=2}^T{b}_{i_t}(o_t)\cdot {\alpha }_{i_{t-1}{i}_t}]\cdot P(o_1\mid i_1,\lambda )\cdot P(i_1\mid \lambda )\\ & =[\prod _{t=2}^T{b}_{i_t}(o_t)\cdot {\alpha }_{i_{t-1}{i}_t}]\cdot b_{i_1}(o_1)\cdot {\pi }_{i_1}\\ & =[\prod _{t=1}^T{b}_{i_t}(o_t)]\cdot [\prod _{t=2}^T{\alpha }_{i_{t-1}{i}_t}]\cdot {\pi }_{i_1}\end{array}$$
现在我们已经得到
$$P(O,I\mid \lambda )=[\prod _{t=1}^T{b}_{i_t}(o_t)]\cdot [\prod _{t=2}^T{\alpha }_{i_{t-1}{i}_t}]\cdot {\pi }_{i_1}$$
然后
$$\begin{array}{rl}{\lambda }^{t+1}& =\underset{\theta }{\mathrm{argmax}}\sum _{I}logP(O,I|\lambda )\cdot P(O,I|{\lambda }^t)\\ & =\underset{\theta }{\mathrm{argmax}}\sum _I(\sum _{t=1}^{T}log(b_{i_t}(o_t))+\sum _{t=2}^{T}log({\alpha }_{i_{t-1}{i}_t})+log({\pi }_{i_1}))\cdot P(O,I|{\lambda }^t)\end{array}$$
记
$$Q(\lambda ,{\lambda }^t)=\underset{\theta }{\mathrm{argmax}}\sum _I(\sum _{t=1}^{T}log(b_{i_t}(o_t))+\sum _{t=2}^{T}log({\alpha }_{i_{t-1}{i}_t})+log({\pi }_{i_1}))\cdot P(O,I|{\lambda }^t)$$
为了求${\pi }^{t+1}$，我们需要求${\pi }_i^{t+1}$

在${\sum }_{i=1}^N{\pi }_i=1$的约束下，

记拉格朗日函数
$$\begin{array}{rl}L(\pi ,\eta )& =\sum _I(log({\pi }_{i_1})\cdot P(O|I,{\lambda }^t))+\eta \cdot (\sum _{i=1}^N{\pi }_i-1)\\ & =\sum _{i_1}...\sum _{i_T}(log({\pi }_{i_1})\cdot P(O|I,{\lambda }^t))+\eta \cdot (\sum _{i=1}^N{\pi }_i-1)\\ & =\sum _{i_1}log({\pi }_{i_1})\cdot P(O|i_1,{\lambda }^t)+\eta \cdot (\sum _{i=1}^N{\pi }_i-1)\\ & =\sum _{i=1}^{N}log({\pi }_i)\cdot P(O|i_1=q_i,{\lambda }^t)+\eta \cdot (\sum _{i=1}^N{\pi }_i-1)\end{array}$$

$$\begin{array}{rl}\frac{\partial L}{\partial {\pi }_i}& =\frac{P(O|i_1=q_i,{\lambda }^t)}{{\pi }_i}+\eta =0\end{array}$$

我们对得到的N个式子求和即可求出$\eta$
$$\begin{array}{r}\sum _{i=1}^N(P(O|i_1=q_i,{\lambda }^t)+\eta \cdot {\pi }_i)\\ =\sum _{i=1}^{N}P(O|i_1=q_i,{\lambda }^t)+\eta =0\end{array}$$
得
$$\eta =-\sum _{i=1}^{N}P(O|i_1=q_i,{\lambda }^t)=-P(O|{\lambda }^t)$$
得
$${\pi }_i^{t+1}=-\frac{P(O|i_1=q_i,{\lambda }^t)}{\eta }=\frac{P(O|i_1=q_i,{\lambda }^t)}{P(O|{\lambda }^t)}$$
为了求$A^{t+1}$，我们即需要求${\alpha }_{ij}^{t+1}$

在状态转移矩阵的每一行和为1的约束下，定义拉格朗日函数
$$\begin{array}{rl}L(\alpha ,\eta )& =\sum _I(\sum _{t=2}^{T}log({\alpha }_{i_{t-1}{i}_t})\cdot P(O|I,{\lambda }^t))+\sum _{i=1}^N{\eta }_i\cdot (\sum _{j=1}^N{\alpha }_{ij}-1)\\ & =\sum _{i=1}^N\sum _{j=1}^N\sum _{t=2}^T\log a_{ij}P\left(O,i_{t-1}=i,i_t=j\mid {\lambda }^t\right)+\sum _{i=1}^N{\eta }_i\cdot (\sum _{j=1}^N{\alpha }_{ij}-1)\end{array}$$
然后对${\alpha }_{ij}$求偏导
$$\frac{\partial L}{\partial {\alpha }_{ij}}=\sum _{t=2}^T\frac{P\left(O,i_{t-1}=i,i_t=j\mid {\lambda }^t\right)}{{\alpha }_{ij}}+{\eta }_i=0$$

$$\sum _{t=2}^{T}P\left(O,i_{t-1}=i,i_t=j\mid {\lambda }^t\right)+{\eta }_i\cdot {\alpha }_{ij}=0$$

我们对j求和
$$\sum _{j=1}^N\sum _{t=2}^{T}P\left(O,i_{t-1}=i,i_t=j\mid {\lambda }^t\right)+{\eta }_i\cdot \sum _{j=1}^N{\alpha }_{ij}=0$$
得到
$${\eta }_i=-\sum _{j=1}^N\sum _{t=2}^{T}P\left(O,i_{t-1}=q_i,i_t=q_j\mid {\lambda }^t\right)=\sum _{t=2}^{T}P(O,i_{t-1}=q_i|{\lambda }^t)$$
故
$$\begin{array}{rl}{\alpha }_{ij}^{t+1}& =-\frac{{\sum }_{t=2}^{T}P\left(O,i_{t-1}=i,i_t=j\mid {\lambda }^t\right)}{{\eta }_i}\\ & =\frac{{\sum }_{t=2}^{T}P\left(O,i_{t-1}=i,i_t=j\mid {\lambda }^t\right)}{{\sum }_{t=2}^{T}P(O,i_{t-1}=q_i|{\lambda }^t)}\end{array}$$
故
$$A^{t+1}=[{\alpha }_{ij}^{t+1}{]}_{N\times N}$$
为了求$B^{t+1}$，我们需要求$b_j(k{)}^{t+1}$

在每行和为1的约束下，得到拉格朗日函数如下
$$\begin{array}{rl}L(b,\eta )& =\sum _I[\sum _{t=1}^{T}log(b_{i_t}(o_t))]\cdot P(O,I|{\lambda }^t)+\sum _{j=1}^N{\eta }_j\cdot [\sum _{k=1}^M{b}_j(k)-1]\\ & =\sum _{j=1}^N\sum _{t=1}^T\log b_j\left(o_t\right)P(O,i_t=q_j\mid {\lambda }^t)+\sum _{j=1}^N{\eta }_j\cdot [\sum _{k=1}^M{b}_j(k)-1]\end{array}$$
然后对$b_j(k)$求偏导
$$\begin{array}{r}\frac{\partial L}{\partial b_j(k)}=\sum _{t=1}^T\frac{P(O,i_t=q_j\mid {\lambda }^t)\cdot I(o_t=v_k)}{b_j(k)}+{\eta }_j=0\end{array}$$

$$\sum _{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\cdot I(o_t=v_k)+{\eta }_j\cdot b_j(k)=0$$

对k求和
$$\sum _{k=1}^M\sum _{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\cdot I(o_t=v_k)+{\eta }_j=0$$
得
$$\begin{array}{rl}{\eta }_j& =-\sum _{k=1}^M\sum _{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\cdot I(o_t=v_k)\\ & =\sum _{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\sum _{k=1}^{M}I(o_t=v_k)\\ & =\sum _{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\end{array}$$
从而得
$$\begin{array}{rl}{b}_j(k{)}^{t+1}& =-\frac{{\sum }_{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\cdot I(o_t=v_k)}{{\eta }_j}\\ & =\frac{{\sum }_{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\cdot I(o_t=v_k)}{{\sum }_{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)}\end{array}$$
综上，得到递推式如下：
$$\begin{array}{rl}& {\pi }_i^{t+1}=\frac{P(O|i_1=q_i,{\lambda }^t)}{P(O|{\lambda }^t)}\\ & \\ & {\alpha }_{ij}^{t+1}=\frac{{\sum }_{t=2}^{T}P\left(O,i_{t-1}=i,i_t=j\mid {\lambda }^t\right)}{{\sum }_{t=2}^{T}P(O,i_{t-1}=q_i|{\lambda }^t)}\\ & \\ & b_j(k{)}^{t+1}=\frac{{\sum }_{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)\cdot I(o_t=v_k)}{{\sum }_{t=1}^{T}P(O,i_t=q_j\mid {\lambda }^t)}\end{array}$$

定义
$$\begin{array}{r}{\xi }_t(i,j)=\frac{p(O,i_t=q_i,i_{t+1}=q_j|{\lambda }^t)}{P(O|{\lambda }^t)}\end{array}$$

$${\gamma }_t(i)=\frac{P(O|i_t=q_i,{\lambda }^t)}{P(O|{\lambda }^t)}$$

于是
$$\begin{array}{rl}& {\pi }_i^{t+1}={\gamma }_1(i)\\ & \\ & {\alpha }_{ij}^{t+1}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_1^{T-1}{\gamma }_t(i)}\\ & \\ & b_j(k{)}^{t+1}=\frac{{\sum }_{t=1}^T{\gamma }_t(j)\cdot I(o_t=v_k)}{{\sum }_{t=1}^T{\gamma }_t(j)}\end{array}$$

算法过程总结

（1）初始化

选取${\alpha }_{ij}^0,b_j(k{)}^0,{\pi }_i^0$

（2）递推
$$\begin{array}{rl}& {\pi }_i^{t+1}={\gamma }_1(i)\\ & \\ & {\alpha }_{ij}^{t+1}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_1^{T-1}{\gamma }_t(i)}\\ & \\ & b_j(k{)}^{t+1}=\frac{{\sum }_{t=1}^T{\gamma }_t(j)\cdot I(o_t=v_k)}{{\sum }_{t=1}^T{\gamma }_t(j)}\end{array}$$
（3）终止

得到模型参数${\lambda }^T=({\pi }_i^T,A^T,B^T)$

补充：忘了讲$\gamma$和$\xi$的求法，读者可尝试自行推导，或者查看《统计学习方法》书中hmm本节内容，然后可与下方_gamma函数和_xi对照

算法实现

我们还是以盒子和球模型为例

生成观测序列

def generate(T):
    boxes = [
        [0,0,0,0,0,1,1,1,1,1],
        [0,0,0,1,1,1,1,1,1,1],
        [0,0,0,0,0,0,1,1,1,1],
        [0,0,0,0,0,0,0,0,1,1]
    ]
    I = []
    O = []
    # 第一步，选择一个盒子
    index = np.random.choice(len(boxes))
    for t in range(T):
        I.append(index)
        # 第二步，抽球
        O.append(np.random.choice(boxes[index]))
        # 第三步，重新选择盒子
        if index==0:
            index = 1
        elif index==1 or index==2:
            index = np.random.choice([index-1,index-1,index+1,index+1,index+1])
        elif index==3:
            index = np.random.choice([index,index-1])
    return I,O

定义模型

class Model:
    def __init__(self,M,N) -> None:
        # 状态集合
        self.N = N
        # 观测集合
        self.M = M

    def _gamma(self,t,i):
        return self.alpha[t][i]*self.beta[t][i]/np.sum(self.alpha[t]*self.beta[t])

    def _xi(self,t,i,j):
        fenmu = 0
        for i1 in range(self.N):
            for j1 in range(self.N):
                fenmu += self.alpha[t][i1]*self.A[i1][j1]*self.B[j1][self.O[t+1]]*self.beta[t+1][j1]
        return (self.alpha[t][i]*self.A[i][j]*self.B[j][self.O[t+1]]*self.beta[t+1][j])/fenmu

    def backward(self):
        # 初值
        self.beta = np.ones(shape=(self.T,self.N))
        # 递推
        for t in reversed(range(1,self.T)):
            for i in range(self.N):
                for j in range(self.N):
                    self.beta[t-1][i] += self.beta[t][j]*self.B[j][self.O[t]]*self.A[i][j]

    def forward(self):
        self.alpha = np.zeros(shape=(self.T,self.N),dtype=float)
        # 初值
        for i in range(self.N):
            self.alpha[0][i] = self.pi[i]*self.B[i][self.O[0]]
        # 递推
        for t in range(self.T-1):
            for i in range(self.N):
                for j in range(self.N):
                    self.alpha[t+1][i] += self.alpha[t][j]*self.A[j][i]
                self.alpha[t+1][i] = self.alpha[t+1][i]*self.B[i][self.O[t]]

    def train(self,O,max_iter=100,toler=1e-5):
        self.O = O
        self.T = len(O)
        # 初始化

        # pi是一个N维的向量
        self.pi = np.ones(shape=(self.N,))/self.N
        # A是一个NxN的矩阵
        self.A = np.ones(shape=(self.N,self.N))/self.N
        # B是一个NxM的矩阵
        self.B = np.ones(shape=(self.N,self.M))/self.M

        # pi是一个N维的向量
        pi = np.ones(shape=(self.N,))/self.N
        # A是一个NxN的矩阵
        A = np.ones(shape=(self.N,self.N))/self.N
        # B是一个NxM的矩阵
        B = np.ones(shape=(self.N,self.M))/self.M


        ## 递推
        for epoch in range(max_iter):
            # 计算前向概率和后向概率
            self.forward()
            self.backward()

            # 计算gamma
            for i in range(self.N):
                pi[i] = self._gamma(1,i)
                for j in range(self.N):
                    fenzi = 0
                    fenmu = 0
                    for t2 in range(self.T-1):
                        fenzi += self._xi(t2,i,j)
                        fenmu += self._gamma(t2,j)
                    A[i][j] = fenzi/fenmu
            for j in range(self.N):
                for k in range(self.M):
                    fenzi = 0
                    fenmu = 0
                    for t2 in range(self.T):
                        fenzi += self._gamma(t2,j)*(self.O[t2]==k)
                        fenmu += self._gamma(t2,j)
                    B[j][k] = fenzi/fenmu
            if np.max(abs(self.pi - pi)) < toler and \
                            np.max(abs(self.A - A)) < toler and \
                            np.max(abs(self.B - B)) < toler:
                self.pi = pi
                self.A = A
                self.B = B
                return pi,A,B
            self.pi = pi.copy()
            self.A = A.copy()
            self.B = B.copy()

        return self.pi,self.A,self.B

    def predict(self,O):
        self.O = O
        self.T = len(O)
        self.forward()
        return np.sum(self.alpha[self.T-1])

训练

model = Model(2,4)
model.train(O,toler=1e-10)

输出如下

(array([0.25, 0.25, 0.25, 0.25]),
 array([[0.25, 0.25, 0.25, 0.25],
        [0.25, 0.25, 0.25, 0.25],
        [0.25, 0.25, 0.25, 0.25],
        [0.25, 0.25, 0.25, 0.25]]),
 array([[0.625, 0.375],
        [0.625, 0.375],
        [0.625, 0.375],
        [0.625, 0.375]]))

完全不对劲，但是可以理解。

监督式学习算法-极大似然估计

推导过Baum-Welch算法之后推导极大似然估计就轻松很多了

假设我们的样本如下
{ ( O 1 , I 1 ) , ( O 2 , I 2 ) , ⋯   , ( O S , I S ) } \left\{\left(O_{1}, I_{1}\right),\left(O_{2}, I_{2}\right), \cdots,\left(O_{S}, I_{S}\right)\right\} {(O1​,I1​),(O2​,I2​),⋯,(OS​,IS​)}
根据极大似然估计则有
$$\begin{array}{rl}\hat{\lambda }& =\underset{\lambda }{\mathrm{argmax}}P(O,I|\lambda )\\ & =\underset{\lambda }{\mathrm{argmax}}\prod _{j=1}^{S}P(O_j,I_j|\lambda )\\ & =\underset{\lambda }{\mathrm{argmax}}\sum _{j=1}^{S}logP(O_j,I_j|\lambda )\\ & =\underset{\lambda }{\mathrm{argmax}}\sum _{j=1}^{S}log\left(\prod _{t=1}^T{b}_{i_t}^j(o_t)\cdot \prod _{t=2}^T{\alpha }_{i_{t-1}{t}_t}^j\cdot {\pi }_{i_1}^j\right)\\ & =\underset{\lambda }{\mathrm{argmax}}\sum _{j=1}^S\left[\sum _{t=1}^{T}log(b_{i_t}^j(o_t))+\sum _{t=2}^{T}log({\alpha }_{i_{t-1}{t}_t}^j)+log({\pi }_{i_1}^j)\right]\end{array}$$
我们先估计$\pi$（i1省略不写），与$\pi$相关的只有最后一项
$$\hat{\pi }=\underset{\pi }{\mathrm{argmax}}\sum _{j=1}^S\log {\pi }^j$$
对$\pi$有约束
$$\sum _{i=1}^N{\pi }_i=1$$
记拉格朗日函数
$$L(\pi ,\eta )=\sum _{j=1}^S\log {\pi }^j+\eta \cdot \left(\sum _{i=1}^N{\pi }_i-1\right)$$
对${\pi }_i$求偏导
$$\frac{\partial L}{\partial {\pi }_i}=\frac{{\sum }_{j=1}^{S}I\left(I_{1j}=i\right)}{{\pi }_i}+\eta =0$$
即
$$\sum _{j=1}^{S}I\left(I_{1j}=i\right)+\eta \cdot {\pi }_i=0$$
对i求和
$$\sum _{i=1}^N\sum _{j=1}^{S}I(I_{1j}=i)+\eta \cdot \sum _{i=1}^N{\pi }_i=0$$
即
$$\begin{array}{rl}& \sum _{j=1}^S\sum _{i=1}^{N}I(I_{1j}=i)+\eta =0\\ & \sum _{j=1}^S\sum _{i=1}^{N}I(I_{1j}=i)+\eta =0\\ & \sum _{j=1}^{S}1+\eta =0\\ & \eta =-S\end{array}$$
所以
$${\hat{\pi }}_i=-\frac{{\sum }_{j=1}^{S}I\left(I_{1j}=i\right)}{\eta }=\frac{{\sum }_{j=1}^{S}I\left(I_{1j}=i\right)}{S}$$
剩下A，B的推导过程，请允许我贴两张图，打公式打烦了

图里面写的有一些问题，希望读者阅读过程中可以自行发现，如果阅读的时候没发现，你看到下面也会发现的

于是我们得到最终估计为
$$\begin{array}{rl}& {\hat{\pi }}_i=\frac{{\sum }_{j=1}^{S}I\left(i_1^j=q_i\right)}{S}\\ & {\hat{\alpha }}_{ij}=\frac{{\sum }_{k=1}^S{\sum }_{t=1}^{T-1}I\left(i_t^k=q_i,i_{t+1}^k=q_j\right)}{{\sum }_{k=1}^S{\sum }_{t=1}^{T-1}I\left(i_t^k=q_i\right)}\\ & {\hat{b}}_j(k)=\frac{{\sum }_{i=1}^S{\sum }_{t=1}^{T}I\left(i_t^i=q_j,o_t^i=v_k\right)}{{\sum }_{i=1}^S{\sum }_{t=1}^{T}I\left(i_t^i=q_j\right)}\end{array}$$

算法实现

生成数据

data = []
for i in range(100):
    I,O = generate(100)
    data.append([I,O])

定义模型

class SupervisedModel:
    def __init__(self,M,N) -> None:
        self.M = M
        self.N = N

    def train(self,data):
        # data:Sx2xT
        # [[[i1,i2,...],[o1,o2,...]],[],[]]
        self.pi = np.zeros(shape=(self.N,))
        self.A = np.zeros(shape=(self.N,self.N))
        self.B = np.zeros(shape=(self.N,self.M))

        S = len(data)
        self.T = len(data[0][0])
        # 求 pi
        for i in range(self.N):
            for j in range(S):
                self.pi[i] += data[j][0][0]==i
            self.pi[i] = self.pi[i]/S
        # 求 A
        for i in range(self.N):
            for j in range(self.N):
                fenzi = 0
                fenmu = 0
                for k in range(S):
                    for t in range(self.T-1):
                        fenzi += data[k][0][t]==i and data[k][0][t+1]==j
                        fenmu += data[k][0][t]==i
                self.A[i][j] = fenzi/fenmu
        
        # 求 B
        for j in range(self.N):
            for k in range(self.M):
                fenzi = 0
                fenmu = 0
                for i in range(S):
                    for t in range(self.T):
                        fenzi += data[i][0][t]==j and data[i][1][t]==k
                        fenmu += data[i][0][t]==j
                self.B[j][k] = fenzi/fenmu
        return self.pi,self.A,self.B

    def predict(self,O):
        # 初值
        beta = np.ones(shape=(self.T,))
        # 递推
        for o in reversed(O):
            temp = np.zeros_like(beta)
            for i in range(self.T):
                for j in range(self.T):
                    temp[i] += beta[j]*self.B[j][o]*self.A[i][j]
            beta = temp
        # 终止
        res = 0
        for i in range(self.T):
            res += self.B[i][O[0]]*beta[i]*self.pi[i]
        return res

训练模型

model = SupervisedModel(2,4)
model.train(data)

(array([0.19, 0.25, 0.2 , 0.36]),
 array([[0.        , 1.        , 0.        , 0.        ],
        [0.38930233, 0.        , 0.61069767, 0.        ],
        [0.        , 0.4073957 , 0.        , 0.5926043 ],
        [0.        , 0.        , 0.49933066, 0.50066934]]),
 array([[0.47313084, 0.52686916],
        [0.30207852, 0.69792148],
        [0.60162602, 0.39837398],
        [0.80851627, 0.19148373]]))

可以看到监督学习的效果非常好