HMM-4-学习模型参数(Baum-Welch)

4 从观测序列中学习模型参数

4.1 最大期望对数似然

设有\(K\)条观测序列\( ^{[5]} \):

\[ X=\{O^k\}_{k=1}^K \]

模型参数\( ^{[5]} \):

\[ \lambda = (A,B,\Pi) \]

调整模型参数使得观测序列概率最大化,即可完成参数学习过程\( ^{[1]} \).

\[ \begin{align}
P(X|\lambda) &= \prod_{k=1}^K P(O^k | \lambda) \\
&= \prod_{k=1}^K \sum_Q P(Q, O^k | \lambda) \\
\log P(X|\lambda) &= \sum_{k=1}^K \log P(O^k | \lambda) \\
&= \sum_{k=1}^K \log \sum_Q P(Q, O^k | \lambda)
\end{align} \]

状态序列\( Q \)未知, 如遍历所有状态序列求和, 计算代价太高且太复杂. 改为计算对数期望. 并在后续简化计算.

Baum-Welch 算法:

\[ \begin{align}
Q(\lambda, \bar{\lambda}) &= \sum_{k=1}^K \sum_Q P(Q|O^k, \lambda) \log P(Q, O^k |\bar{\lambda}) \\
\bar{\lambda} &= \arg \max_{\bar{\lambda}} Q(\lambda, \bar{\lambda})
\end{align} \]

\( Q \)状态概率未知(此时还不存在最佳模型参数), 使用交替迭代的方式计算,第一步随机生成模型参数\( \lambda \),使用此参数在观测序列上计算得到期望值,再调整参数\( \bar{\lambda} \)使得\( Q(\lambda, \bar{\lambda}) \)最大化, 反复迭代. 最终收敛到极值.

4.1.1 对数概率期望

\[ \begin{align}
\bar{\lambda} &= \arg \max_{\bar{\lambda}} Q(\lambda, \bar{\lambda}) \\
& \Rightarrow \sum_{k=1}^K \sum_Q P(Q|O^k, \lambda) \log P(Q, O^k |\bar{\lambda}) \\
& \Rightarrow \sum_{k=1}^K \sum_Q P(Q|O^k, \lambda) \log \left[ \bar{\pi}_{q_1} \bar{b}_{q_1}(O_1^k) \bar{a}_{q_1 q_2} \bar{b}_{q_2}(O_2^k) \cdots \bar{a}_{q_{t-1} q_T} \bar{b}_{q_T}(O_T^k) \right] \\
& \Rightarrow \left[ \sum_{k=1}^K \sum_Q P(Q|O^k, \lambda) \log \bar{\pi}_{q_1} \right] + \left[ \sum_{k=1}^K \sum_Q \sum_{t=1}^{T_{O^k}-1} P(Q|O^k, \lambda) \log \bar{a}_{q_t q_{t+1}} \right] + \left[ \sum_{k=1}^K \sum_Q \sum_{t=1}^{T_{O^k}} P(Q|O^k, \lambda) \log \bar{b}_{q_t}(O_t^k) \right]
\end{align} \]

于是,在这里\( \bar{A}, \bar{B}, \bar{\Pi} \)参数被分为三组互相无关的求和,可分别计算最大化\( ^{[3][4]} \).

4.1.2 最大化

