11.1 概率无向图模型

【补充解释】成对马尔可夫性可以概括为:概率无向图模型中,任意两个没有边连接的结点是相互独立的。

【补充解释】局部马尔可夫性可以概括为:概率无向图模型中,任意两个没有边直接相连的结点是相互独立的。

【补充解释】全局马尔可夫性可以概括为:概率无向图模型中,任意两个没有边直接相连的集合是相互独立的。

11.2 条件随机场的定义与形式

【补充解释】条件随机场的定义中,因为具有局部马尔可夫性,所以结点 y v y_v yv​与所有 Y w ( w ∼ v ) Y_w (w \sim v) Yw​(w∼v)以外的结点相互独立,于是有 P ( Y v ∣ X , Y w , w ≠ v ) = P ( Y v ∣ X , Y w , w ∼ v ) P(Y_v|X,Y_w,w \ne v)=P(Y_v|X,Y_w,w \sim v) P(Yv​∣X,Yw​,w​=v)=P(Yv​∣X,Yw​,w∼v)。

【重点】线性链随机场中的最大团是相邻两个结点的集合。每个状态只受对应观测及相邻状态的影响。

【补充解释】线性链条件随机场的参数化形式中: t k t_k tk​为第 k k k个转移特征; s l s_l sl​为第 l l l个状态特征。

【补充解释】 M i ( x ) M_i(x) Mi​(x)表示矩阵, M i ( y i − 1 , y i ∣ x ) M_i(y_{i-1},y_i|x) Mi​(yi−1​,yi​∣x)表示矩阵总的一个元素。

【重点】规范化因子 Z w ( x ) Z_w(x) Zw​(x)是以start为起点,以stop为重点通过状态的所有路径 y 1 y 2 ⋯ y n y_1 y_2 \cdots y_n y1​y2​⋯yn​的非规范化概率 ∏ i = 1 n + 1 M i ( y i − 1 , y i ∣ x ) \prod_{i=1}^{n+1} M_i(y_{i-1},y_i|x) ∏i=1n+1​Mi​(yi−1​,yi​∣x)之和。

【补充解释】例11.2中状态图的理解:【机器学习】【条件随机场CRF-1】CRF的矩阵形式表示的示例讲解 + Python实现 - CV_ML_DP的文章

条件随机场中特征函数如何理解?

【阅读文献】条件随机场中特征函数如何理解? - 隔壁大王的回答 - 知乎

条件随机场矩阵形式的理解

根据矩阵的定义式(11.21),以及矩阵中元素 W i ( y i − 1 , y i ∣ x ) W_i(y_{i-1},y_i|x) Wi​(yi−1​,yi​∣x)的定义式(11.23)和 M i ( y i − 1 , y i ∣ x ) M_i(y_{i-1},y_i|x) Mi​(yi−1​,yi​∣x)的定义式(11.22)。对矩阵中第 m m m行第 n n n列的元素,我们有
M i ( x ) m n = [ M i ( y i − 1 , y i ∣ x ) ] m n = e x p ( ∑ k = 1 K w k f k ( y i − 1 = Y m , y i = Y n , x , i ) ) M_i(x)_{mn} = \begin{bmatrix}M_i(y_{i-1},y_i|x)\end{bmatrix}_{mn} = exp(\sum_{k=1}^K w_k f_k(y_{i-1}=\mathcal{Y}_m,y_i=\mathcal{Y}_n,x,i)) Mi​(x)mn​=[Mi​(yi−1​,yi​∣x)​]mn​=exp(k=1K​wk​fk​(yi−1​=Ym​,yi​=Yn​,x,i))

其中 x x x表示观测序列; Y \mathcal{Y} Y表示标记序列所有可能取值的集合, Y m \mathcal{Y}_m Ym​表示集合中的第 m m m个标记, Y n \mathcal{Y}_n Yn​表示集合中的第 n n n个标记。

矩阵元素的理解

矩阵 M i ( x ) M_i(x) Mi​(x)中第 m m m行第 n n n列的元素,可以理解为:在观测序列 x x x固定的情况下,当第 i − 1 i-1 i−1个标记为 Y m \mathcal{Y}_m Ym​,第 i i i个标记为 Y n \mathcal{Y}_n Yn​时,所有转移与状态特征在位置 i i i的加权和。为方便理解,或也可以不严谨地将其理解为在观测序列 x x x固定的情况下,从标记 y i − 1 = Y m y_{i-1}=\mathcal{Y}_m yi−1​=Ym​转移到标记 y i = Y n y_i=\mathcal{Y}_n yi​=Yn​的概率。

矩阵相乘的理解

我们不妨令矩阵 M i − 1 ( x ) M_{i-1}(x) Mi−1​(x)与矩阵 M i ( x ) M_i(x) Mi​(x)相乘,并将结果称为 M i − 1 , i ( x ) M_{i-1,i}(x) Mi−1,i​(x)
$$
M_{i-1}(x) M_{i}(x) =

\begin{bmatrix}
a_{11} & a_{12} \
a_{21} & a_{22}
\end{bmatrix}

\begin{bmatrix}
b_{11} & b_{12} \
b_{21} & b_{22}
\end{bmatrix}

=

\begin{bmatrix}
a_{11}b_{11} + a_{12}b_{21} & a_{11}b_{12} + a_{12}b_{22} \
a_{21}b_{11} + a_{22}b_{21} & a_{21}b_{12} + a_{22}b_{22}
\end{bmatrix}

= M_{i-1,i}(x)
于 是 有 于是有 于是有
\begin{align*}
M_{i-1,i}(x){mn}
& = \begin{bmatrix}M
{i-1,i}(y_{i-1},y_i|x)\end{bmatrix}{mn} \
& = \sum
{{y|y_{i-2}=\mathcal{Y}m,y_i=\mathcal{Y}n}} \prod{j=i-1}^i M_j(y{j-1},y_j|x)
\end{align*}
$$
为方便理解,也可以不严谨地将其理解为在观测序列 x x x固定的情况下,从标记 y i − 2 = Y m y_{i-2}=\mathcal{Y}_m yi−2​=Ym​转移到标记 y i = Y n y_i=\mathcal{Y}_n yi​=Yn​的所有可能路径的概率和。

式11.24的推导(条件随机场的矩阵形式)

