保姆韦尔奇方法Baum-Welch

保姆韦尔奇方法Baum-Welch

本例中的代码为用保姆韦尔奇方法实现的机器学习算法，请轻拍

# -*- coding:utf-8 -*-  
import sys

# alpha函数
def alpha(A, B, PI, O): 
    a = [[0 for _ in range(len(A))] for _ in range(len(O))]
    # 设置初始化的值
    for i in range(len(A)):
        a[0][i] = PI[i] * B[i][O[0]]
    # 根据前向算法计算概率
    for t in range(1, len(O)):
        for j in range(len(A)):   # 状态转换到达的下一个状态
            sum = 0
            for i in range(len(A)):
                sum += a[t-1][i] * A[i][j]
            a[t][j] = sum * B[j][O[t]]
    return a


# beta函数
def beta(A, B, PI, O):
    b = [[0 for _ in range(len(A))] for _ in range(len(O))]
    # 初始化
    for j in range(len(A)):
        b[len(O) - 1][j] = 1
    # 根据后向算法公式计算后向概率
    for t in reversed(range(0, len(O) - 1)):
        for i in range(len(A)):
            sum = 0 
            for j in range(len(A)):
                sum += A[i][j] * B[j][O[t + 1]] * b[t + 1][j]
            b[t][i] = sum
    return b


# 计算zeta的值
def zeta(alpha, beta, A, B, O):
    numerator = [[[0 for _ in range(len(A))] for _ in range(len(A))] for _ in range(len(O))]
    # 计算分子
    for t in range(len(O)):
        for i in range(len(A)):
            for j in range(len(A)):
                if t < len(O) - 2:
                    numerator[t][i][j] = alpha[t][i] * A[i][j] * B[j][O[t + 1]] * beta[t + 1][j]
                elif t < len(O) - 1:
                    numerator[t][i][j] = alpha[t][i] * A[i][j] * B[j][O[t + 1]]
                else:
                    numerator[t][i][j] = alpha[t][i]
    # 计算分母
    denominator = [0 for _ in range(len(O))]
    for t in range(len(O)):
        denominator[t] = 0 ;
        for i in range(len(A)):
            for j in range(len(A)):
                denominator[t] += numerator[t][i][j]

    # 计算整个zeta的值
    zeta = [[[0 for _ in range(len(A))] for _ in range(len(A))] for _ in range(len(O))]
    for t in range(len(O)):
        for i in range(len(A)):
            for j in range(len(A)):
                zeta[t][i][j] = numerator[t][i][j] / denominator[t]
    return zeta


def gamma(zeta, O):
    gamma = [[0 for _ in range(len(zeta[0]))] for _ in range(len(O))]
    for t in range(len(O)):
        for i in range(len(zeta[0])):
            gamma[t][i] = 0
            for j in range(len(zeta[0])):
                gamma[t][i] += zeta[t][i][j]
    return gamma

def calc_A(zeta, gamma):
    a = [[0 for _ in range(len(zeta[0]))] for _ in range(len(zeta[0]))]
    for i in range(len(zeta[0])):
        for j in range(len(zeta[0])):
            a[i][j] = sum(zeta[t][i][j] for t in range(len(zeta)-1)) / sum(gamma[t][i] for t in range(len(zeta)-1))
    return a

def calc_B(gamma, O, v):
    b = [[0 for _ in range(len(set(O)))] for _ in range(len(gamma[0]))]
    for j in range(len(gamma[0])):
        for k in range(len(set(O))):
            b[j][k] = sum([gamma[t][j] for t in range(len(O)) if O[t] == v[k] ]) / sum([gamma[t][j] for t in range(len(O))])
    return b

def calc_PI(gamma, A):
    pi = [0 for _ in range(len(A))]
    for i in range(len(A)):
        pi[i] = gamma[0][i]
    return pi

# v = observations
def baum_welch(A, B, PI , O, e, v):
    done = True
    while done:
        a = alpha(A, B, PI, O)
        b = beta(A, B, PI, O)
        z = zeta(a, b, A, B, O)
        g = gamma(z, O)
        newA = calc_A(z, g)
        newB = calc_B(g, O, v)
        newPI = calc_PI(g, A)
        print newPI
        A = newA
        B = newB
        PI = newPI
        # 达到一个阈值就可以停止循环。否则一直下去会等于0.
        if False:
            break



# 把状态转换map转化成为矩阵
def matrix(X, index1, index2):
    matrix = [[0 for _ in range(len(index2))] for _ in range(len(index1))]
    for row in X :
        for col in X[row] :
            matrix[index1.index(row)][index2.index(col)] = X[row][col]
    return matrix


if __name__ == "__main__":
    # 初始化，随机的给参数A、B、PI赋值
    states = ["相处", "拜拜"]   # 状态列表
    # 撒娇，小拳拳捶你胸口 
    # 状态输出概率，状态发射概率。
    observations = ["撒娇", "低头玩手机", "眼神很友好", "主动留下联系方式"] 
    A = {"相处": {"相处": 0.5, "拜拜": 0.5}, "拜拜": {"相处":.5, "拜拜": 0.5}}
    B = {"相处":{"撒娇":0.4,"低头玩手机":0.1,"眼神很友好":0.3,"主动留下联系方式":0.2},
         "拜拜":{"撒娇":0.1,"低头玩手机":0.5,"眼神很友好":0.2,"主动留下联系方式":0.2}
        }
    PI = [0.5, 0.5]                     # 初始状态概率
    O  = [1, 2, 0, 2, 3, 0]             # 时间t范围内的输出
    print len(set(O))
    # sys.exit()
    A = matrix(A, states, states)       # 状态转换概率
    B = matrix(B, states, observations) # 状态发射（输出）概率

    observations = [0, 1, 2, 3]
    print A
    print B
    print observations
    a = alpha(A, B, PI, O)
    b = beta(A, B, PI, O)
    baum_welch(A, B, PI , O, 0.05, observations)
... prompt'''