\[ \begin{align}
\bar{\Pi} &= \arg \max_{\bar{\Pi}} \sum_{k=1}^K \sum_Q P(Q|O^k, \lambda) \log \bar{\pi}_{q_1} \\
& \Rightarrow \sum_{k=1}^K \sum_{q_1} \sum_{q_2 q_3 \cdots q_{T_{O^k}}} P(q_1 q_2 \cdots q_{T_{O^k}} |O^k, \lambda) \log \bar{\pi}_{q_1} \\
& \Rightarrow \sum_{k=1}^K \sum_{i=1}^N P(q_1 = S_i |O^k, \lambda) \log \bar{\pi}_i \\
& s.t. \sum_{i=1}^N \bar{\pi}_i = 1 \text{ Lagrange Multiplier } \\
& \Rightarrow \max \sum_{k=1}^K \sum_{i=1}^N P(q_1 = S_i |O^k, \lambda) \log \bar{\pi}_i + \mu \left[ \sum_{i=1}^N \bar{\pi}_i – 1 \right] \\
& \Rightarrow \frac{\partial}{\partial \bar{\pi}_i} \sum_{k=1}^K \sum_{i=1}^N P(q_1 = S_i |O^k, \lambda) \log \bar{\pi}_i + \mu \left[ \sum_{i=1}^N \bar{\pi}_i – 1 \right] = 0 \\
& \Rightarrow \frac{\sum_{k=1}^K P(q_1 = S_i |O^k, \lambda)}{\bar{\pi}_i} + \mu = 0 \\
& \Rightarrow \sum_{k=1}^K P(q_1 = S_i |O^k, \lambda) + \mu \bar{\pi}_i = 0 \\
& \Rightarrow \sum_{i=1}^N \left[ \sum_{k=1}^K P(q_1 = S_i |O^k, \lambda) + \mu \bar{\pi}_i \right] = 0 \\
\text{note:} \\
& \frac{\partial}{\partial \mu} \sum_{k=1}^K \sum_{i=1}^N P(q_1 = S_i |O^k, \lambda) \log \bar{\pi}_i + \mu \left[ \sum_{i=1}^N \bar{\pi}_i – 1 \right] = \sum_{i=1}^N \bar{\pi}_i – 1 = 0 \\
\text{so:} \\
& \Rightarrow \sum_{k=1}^K \sum_{i=1}^N P(q_1 = S_i | O^k, \lambda) + \mu \sum_{i=1}^N \bar{\pi}_i = 0 \\
& \Rightarrow K + \mu = 0 \\
& \Rightarrow \frac{\sum_{k=1}^K P(q_1 = S_i |O^k, \lambda)}{\bar{\pi}_i} + \mu = K + \mu \\
{\bar{\pi}_i} &= \frac{\sum_{k=1}^K P(q_1 = S_i |O^k, \lambda)}{K} \\
\bar{A} &= \arg \max_{\bar{A}} \sum_{k=1}^K \sum_Q \sum_{t=1}^{T_{O^k}-1} P(Q|O^k, \lambda) \log \bar{a}_{q_t q_{t+1}} \\
& \Rightarrow \sum_{k=1}^K \sum_{q_t q_{t+1}} \sum_{q_1 \cdots q_{t-1}, q_{t+2} \cdots q_{T_{O^k}}} \sum_{t=1}^{T_{O^k}-1} P(q_1 \cdots q_{t-1},q_t,q_{t+1},q_{t+2} \cdots q_{T_{O^k}}|O^k, \lambda) \log \bar{a}_{q_t q_{t+1}} \\
& \Rightarrow \sum_{k=1}^K \sum_{q_t q_{t+1}} \sum_{t=1}^{T_{O^k}-1} P(q_t,q_{t+1}|O^k, \lambda) \log \bar{a}_{q_t q_{t+1}} \\
& \Rightarrow \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda) \log \bar{a}_{ij} \\
& s.t. \sum_{j=1}^N \bar{a}_{ij} = 1 \text{ Lagrange Multiplier } \\
& \Rightarrow \max \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda) \log \bar{a}_{ij} + \sum_{i=1}^N \mu_i \left[ \sum_{j=1}^N \bar{a}_{ij} – 1 \right] \\
& \Rightarrow \frac{\partial}{\partial \bar{a}_{ij}} \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda) \log \bar{a}_{ij} + \sum_{i=1}^N \mu_i \left[ \sum_{j=1}^N \bar{a}_{ij} – 1 \right] = 0 \\
& \Rightarrow \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda)}{\bar{a}_{ij}} + \mu_i = 0 \\
& \Rightarrow \sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda) + \mu_i \bar{a}_{ij} = 0 \\
& \Rightarrow \sum_{j=1}^N \left[ \sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda) + \mu_i \bar{a}_{ij} \right] = 0 \\
\text{note:} \\
& \frac{\partial}{\partial \mu_i } \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda) \log \bar{a}_{ij} + \sum_{i=1}^N \mu_i \left[ \sum_{j=1}^N \bar{a}_{ij} – 1 \right] = \sum_{j=1}^N \bar{a}_{ij} – 1 = 0 \\
\text{so:} \\
& \Rightarrow \sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} \sum_{j=1}^N P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda) + \mu_i \sum_{j=1}^N \bar{a}_{ij} = 0 \\
& \Rightarrow \sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i|O^k, \lambda) + \mu_i = 0 \\
& \Rightarrow \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda)}{\bar{a}_{ij}} + \mu_i = \sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i|O^k, \lambda) + \mu_i \\
\bar{a}_{ij} &= \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda)}{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i|O^k, \lambda)} \\
\bar{B} &= \arg \max_{\bar{B}} \sum_{k=1}^K \sum_Q \sum_{t=1}^{T_{O^k}} P(Q|O^k,\lambda) \log \bar{b}_{q_t}(O_t^k) \\
& \Rightarrow \sum_{k=1}^K \sum_{q_t} \sum_{q_1 \cdots q_{t-1}, q_{t+1} \cdots q_{T_{O^k}}} \sum_{t=1}^{T_{O^k}} P(q_1 \cdots q_{t-1}, q_t, q_{t+1} \cdots q_{T_{O^k}}|O^k,\lambda) \log \bar{b}_{q_t}(O_t^k) \\
& \Rightarrow \sum_{k=1}^K \sum_{i=1}^N \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) \log \bar{b}_i(O_t^k) \\
& \Rightarrow \sum_{k=1}^K \sum_{i=1}^N \sum_{m=1}^M \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) \log \bar{b}_i(m)1(O_t^k=v_m) \\
& s.t. \sum_{m=1}^M \bar{b}_i(m) = 1 \text{ Lagrange Multiplier } \\
& \Rightarrow \max \sum_{k=1}^K \sum_{i=1}^N \sum_{m=1}^M \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) \log \bar{b}_i(m)1(O_t^k=v_m) + \sum_{i=1}^N \mu_i \left[ \sum_{m=1}^M \bar{b}_i(m) – 1 \right] \\
& \Rightarrow \frac{\partial}{\partial \bar{b}_i(m)} \sum_{k=1}^K \sum_{i=1}^N \sum_{m=1}^M \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) \log \bar{b}_i(m)1(O_t^k=v_m) + \sum_{i=1}^N \mu_i \left[ \sum_{m=1}^M \bar{b}_i(m) – 1 \right] = 0 \\
& \Rightarrow \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) 1(O_t^k=v_m)}{\bar{b}_i(m)} + \mu_i = 0 \\
& \Rightarrow \sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) 1(O_t^k=v_m) + \mu_i \bar{b}_i(m) = 0 \\
& \Rightarrow \sum_{m=1}^M \left[ \sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) 1(O_t^k=v_m) + \mu_i \bar{b}_i(m) \right] = 0 \\
\text{note:} \\
& \frac{\partial}{ \partial \mu_i } \sum_{k=1}^K \sum_{i=1}^N \sum_{m=1}^M \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) \log \bar{b}_i(m)1(O_t^k=v_m) + \sum_{i=1}^N \mu_i \left[ \sum_{m=1}^M \bar{b}_i(m) – 1 \right] = \sum_{m=1}^M \bar{b}_i(m) – 1 = 0 \\
\text{so:} \\
& \Rightarrow \sum_{k=1}^K \sum_{t=1}^{T_{O^k}} \sum_{m=1}^M P(q_t=S_i|O^k,\lambda) 1(O_t^k=v_m) + \mu_i \sum_{m=1}^M \bar{b}_i(m) = 0 \\
& \Rightarrow \sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) + \mu_i = 0 \\
& \Rightarrow \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) 1(O_t^k=v_m)}{\bar{b}_i(m)} + \mu_i = \sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) + \mu_i \\
\bar{b}_i(m) &= \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) 1(O_t^k=v_m)}{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda)}
\end{align} \]