根据条件随机场的简化形式,线性链条件随机场可以表示为,即式(11.15)
P ( y ∣ x ) = 1 Z ( x ) e x p ( ∑ k = 1 K w k f k ( y , x ) ) (11.24.1) P(y|x) = \frac{1}{Z(x)} exp \Big( \sum_{k=1}^K w_k f_k(y,x) \Big) \tag{11.24.1} P(y∣x)=Z(x)1​exp(k=1K​wk​fk​(y,x))(11.24.1)
将式(11.13)代入上式(11.24.1),同时考虑到引进了特殊的起点和终点状态标记,需要求和的位置为 i = 1 , 2 , ⋯   , n + 1 i=1,2,\cdots,n+1 i=1,2,⋯,n+1,于是有
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ P(y|x) & = \fr…
将矩阵随机变量的元素 W i ( y i − 1 , y i ∣ x ) W_i(y_{i-1},y_i|x) Wi​(yi−1​,yi​∣x)的定义式(11.23)和 M i ( y i − 1 , y i ∣ x ) M_i(y_{i-1},y_i|x) Mi​(yi−1​,yi​∣x)的定义式(11.22)代入到上式(11.24.2),于是有
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ P(y|x) & = \fr…
上式即式(11.24)。此时有
Z ( x ) = ∑ y ∏ i = 1 n + 1 M i ( y i − 1 , y i ∣ x ) Z(x) = \sum_y \prod_{i=1}^{n+1} M_i(y_{i-1},y_i|x) Z(x)=yi=1n+1​Mi​(yi−1​,yi​∣x)

增加例11.4

(因为书中例子的特征函数的特征条件,均只受状态序列的影响,而不受观测序列的影响;所以提供既受状态序列影响,又受观测序列的影响的例子如下)

设有一标注问题:输入观测序列为 X = ( X 1 , X 2 , X 3 ) X=(X_1,X_2,X_3) X=(X1​,X2​,X3​), X 1 , X 2 , X 3 X_1,X_2,X_3 X1​,X2​,X3​取值于 X = { 0 , 1 } \mathcal{X} = \{0,1\} X={0,1};输出标记序列为 Y = ( Y 1 , Y 2 , Y 3 ) Y=(Y_1,Y_2,Y_3) Y=(Y1​,Y2​,Y3​), Y 1 , Y 2 , Y 3 Y_1,Y_2,Y_3 Y1​,Y2​,Y3​取值于 Y = { 0 , 1 } \mathcal{Y} = \{0,1\} Y={0,1}。

假设特征 t k t_k tk​, s l s_l sl​和对应的权值 λ k \lambda_k λk​, μ l \mu_l μl​如下(下面仅注明特征取值为 1 1 1的条件,取值为 0 0 0的条件省略):

  • t 1 = t 1 ( y i − 1 = 0 , y i = 1 , x ∈ { ( 0 , 1 , 0 ) , ( 0 , 1 , 1 ) , ( 0 , 0 , 0 ) , ( 0 , 0 , 1 ) } , i ∈ { 1 , 2 } ) t_1 = t_1(y_{i-1}=0,y_i=1,x \in \{(0,1,0),(0,1,1),(0,0,0),(0,0,1)\},i \in \{1,2\}) t1​=t1​(yi−1​=0,yi​=1,x∈{(0,1,0),(0,1,1),(0,0,0),(0,0,1)},i∈{1,2}), λ 1 = 1 \lambda_1=1 λ1​=1;
  • t 2 = t 2 ( y i − 1 = 0 , y i = 0 , x ∈ { ( 1 , 1 , 0 ) , ( 1 , 1 , 1 ) , ( 1 , 0 , 0 ) , ( 1 , 0 , 1 ) } , i ∈ { 1 } ) t_2 = t_2(y_{i-1}=0,y_i=0,x \in \{(1,1,0),(1,1,1),(1,0,0),(1,0,1)\},i \in \{1\}) t2​=t2​(yi−1​=0,yi​=0,x∈{(1,1,0),(1,1,1),(1,0,0),(1,0,1)},i∈{1}), λ 2 = 0.6 \lambda_2=0.6 λ2​=0.6;
  • t 3 = t 3 ( y i − 1 = 1 , y i ∈ { 0 , 1 } , x ∈ { ( 0 , 0 , 0 ) , ( 1 , 1 , 1 ) } , i ∈ { 2 } ) t_3 = t_3(y_{i-1}=1,y_i \in \{0,1\},x \in \{(0,0,0),(1,1,1)\},i \in \{2\}) t3​=t3​(yi−1​=1,yi​∈{0,1},x∈{(0,0,0),(1,1,1)},i∈{2}), λ 3 = 1.2 \lambda_3=1.2 λ3​=1.2;
  • t 4 = t 4 ( y i − 1 = 1 , y i = 1 , x ∈ { ( 0 , 1 , 0 ) , ( 0 , 1 , 1 ) , ( 0 , 0 , 0 ) , ( 0 , 0 , 1 ) , ( 1 , 1 , 0 ) , ( 1 , 1 , 1 ) , ( 1 , 0 , 0 ) , ( 1 , 0 , 1 ) } , i ∈ { 2 } ) t_4 = t_4(y_{i-1}=1,y_i=1,x \in \{(0,1,0),(0,1,1),(0,0,0),(0,0,1),(1,1,0),(1,1,1),(1,0,0),(1,0,1)\},i \in \{2\}) t4​=t4​(yi−1​=1,yi​=1,x∈{(0,1,0),(0,1,1),(0,0,0),(0,0,1),(1,1,0),(1,1,1),(1,0,0),(1,0,1)},i∈{2}), λ 4 = 0.2 \lambda_4=0.2 λ4​=0.2;
  • t 5 = t 5 ( y i − 1 = { 0 , 1 } , y i = 0 , x ∈ { ( 0 , 1 , 0 ) , ( 0 , 1 , 1 ) , ( 1 , 1 , 0 ) , ( 1 , 1 , 1 ) } , i ∈ { 1 , 2 } ) t_5 = t_5(y_{i-1}=\{0,1\},y_i =0,x \in \{(0,1,0),(0,1,1),(1,1,0),(1,1,1)\},i \in \{1,2\}) t5​=t5​(yi−1​={0,1},yi​=0,x∈{(0,1,0),(0,1,1),(1,1,0),(1,1,1)},i∈{1,2}), λ 5 = 1.4 \lambda_5=1.4 λ5​=1.4;
  • s 1 = s 1 ( y i = 0 , x ∈ ( 0 , 1 , 1 ) , ( 1 , 1 , 0 ) , ( 1 , 0 , 1 ) , i ∈ { 0 , 1 , 2 } ) s_1 = s_1(y_i = 0,x\in {(0,1,1),(1,1,0),(1,0,1)},i \in \{0,1,2\}) s1​=s1​(yi​=0,x∈(0,1,1),(1,1,0),(1,0,1),i∈{0,1,2}), μ 1 = 1 \mu_1=1 μ1​=1;
  • s 2 = s 2 ( y i = 1 , x ∈ ( 0 , 1 , 0 ) , ( 0 , 1 , 1 ) , ( 0 , 0 , 0 ) , ( 0 , 0 , 1 ) , ( 1 , 1 , 0 ) , ( 1 , 1 , 1 ) , ( 1 , 0 , 0 ) , ( 1 , 0 , 1 ) , i = 0 ) s_2 = s_2(y_i = 1,x\in {(0,1,0),(0,1,1),(0,0,0),(0,0,1),(1,1,0),(1,1,1),(1,0,0),(1,0,1)},i = 0) s2​=s2​(yi​=1,x∈(0,1,0),(0,1,1),(0,0,0),(0,0,1),(1,1,0),(1,1,1),(1,0,0),(1,0,1),i=0), μ 2 = 0.2 \mu_2=0.2 μ2​=0.2;
  • s 3 = s 3 ( y i = 0 , x ∈ ( 1 , 0 , 0 ) , ( 0 , 1 , 0 ) , ( 0 , 0 , 1 ) , i ∈ { 0 , 1 } ) s_3 = s_3(y_i = 0,x\in {(1,0,0),(0,1,0),(0,0,1)},i \in \{0,1\}) s3​=s3​(yi​=0,x∈(1,0,0),(0,1,0),(0,0,1),i∈{0,1}), μ 3 = 0.8 \mu_3=0.8 μ3​=0.8;
  • s 4 = s 4 ( y i = 1 , x ∈ ( 1 , 0 , 1 ) , ( 0 , 1 , 0 ) , i ∈ { 0 , 2 } ) s_4 = s_4(y_i = 1,x\in {(1,0,1),(0,1,0)},i \in \{0,2\}) s4​=s4​(yi​=1,x∈(1,0,1),(0,1,0),i∈{0,2}), μ 4 = 0.5 \mu_4=0.5 μ4​=0.5。

