### 使用EM算法解决HMM问题
#### EM算法简介
EM（Expectation-Maximization）算法是一种迭代优化方法，适用于含有隐含变量的概率模型参数的最大似然估计或最大后验概率估计。该算法通过交替执行两个步骤——E步（期望步）和M步（最大化步），逐步逼近最优解[^4]。
#### HMM中的EM算法应用
对于隐马尔可夫模型(Hidden Markov Model)，当已知观察序列而未知隐藏状态序列时，可以通过Baum-Welch算法来利用EM框架进行训练。此过程不需要预先知道确切的状态转移情况；相反，在每次迭代过程中更新这些参数直到收敛为止。具体来说：
- **初始化**：随机设定初始发射矩阵A、转移矩阵B和平滑因子π。
- **E-step (期望步)**：基于当前参数计算每个时间点上各个潜在状态发生的概率分布γ(t,i)=P(qt=i|O,λ) 和双线性函数ξ(t,i,j)= P(qt=i,qt+1=j | O, λ)，其中q表示时刻t下的真实但未见的状态，O代表整个观测序列，λ=(A,B,π) 表示待估参数集。
- **M-step (最大化步)**：重新评估并调整上述三个核心组件(A,B,π),使得下一次循环中得到更高的似然度L(θ|X)[^2]。
以下是Python版本的一个简化版实现案例:
```python
import numpy as np
def forward_backward(obs_seq, A, B, pi):
T = len(obs_seq)
N = A.shape[0]
alpha = np.zeros((T,N))
beta = np.ones((T,N))
# Forward algorithm
for t in range(T):
if t==0:
alpha[t,:] = pi * B[:,obs_seq[t]]
else:
alpha[t,:]=np.dot(alpha[t-1],A)*B[:,obs_seq[t]]
scale_factor=alpha.sum(axis=-1).reshape(-1,1)+1e-8
alpha /=scale_factor
# Backward algorithm
for t in reversed(range(T)):
if t<T-1:
beta[t,:]=(beta[t+1]*B[:,obs_seq[t+1]]@A.T)/scale_factor[t+1]
gamma=np.multiply(alpha,beta)
xi_sum=np.zeros_like(A,dtype=float)
for t in range(T-1):
partial_xi=B[:,obs_seq[t+1]].reshape((-1,1)) @ \
(alpha[t].reshape((1,-1))*A*beta[t+1])
xi_sum+=partial_xi/scale_factor[t+1]
return gamma,xi_sum
def em_hmm_train(obs_sequences,max_iter=50,tol=1e-6):
num_states=len(pi)
obs_symbols=max([max(seq)for seq in obs_sequences])+1
old_log_likelihood=None
for iteration in range(max_iter):
new_A=np.zeros((num_states,num_states))
new_B=np.zeros((num_states,obs_symbols))
log_likelihood=0
for sequence in obs_sequences:
gamma,xi_sum=forward_backward(sequence,A,B,pi)
# Update transition probabilities
start_probabilities=sum(gamma[0])/len(obs_sequences)
end_probabilities=sum(gamma[-1])
trans_counts=np.sum(xi_sum,axis=0)
state_occupancies=np.sum(gamma[:-1],axis=0)+end_probabilities
new_A[:]+=trans_counts/(state_occupancies.reshape((-1,1))+1e-8)
# Update emission probabilities
emit_counts=np.bincount(
[s*num_states+i
for i,symbol in enumerate(sequence)
for s in range(num_states)],
weights=[gamma[i][s]
for i,_ in enumerate(sequence)
for s in range(num_states)]
).reshape((num_states,-1)).sum(axis=0)
symbol_totals=new_B.sum(axis=1).reshape((-1,1))+1e-8
new_B[:,:]+=(emit_counts.reshape((num_states,-1))/symbol_totals)
log_likelihood += sum(np.log(scale_factors.flatten()))
A,new_B=A.copy(),new_B.copy()
change_in_ll=log_likelihood-old_log_likelihood if old_log_likelihood is not None else float('inf')
print(f"Iter {iteration}: Log Likelihood={log_likelihood:.3f}, Change={change_in_ll:.3f}")
if abs(change_in_ll)<tol and iteration>0: break
old_log_likelihood=log_likelihood
return {'A':A,'B':B,'pi':start_probabilities}
```
这段代码展示了如何使用前向-后向算法作为内部工具来进行 Baum–Welch 迭代以改进HMM 参数估计的过程[^3]。
扫一扫
1706
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