于是最终得到了参数更新:

\[\begin{align}
{\bar{\pi}_i} &= \frac{\sum_{k=1}^K P(q_1 = S_i |O^k, \lambda)}{K} \\
\bar{a}_{ij} &= \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i,q_{t+1}=S_j|O^k, \lambda)}{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} P(q_t=S_i|O^k, \lambda)} \\
\bar{b}_i(m) &= \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda) 1(O_t^k=v_m)}{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} P(q_t=S_i|O^k,\lambda)}
\end{align}\]

定义:

\[ \begin{align}
\xi_t(i, j) & \equiv P(q_t=S_i, q_{t+1}=S_j | O, \lambda) \\
&= \frac{P(q_t=S_i, q_{t+1}=S_j, O | \lambda)}{P(O | \lambda)} \\
&= \frac{P(O | q_t=S_i, q_{t+1}=S_j, \lambda)P(q_t=S_i, q_{t+1}=S_j | \lambda)}{P(O | \lambda)} \\
&= \frac{P(O | q_t=S_i, q_{t+1}=S_j, \lambda)P(q_{t+1}=S_j | q_t=S_i, \lambda)P(q_t=S_i | \lambda)}{P(O | \lambda)} \\
&= \frac{P(O_1 \cdots O_t | q_t=S_i, \lambda)P(O_{t+1} | q_{t+1}=S_j, \lambda)P(O_{t+2} \cdots O_T | q_{t+1}=S_j, \lambda)P(q_{t+1}=S_j | q_t=S_i, \lambda)P(q_t=S_i | \lambda)}{P(O | \lambda)} \\
&= \frac{P(O_1 \cdots O_t, q_t=S_i | \lambda)P(O_{t+1} | q_{t+1}=S_j, \lambda)P(O_{t+2} \cdots O_T | q_{t+1}=S_j, \lambda)P(q_{t+1}=S_j | q_t=S_i, \lambda)}{P(O | \lambda)} \\
&= \frac{ \alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j) }{ \sum_k^N \sum_l^N \alpha_t(k) a_{kl} b_l(O_{t+1}) \beta_{t+1}(l) } \\
\gamma_t(i) & \equiv P(q_t=S_i | O, \lambda) = \sum_{j=1}^N \xi_t(i,j)
\end{align} \]