已知模型计算条件概率(原生Python实现)

def count_conditional_probability(w1, t, w2, s, x, y):
    """已知条件随机场模型计算状态序列关于观测序列的非规范化条件概率

    :param w1: 模型的转移特征权重
    :param t: 模型的转移特征函数
    :param w2: 模型的状态特征权重
    :param s: 模型的状态特征函数
    :param x: 需要计算的观测序列
    :param y: 需要计算的状态序列
    :return: 状态序列关于观测序列的条件概率
    """
    n_features_1 = len(w1)  # 转移特征数
    n_features_2 = len(w2)  # 状态特征数
    n_position = len(x)  # 序列中的位置数

    res = 0
    for k in range(n_features_1):
        for i in range(1, n_position):
            res += w1[k] * t[k](y[i - 1], y[i], x, i)
    for k in range(n_features_2):
        for i in range(n_position):
            res += w2[k] * s[k](y[i], x, i)
    return pow(math.e, res)

【测试】例11.4(下文省略特征函数及参数部分)

if __name__ == "__main__":
    def t1(y0, y1, x, i):
        return int(y0 in {0} and y1 in {1} and x in {(0, 1, 0), (0, 1, 1), (0, 0, 0), (0, 0, 1)} and i in {1, 2})
    def t2(y0, y1, x, i):
        return int(y0 in {0} and y1 in {0} and x in {(1, 1, 0), (1, 1, 1), (1, 0, 0), (1, 0, 1)} and i in {1})
    def t3(y0, y1, x, i):
        return int(y0 in {1} and y1 in {0, 1} and x in {(0, 0, 0), (1, 1, 1)} and i in {2})
    def t4(y0, y1, x, i):
        return int(y0 in {1} and y1 in {1} and x in {(0, 1, 0), (0, 1, 1), (0, 0, 0), (0, 0, 1), (1, 1, 0), (1, 1, 1), (1, 0, 0), (1, 0, 1)} and i in {2})
    def t5(y0, y1, x, i):
        return int(y0 in {0, 1} and y1 in {0} and x in {(0, 1, 0), (0, 1, 1), (1, 1, 0), (1, 1, 1)} and i in {1, 2})
    def s1(y0, x, i):
        return int(y0 in {0} and x in {(0, 1, 1), (1, 1, 0), (1, 0, 1)} and i in {0, 1, 2})
    def s2(y0, x, i):
        return int(y0 in {1} and x in {(0, 1, 0), (0, 1, 1), (0, 0, 0), (0, 0, 1), (1, 1, 0), (1, 1, 1), (1, 0, 0), (1, 0, 1)} and i in {0})
    def s3(y0, x, i):
        return int(y0 in {0} and x in {(0, 1, 1), (1, 1, 0), (1, 0, 1)} and i in {0, 1})
    def s4(y0, x, i):
        return int(y0 in {1} and x in {(1, 0, 1), (0, 1, 0)} and i in {0, 2})

    w1 = [1, 0.6, 1.2, 0.2, 1.4]
    t = [t1, t2, t3, t4, t5]
    w2 = [1, 0.2, 0.8, 0.5]
    s = [s1, s2, s3, s4]

    for x in {(0, 1, 0), (0, 1, 1), (0, 0, 0), (0, 0, 1), (1, 1, 0), (1, 1, 1), (1, 0, 0), (1, 0, 1)}:
        for y in {(0, 1, 0), (0, 1, 1), (0, 0, 0), (0, 0, 1), (1, 1, 0), (1, 1, 1), (1, 0, 0), (1, 0, 1)}:
            print(x, "->", y, ":", count_conditional_probability(w1, t, w2, s, x, y))  

# 输出摘要:
# (1, 0, 1) -> (1, 0, 1) : 20.085536923187664
# (1, 0, 1) -> (1, 1, 0) : 5.473947391727199
# (1, 0, 1) -> (0, 1, 0) : 16.444646771097045
# ......

已知模型构造随机样本集(原生Python实现)

import bisect
import random

