元理系院生の新入社員がPythonとJavaで色々頑張るブログ

プログラミングや機械学習について調べた事を書いていきます

PythonでLeft to Right型隠れマルコフモデルとBaumwelchアルゴリズム

この記事に書かれていること

Pythonによる隠れマルコフモデル(Hidden Markov Model)とBaumwelchの実装
http://unicorn.ike.tottori-u.ac.jp/murakami/paper/STUDY/IEICE_2010_07/main.pdf

隠れマルコフモデルとBaumwelchアルゴリズムのとても親切な資料が公開されて居たので、Pythonで実装してみました。

時間が有る時に、隠れマルコフモデルの記事も書きたいですね。

また、Androidにおいて手書き数字文字のstroke情報をサンプリングし、strokeの移動角度を16分割にコードブック化し文字認識をした所、状態数1の時が最も識別率が高くなりました。(90%くらい)

Pythonによる隠れマルコフモデルの実装

__author__ = 'emoson'


def forward(A, B, O):
    """
    forwardアルゴリズム
    """
    alpha = [[0.0 for j in range(len(A)+1)] for t in range(len(O)+1)]
    alpha[0][0] = 1.0
    for t in range(1, len(alpha)):
        for j in range(len(alpha[t])):
            if j < len(alpha[t])-1:
                alpha[t][j] += alpha[t-1][j] * A[j][j] * B[j][O[t-1]]
            if j > 0:
                alpha[t][j] += alpha[t-1][j-1] * A[j-1][j] * B[j-1][O[t-1]]
    return alpha


def backward(A, B, O):
    """
    backwardアルゴリズム
    """
    beta = [[0.0 for j in range(len(A)+1)] for t in range(len(O)+1)]
    beta[-1][-1] = 1.0
    for t in range(len(beta)-2, -1, -1):
        for j in range(len(beta)-1, -1, -1):
            if j < len(beta[t])-1:
                beta[t][j] = beta[t+1][j+1] * A[j][j+1] * B[j][O[t]] + beta[t+1][j] * A[j][j] * B[j][O[t]]
    return beta


def expectation(A, B, O, alpha, beta):
    """
    期待値ステップ
    """
    gamma = [[[0.0 for j in range(len(A)+1)] for i in range(len(A))] for t in range(len(O))]
    for t in range(len(gamma)):
        for i in range(len(gamma[t])):
            for j in range(len(gamma[t][i])):
                gamma[t][i][j] = alpha[t][i]*A[i][j]*B[i][O[t]]*beta[t+1][j] / alpha[-1][-1]
    return gamma


def maximization(A, B, O, gamma):
    """
    最大化ステップ
    """
    #状態遷移確率の再計算
    for i in range(len(A)):
        for j in range(len(A[i])):
            a, b = 0, 0
            for t in range(len(gamma)):
                a += gamma[t][i][j]
                for _j in range(len(A[i])):
                    b += gamma[t][i][_j]
            A[i][j] = a / b

    #出力確率の再計算
    for j in range(len(B)):
        for o in range(len(B[j])):
            a, b = 0, 0
            for t in range(len(gamma)):
                for k in range(len(gamma[t][j])):
                    if o == O[t]:
                        a += gamma[t][j][k]
                    b += gamma[t][j][k]
            B[j][o] = a / b


def learning(N, K, O, max_step=None, A=None, B=None):
    """
    Left to Right HiddenMarcovModelをbaum-welchにより学習
    :param N:           状態数
    :param K:           シンボル数
    :param O:           時系列データ
    :param max_step:    最大学習回数
    :return:
    """
    #初期状態確率
    Pi = [1.0 if n == 0 else 0.0 for n in range(N)]
    #状態遷移確率
    if A is None:
        A = [[1.0/(N+1) for j in range(N+1)] for i in range(N)]
    #シンボル出力確率
    if B is None:
        B = [[1.0/K for j in range(K)] for i in range(N)]

    prev_likelihood = 0
    step = 0
    while True:
        alpha = forward(A, B, O)
        beta = backward(A, B, O)
        if not prev_likelihood < beta[0][0]:
            break
        if max_step is not None and step > max_step:
            break

        prev_likelihood = beta[0][0]
        gamma = expectation(A, B, O, alpha, beta)
        maximization(A, B, O, gamma)
        step += 1
    return A, B