\( \gamma_T(i) \)计算时会超出序列长度,可以更改为:

\[ \gamma_t(i) = \frac{ \alpha_t(i) \beta_t(i) }{ \sum_l^N \alpha_t(l) \beta_t(l) } \]

于是更新右侧可以改写为\( ^{[1]} \):

\[\begin{align}
{\bar{\pi}_i} &= \frac{\sum_{k=1}^K \gamma_1^k(i)}{K} \\
\bar{a}_{ij} &= \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} \xi_t^k(i, j)}{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}-1} \gamma_t^k(i)} \\
\bar{b}_i(m) &= \frac{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} \gamma_t^k(i) 1(O_t^k=v_m)}{\sum_{k=1}^K \sum_{t=1}^{T_{O^k}} \gamma_t^k(i)}
\end{align}\]

4.1.3 停止迭代

当更新前后参数差小于特定值时停止迭代,模型参数学习完成.

\[ || \lambda – \bar{\lambda} ||_1 < \epsilon \text{ or } || \lambda – \bar{\lambda} ||_2 < \epsilon \]

4.1.4 局部最优

Baum-Welch算法得到的参数无法保证是全局最优,其效果严重依赖随机初始化参数值,为尽量提高参数优化效果,工程中,需要使用多组随机参数作为Baum-Welch算法初值参数,计算完毕后再择优选择最终得到的最好的参数值.

4.2 数值稳定性

使用\( \alpha_t(i), \beta_t(i) \)时,当观测序列过长,计算过程会超出机器精度,这里需要使用 \( \hat{\alpha}_t(i), \hat{\beta}_t(i) \)替代\( ^{[2]} \).