def make_hidden_sequence(w1, t, w2, s, x_range, y_range, n_samples=1000, random_state=0):
    """已知模型构造随机样本集

    :param w1: 模型的转移特征权重
    :param t: 模型的转移特征函数
    :param w2: 模型的状态特征权重
    :param s: 模型的状态特征函数
    :param x_range: 观测序列的可能取值
    :param y_range: 状态序列的可能取值
    :param n_samples: 生成样本集样本数(近似)
    :return: 状态序列关于观测序列的条件概率
    """
    P = [[0.0] * len(y_range) for _ in range(len(x_range))]  # 条件概率分布

    lst = []
    sum_ = 0
    for i, x in enumerate(x_range):
        for j, y in enumerate(y_range):
            P[i][j] = round(count_conditional_probability(w1, t, w2, s, x, y), 1)
            sum_ += P[i][j]
            lst.append(sum_)

    X, Y = [], []

    random.seed(random_state)
    for _ in range(n_samples):
        r = random.uniform(0, sum_)
        idx = bisect.bisect_left(lst, r)
        i, j = divmod(idx, len(y_range))
        X.append(x_range[i])
        Y.append(y_range[j])

    return X, Y

【测试】例11.4(省略特征函数及参数部分)

if __name__ == "__main__":
    X, Y = make_hidden_sequence(
        w1, t, w2, s,
        [(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)],
        [(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)])

    for row in zip(X, Y):
        print(row)

# 输出摘要:
# ((1, 1, 0), (0, 0, 0))
# ((1, 1, 0), (0, 0, 0))
# ((1, 0, 1), (0, 0, 0))
# ......

11.3 条件随机场的概率计算问题

【补充解释】 α i ( x ) \alpha_i(x) αi​(x)每次均以行向量的形式,与各个位置的随机矩阵相乘时; β i ( x ) \beta_i(x) βi​(x)每次均以列向量的形式,与各个位置的随机矩阵相乘。

式11.34的推导(条件随机场特征函数关于条件分布的数学期望)

显然地,有特征函数 f k f_k fk​关于条件分布 P ( Y ∣ X ) P(Y|X) P(Y∣X)的数学期望,即式(11.34)第一行,其中 x x x为输入观测序列的取值, y y y为对应状态序列的取值:
E P ( Y ∣ X ) [ f k ] = ∑ y P ( y ∣ x ) f k ( y , x ) (11.34.1) E_{P(Y|X)}[f_k] = \sum_y P(y|x) f_k(y,x) \tag{11.34.1} EP(Y∣X)​[fk​]=y​P(y∣x)fk​(y,x)(11.34.1)

其中 k = 1 , 2 , ⋯   , K k=1,2,\cdots,K k=1,2,⋯,K,下文省略 k k k的取值范围。

根据特征函数 f k ( y , x ) f_k(y,x) fk​(y,x)的定义式(11.13),同时考虑到引进了特殊的起点和终点状态标记,需要求和的位置为 i = 1 , 2 , ⋯   , n + 1 i=1,2,\cdots,n+1 i=1,2,⋯,n+1,于是有
f k ( y , x ) = ∑ i = 1 n + 1 f k ( y i − 1 , y i , x , i ) (11.34.2) f_k(y,x) = \sum_{i=1}^{n+1} f_k(y_{i-1},y_i,x,i) \tag{11.34.2} fk​(y,x)=i=1n+1​fk​(yi−1​,yi​,x,i)(11.34.2)

