本来是想讲BERT的，但是BERT的重点是部署应用，而且用BERT跑一些NLP领域的很多任务时，一般做法都是BERT后面再串一个概率模型来约束输出，比如串联一个条件随机场CRF模型。而我们还没讲CRF呢，而且要了解CRF需要首先了解隐马尔可夫模型HMM(Hidden Markov Model)，但HMM又牵扯了好几个比较晦涩的算法，比如前向算法、后向算法、维特比算法(viterbi algorothm)、极大似然估计、Baum-Welch算法、EM算法 等，内容有点难也有点多，所以这里我单独针对HMM和CRF写一两个章节，为BERT铺个垫。
一、HMM的原理
隐马尔可夫模型，HMM (Hidden Markov Model)，是一个拟合序列数据的概率模型，也叫统计模型，是用来描述马尔可夫过程的模型。还是老套路，从原理、手动计算、调包实现三个角度来一次性说清楚HMM。
网上搜HMM，大部分文章一上来就整篇整篇的数学推导，随意晦涩大量的符号和脚标，让大部分读者上来就一脸懵。我最讨厌这种卖弄了，其实底层的数学没那么难理解，搞一堆数学符号让人抓狂。本篇通过一个最通俗的例子，来抽取HMM的底层原理。
1、做一个简单的小概率题
这个概率题就是HMM的概型，理解了这道题你就理解了HMM。我拿这道题来对比理解HMM。
这道题的解答案其实不难，就是一系列概率的相乘相加，但我们的目的是从它抽取HMM的原理。所以我们先不急着解题。
这道题最大的难点可能就是有的人不理解矩阵A和矩阵B。这里我先解释一下：
（1）我们先按行看矩阵A：
第一行的0.5表示：我前一次是从盒子1中抽球的，那我这一次还从盒子1中抽球的概率。
第一行的0.4表示：我前一次是从盒子1中抽球的，那我这一次从盒子2中抽球的概率。
第一行的0.1表示：我前一次是从盒子1中抽球的，那我这一次从盒子3中抽球的概率。同理A矩阵的第二行、第三行：
第二行就表示：我前一次是从盒子2中抽球的，那我这一次去盒子1、盒子2、盒子3中抽取的概率。
第三行就表示：我前一次是从盒子3中抽球的，那我这一次去盒子1、盒子2、盒子3中抽取的概率。所以，矩阵A的每行3个数字之和都是1。
而且A矩阵是一个方阵，就是行数=列数=盒子的种类数
（2）矩阵B就更好理解了，矩阵B的第一行表示我已经抽取的是盒子1，那我抽取白球的概率就是0.4，抽取黑球的概率就是0.6了，因为第一个盒子只有4白6黑嘛。同理矩阵B的第二行和第三行。所以矩阵B的每行数字之和也是1。但是B矩阵的行数=盒子的种类数，列数=小球的种类数。
2、从上面的概率题，抽出HMM概型
（1）抽取HMM中的相关概念：
一是，上题中的盒1盒2盒3，对应在HMM中就是隐含状态(states)。对应下图A处。
二是，上题中的最后抽取的序列是“白黑白白黑”，这个序列叫观测序列(observations)，也叫可见状态链，对应下图C处。
三是，由于每次抽取的小球的颜色的概率不仅取决于，小球所在的盒子中的黑白球的比例，还取决于你是从哪个盒子中抽的，所以其实观测序列还对应着一个隐含状态序列，就是下图的B处（说明：B只是我胡乱做的一个隐藏状态序列，真实的B可能不是这个序列）。这个隐藏状态序列叫隐藏状态链。
意思就是我们是可以观测到小球颜色序列，但是观测不到每次小球都是从哪个盒子中取的，也就是观测不到盒子序列。但是我们可以肯定的是：当我们得到“白黑白白黑”这个序列时，一定是存在一个对应的盒子序列的，虽然我们观测不到这个盒子序列，但它是真实存在的。这就好比我们只能看到股票的涨涨跌跌，但我们看不到涨跌背后的原因，但是我们笃定涨跌背后肯定是有原因的。
四是，上题的观测序列长度=隐藏状态序列长度=5。其实现实很多问题中，二者不一定非得相等。
可见：一个隐马尔科夫链，其实是一个双链条！一条C链条，另一条B链条。
而且B链和C链都是一个随机过程，所以隐马尔可夫链也是一个双重随机过程。也就是一个HMM有两个随机序列。
在很多现实问题中，有些现象不符合隐马尔可夫过程，但人们会通过扩展"可观测"和"隐藏"状态的概念来构造一个HMM过程，然后用HMM来拟合实际问题。一般我们的做法是：把C当成"果"，把对应的B当成"因"。先寻找多种"因"，并总结出因与因之间的关系(就是矩阵A)、因与果之间的关系(就是矩阵B),C又是可观测到的，就可以用HMM进行拟合了。
（2）抽取HMM中的两个假设
隐马尔科夫模型有两个基本假设：
一是，有限状态马尔科夫性假设（First-Order Markov Assumption）：
这个假设表明，系统在某个时间点的状态只依赖于前一个时间点的状态，与更早的状态无关。
数学表示为：P(qt+1=sj∣qt=si,qt−1,…,q1)=P(qt+1=sj∣qt=si)
这个假设简化了状态转移的计算，使得计算量大大降低。
二是，观测独立性假设（Observation Independence Assumption）：
这个假设表明，给定当前的隐藏状态，观测值与之前或之后的观测值是相互独立的。
数学表示为：P(Ot=vk∣qt=sj,Ot−1,Ot−2,…,O1)=P(Ot=vk∣qt=sj)
这一假设使得观测序列的概率计算更为简便，只需考虑当前隐藏状态下的观测概率。
如果上面的官话和数学表示看不懂，我给大家翻译成大白话，就是下图：
这两个假设也不难理解，第一个假设就是说隐藏状态只与其前一个隐藏状态相关，与其他无关。第二个假设就是说每个可观测状态都是由它对应的隐藏状态生成的，与其他可观测状态无关。
又是车轱辘话，这不就是矩阵A和矩阵B嘛。矩阵A就是假设1的数学语言，矩阵B就是假设2的数学语言。也所以A矩阵我们叫状态转移概率矩阵。B矩阵叫发射概率矩阵，就是从隐藏状态生成观测结果的过程。矩阵A是因与因之间的关系，矩阵B是每个因与果之间的关系。
（3）抽取HMM中的三元组：λ = (π, A, B)
一个隐马尔科夫模型可以用一个三元组：λ=(π,A,B)来表示。其中：π是初始状态分布向量；A是状态转移概率矩阵；B是观测概率矩阵。
还是车轱辘话，三元组里面的π,A,B，前面解释多次了。但是这里要说明的是：我们把(π,A,B)用λ来表示，λ是参数的意思，也就是说π,A,B都是模型的参数。也就是一旦π,A,B定下来后，那HMM概型就确定了。也就是说一组λ=(π,A,B)代表一个HMM。如果π,A,B中的任何一个或多个改变了，那就变成另外一个HMM了。
3、HMM的三个基本问题
隐马尔科夫模型(HMM)涉及三个基本问题：评估问题（Evaluation Problem）、解码问题（Decoding Problem）和学习问题（Learning Problem）。每个问题都有其特定的数学表达方式。
很多人看到这里都会冒出一个疑问，为什么会有这三个问题？！很多资料也不解释就直接堆数学公式，真是云里雾里。这里我解释一下：
前面写的概率题，也就是小标题1，其实就是一个HMM问题。
前面写的小标题2，其实就是一个完整的HMM概型包括的所有东西：观测序列Observations、掩藏序列States、初始状态概率π、隐藏状态转换矩阵A、发射矩阵B。就这五大要素。
这就类似于，我告诉了你"2+3=5",也类似于我告诉了你HMM概型有五要素：
（1） 你可以求2+3=？这就是评估问题，有的资料上叫计算问题。我们小标题1的概率题其实就是这个问题，就是计算问题，让你算一下得到“白黑白白黑”这个观测序列的概率。就是已知观测序列Observations、初始状态概率π、隐藏状态转换矩阵A、发射矩阵B，计算观测序列的概率。而解决这个问题用的是前向算法。后面我会讲这个算法。
（2） 当你知道了2和5，你可以反推2+？=5，这就是解码问题，有的地方也叫预测问题。就是已知观测序列Observations、初始状态概率π、隐藏状态转换矩阵A、发射矩阵B，我反推一下隐藏序列States。反推最大概率的隐藏序列使用的方法是维特比算法(Viterbi Algorithm)。后面也会讲该算法。
这里面还有一层逻辑，就是我们为什么要反推最大概率的隐藏序列？一是为了预测喽。当你知道了隐藏序列，那是不是就知道了下一个最大概率的隐藏状态-->知道了下一个隐藏状态，是不是就知道了这个隐藏状态对应的最大概率观测值了，是不是就做了预测了？！是不是此时就是一个生成模型了？！二是你就可以进行自动标注了，很多NLP任务就是这个问题，后面有时间我会单独写一个篇章介绍一个小案例。
（3） 当你知道了答案是5，你推测一下？+ ？=5，这就是学习问题。就是我们试着去学一个最优的(π, A, B)，让在这个(π, A, B)下，观测序列Observations的概率最大。这种反推模型参数的任务就是学习问题。但是学习问题又分两种情况：
情况1：已知观测序列Observations和隐藏序列States的情况下，来反推(π, A, B)，此时就是一个有监督的学习任务。那根据参数的定义，进行频数统计来算参数(π, A, B)即可。
情况2：我们仅仅只知道观测序列Observations，不知道隐藏序列States，来反推(π, A, B)，此时是一个无监督的学习任务。如果你的观测序列是离散型数据，那此时你得使用Baum-Welch算法(是EM算法的离散形式)来估计参数。如果你的观测序列是连续型数据，那你就直接使用EM算法来估计模型参数。
但是，情况2的两种算法非常难以理解！！！Baum-Welch算法会牵扯到前向算法和后向算法；EM算法还会牵涉到极大似然估计等更多内容。
二、HMM的相关算法：前向算法、维特比算法¶
与HMM相关的计算就是上面三大问题的计算。本部分给大家展示评估问题和解码问题的手动计算过程，也就是展示前向算法和维特比算法的计算过程。
1、HMM计算问题之：前向算法
前向算法就是计算观测序列出现的概率。
我们的小概率题就是这个问题。小概率题就是已知观测序列Observations、初始状态概率π、隐藏状态转换矩阵A、发射矩阵B，求观测序列的概率。
上图是得到"白黑白白黑"序列的所有路径，这些路径的概率之和就是序列"白黑白白黑"出现的概率。
（1）先手动算一下：
上面的手动计算过程是基于最基本的概率常识来计算的。计算过程和答案肯定都不会错。
我们回看这个计算过程，其实这个计算过程就是一个递推的过程，直到序列结束，递推也结束，最后加和所有状态下的概率即可。
（2）李航课本上的前向算法是这样滴：
公式看着真是有些不爽，其实这个公式对应的就是我上面的手写过程。
（3）下面我用代码把算法写一下：
上述代码的计算结果和手动计算的结果是一致的。
（4）把前向算法的代码给大家
# 初始化观测序列、隐藏状态集合
Observations = ['白球', '黑球', '白球', '白球', '黑球']
Obs_set = ['白球', '黑球']
State_set = ['盒子1', '盒子2', '盒子3']
# 将观测序列、状态集合 转换为索引
obs_seq = torch.tensor([0, 1, 0, 0, 1]) # 白黑白白黑
state_set = torch.tensor([0, 1, 2]) #盒1盒2盒3
# 初始化参数
Start_prob = torch.tensor([0.2, 0.5, 0.3])
Trans_prob = torch.tensor([[0.5, 0.4, 0.1], [0.2, 0.2, 0.6], [0.2, 0.5, 0.3]])
Emission_prob = torch.tensor([[0.4, 0.6],[0.8, 0.2],[0.5, 0.5]])
# 前向算法
def forward(obs_seq, state_set, start_prob, trans_prob, emission_prob):
T = len(obs_seq) #观测序列的长度
N = len(state_set) #状态种类个数
alpha = torch.zeros((T, N))
alpha[0, :] = start_prob * emission_prob[:, obs_seq[0]]
for i in range(1, T):
alpha[i] = torch.mm(alpha[i-1].unsqueeze(0), trans_prob) * emission_prob[:, obs_seq[i]]
forward_sum = alpha[-1].sum()
return alpha, forward_sum
#测试结果
alpha, forward_sum = forward(obs_seq=obs_seq, state_set=state_set,
start_prob=Start_prob, trans_prob=Trans_prob, emission_prob=Emission_prob)
alpha, forward_sum(tensor([[0.0800, 0.4000, 0.1500],
[0.0900, 0.0374, 0.1465],
[0.0327, 0.0934, 0.0377],
[0.0170, 0.0405, 0.0353],
[0.0142, 0.0065, 0.0183]]),
tensor(0.0390))
（5）我再拿上面的代码跑一下李航的《统计学习方法》里面的例子，看和书里的计算结果是否一致：
2、HMM解码问题之：维特比算法(Viterbi Algorithm)
当给定了HMM模型的所有参数和一个观测序列，维特比算法是用来反推那个最大概率的隐藏序列的。在实际很多NLP任务中都会用到这个算法，比如词性标注、命名实体识别、语音识别、机器翻译、拼音转汉字、分词、句法分析等任务。所以这个方法是必须要想得非常清楚的！限于文章篇幅和整体逻辑，这里只讲维特比算法的原理，不讲案例了，后面有机会单独开一个新篇章展示案例。
（1）前向算法是不是已经推算出最大概率的隐藏序列了？
上图是我们前向计算时的手动计算过程，当时已经给大家标出了，每步摸出的小球最大可能性来自哪个盒子。所以我们轻松就得到隐藏序列是:盒2-盒3-盒2-盒2-盒3。
但是其实仔细想一下，这个方法是不严谨的。因为它没有考虑每个时间步上的最大值可能来自的不同前驱状态。它计算的只是每个时间步上所有可能的前驱状态之和对应的最大值。而我们要寻找的是所有时间步的最大联合概率，所以前驱状态是必须要考虑的。
其实对于这个问题，在李航的《统计学习方法》中是给出了两种解法，一种就是我们上面的近似算法，另一种就是维特比算法。所以至此你就知道，为什么会有维特比算法，这个算法是干什么的，这些问题了。
（2）理解维特比算法的关键点
维特比算法还是比较难以理解的，我也是参考了很多资料和例子，最后发现豁然开朗的是下面这句话，建议先反复多理解几遍这句话：
如果最优路径在某时刻t通过节点i，那么这条路径从节点i到终点的部分路径，在节点i到终点的路径中，必须是最优的。这是用动态规划求解概率最大路径的原理。那么，我们根据这个原理，就可以从t=1时刻开始，不断向后递推到下一个状态的路径的最大概率，直到在最后到达最终的最优路径终点，然后依据终点回溯到起始点，这样就能得到最优路径了。
所以为了找到最优路径的各个节点，第一步我们先从前往后递推计算每个状态节点之间的路径的最大概率；第二步从后往前回溯最优路径。这就是维特比算法。
下面先展示一下李航课本里面的维特比算法的数学推导：
（3）手动求解
还是老习惯，先手动求解一番，因为容易理解，不费脑细胞嘛。
第一步：从前往后计算每个状态节点的最大概率，也就是对每个节点进行筛选，去掉不用考虑的路径：
第一步我们对每个step的每个状态求max，就是筛选每个节点的不需要考虑的路径，这样就减少了很多待选路径。
第一步我主要想得到的就是图4，图4的路径就已经非常少了，更加方便我们寻找最优路径。
第一步我们求的局部最优路径(就是上图图4中画的紫色线)其实是毫无意义的，是用不上的。它只是在每个时间步上的最大的隐藏状态，但它并不一定是最优的隐藏序列。因为每个时间步上的最大值可能来自的是不同的前驱状态，我们要寻找的是所有时间步的最大联合概率而不仅仅是每个局部时间步的最大概率。所以简单粗暴的选择每个时间步上的最大值，可能会出现路径不连贯或者路径中有实际中不发生的部分或者可能会出现整体概率下降等情况。从上图也可见；'step3的盒2'转化为'step4的盒3'概率最高，但是step5却是从'step4的盒2'到'step5的盒3'概率最高。所以按照局部最优路径，这个路径是无法实现的。所以我们要开启第二步：根据前驱状态，倒着回溯最优路径。
第二步：从后往前回溯最优路径
为什么要从后往前：如果最优路径在某时刻t通过节点i，那么这条路径从节点i到终点的部分路径，在节点i到终点的路径中，必须是最优的。所以我们要从后往前！
回溯的路径就是上图5中的红色箭头。
我们先从step5看，看step5中谁的概率最大？从上图B处看，B比盒1盒2的概率都大，所以step5的隐状态就是盒3。
然后考虑前驱、开始回溯：'step5中的盒3'是从哪里转化来的呢？看C处，是不是就是从'step4的盒子2'转换而来的。那step4的最优隐状态就是盒2。
继续倒推：'step4中的盒2'从step3的谁转化而来？'step4的盒2'从'step3的盒2'转化而来的概率最大，所以step3的最优隐状态就是盒2。
同理，'step2中的盒3'到'step3的盒2'概率最大，'step1中的盒2'到'step2中的盒3'概率最大。
所以我们最优的因状态序列就是：盒2-->盒3-->盒2-->盒2-->盒3。而且这个路径都是考虑前驱下得来的路径，所以都是连续的，也是最优的。
这就是维特比算法的全部过程。下面展示一下回溯过程：
代码实现维特比算法的核心就是计算上图右边的两个矩阵：
第一个矩阵是记录每个时间步上，每种状态和它之前所有可能状态中的最大路径。虽然这个矩阵意义不大，但它的意义都在计算的过程中，所以我们在计算这个矩阵时，重点要记录第二个矩阵。
第二个矩阵是记录每个时间步上，每种状态的最大路径是由它前面的哪个状态引起的。这个矩阵是我们回溯路径的参考依据！
用一句话概括维特比算法就是：在每一时刻，计算当前时刻落在每种隐状态的最大概率，并记录这个最大概率是从前一时刻哪一个隐状态转移过来的，最后再从结尾回溯最优路径，这个路径就是最优的最大概率路径。
（4）代码实现维特比算法
#为了一致性，这里放一个torch版本的，上面截图是np版本的。
#初始化观测序列、隐藏状态集合
Observations = ['白球', '黑球', '白球', '白球', '黑球']
Obs_set = ['白球', '黑球'] #这个变量可以不要
State_set = ['盒子1', '盒子2', '盒子3']
# 将观测序列、状态集合 转换为索引
obs_seq = torch.tensor([0, 1, 0, 0, 1]) # 白0黑1
state_set = torch.tensor([0, 1, 2]) #0表示盒1， 1表示盒2， 2表示盒3
# 初始化参数
Start_prob = torch.tensor([0.2, 0.5, 0.3])
Trans_prob = torch.tensor([[0.5, 0.4, 0.1], [0.2, 0.2, 0.6], [0.2, 0.5, 0.3]])
Emission_prob = torch.tensor([[0.4, 0.6],[0.8, 0.2],[0.5, 0.5]])
#维特比算法
def viterbi(obs_seq, start_prob, trans_prob, emission_prob):
N = len(State_set) #3种状态
T = len(obs_seq) #5个时间步
delta = torch.zeros((T, N)) #5行3列
psi = torch.zeros((T, N)) #5行3列
delta[0] = start_prob * emission_prob[:,obs_seq[0]]
for i in range(1, T): #i=1,2,3,4 时间步
for j in range(N): #j=0,1,2 隐状态
delta[i,j] = torch.max(delta[i-1] * Trans_prob[:,j]) * emission_prob[j, obs_seq[i]]
psi[i,j] = torch.argmax(delta[i-1] * Trans_prob[:,j])
path = torch.zeros(T, dtype=int)
path[-1] = torch.argmax(delta[-1])
for t in range(T-2, -1, -1): # t=3,2, 1,0
path[t] = psi[t+1, path[t+1]]
return delta, psi, path
#测试算法的计算结果
viterbi(obs_seq=obs_seq, start_prob=Start_prob, trans_prob=Trans_prob, emission_prob=Emission_prob)(tensor([[0.0800, 0.4000, 0.1500],
[0.0480, 0.0160, 0.1200],
[0.0096, 0.0480, 0.0180],
[0.0038, 0.0077, 0.0144],
[0.0017, 0.0014, 0.0023]]),
tensor([[0., 0., 0.],
[1., 1., 1.],
[0., 2., 2.],
[1., 1., 1.],
[2., 2., 1.]]),
tensor([1, 2, 1, 1, 2]))
（5）写在最后
现实世界充满了我们可以看到最终结果的现象，但实际上却无法观察到产生这些结果的潜在因素。典型的例子就是股票，我们只看到了股票价格的涨涨跌跌，我们无法准确的知道这些涨跌背后的原因。然而，虽然我们无法观察到背后的驱动因素，但我们可以将这些驱动因素罗列出来，作为隐状态，就可以使用隐马尔可夫模型，将这些驱动因素和涨跌现象进行建模，将这些现象建模为概率系统，进而进行股票涨跌的预测。
在NLP领域中，使用维特比算法的任务也是非常非常的多，但是过程和细节也同样非常多。所以为了保持文章思路的一致性，关于案例，我后面打算单开一个篇章来写。下面我继续写HMM的其他算法。
三、HMM参数估计(监督学习)————频数统计法¶
监督学习的意思就是我们的观测序列是有对应的状态序列的，或者说是观测序列有对应的标注数据的，标注数据就是隐藏状态。我们的目标是找到一个最优的HMM，也就是找一组最优的参数，让这个HMM可以最好的拟合这两个序列。那拟合完毕后我们就可以用这个模型进行自动标注任务了。
当你有大量的带标注的数据，那你可以使用频数统计来推断HMM的参数。所以：
- 一是我们先创建数据，这个数据要有尽量多的观测序列和对应的标注序列样本。因为尽量多时，频率才无限接近概率嘛。
- 二是根据参数的定义进行频数统计。大数定理告诉我们“频率的极限是概率”，当我们无法统计无数次而得到概率时，我们是可以把足够多的有限次的频率当作概率的。
下面我们开始展示频数统计法的计算过程。还以上面的小球为例来制作数据：
我就根据这组参数来随机生成一个100个样本的(观测值，隐藏值)序列和一个10万条样本的序列，然后再根据这两个序列来统计模型参数，看谁的统计结果最接近真实参数：
# 频数统计法的代码：
list_hz = ['盒子1','盒子1','盒子2','盒子2','盒子2','盒子2','盒子2','盒子3','盒子3','盒子3'] #0.2-0.5-0.3 初始状态概率
hz1_ball = ['白球','白球','白球','白球','黑球','黑球','黑球','黑球','黑球','黑球'] #0.4--0.6 发射矩阵
hz2_ball = ['白球','白球','白球','白球','白球','白球','白球','白球','黑球','黑球'] #0.8--0.2
hz3_ball = ['白球','白球','白球','白球','白球','黑球','黑球','黑球','黑球','黑球'] #0.5--0.5
hz1_hz = ['盒子1','盒子1','盒子1','盒子1','盒子1','盒子2','盒子2','盒子2','盒子2','盒子3'] #0.5-0.4-0.1 转换矩阵
hz2_hz = ['盒子1','盒子1','盒子2','盒子2','盒子3','盒子3','盒子3','盒子3','盒子3','盒子3'] #0.2-0.2-0.6
hz3_hz = ['盒子1','盒子1','盒子2','盒子2','盒子2','盒子2','盒子2','盒子3','盒子3','盒子3'] #0.2-0.5-0.3
#制作数据
def create_data(list_hz, hz1_ball, hz2_ball, hz3_ball, hz1_hz, hz2_hz, hz3_hz):
random.seed(1)
data = []
hezi1 = random.choice(list_hz)
if hezi1 == '盒子1': data.append([hezi1, random.choice(hz1_ball)])
if hezi1 == '盒子2': data.append([hezi1, random.choice(hz2_ball)])
if hezi1 == '盒子3': data.append([hezi1, random.choice(hz3_ball)])
pre_hezi = hezi1
for _ in range(1, 100000): #你想生成多少个数据，你就调整这个参数
if pre_hezi == '盒子1':
box = random.choice(hz1_hz)
if box == '盒子1': data.append([box, random.choice(hz1_ball)])
if box == '盒子2': data.append([box, random.choice(hz2_ball)])
if box == '盒子3': data.append([box, random.choice(hz3_ball)])
pre_hezi = box
continue
if pre_hezi == '盒子2':
box = random.choice(hz2_hz)
if box == '盒子1': data.append([box, random.choice(hz1_ball)])
if box == '盒子2': data.append([box, random.choice(hz2_ball)])
if box == '盒子3': data.append([box, random.choice(hz3_ball)])
pre_hezi = box
continue
if pre_hezi == '盒子3':
box = random.choice(hz3_hz)
if box == '盒子1': data.append([box, random.choice(hz1_ball)])
if box == '盒子2': data.append([box, random.choice(hz2_ball)])
if box == '盒子3': data.append([box, random.choice(hz3_ball)])
pre_hezi = box
continue
return data
#生成数据
data = create_data(list_hz, hz1_ball, hz2_ball, hz3_ball, hz1_hz, hz2_hz, hz3_hz)
#计算π #0.2-0.5-0.3
pai = np.zeros([3])
for i in data:
if i[0]=='盒子1': pai[0]+=1
if i[0]=='盒子2': pai[1]+=1
if i[0]=='盒子3': pai[2]+=1
π = (pai/len(data)).round(3)
π
#计算A 转换矩阵A
trans = np.zeros([3,3])
prebox = data[0][0]
for i in range(1, len(data)):
if prebox == '盒子1' and data[i][0]=='盒子1':
trans[0,0] += 1
prebox = data[i][0]
continue
if prebox == '盒子1' and data[i][0]=='盒子2':
trans[0,1] += 1
prebox = data[i][0]
continue
if prebox == '盒子1' and data[i][0]=='盒子3':
trans[0,2] += 1
prebox = data[i][0]
continue
if prebox == '盒子2' and data[i][0]=='盒子1':
trans[1,0] += 1
prebox = data[i][0]
continue
if prebox == '盒子2' and data[i][0]=='盒子2':
trans[1,1] += 1
prebox = data[i][0]
continue
if prebox == '盒子2' and data[i][0]=='盒子3':
trans[1,2] += 1
prebox = data[i][0]
continue
if prebox == '盒子3' and data[i][0]=='盒子1':
trans[2,0] += 1
prebox = data[i][0]
continue
if prebox == '盒子3' and data[i][0]=='盒子2':
trans[2,1] += 1
prebox = data[i][0]
continue
if prebox == '盒子3' and data[i][0]=='盒子3':
trans[2,2] += 1
prebox = data[i][0]
continue
A = (trans/trans.sum(1).reshape(-1,1)).round(3)
A
#计算B 发射矩阵B
emit = np.zeros([3,2])
for i in data:
if i[0] == '盒子1' and i[1] == '白球':
emit[0,0] +=1
continue
if i[0] == '盒子1' and i[1] == '黑球':
emit[0,1] +=1
continue
if i[0] == '盒子2' and i[1] == '白球':
emit[1,0] +=1
continue
if i[0] == '盒子2' and i[1] == '黑球':
emit[1,1] +=1
continue
if i[0] == '盒子3' and i[1] == '白球':
emit[2,0] +=1
continue
if i[0] == '盒子3' and i[1] == '黑球':
emit[2,1] +=1
continue
B = (emit/emit.sum(1).reshape(-1, 1)).round(3)
Barray([0.285, 0.362, 0.353])
array([[0.501, 0.402, 0.098],
[0.199, 0.202, 0.6 ],
[0.199, 0.495, 0.306]])
array([[0.4 , 0.6 ],
[0.8 , 0.2 ],
[0.501, 0.499]])
这是计算过程的数学公式：
小结：频数统计比较简单易理解，但前提是你的数据量要足够大，你才能比较准确的统计出隐藏在事物背后的规律。如果你的数据量比较小，结果就不会很准确。
四、HMM参数估计(非监督学习)————后向算法、Baum-Welch算法¶
1、后向算法
HMM的非监督学习参数估计，需要用到后向算法。所以这里先提前介绍一下后向算法。
其实后向算法还是非常难理解的。网上很多资料要么是模模糊糊一带而过，要么就是不解释为什么要计算后向概率？后向概率又是个啥？后向算法是干什么用的？计算它的目的是什么？等等一些列的疑问都没有说明，上来就开始堆砌一大片一大片的数学推导，真是相当崩溃..... 我是考证很多博文后，已经把后向算法的计算过程都手动计算出来了，也还是对这个算法的目的迷迷糊糊。直到我写完Viterbi算法，开始探索Baum-Welch和EM时，又遇到后向算法，才算对这个算法有了一个新的认识。
（1）我们先看看后向算法的计算过程：
参考很多资料后，我发现后向算法就是，已知模型参数(π, A, B)和一个观测序列O，然后从后往前、计算每个step之后的序列、出现的原因的概率。
依旧使用前向算法用的小案例，下面我把这个小案例中的所有可能路径再画一遍：
前向算法是根据状态概率、一步步从前往后求概率。但是后向算法是从结果一步步倒退所有可能的因！理解起来有些费劲，所以我把所有的可能路径都给大家标出来了，这样方便倒着推算。
下面开始手动倒推一下：
根据手动推导过程，下面写后向算法的代码实现：
# 初始化观测序列、隐藏状态集合
Observations = ['白球', '黑球', '白球', '白球', '黑球']
Obs_set = ['白球', '黑球']
State_set = ['盒子1', '盒子2', '盒子3']
# 将观测序列、状态集合 转换为索引
obs_seq = torch.tensor([0, 1, 0, 0, 1]) # 白黑白白黑
state_set = torch.tensor([0, 1, 2]) #盒1盒2盒3
# 初始化参数
Start_prob = torch.tensor([0.2, 0.5, 0.3])
Trans_prob = torch.tensor([[0.5, 0.4, 0.1], [0.2, 0.2, 0.6], [0.2, 0.5, 0.3]])
Emission_prob = torch.tensor([[0.4, 0.6],[0.8, 0.2],[0.5, 0.5]])
# 后向算法
def backward(obs_seq, state_set, start_prob, trans_prob, emission_prob):
T = len(obs_seq) #观测序列的长度
N = len(state_set) #状态种类个数
beta = torch.ones((T, N))
beta[-2] = torch.mm(Trans_prob,Emission_prob[:, obs_seq[-1]].unsqueeze(1)).reshape(1,-1)
for i in range(T-3, -1, -1): #i=210
beta[i] = (Trans_prob * Emission_prob[:, obs_seq[i+1]] *beta[i+1]).sum(1)
backward_sum = (Start_prob * Emission_prob[:,obs_seq[0]] * beta[0]).sum()
return beta, backward_sum
# 测试结果
beta, backward_sum = backward(obs_seq=obs_seq, state_set=state_set,
start_prob=Start_prob, trans_prob=Trans_prob, emission_prob=Emission_prob)
beta.round(decimals=3), backward_sum(tensor([[0.0590, 0.0660, 0.0520],
[0.1340, 0.1370, 0.1490],
[0.2520, 0.2190, 0.2740],
[0.4300, 0.4600, 0.3700],
[1.0000, 1.0000, 1.0000]]),
tensor(0.0390))
不放心的话，我用这个代码也跑一下李航的那个例子：
网上这篇博文也针对这个例子进行了计算： 【大道至简】机器学习算法之隐马尔科夫模型(Hidden Markov Model, HMM)详解（2）---计算问题：前向算法和后向算法原理详解公式推导及Python实现_pycharm中hmm算法-CSDN博客 ，感兴趣的可以参考，写得非常不错，我一开始也是看这篇博文才得到的启发。
2、为什么需要后向算法，后向算法是干啥用的？
我们前向算法得到的是alpha矩阵，后向算法得到的是beta矩阵，下面我们看看这两个矩阵的意思：
可见：
alpha的step1(第1列)*beta的第1列就是序列【白黑白白黑】的概率
alpha的step2(第2列)*beta的第2列也是序列【白黑白白黑】的概率
alpha的step3(第3列)*beta的第3列也是序列【白黑白白黑】的概率
alpha的step4(第4列)*beta的第4列也是序列【白黑白白黑】的概率
alpha的step5(第5列)*beta的第5列也是序列【白黑白白黑】的概率
所以，alpha的第1列的第1个数字*beta的第1列的第1个数字就是：step1中状态1(盒子1)的概率，也就是P(S1|O,λ)的概率
也就是，上图的gamma矩阵，就是各个step中所有状态的后验概率！也就是根据前向结果和后向结果，我们可以计算出任意节点上的k种隐藏状态的概率是多少。
这个概率正是Baum-Welch算法中最关键的一环。
3、Baum-Welch算法
虽然很多资料上说Baum-Welch算法是EM算法的离散形式，但是我个人感觉还是差别挺大。首先EM不是像HMM那样的双随机链，而且当前隐藏状态仅和前一个隐藏状态相关这样的假设。所以二者还是有很大的区别的。
参考很多资料后，我大概总结Baum-Welch的思路是这样的：先随机初始化一组参数，然后根据这组参数，就像上图那样计算出各个step的各个状态的后验概率，然后根据这个后验概率推出一个组新的参数，然后用这组新参数再算出各个状态的后验概率-->更新参数-->计算后验概率-->更新参数。。。。如此迭代。在这个迭代过程中P(O)是不断提高的过程，也就是新参数更能式序列O出现的概率最大。
下面还是针对我们的小概率题【白黑白白黑】使用Baum-Welch算法，看看结果：
import numpy as np
# 随机初始化模型参数
np.random.seed(1) # 确保示例的可重复性
A_init = np.random.rand(3, 3)
B_init = np.random.rand(3, 2)
pi_init = np.random.rand(3)
# 归一化
A_init /= A_init.sum(axis=1)[:, np.newaxis] #相当于reshape(-1, 1)
B_init /= B_init.sum(axis=1)[:, np.newaxis]
pi_init /= pi_init.sum()
# 鲍姆-韦尔奇算法（Baum-Welch Algorithm）
def baum_welch(observations, A_init, B_init, pi_init, iterations=15):
N = A_init.shape[0] # 状态数量：盒1盒2盒3
M = B_init.shape[1] # 观察符号数量:白、黑
T = len(observations) # 观察序列长度
for _ in range(iterations):
# E步骤：使用当前参数计算前向概率和后向概率
alpha = np.zeros((N, T))
beta = np.zeros((N, T))
alpha[:, 0] = pi_init * B_init[:, observations[0]]
beta[:, -1] = 1
# 计算alpha
for t in range(1, T):
for j in range(N):
alpha[j, t] = np.sum(alpha[:, t-1] * A_init[:, j]) * B_init[j, observations[t]]
print(sum(alpha[:,-1])) ###每次迭代一轮就查看一下序列出现的概率
# 计算beta
for t in range(T-2, -1, -1):
for i in range(N):
beta[i, t] = np.sum(A_init[i, :] * B_init[:, observations[t+1]] * beta[:, t+1])
# 计算gamma和xi（辅助变量）
gamma = np.zeros((N, T)) #3行5列
xi = np.zeros((N, N, T-1)) # 3个 3行4列
for t in range(T-1): #t=0 1 2, 3
denom = np.sum([alpha[i, t] * beta[i, t] for i in range(N)])
for i in range(N):
gamma[i, t] = (alpha[i, t] * beta[i, t]) / denom
for j in range(N):
xi[i, j, t] = (alpha[i, t] * A_init[i, j] * B_init[j, observations[t+1]] * beta[j, t+1]) / denom
gamma[:, -1] = alpha[:, -1] * beta[:, -1] / np.sum(alpha[:, -1] * beta[:, -1])
# M步骤：更新参数
pi_init = gamma[:, 0]
A_init = np.sum(xi, 2) / np.sum(gamma[:, :-1], axis=1)[:, np.newaxis]
for k in range(M):
mask = observations == k
B_init[:, k] = np.sum(gamma[:, mask], axis=1) / np.sum(gamma, axis=1)
return A_init.round(2), B_init.round(2), pi_init.round(2)
# 应用鲍姆-韦尔奇算法
A_final, B_final, pi_final = baum_welch(observations, A_init, B_init, pi_init)
A_final, B_final, pi_final0.02687165892801032 0.03806650915711237 0.042911050341567454 0.05018091281686047 0.06086689214550695 0.07905763580213726 0.11846425975523714 0.2172489684996949 0.44696281607804095 0.7528413601845925 0.9578196558943375 0.9990559100140414 0.9999995537535143 0.9999999999999005 1.0
(array([[0., 1., 0.],
[0., 0., 1.],
[1., 0., 0.]]),
array([[0., 1.],
[1., 0.],
[1., 0.]]),
array([0., 0., 1.]))
可见经过十几轮迭代，序列出现的概率就逼近1了，而且最终的λ也很奇葩，毫无疑问，这是完全的过拟合了！可见，序列数量太少，这种方法非常不合适。
那我们再用前面监督学习的小例子，生成100个数据的序列，看看能不能估计到更加真实的参数：
list_hz = ['盒子1','盒子1','盒子2','盒子2','盒子2','盒子2','盒子2','盒子3','盒子3','盒子3'] #0.2-0.5-0.3 初始状态概率
hz1_ball = ['白球','白球','白球','白球','黑球','黑球','黑球','黑球','黑球','黑球'] #0.4--0.6 发射矩阵
hz2_ball = ['白球','白球','白球','白球','白球','白球','白球','白球','黑球','黑球'] #0.8--0.2
hz3_ball = ['白球','白球','白球','白球','白球','黑球','黑球','黑球','黑球','黑球'] #0.5--0.5
hz1_hz = ['盒子1','盒子1','盒子1','盒子1','盒子1','盒子2','盒子2','盒子2','盒子2','盒子3'] #0.5-0.4-0.1 转换矩阵
hz2_hz = ['盒子1','盒子1','盒子2','盒子2','盒子3','盒子3','盒子3','盒子3','盒子3','盒子3'] #0.2-0.2-0.6
hz3_hz = ['盒子1','盒子1','盒子2','盒子2','盒子2','盒子2','盒子2','盒子3','盒子3','盒子3'] #0.2-0.5-0.3
def create_data(list_hz, hz1_ball, hz2_ball, hz3_ball, hz1_hz, hz2_hz, hz3_hz):
random.seed(1)
data = []
hezi1 = random.choice(list_hz)
if hezi1 == '盒子1': data.append([hezi1, random.choice(hz1_ball)])
if hezi1 == '盒子2': data.append([hezi1, random.choice(hz2_ball)])
if hezi1 == '盒子3': data.append([hezi1, random.choice(hz3_ball)])
pre_hezi = hezi1
for _ in range(1, 100):
if pre_hezi == '盒子1':
box = random.choice(hz1_hz)
if box == '盒子1': data.append([box, random.choice(hz1_ball)])
if box == '盒子2': data.append([box, random.choice(hz2_ball)])
if box == '盒子3': data.append([box, random.choice(hz3_ball)])
pre_hezi = box
continue
if pre_hezi == '盒子2':
box = random.choice(hz2_hz)
if box == '盒子1': data.append([box, random.choice(hz1_ball)])
if box == '盒子2': data.append([box, random.choice(hz2_ball)])
if box == '盒子3': data.append([box, random.choice(hz3_ball)])
pre_hezi = box
continue
if pre_hezi == '盒子3':
box = random.choice(hz3_hz)
if box == '盒子1': data.append([box, random.choice(hz1_ball)])
if box == '盒子2': data.append([box, random.choice(hz2_ball)])
if box == '盒子3': data.append([box, random.choice(hz3_ball)])
pre_hezi = box
continue
return data
data = create_data(list_hz, hz1_ball, hz2_ball, hz3_ball, hz1_hz, hz2_hz, hz3_hz)
len(data)
data_index = []
for i in data:
if i[1] == '白球': data_index.append(0)
if i[1] == '黑球': data_index.append(1)
observations = np.array(data_index)
observations
import numpy as np
# 有资料说初始化也是影响模型效果的重要因素，那我们就直接从上帝之手那个点开始初始化
pi_init = np.array([0.2, 0.5, 0.3])
A_init = np.array([[0.5, 0.4, 0.1], [0.2, 0.2, 0.6], [0.2, 0.5, 0.3]])
B_init = np.array([[0.4, 0.6],[0.8, 0.2],[0.5, 0.5]])
# 鲍姆-韦尔奇算法（Baum-Welch Algorithm）
def baum_welch(observations, A_init, B_init, pi_init, iterations=1000):
N = A_init.shape[0] # 状态数量：盒1盒2盒3
M = B_init.shape[1] # 观察符号数量:白、黑
T = len(observations) # 观察序列长度
for _ in range(iterations):
# E步骤：使用当前参数计算前向概率和后向概率
alpha = np.zeros((N, T))
beta = np.zeros((N, T))
alpha[:, 0] = pi_init * B_init[:, observations[0]]
beta[:, -1] = 1
# 计算alpha
for t in range(1, T):
for j in range(N):
alpha[j, t] = np.sum(alpha[:, t-1] * A_init[:, j]) * B_init[j, observations[t]]
print(sum(alpha[:,-1]))
# 计算beta
for t in range(T-2, -1, -1):
for i in range(N):
beta[i, t] = np.sum(A_init[i, :] * B_init[:, observations[t+1]] * beta[:, t+1])
# 计算gamma和xi（辅助变量）
gamma = np.zeros((N, T)) #3行5列
xi = np.zeros((N, N, T-1)) # 3个 3行4列
for t in range(T-1): #t=0 1 2, 3
denom = np.sum([alpha[i, t] * beta[i, t] for i in range(N)])
for i in range(N):
gamma[i, t] = (alpha[i, t] * beta[i, t]) / denom
for j in range(N):
xi[i, j, t] = (alpha[i, t] * A_init[i, j] * B_init[j, observations[t+1]] * beta[j, t+1]) / denom
gamma[:, -1] = alpha[:, -1] * beta[:, -1] / np.sum(alpha[:, -1] * beta[:, -1])
# M步骤：更新参数
pi_init = gamma[:, 0]
A_init = np.sum(xi, 2) / np.sum(gamma[:, :-1], axis=1)[:, np.newaxis]
for k in range(M):
mask = observations == k
B_init[:, k] = np.sum(gamma[:, mask], axis=1) / np.sum(gamma, axis=1)
return A_init.round(2), B_init.round(2), pi_init.round(2)
# 应用鲍姆-韦尔奇算法
A_final, B_final, pi_final = baum_welch(observations, A_init, B_init, pi_init)
A_final, B_final, pi_final效果还是差强人意，还没有频数统计更接近真实情况。所以大家使用Baum-Welch算法时，还是得考虑清楚，也许结果相差甚远。
网上但凡是介绍Baum-Welch算法的，都是罗列一大堆一大堆的数学公式，实在是头疼。其实看完上面的代码，我们可以反推Baum-Welch算法的思路：
结合左边的所有路径图，就非常直观了。
五、使用hmmlearn库中的HMM
1、安装hmmlearn库
有一个专门实现隐马尔可夫模型的python库：hmmlearn，这个库提供了多种类型的HMM模型以适应不同的数据类型和应用场景。所以需要提前安装一下这个库：
2、Python中的实现示例：处理离散型观测数据
from hmmlearn import hmm
#定义模型参数
pi_init = np.array([0.2, 0.5, 0.3])
A_init = np.array([[0.5, 0.4, 0.1], [0.2, 0.2, 0.6], [0.2, 0.5, 0.3]])
B_init = np.array([[0.4, 0.6],[0.8, 0.2],[0.5, 0.5]])
#observations = np.array([0,1,0,0,1])[:, np.newaxis] #这个序列数据太少，模型就老报错！
#创建hmm实例，并设置参数
model = hmm.CategoricalHMM(n_components=3, n_iter=100)
model.startprob_prior = pi_init
model.transmat_prior = A_init
model.emissionprob_prior = B_init
#用观测序列训练模型
model.fit(observations)
#通过predict接口，预测最大概率的隐藏序列
hidden_states1 = model.predict(observations)
hidden_states1
#也可以通过decode接口，获取观测序列的logp、以及最大概率的隐藏序列
log_probability, hidden_states2 = model.decode(observations,
lengths = len(observations),
algorithm ='viterbi' )
log_probability
hidden_states2 #hidden_states1和hidden_states2的结果一样！3、Python中的实现示例：处理连续型观测数据
如果每个隐藏状态都对应一个高斯分布（或者多元高斯分布），你可以使用hmm.GaussianHMM(n_components, covariance_type, n_iter, tol, init_params)类。
n_components: 隐状态的数量（即 HMM 中的状态数）。必须为正整数。
covariance_type: 指定协方差矩阵的类型，可以为：
'full': 每个高斯分布有其自己的协方差矩阵。
'tied': 所有高斯分布共享相同的协方差矩阵。
'diag': 每个高斯分布有对角协方差矩阵。
'spherical': 每个高斯分布有相同的方差。
n_iter: 最大迭代次数，用于训练模型。
tol: 停止训练的阈值，决定何时停止优化。
init_params: 指定用于初始化模型的参数（如‘s’表示起始状态概率，‘t’表示转移概率，‘c’表示模型的协方差）。
如果你的隐藏状态比较复杂，单一高斯分布不能很好的拟合你的数据，那你可以使用hmm.GMMHMM(n_components, n_mix, covariance_type, n_iter)类创建一个高斯混合模型实例：
n_mix：指定每个隐状态下的高斯成分数量（即每个隐状态由n_mix个高斯分布组成）。
import numpy as np
from hmmlearn import hmm
import matplotlib.pyplot as plt
# 生成模拟数据
np.random.seed(42) # 为了结果的可重现性
n_samples = 200
mean_sunny = 30 # 假设晴天下的温度均值为30°C，
cov_sunny = 5 # 假设晴天温度波动较大
mean_rainy = 15 #假设雨天的均值为15°C
cov_rainy = 2 # 雨天温度波动较小
# 生成隐状态序列（0: 晴天, 1: 雨天）
states = np.array([0] * 100 + [1] * 100) # 模拟前100个观测为晴天，后100个为雨天
np.random.shuffle(states) # 随机打乱状态顺序
# 生成观测数据
observations = []
for state in states:
if state == 0: # 晴天
temperature = np.random.normal(mean_sunny, cov_sunny)
else: # 雨天
temperature = np.random.normal(mean_rainy, cov_rainy)
observations.append(temperature)
observations = np.array(observations).reshape(-1, 1) # 转换为列向量
# 创建 GaussianHMM，设置隐状态数量为2
model = hmm.GaussianHMM(n_components=2, covariance_type="full", n_iter=100)
# 或者你可以创建高斯混合模型：
#model = hmm.GMMHMM(n_components=2, n_mix=2, covariance_type="full", n_iter=100)
# 训练模型
model.fit(observations)
# 预测隐状态
hidden_states = model.predict(observations)
# 绘制结果
plt.figure(figsize=(12, 6))
plt.plot(observations, color='lightgray', label='Temperature Observations')
plt.scatter(range(n_samples), observations, c=hidden_states, s=30, cmap='viridis', label='Predicted Weather States')
plt.title('Gaussian HMM for Sunny and Rainy Days')
plt.xlabel('Sample Index')
plt.ylabel('Temperature (°C)')
plt.legend()
plt.show()
# 输出模型参数
print("Estimated Means:\n", model.means_)
print("Estimated Covariances:\n", model.covars_)
Estimated Means: [[31.07484444] [15.2120914 ]] Estimated Covariances: [[[25.25691875]] [[ 3.5806111 ]]]
从上图中可以看出温度数据与隐状态（晴天和雨天）之间的关系。黑色代表的隐状态是晴天，黄色代表的隐状态是雨天，我们可以根据模型预测的隐状态来分析天气变化。
完结。
------------------------------------------------------------------------------------------------
参考资料：HMM参数估计(监督学习)————极大似然估计
极大似然估计就要牵扯到一些概率的底层的东西了，所以必须得把底层的东西也说清楚。
（1）什么是极大似然估计法maximum-likelihood-estimation(MLE)？
极大似然估计的底层是大数定理和数学期望：
什么是大数定理？
大数定理的结论是频率的极限是概率。你可以通俗的理解为某件事情发生了无数多次后，人们总结出来的规律。比如抛硬币，我抛了无数多次后统计发现，正面向上的概率和反面向上的概率基本都是0.5，那此后即使我只抛一次，我也可以认为：P(硬币正面向上)=0.5，P(硬币正面向下)=0.5。所以大数定理是概率的理论基础。
然而在现实生活中，一些事情是不会出现无数次让我们观察个够，我们经常遇到的都是非常有限的数据，还让我们去预测事情发生的概率！那此时我们就得反过来想，假设一个事件的发生肯定是符合某种规律(隐藏在现象背后的原因)的，那我们观测到的这有限次现象肯定也是在这个规律下发生的！所以我们是不是可以从现象倒退其背后的原因？！答案是可以的。目前最常见的倒推方法就是最大似然估计法。那如何倒推呢？首先你得知道什么是数学期望：
期望的公式是E(X)=的∑x*p。你可以这样理解：期望就是一个事件可能获得的所有结果和其获得相应结果的概率加权求和，作为该事件的整体期望结果。所以期望是建立的概率的基础上的。比如一个同学抛硬币，如果硬币向上，他可以得到10元，如果硬币向下，他可以得到5元。前提是我们已经知道P(硬币正面向上)=0.5，P(硬币正面向下)=0.5，那这个同学可以得到money的期望就是10*0.5+5*0.5=7.5。这就是期望。其实现实中他要么是得到5元，要么就是得到10元。但是就是有了这有限次的观察，他就有了一个期望的"得到的钱数"。如果这位同学可以抛无数次，那随着他抛的次数越多，他平均每次得到的钱数就越接近7.5元。所以我们是不是可以通过有限次的期望去逼近他无限次的规律？！其实这就已经是极大似然的底层逻辑了。什么是极大似然估计？
似然估计是计算一件只重复有限次的事件的期望。极大似然估计就是求一件事情的期望的最大值，也就是最逼近无限次重复时的这件事情的期望，也就是这件事情背后的最接近真相的规律。所以这个逻辑链条是：大数定理-->概率-->期望-->极大似然估计
所以，极大似然估计法的步骤一般是：
写出似然函数；
对似然函数取对数，并整理；
求导数，令导数为0，得到似然方程；
解似然方程，得到的参数即为所求。
所以，使用极大似然估计法，第一步你必须得写出似然函数。第二步取对数是因为，似然函数一般都是概率连乘的形式，而概率都是0-1之间的数字，又连乘，岂不更小，那我取对数至少就变成了连加的形式，方便计算分析。最后对似然方程求导解参数即可。
（2）用上述案例中的100个样本，使用极大似然估计，来估计HMM的参数
我们以上帝之眼已经生成了两个对应的序列：观测序列和隐状态序列。现在我假装不知道参数，就是不知道(π, A, B)，我用统计方法来计算一下(π', A', B'):
扫一扫
7362
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