\[\begin{align}
\hat{\xi}_t(i, j) &= \frac{1}{C_{t+1}} \hat{\alpha}_t(i) a_{ij} b_j(O_{t+1}) \hat{\beta}_{t+1}(j) \\
\hat{\gamma}_t(i) &= \hat{\alpha}_t(i) \hat{\beta}_t(i)
\end{align}\]

4.3 代码试验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
import numpy as np
from matplotlib import pyplot as plt
 
np.random.seed(42)
 
# normal forward algorithm
class Alpha:
  def __init__(self, O, Lambda):
    self.O = O
    self.A, self.B, self.Pi = Lambda
    self.T = len(O)
    self.N = self.A.shape[0]
    self.M = self.B.shape[1]
    assert( self.A.shape[0] == self.A.shape[1] == self.B.shape[0] == self.Pi.shape[0] )
    assert( min(self.O) >= 0 and max(self.O) < self.B.shape[1] )
    '''
    \alpha_t(i) \equiv P(O_1 \cdots O_t, q_t = S_i | \lambda)
    notice: param t: 0=>t1, 1=>t2, ..., T-1=>T
    '''
    self.alpha = np.zeros(shape=(self.T, self.N))
    for t in range(self.T):
      if t == 0:
        self.alpha[t] = [self.Pi[i] * self.B[i, self.O[t]] for i in range(self.N)]
      else:
        self.alpha[t] = [sum([self.alpha[t-1, h] * self.A[h, i] for h in range(self.N)]) * self.B[i, self.O[t]] for i in range(self.N)]
 
  def __call__(self, t, i):
    return self.alpha[t, i]
 
# scaling forward algorithm
class Alpha_Hat:
  def __init__(self, O, Lambda):
    self.O = O
    self.A, self.B, self.Pi = Lambda
    self.T = len(O)
    self.N = self.A.shape[0]
    self.M = self.B.shape[1]
    assert( self.A.shape[0] == self.A.shape[1] == self.B.shape[0] == self.Pi.shape[0] )
    assert( min(self.O) >= 0 and max(self.O) < self.B.shape[1] )
    '''
    \hat{\alpha}_t(i) & \equiv P(q_t = S_i | O_1 \cdots O_t, \lambda)
    notice: param t: 0=>t1, 1=>t2, ..., T-1=>T
    '''
    self.C = np.zeros(shape=(self.T,))
    self.alpha_hat = np.zeros(shape=(self.T, self.N))
    for t in range(self.T):
      if t==0:
        self.alpha_hat[t] = [self.Pi[i] * self.B[i, self.O[t]] for i in range(self.N)]
      else:
        self.alpha_hat[t] = [sum([self.alpha_hat[t-1, h] * self.A[h, i] for h in range(self.N)]) * self.B[i, self.O[t]] for i in range(self.N)]
      self.C[t] = self.alpha_hat[t].sum()
      self.alpha_hat[t] /= self.C[t]
 
  def __call__(self, t, i):
    return self.alpha_hat[t,i]
 
# normal backword algorithm
class Beta:
  def __init__(self, O, Lambda):
    self.O = O
    self.A, self.B, self.Pi = Lambda
    self.T = len(O)
    self.N = self.A.shape[0]
    self.M = self.B.shape[1]
    assert( self.A.shape[0] == self.A.shape[1] == self.B.shape[0] == self.Pi.shape[0] )
    assert( min(self.O) >= 0 and max(self.O) < self.B.shape[1] )
    '''
    \beta_t(i) \equiv P(O_{t+1} \cdots O_T | q_t = S_i, \lambda)
    notice: param t: T=>T,T-1=>T-1,...,1=>t1,0=>t0
    '''
    self.beta = np.ones(shape=(self.T+1, self.N))
    for t in range(self.T-1, -1, -1):
      if t == 0:
        self.beta[t] = np.array([ self.Pi[j] * self.B[j, self.O[t]] * self.beta[t+1, j] for j in range(self.N)]).sum()
      else:
        self.beta[t] = np.array([ self.A[:,j] * self.B[j, self.O[t]] * self.beta[t+1, j] for j in range(self.N)]).sum(axis=0)
 
  def __call__(self, t, i):
    return self.beta[t, i]
 
