用隐马尔可夫模型推测女朋友心情(Ⅳ)：训练样本的生成

曾在第一期提到过，王同学没有女朋友心情序列与事件序列的样本，他只凭着对女友的了解，预估了三个参数矩阵的概率分布。
这一期我们要做一件有趣的事，就是用王同学给的模型参数来生成训练样本。下一期，我们再用所生成的样本训练一个新的模型，看看这个新模型的参数值是否与王同学给的参数值一致（或者说十分接近）。

1 样本生成算法

图-1 隐马尔可夫模型

从图-1我们可以看出，生成样本的过程就是从第一个时刻开始，逐步得到每一时刻的隐状态与观测值的过程。
在第一个时刻，从初始概率向量$\mathbf{\pi }$中抽取一个隐状态，再利用抽取到的隐状态，从发射概率矩阵的概率分布中抽取一个观测值。
在第二个时刻，利用上一时刻的隐状态从转移概率矩阵$\mathbf{A}$的概率分布中抽取一个隐状态，再利用抽取到的隐状态，从发射概率矩阵$\mathbf{B}$的概率分布中抽取一个观测值。依此类推，得到给定长度的隐状态序列与观测序列。
当然，也可以先依次采样每个时刻的隐状态，再依次采样每个时刻的观测值。下面用数学语言描述这个过程：

假设序列长度为$T$，生成步骤如下：
(1) 根据初始状态概率向量$\mathbf{\pi }$采样第一个时刻的状态$y_1=s_i$，即$y_1\sim \mathbf{\pi }$；
(2) $y_t$采样结束得到$s_i$后，根据状态转移概率矩阵第$i$行的概率向量${\mathbf{A}}_{\mathbf{i},:}$采样下一时刻的状态$y_{t+1}$，即$y_{t+1}\sim {\mathbf{A}}_{\mathbf{i},:}$；
(3) 对每个$y_t=s_i$，根据发射概率矩阵的第$i$行${\mathbf{B}}_{\mathbf{i},:}$采样$x_t$，即$x_t\sim {\mathbf{B}}_{\mathbf{i},:}$；
(4) 重复步骤(2)共计$T-1$次，重复步骤(3)共计$T$次，输出序列$\mathbf{x}$与$\mathbf{y}$；

2 代码实现

下面我们根据上面的算法描述，用Python语言实现这个算法。

(1) 采样函数

上述算法中多次出现“根据概率矩阵的某行向量来采样“的文字表述，也就是说，代码要实现依概率分布来采样的功能。依概率分布采样功能的实现如下方代码所示：

def sort_index(lst):
    lst_idx = [i for i in range(len(lst))]
    idx_sorted = []
    lst_sorted = []

    for i in range(len(lst) - 1):
        min_idx = 0
        for j in range(1, len(lst_idx)):
            if lst[lst_idx[min_idx]] > lst[lst_idx[j]]:
                min_idx = j

        idx = lst_idx.pop(min_idx)
        idx_sorted.append(idx)
        lst_sorted.append(lst[idx])

    idx_sorted.append(lst_idx[0])
    lst_sorted.append(lst[lst_idx[0]])

    return lst_sorted, idx_sorted

def sample_from(p):
    p_sorted, idx_sorted = sort_index(p)
    p2dis = []
    p2dis.append(p_sorted[0])
    sum = p_sorted[0]
    for i in range(1, len(p_sorted)):
        sum += p_sorted[i]
        p2dis.append(sum)

    a, b = 0, 1
    n = random.uniform(a, b)
    i = bisect.bisect_left(p2dis, n)
    val = idx_sorted[i]

    return val

依概率采样的基本思想是：依均匀分布随机生成一个[0, 1]的数，因为是均匀分布，所以这个数落在某个区间的概率只与该区间长度有关。因此，我们只要根据概率分布数组元素的概率值构建好大小不同的区间，区间的长度等于概率分布数组的元素值（概率值）。那么依均分布生成的随机数$x$落在不同区间的次数自然也就依概率分布了。

(2) 单个样本生成函数

# 初始状态概率数组
PI = [0.7, 0.2, 0.1]

# 状态转移概率数组
A = [
    [0.5, 0.3, 0.2],
    [0.4, 0.5, 0.1],
    [0.6, 0.1, 0.3],
]

# 发射概率数组
B = [
    [0.4, 0.1, 0.2, 0.3],
    [0.2, 0.5, 0.1, 0.2],
    [0.1, 0.1, 0.6, 0.2],
]

def generate_one(length):
    y_x = [[-1] * length for i in range(2)]
    y_x[0][0] = sample_from(PI)  # 采样首个隐状态
    y_x[1][0] = sample_from(B[y_x[0][0]])  # 根据首个隐状态采样观测值
    for t in range(1, length):
        y_x[0][t] = sample_from(A[y_x[0][t - 1]])
        y_x[1][t] = sample_from(B[y_x[0][t]])

    return y_x

(3) 多个样本生成函数

def generate(minLen, maxLen, num):
    """
    :param minLen: 样本序列最小长度
    :param maxLen: 样本序列最大长度
    :param num: 样本个数
    :return: 样本
    """
    seqs = []
    for i in range(num):
        length = random.randint(minLen, maxLen)
        y_x = generate_one(length)
        seqs.append(y_x)

    return seqs

generate()函数所生成样本序列的最小长度为minLen，最大长度为maxLen，共生成num个样本。

(4) 将生成结果转换成文字，看看都生成些了什么。

# 隐状态字典
MAP_STATE = {
    0: 'neutral',
    1: 'happy',
    2: 'unhappy',
}

# 事件字典
MAP_ACTIVITY = {
    0: 'reading',
    1: 'shopping',
    2: 'running',
    3: 'watching',
}

def index2activity(idx):
    activity = []
    for i in idx:
        act = MAP_ACTIVITY[i]
        activity.append(act)

    return activity


def index2state(index):
    emo_list = []
    for i in index:
        emo_list.append(MAP_STATE[i])

    return emo_list
    
def convert_index_to_word(mat):
    for yx in mat:
        emo = index2state(yx[0])
        act = index2activity(yx[1])
        print(" ".join(a + '/' + e for a, e in zip(act, emo)))
        

(5) main函数：生成样本，打印结果。

if __name__ == '__main__':
    mat = generate(5, 7, 2)
    convert_index_to_word(mat)    

(6) 程序运行结果