将上式(11.34.2)代入式(11.34.1)得
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ E_{P(Y|X)}[f_k…
考虑到对于特征函数,只有当满足特征条件时取值为1,否则为0。所以,在上式(11.34.3)中,只需要考虑满足当前特征函数的特征条件的状态序列。

例如,对于 f k ( y i − 1 , y i , x , i ) f_k(y_{i-1},y_i,x,i) fk​(yi−1​,yi​,x,i),只需要考虑状态 y i − 1 y_{i-1} yi−1​和状态 y i y_i yi​满足特征条件的状态序列,记作 { y ′ ∣ y i − 1 ′ = y i − 1 , y i ′ = y i } \{y'|y'_{i-1}=y_{i-1},y'_i=y_i\} {y∣yi−1​=yi−1​,yi​=yi​}。

于是,上式可以写成
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ E_{P(Y|X)}[f_k…
按照书中写法,将上式(11.34.4)的求和符号中的 { y ′ ∣ y i − 1 ′ = y i − 1 , y i ′ = y i } \{y'|y'_{i-1}=y_{i-1},y'_i=y_i\} {y∣yi−1​=yi−1​,yi​=yi​}简写为 y i − 1 y i y_{i-1} y_i yi−1​yi
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ E_{P(Y|X)}[f_k…
因为 P ( y ′ ∣ x ) = P ( Y i − 1 = y i − 1 , Y i = y i ∣ x ) P(y'|x) = P(Y_{i-1}=y_{i-1},Y_i=y_i|x) P(y∣x)=P(Yi−1​=yi−1​,Yi​=yi​∣x),所以将式(11.33)代入上式,得
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ E_{P(Y|X)}[f_k…
其中, Z ( x ) Z(x) Z(x)由式(11.33)定义,有
Z ( x ) = α n T ( x ) 1 = 1 β 1 ( x ) Z(x) = \alpha_n^T(x) 1 = 1 \beta_1(x) Z(x)=αnT​(x)1=1β1​(x)

11.4 条件随机场的学习算法

【补充解释】对数似然函数中的 N N N,是指训练数据集的样本数。在迭代尺度法的条件中,已知训练数据集和经验概率分布 P ^ ( X , Y ) \hat{P}(X,Y) P^(X,Y)。

【补充说明】 E p ( t k ) E_p(t_k) Ep​(tk​)和 E p ( s l ) E_p(s_l) Ep​(sl​)即 E P ( X , Y ) [ f k ] E_{P(X,Y)}[f_k] EP(X,Y)​[fk​](式11.35)。

【勘误补充】P. 228 倒数第9行的 P ~ ( x , y ) \tilde{P}(x,y) P~(x,y)应写作 P ~ ( X , Y ) \tilde{P}(X,Y) P~(X,Y)。

【勘误补充】P. 230 第8行的“期待值”应写作“期望值”。

改进的迭代尺度法求条件随机场推导

(一) 输入

给定训练数据集
T = { ( x ( 1 ) , y ( 1 ) ) , ( x ( 2 ) , y ( 2 ) ) , ⋯   , ( x ( N ) , y ( N ) ) } (11.A.1) T = \{(x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\cdots,(x^{(N)},y^{(N)})\} \tag{11.A.1} T={(x(1),y(1)),(x(2),y(2)),⋯,(x(N),y(N))}(11.A.1)
其中 x ( i ) = ( x 1 ( i ) , x 2 ( i ) , ⋯   , x n ( i ) ) x^{(i)} = (x^{(i)}_1,x^{(i)}_2,\cdots,x^{(i)}_n) x(i)=(x1(i)​,x2(i)​,⋯,xn(i)​)表示第 i i i个样本的观测序列; x j ( i ) x^{(i)}_j xj(i)​表示观测序列 x ( i ) x^{(i)} x(i)中第 j j j个位置的值。同理 y ( i ) = ( y 1 ( i ) , y 2 ( i ) , ⋯   , y n ( i ) ) y^{(i)} = (y^{(i)}_1,y^{(i)}_2,\cdots,y^{(i)}_n) y(i)=(y1(i)​,y2(i)​,⋯,yn(i)​)表示第 i i i个样本的状态序列, y j ( i ) y^{(i)}_j yj(i)​表示状态序列 y ( i ) y^{(i)} y(i)中第 j j j个位置的值。

给定 K 1 K_1 K1​个转移特征函数
t k ( x , y ) = ∑ i = 1 n t k ( y i − 1 , y i , x , i ) , k = 1 , 2 , ⋯   , K 1 (11.A.2) t_k(x,y) = \sum_{i=1}^n t_k(y_{i-1},y_i,x,i), \hspace{1em} k=1,2,\cdots,K_1 \tag{11.A.2} tk​(x,y)=i=1n​tk​(yi−1​,yi​,x,i),k=1,2,⋯,K1​(11.A.2)
给定 K 2 K_2 K2​个状态特征函数
s l ( x , y ) = ∑ i = 1 n s l ( y i , x , i ) , k = 1 , 2 , ⋯   , K 2 (11.A.3) s_l(x,y) = \sum_{i=1}^n s_l(y_i,x,i), \hspace{1em} k=1,2,\cdots,K_2 \tag{11.A.3} sl​(x,y)=i=1n​sl​(yi​,x,i),k=1,2,⋯,K2​(11.A.3)

(二) 模型

根据式(11.15)和式(11.16),条件随机场模型可表示为

P ( y ∣ x ) = 1 Z ( x ) e x p ∑ k = 1 K w k f k ( y , x ) (11.15) P(y|x) = \frac{1}{Z(x)} exp \sum_{k=1}^K w_k f_k(y,x) \tag{11.15} P(y∣x)=Z(x)1​expk=1K​wk​fk​(y,x)(11.15)

Z ( x ) = ∑ y e x p ∑ k = 1 K w k f k ( y , x ) (11.16) Z(x) = \sum_y exp \sum_{k=1}^K w_k f_k(y,x) \tag{11.16} Z(x)=y​expk=1K​wk​fk​(y,x)(11.16)

因为 P ( y ∣ x ) P(y|x) P(y∣x)和 Z ( x ) Z(x) Z(x)中均包含参数 w k w_k wk​,于是将 P ( y ∣ x ) P(y|x) P(y∣x)记作 P w ( y ∣ x ) P_w(y|x) Pw​(y∣x),将 Z ( x ) Z(x) Z(x)记作 Z w ( x ) Z_w(x) Zw​(x)。于是条件随机场模型写作
P w ( y ∣ x ) = 1 Z ( x ) e x p ∑ k = 1 K w k f k ( y , x ) (11.A.4) P_w(y|x) = \frac{1}{Z(x)} exp \sum_{k=1}^K w_k f_k(y,x) \tag{11.A.4} Pw​(y∣x)=Z(x)1​expk=1K​wk​fk​(y,x)(11.A.4)

Z w ( x ) = ∑ y e x p ∑ k = 1 K w k f k ( y , x ) (11.A.5) Z_w(x) = \sum_y exp \sum_{k=1}^K w_k f_k(y,x) \tag{11.A.5} Zw​(x)=y​expk=1K​wk​fk​(y,x)(11.A.5)

(三) 推导

已知条件随机场模型的对数似然函数为
L ( w ) = ∑ x , y P ~ ( x , y ) l o g P w ( y ∣ x ) (11.A.6) L(w) = \sum_{x,y} \tilde{P}(x,y) log P_w(y|x) \tag{11.A.6} L(w)=x,yP~(x,y)logPw​(y∣x)(11.A.6)
将式(11.A.4)和式(11.A.5)代入上式(11.A.6)
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ L(w) & = \sum…
考虑到 Z w ( x ) Z_w(x) Zw​(x)不受变量 y y y的影响,于是上式(11.A.7)可以写成
L ( w ) = ∑ x , y P ~ ( x , y ) ∑ k = 1 K w k f k ( y , x ) − ∑ x P ~ ( x ) l o g Z w ( x ) (11.A.8) L(w) = \sum_{x,y} \tilde{P}(x,y) \sum_{k=1}^K w_k f_k(y,x) - \sum_x \tilde{P}(x) log Z_w(x) \tag{11.A.8} L(w)=x,yP~(x,y)k=1K​wk​fk​(y,x)−xP~(x)logZw​(x)(11.A.8)
假设条件随机场模型当前的参数向量是 w = ( w 1 , w 2 , ⋯   , w n ) T w = (w_1,w_2,\cdots,w_n)^T w=(w1​,w2​,⋯,wn​)T,那么我们希望找到一个新的参数向量 w + δ = ( w 1 + δ 1 , w 2 + δ 2 , ⋯   , w n + δ n ) T w+\delta = (w_1+\delta_1,w_2+\delta_2,\cdots,w_n+\delta_n)^T w+δ=(w1​+δ1​,w2​+δ2​,⋯,wn​+δn​)T,使得模型的对数似然函数值增大。

对于给定的经验分布 P ~ ( x , y ) \tilde{P}(x,y) P~(x,y),模型参数从 w w w到 w + δ w+\delta w+δ,对数似然函数的改变量是
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ L(w+\delta) - …
上式(11.A.9)与《统计学习方法》6.3.1中对数似然函数的改变量式相同,参考书中即可推导出书中式(11.36)、式(11.37)和式(11.38)。

条件随机场模型学习的改进的迭代尺度法(原生Python实现)

import math
from code.maximum_entropy_model import newton_method_linear  # 参见最大熵模型部分

def improved_iterative_scaling(x, y, transfer_features, state_features, tol=1e-4, max_iter=1000):
    """改进的迭代尺度法学习条件随机场模型

    :param x: 输入变量
    :param y: 输出变量
    :param transfer_features: 转移特征函数
    :param state_features: 状态特征函数
    :param tol: 容差
    :param max_iter: 最大迭代次数
    :return: 条件随机场模型
    """
    n_samples = len(x)  # 样本数
    n_transfer_features = len(transfer_features)  # 转移特征数
    n_state_features = len(state_features)  # 状态特征数

    # 坐标压缩
    y_list = list(set(y))
    y_mapping = {c: i for i, c in enumerate(y_list)}
    x_list = list(set(tuple(x[i]) for i in range(n_samples)))
    x_mapping = {c: i for i, c in enumerate(x_list)}

    n_x = len(x_list)  # 观测序列可能取值数
    n_y = len(y_list)  # 状态序列可能取值数
    n_total = n_x * n_y  # 观测序列可能取值和状态序列可能取值的笛卡尔积

    print(x_list, x_mapping)
    print(y_list, y_mapping)

    # 计算联合分布的经验分布:P(X,Y) (empirical_joint_distribution)
    d1 = [[0.0] * n_y for _ in range(n_x)]  # empirical_joint_distribution
    for i in range(n_samples):
        d1[x_mapping[tuple(x[i])]][y_mapping[y[i]]] += 1 / n_samples
    # print("联合分布的经验分布:", d1)

    # 计算边缘分布的经验分布:P(X) (empirical_marginal_distribution)
    d2 = [0.0] * n_x  # empirical_marginal_distribution
    for i in range(n_samples):
        d2[x_mapping[tuple(x[i])]] += 1 / n_samples
    # print("边缘分布的经验分布", d2)

    # 所有特征在(x,y)出现的次数:f#(x,y) (samples_n_features)
    nn = [[0.0] * n_y for _ in range(n_x)]  # samples_n_features

    # 计算转移特征函数关于经验分布的期望值:EP(tk) (empirical_joint_distribution_each_feature)
    d3 = [0.0] * n_transfer_features  # empirical_joint_distribution_each_feature
    for k in range(n_transfer_features):
        for xi in range(n_x):
            for yi in range(n_y):
                x = x_list[xi]
                y = y_list[yi]
                n_position = len(x_list[xi])  # 序列中的位置数
                for i in range(1, n_position):
                    if transfer_features[k](y[i - 1], y[i], x, i):
                        d3[k] += d1[xi][yi]
                        nn[xi][yi] += 1

    # 计算状态特征函数关于经验分布的期望值:EP(sl) (empirical_joint_distribution_each_feature)
    d4 = [0.0] * n_state_features  # empirical_joint_distribution_each_feature
    for l in range(n_state_features):
        for xi in range(n_x):
            for yi in range(n_y):
                x = x_list[xi]
                y = y_list[yi]
                n_position = len(x_list[xi])  # 序列中的位置数
                for i in range(n_position):
                    if state_features[l](y[i], x, i):
                        d4[l] += d1[xi][yi]
                        nn[xi][yi] += 1

    # print("转移特征函数关于经验分布的期望值:", d3)
    # print("状态特征函数关于经验分布的期望值:", d4)
    # print("所有特征在(x,y)出现的次数:", nn)

    # 定义w的初值和模型P(Y|X)的初值
    w1 = [0] * n_transfer_features  # w的初值:wi=0
    w2 = [0] * n_state_features  # w的初值:wi=0
    p0 = [[1 / n_total] * n_y for _ in range(n_x)]  # 当wi=0时,P(Y|X)的值

    for _ in range(max_iter):
        change = False

        # 遍历所有转移特征以更新w
        for k in range(n_transfer_features):
            def func(d, kk):
                """目标方程"""
                res = 0
                for xxi in range(n_x):
                    for yyi in range(n_y):
                        xx = x_list[xxi]
                        yy = y_list[yyi]
                        n_position = len(x_list[xxi])  # 序列中的位置数
                        val = 0
                        for i in range(1, n_position):
                            val += transfer_features[kk](yy[i - 1], yy[i], xx, i)
                        val *= d2[xxi] * p0[xxi][yyi] * pow(math.e, d * nn[xxi][yyi])
                        res += val
                res -= d3[kk]
                return res

            # 牛顿法求解目标方程的根
            dj = newton_method_linear(func, args=(k,))

            # 更新wi的值
            w1[k] += dj
            if abs(dj) >= tol:
                change = True

        for l in range(n_state_features):
            def func(d, ll):
                """目标方程"""
                res = 0
                for xxi in range(n_x):
                    for yyi in range(n_y):
                        xx = x_list[xxi]
                        yy = y_list[yyi]
                        n_position = len(x_list[xxi])  # 序列中的位置数
                        val = 0
                        for i in range(n_position):
                            val += state_features[ll](yy[i], xx, i)
                        val *= d2[xxi] * p0[xxi][yyi] * pow(math.e, d * nn[xxi][yyi])
                        res += val
                res -= d4[ll]
                return res

            # 牛顿法求解目标方程的根
            dj = newton_method_linear(func, args=(l,))

            # 更新wi的值
            w2[l] += dj
            if abs(dj) >= tol:
                change = True

        # 计算新的模型
        p1 = [[0.0] * n_y for _ in range(n_x)]
        for xi in range(n_x):
            for yi in range(n_y):
                x = x_list[xi]
                y = y_list[yi]
                n_position = len(x_list[xi])  # 序列中的位置数
                res = 0
                for k in range(n_transfer_features):
                    for i in range(1, n_position):
                        res += w1[k] * t[k](y[i - 1], y[i], x, i)
                for l in range(n_state_features):
                    for i in range(n_position):
                        res += w2[l] * s[l](y[i], x, i)
                p1[xi][yi] = pow(math.e, res)
            total = sum(p1[xi][yi] for yi in range(n_y))
            if total > 0:
                for yi in range(n_y):
                    p1[xi][yi] /= total

        if not change:
            break

        p0 = p1

    ans = {}
    for xi in range(n_x):
        for yi in range(n_y):
            ans[(tuple(x_list[xi]), y_list[yi])] = p0[xi][yi]
    return w1 + w2, ans

【测试】例11.4(省略特征函数及参数部分)

if __name__ == "__main__":
    # 生成随机模型
    X, Y = make_hidden_sequence(
        w1, t, w2, s,
        [(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)],
        [(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)])

    w, P = improved_iterative_scaling(X, Y, t, s)
    print("学习结果:", [round(elem, 2) for elem in w])
    # 生成训练数据集权重: [1, 0.6, 1.2, 0.2, 1.4, 1, 0.2, 0.8, 0.5]
    # 学习结果: [1.07, 0.75, 0.75, 0.35, 1.38, 1.04, 0.22, 0.67, 0.4] (迭代次数:613)

S算法和T算法的时间复杂度

朴素的改进的迭代尺度法在求解 δ k \delta_k δk​的方程时,需要通过一维搜索求解,且每次迭代中均需要计算方程的值。其时间复杂度为 O ( x y n k ) O(xynk) O(xynk),其中 x x x为观测序列的可能取值数, y y y为状态序列的可能取值数, n n n为序列中的元素数, k k k为迭代次数。

S算法在求解 δ k \delta_k δk​时,可以直接通过解析解求解。其时间复杂度为 O ( x y n ) O(xyn) O(xyn)。

T算法在求解 δ k \delta_k δk​时,需要通过一维搜索求解,但在计算方程的值时,需要遍历 y y y和 n n n的部分成为了定值,只需要计算一次。其时间复杂度为 O ( x y n + x k ) O(xyn+xk) O(xyn+xk)。

条件随机场模型学习的BFGS算法(原生Python实现)

import math
import numpy as np
from code.gradient_descent import golden_section_for_line_search  # 详见梯度下降法部分
from code.gradient_descent import partial_derivative  # 详见梯度下降法部分

def bfgs_algorithm(x, y, transfer_features, state_features, tol=1e-4, distance=20, max_iter=100):
    """拟牛顿法学习条件随机场模型

    :param x: 输入变量
    :param y: 输出变量
    :param transfer_features: 转移特征函数
    :param state_features: 状态特征函数
    :param tol: 容差
    :param distance: 一维搜索倍率
    :param max_iter: 最大迭代次数
    :return: 条件随机场模型
    """
    n_samples = len(x)  # 样本数
    n_transfer_features = len(transfer_features)  # 转移特征数
    n_state_features = len(state_features)  # 状态特征数
    n_features = n_transfer_features + n_state_features  # 特征数

    # 坐标压缩
    y_list = list(set(y))
    y_mapping = {c: i for i, c in enumerate(y_list)}
    x_list = list(set(tuple(x[i]) for i in range(n_samples)))
    x_mapping = {c: i for i, c in enumerate(x_list)}

    n_x = len(x_list)  # 观测序列可能取值数
    n_y = len(y_list)  # 状态序列可能取值数
    n_total = n_x * n_y  # 观测序列可能取值和状态序列可能取值的笛卡尔积

    # 计算联合分布的经验分布:P(X,Y) (empirical_joint_distribution)
    d1 = [[0.0] * n_y for _ in range(n_x)]  # empirical_joint_distribution
    for i in range(n_samples):
        d1[x_mapping[tuple(x[i])]][y_mapping[y[i]]] += 1 / n_samples
    # print("联合分布的经验分布:", d1)

    # 计算边缘分布的经验分布:P(X) (empirical_marginal_distribution)
    d2 = [0.0] * n_x  # empirical_marginal_distribution
    for i in range(n_samples):
        d2[x_mapping[tuple(x[i])]] += 1 / n_samples
    # print("边缘分布的经验分布", d2)

    # 所有特征在(x,y)出现的次数:f#(x,y) (samples_n_features)
    nn = [[0.0] * n_y for _ in range(n_x)]  # samples_n_features

    # 计算转移特征函数关于经验分布的期望值:EP(tk) (empirical_joint_distribution_each_feature)
    d3 = [0.0] * n_transfer_features  # empirical_joint_distribution_each_feature
    for k in range(n_transfer_features):
        for xi in range(n_x):
            for yi in range(n_y):
                x = x_list[xi]
                y = y_list[yi]
                n_position = len(x_list[xi])  # 序列中的位置数
                for i in range(1, n_position):
                    if transfer_features[k](y[i - 1], y[i], x, i):
                        d3[k] += d1[xi][yi]
                        nn[xi][yi] += 1

    # 计算状态特征函数关于经验分布的期望值:EP(sl) (empirical_joint_distribution_each_feature)
    d4 = [0.0] * n_state_features  # empirical_joint_distribution_each_feature
    for l in range(n_state_features):
        for xi in range(n_x):
            for yi in range(n_y):
                x = x_list[xi]
                y = y_list[yi]
                n_position = len(x_list[xi])  # 序列中的位置数
                for i in range(n_position):
                    if state_features[l](y[i], x, i):
                        d4[l] += d1[xi][yi]
                        nn[xi][yi] += 1

    # print("转移特征函数关于经验分布的期望值:", d3)
    # print("状态特征函数关于经验分布的期望值:", d4)
    # print("所有特征在(x,y)出现的次数:", nn)

    def func(ww):
        """目标函数"""
        res = 0
        for xxi in range(n_x):
            xx = x_list[xxi]
            n_position = len(xx)
            t1 = 0
            for yyi in range(n_y):
                yy = y_list[yyi]
                t2 = 0
                for kk in range(n_transfer_features):
                    for i in range(1, n_position):
                        if transfer_features[kk](yy[i - 1], yy[i], xx, i):
                            t2 += ww[kk]
                for ll in range(n_state_features):
                    for i in range(n_position):
                        if state_features[ll](yy[i], xx, i):
                            t2 += ww[ll + n_transfer_features]
                t1 += pow(math.e, t2)
            res += d2[xxi] * math.log(t1, math.e)

        for xxi in range(n_x):
            xx = x_list[xxi]
            n_position = len(xx)
            for yyi in range(n_y):
                yy = y_list[yyi]
                t3 = 0
                for kk in range(n_transfer_features):
                    for i in range(1, n_position):
                        if transfer_features[kk](yy[i - 1], yy[i], xx, i):
                            t3 += ww[kk]
                for ll in range(n_state_features):
                    for i in range(n_position):
                        if state_features[ll](yy[i], xx, i):
                            t3 += ww[ll + n_transfer_features]
                res -= d1[xxi][yyi] * t3

        return res

    # 定义w的初值和B0的初值
    w0 = [0] * n_features  # w的初值:wi=0
    B0 = np.identity(n_features)  # 构造初始矩阵G0(单位矩阵)

    for k in range(max_iter):
        print("迭代次数:", k, "->", [round(elem, 2) for elem in w0])

        # 计算梯度 gk
        nabla = partial_derivative(func, w0)

        g0 = np.matrix([nabla]).T  # g0 = g_k

        # 当梯度的模长小于精度要求时,停止迭代
        if pow(sum([nabla[i] ** 2 for i in range(n_features)]), 0.5) < tol:
            break

        # 计算pk
        if k == 0:
            pk = - B0 * g0  # 若numpy计算逆矩阵时有0,则对应位置会变为inf
        else:
            pk = - (B0 ** -1) * g0

        # 一维搜索求lambda_k
        def f(xx):
            """pk 方向的一维函数"""
            x2 = [w0[jj] + xx * float(pk[jj][0]) for jj in range(n_features)]
            return func(x2)

        lk = golden_section_for_line_search(f, 0, distance, epsilon=1e-6)  # lk = lambda_k

        # print(k, "lk:", lk)

        # 更新当前点坐标
        w1 = [w0[j] + lk * float(pk[j][0]) for j in range(n_features)]

        # print(k, "w1:", w1)

        # 计算g_{k+1},若模长小于精度要求时,则停止迭代
        # 计算新的模型

        nabla = partial_derivative(func, w1)

        g1 = np.matrix([nabla]).T  # g0 = g_{k+1}

        # 当梯度的模长小于精度要求时,停止迭代
        if pow(sum([nabla[i] ** 2 for i in range(n_features)]), 0.5) < tol:
            w0 = w1
            break

        # 计算G_{k+1}
        yk = g1 - g0
        dk = np.matrix([[lk * float(pk[j][0]) for j in range(n_features)]]).T

        # 当更新的距离小于精度要求时,停止迭代
        if pow(sum([lk * float(pk[j][0]) ** 2 for j in range(n_features)]), 0.5) < tol:
            w0 = w1
            break

        B1 = B0 + (yk * yk.T) / (yk.T * dk) + (B0 * dk * dk.T * B0) / (dk.T * B0 * dk)

        B0 = B1
        w0 = w1

    # 计算模型结果
    p0 = [[0.0] * n_y for _ in range(n_x)]
    for xi in range(n_x):
        for yi in range(n_y):
            x = x_list[xi]
            y = y_list[yi]
            n_position = len(x_list[xi])  # 序列中的位置数
            res = 0
            for k in range(n_transfer_features):
                for i in range(1, n_position):
                    res += w0[k] * t[k](y[i - 1], y[i], x, i)
            for l in range(n_state_features):
                for i in range(n_position):
                    res += w0[n_transfer_features + l] * s[l](y[i], x, i)
            p0[xi][yi] = pow(math.e, res)
        total = sum(p0[xi][yi] for yi in range(n_y))
        if total > 0:
            for yi in range(n_y):
                p0[xi][yi] /= total

    # 返回最终结果
    ans = {}
    for xi in range(n_x):
        for yi in range(n_y):
            ans[(tuple(x_list[xi]), y_list[yi])] = p0[xi][yi]
    return w0, ans

【测试】例11.4(省略特征函数及参数部分)

if __name__ == "__main__":
    # 生成随机模型
    w1, t, w2, s = crf_model()
    X, Y = make_hidden_sequence(
        w1, t, w2, s,
        [(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)],
        [(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)])

    w, P = bfgs_algorithm(X, Y, t, s)
    print("学习结果:", [round(elem, 2) for elem in w])
    # 生成训练数据集权重: [1, 0.6, 1.2, 0.2, 1.4, 1, 0.2, 0.8, 0.5]
    # 学习结果: [0.82, 0.67, 0.2, 0.04, 1.21, 1.04, 0.11, 0.63, 0.35] (迭代次数:15)

11.5 条件随机场的预测算法

条件随机场预测的维特比算法

def viterbi_algorithm(w1, transfer_features, w2, state_features, x, n_state):
    """维特比算法预测状态序列

    :param w1: 模型的转移特征权重
    :param t: 模型的转移特征函数
    :param w2: 模型的状态特征权重
    :param s: 模型的状态特征函数
    :param x: 需要计算的观测序列
    :param n_state: 状态的可能取值数
    :return: 最优可能的状态序列
    """
    n_transfer_features = len(transfer_features)  # 转移特征数
    n_state_features = len(state_features)  # 状态特征数
    n_position = len(x)  # 序列中的位置数

    # 定义状态矩阵
    dp = [[0.0] * n_state for _ in range(n_position)]  # 概率最大值
    last = [[-1] * n_state for _ in range(n_position)]  # 上一个结点

    # 处理t=0的情况
    for i in range(n_state):
        for l in range(n_state_features):
            dp[0][i] += w2[l] * state_features[l](y0=i, x=x, i=0)

    # 处理t>0的情况
    for t in range(1, n_position):
        for i in range(n_state):
            for j in range(n_state):
                d = dp[t - 1][i]
                for k in range(n_transfer_features):
                    d += w1[k] * transfer_features[k](y0=i, y1=j, x=x, i=t)
                for l in range(n_state_features):
                    d += w2[l] * state_features[l](y0=j, x=x, i=t)
                # print((i, j), "=", d)
                if d >= dp[t][j]:
                    dp[t][j] = d
                    last[t][j] = i
        # print(dp[t], last[t])

    # 计算最优路径的终点
    best_end, best_gamma = 0, 0
    for i in range(n_state):
        if dp[-1][i] > best_gamma:
            best_end, best_gamma = i, dp[-1][i]

    # 计算最优路径
    ans = [0] * (n_position - 1) + [best_end]
    for t in range(n_position - 1, 0, -1):
        ans[t - 1] = last[t][ans[t]]
    return ans

【测试】例11.3(将书中例子从1开始的数组改为从0开始)

if __name__ == "__main__":
    def t1(y0, y1, x, i):
        return int(y0 == 0 and y1 == 1 and i in {1, 2})
    def t2(y0, y1, x, i):
        return int(y0 == 0 and y1 == 0 and i in {1})
    def t3(y0, y1, x, i):
        return int(y0 == 1 and y1 == 0 and i in {2})
    def t4(y0, y1, x, i):
        return int(y0 == 1 and y1 == 0 and i in {1})
    def t5(y0, y1, x, i):
        return int(y0 == 1 and y1 == 1 and i in {2})
    def s1(y0, x, i):
        return int(y0 == 0 and i in {0})
    def s2(y0, x, i):
        return int(y0 == 1 and i in {0, 1})
    def s3(y0, x, i):
        return int(y0 == 0 and i in {1, 2})
    def s4(y0, x, i):
        return int(y0 == 1 and i in {2})
    w1 = [1, 0.6, 1, 1, 0.2]
    t = [t1, t2, t3, t4, t5]
    w2 = [1, 0.5, 0.8, 0.5]
    s = [s1, s2, s3, s4]
	
    import random
    print(viterbi_algorithm(w1, t, w2, s, [random.randint(0, 1) for _ in range(3)], 2))  # [0, 1, 0]