# scaling backward algorithm
class Beta_Hat:
  def __init__(self, O, Lambda, C):
    self.O = O
    self.A, self.B, self.Pi = Lambda
    self.T = len(O)
    self.N = self.A.shape[0]
    self.M = self.B.shape[1]
    self.C = C
    assert( len(self.O) == len(self.C) )
    assert( self.A.shape[0] == self.A.shape[1] == self.B.shape[0] == self.Pi.shape[0] )
    assert( min(self.O) >= 0 and max(self.O) < self.B.shape[1] )
    '''
    \hat{\beta}_t(i) = \frac{ \beta_t(i) }{ \prod_{u=t+1}^T C_u } 
    notice: param t: T=>T,T-1=>T-1,...,1=>t1,0=>t0
    '''
    self.beta_hat = np.ones(shape=(self.T+1, self.N))
    for t in range(self.T-1, -1, -1):
      if t == 0:
        self.beta_hat[t] = np.array([ self.Pi[j] * self.B[j, self.O[t]] * self.beta_hat[t+1, j] for j in range(self.N)]).sum()
      else:
        self.beta_hat[t] = np.array([ self.A[:,j] * self.B[j, self.O[t]] * self.beta_hat[t+1, j] for j in range(self.N)]).sum(axis=0)
      self.beta_hat[t] /= self.C[t]
 
  def __call__(self, t, i):
    return self.beta_hat[t,i]
 
def Xi(t, i, j, O, Lambda, alpha, beta):
  # \xi_t(i, j) = \frac{ \alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j) }{ \sum_k^N \sum_l^N \alpha_t(k) a_{kl} b_l(O_{t+1}) \beta_{t+1}(l) }
  A, B, Pi = Lambda
  N, M = B.shape
  numerator = alpha(t, i) * A[i, j] * B[j, O[t+1]] * beta(t+2, j)
  denominator = sum( [alpha(t, i) * A[i, j] * B[j, O[t+1]] * beta(t+2, j) for i in range(N) for j in range(N)] )
  return numerator/denominator
 
def Xi_Hat(t, i, j, O, Lambda, alpha_hat, beta_hat):
  # \hat{\xi}_t(i, j) &= \frac{1}{C_{t+1}} \hat{\alpha}_t(i) a_{ij} b_j(O_{t+1}) \hat{\beta}_{t+1}(j)
  A, B, Pi = Lambda
  N, M = B.shape
  numerator = alpha_hat(t, i) * A[i, j] * B[j, O[t+1]] * beta_hat(t+2, j)
  denominator = alpha_hat.C[t+1]
  return numerator/denominator
 
def Gamma(t, i, O, Lambda, alpha, beta):
  # \gamma_t(i) = \frac{ \alpha_t(i) \beta_t(i) }{ \sum_l^N \alpha_t(l) \beta_t(l) }
  A, B, Pi = Lambda
  N, M = B.shape
  numerator = alpha(t, i) * beta(t+1, i)
  denominator = sum( [alpha(t, i) * beta(t+1, i) for i in range(N)] )
  return numerator / denominator
 
def Gamma_Hat(t, i, O, Lambda, alpha_hat, beta_hat):
  # \hat{\gamma}_t(i) &= \hat{\alpha}_t(i) \hat{\beta}_t(i)
  return alpha_hat(t, i) * beta_hat(t+1, i)
 
def Prepare_X(X, Lambda):
  from collections import namedtuple
  PreX = namedtuple('PreX', 'O T alpha beta alpha_hat beta_hat')
  ret = []
  for k in range(len(X)):
    O = X[k]
    T = len(O)
    alpha = None # Alpha(O, Lambda)
    beta = None # Beta(O, Lambda)
    alpha_hat = Alpha_Hat(O, Lambda)
    beta_hat = Beta_Hat(O, Lambda, alpha_hat.C)
    ret.append(PreX(O=O,T=T,alpha=alpha,beta=beta,alpha_hat=alpha_hat,beta_hat=beta_hat))
  return ret
 
def Baum_Welch(X, M, N, max_iter=30, epsilon=1e-8):
  # initial state probabilities
  Pi = np.random.dirichlet(np.ones(N), size=1).flatten() # \Pi = [\pi_i]
  # state transition probabilities
  A = np.random.dirichlet(np.ones(N), size=N) # A = [a_{ij}]
  # observation emission probabilities
  B = np.random.dirichlet(np.ones(M), size=N) # B = [b_j(m)]
  # the parameters of the model
  Lambda = (A, B, Pi) # \lambda = (A,B,\Pi)
  # total of observations
  K = len(X)
 
  pltx, plty, pltz = [], [], []
 
  for it in range(max_iter):
    pre_X = Prepare_X(X, Lambda)
 
    #PX = np.prod( [sum([ pre_X[k].alpha(pre_X[k].T-1, i) for i in range(N)]) for k in range(K)] )
    logPX = sum( [sum([ np.log( pre_X[k].alpha_hat.C[t] ).sum() for t in range(pre_X[k].T)]) for k in range(K)] )
    #np.testing.assert_allclose(PX, np.exp(logPX))
 
    pltx.append(it)
    plty.append(logPX)
 
    # Pi_bar = np.zeros(shape=Pi.shape)
    Pi_bar_hat = np.zeros(shape=Pi.shape)
    for i in range(N):
      # Pi_bar[i] = sum([Gamma(0, i, pre_X[k].O, Lambda, pre_X[k].alpha, pre_X[k].beta) for k in range(K)]) / K
      Pi_bar_hat[i] = sum([Gamma_Hat(0, i, pre_X[k].O, Lambda, pre_X[k].alpha_hat, pre_X[k].beta_hat) for k in range(K)]) / K
    # np.testing.assert_allclose(Pi_bar, Pi_bar_hat)
    np.testing.assert_allclose(Pi_bar_hat.sum(), 1.0)
 
    # A_bar = np.zeros(shape=A.shape)
    A_bar_hat = np.zeros(shape=A.shape)
    for i in range(N):
      for j in range(N):
        # A_bar[i, j] = sum([Xi(t, i, j, pre_X[k].O, Lambda, pre_X[k].alpha, pre_X[k].beta) for k in range(K) for t in range(pre_X[k].T-1)]) / sum([Gamma(t, i, pre_X[k].O, Lambda, pre_X[k].alpha, pre_X[k].beta) for k in range(K) for t in range(pre_X[k].T-1)])
        A_bar_hat[i, j] = sum([Xi_Hat(t, i, j, pre_X[k].O, Lambda, pre_X[k].alpha_hat, pre_X[k].beta_hat) for k in range(K) for t in range(pre_X[k].T-1)]) / sum([Gamma_Hat(t, i, pre_X[k].O, Lambda, pre_X[k].alpha_hat, pre_X[k].beta_hat) for k in range(K) for t in range(pre_X[k].T-1)])
    # np.testing.assert_allclose(A_bar, A_bar_hat)
    np.testing.assert_allclose(A_bar_hat.sum(axis=1), np.ones(N))
 
    # B_bar = np.zeros(shape=B.shape)
    B_bar_hat = np.zeros(shape=B.shape)
    for i in range(N):
      for m in range(M):
        # B_bar[i, m] = sum([Gamma(t, i, pre_X[k].O, Lambda, pre_X[k].alpha, pre_X[k].beta) * (m == pre_X[k].O[t]) for k in range(K) for t in range(pre_X[k].T)]) / sum([Gamma(t, i, pre_X[k].O, Lambda, pre_X[k].alpha, pre_X[k].beta) for k in range(K) for t in range(pre_X[k].T)])
        B_bar_hat[i, m] = sum([Gamma_Hat(t, i, pre_X[k].O, Lambda, pre_X[k].alpha_hat, pre_X[k].beta_hat) * (m == pre_X[k].O[t]) for k in range(K) for t in range(pre_X[k].T)]) / sum([Gamma_Hat(t, i, pre_X[k].O, Lambda, pre_X[k].alpha_hat, pre_X[k].beta_hat) for k in range(K) for t in range(pre_X[k].T)])
    # np.testing.assert_allclose(B_bar, B_bar_hat)
    np.testing.assert_allclose(B_bar_hat.sum(axis=1), np.ones(N))
 
    diff = np.sqrt( np.square(A_bar_hat - Lambda[0]).sum() + np.square(B_bar_hat - Lambda[1]).sum() + np.square(Pi_bar_hat - Lambda[2]).sum() )
 
    pltz.append(diff)
 
    Lambda = (A_bar_hat, B_bar_hat, Pi_bar_hat)
 
    print(it, logPX, diff)
 
    if diff < epsilon:
      break
 
  return Lambda, (pltx, plty, pltz)
 
# observation sequence
# O = { O_1,O_2,\cdots,O_T } 
# X=\{O^k\}_{k=1}^K
g_X = [[1,2,3,4,5,6],
        [1,2,3],
        [2,3,0,5],
        [5,4,3],
        [5,2]]
# M: number of distinct observation symbols
# N: number of statesin the model
g_M, g_N = 7, 3
 
Lambda_hat, pltxyz = Baum_Welch(g_X, g_M, g_N, 30, 1e-8)
 
plt.subplot(1,2,1)
plt.plot( pltxyz[0], pltxyz[1] )
plt.title('log P(X)')
plt.xlabel('iter')
plt.ylabel('log P(X)')
plt.subplot(1,2,2)
plt.plot( pltxyz[0], pltxyz[2] )
plt.title('diff')
plt.xlabel('iter')
plt.ylabel('diff')
plt.show()

输出:

0 -40.62497815989722 0.721147173024104
1 -30.508399927426787 0.354986113329732
2 -28.247689077612087 0.4567153042700449
3 -25.137580719169478 0.34841771603747757
4 -22.822641883753615 0.20386784210766162
5 -21.789989273643144 0.16405469446060156
6 -21.04255296942742 0.07327578135195437
7 -20.671984408468745 0.008848699208315375
8 -20.62369397471953 0.001030852735906072
9 -20.621511230954965 0.00046458927822352956
10 -20.62076650836539 0.00020925707720285193
11 -20.620433631826145 9.418600277596499e-05
12 -20.620284267984392 4.238121825715551e-05
13 -20.620217161942527 1.906868108141107e-05
14 -20.62018699705531 8.579476653362092e-06
15 -20.620173433971743 3.860169682862392e-06
16 -20.620167334511994 1.736846135801057e-06
17 -20.62016459116371 7.814946589306341e-07
18 -20.6201633571643 3.5164095716767247e-07
19 -20.62016280204743 1.5822691384362145e-07
20 -20.620162552310738 7.119794001774433e-08
21 -20.620162439952967 3.2037566706928046e-08
22 -20.620162389400523 1.4416361746859321e-08
23 -20.620162366655006 6.487167112221767e-09

hmm-baum-welch-result

References:
[1] Rabiner, L. R. 1989. “A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition.” Proceedings of the IEEE 77:257–286.
[2] Christopher M. Bishop “Pattern Recognition and Machine Learning” 2006 Springer 13.2.4 Scaling factors p627-629.
[3] Advanced Stochastic Modeling Course material for STA 531
[4] Stephen Tu, Derivation of Baum-Welch Algorithm for Hidden Markov Models
[5] Ethem Alpaydın,Introduction to Machine Learning Third Edition p418-432

Leave a Reply

Your email address will not be published. Required fields are marked *

Time limit is exhausted. Please reload the CAPTCHA.

Proudly powered by WordPress   Premium Style Theme by www.gopiplus.com