西尾泰和のブログ

とりあえずHidden Markov model - Wikipedia, the free encyclopediaに載っていたPythonコードをきちんと動くようにしてみました。

def choise(dic):
    "choise key from dictionary using value(probability)"
    from random import random
    rest_prob = 1.0
    for k in dic:
        p = dic[k]
        if random() * rest_prob < p:
            return k
        rest_prob -= p


class StateMachine:
    states = ('Rainy', 'Sunny')

    observations = ('walk', 'shop', 'clean')

    start_probability = {'Rainy': 0.6, 'Sunny': 0.4}

    transition_probability = {
       'Rainy' : {'Rainy': 0.7, 'Sunny': 0.3},
       'Sunny' : {'Rainy': 0.4, 'Sunny': 0.6},
    }

    emission_probability = {
       'Rainy' : {'walk': 0.1, 'shop': 0.4, 'clean': 0.5},
       'Sunny' : {'walk': 0.6, 'shop': 0.3, 'clean': 0.1},
    }

    def __init__(self):
        sp = self.start_probability
        self.currentState = choise(sp)

    def emit(self):
        ep = self.emission_probability
        s = self.currentState
        return choise(ep[s])

    def transit(self):
        tp = self.transition_probability
        s = self.currentState
        self.currentState = choise(tp[s])

でもRainyとかが何度も出てくるし、「states = ('Rainy', 'Sunny')」や「observations = ('walk', 'shop', 'clean')」が全く無意味で面白くないですね。

def generate_dict(keys, values):
    d = {}
    for (i, k) in enumerate(keys):
        d[k] = values[i]
    return d

start_probability = generate_dict(states, (0.6, 0.4))

transition_probability = generate_dict(
    states, [generate_dict(states, x) for x in
             ((0.7, 0.3),
              (0.4, 0.6))])

emission_probability = generate_dict(
    states, [generate_dict(observations, x) for x in
             ((0.1, 0.4, 0.5),
              (0.6, 0.3, 0.1))])

うわー。

def build_prob_table(keysList, valueTable):
    if len(keysList) == 1:
        d = {}
        for (i, k) in enumerate(keysList[0]):
            d[k] = valueTable[i]
        return d
    result = build_prob_table(
        [keysList[0]],
        [build_prob_table(keysList[1:], x) for x in valueTable])
    return result

start_probability = build_prob_table([states], (0.6, 0.4))

transition_probability = build_prob_table(
    [states, states],
    ((0.7, 0.3),
     (0.4, 0.6)))

emission_probability = build_prob_table(
    [states, observations],
    ((0.1, 0.4, 0.5),
     (0.6, 0.3, 0.1)))

これで2次元に限らず何次元でも行けますよ！(でも仕様上2次元までしかないのだけど。過度の抽象化。)

とりあえずViterbiアルゴリズムを実装してみました。'walk', 'walk', 'shop', 'walk', 'walk', 'clean', 'clean', 'clean', 'clean', 'clean'を入力すると'Sunny', 'Sunny', 'Sunny', 'Sunny', 'Sunny', 'Rainy', 'Rainy', 'Rainy', 'Rainy', 'Rainy'が帰ってきます。つまり、shop単体ではその日がRainyであった尤度が高いのだけども、前後をSunnyである尤度の高いwalkで挟まれているので全体としては「shopの日もSunnyであろう」という推定がもっとも尤度が高くなるわけですね。一方、walkとcleanの間にshopを挟んで 'walk', 'walk', 'shop', 'walk', 'shop', 'clean', 'clean', 'clean', 'clean', 'clean'にした場合は、'Sunny', 'Sunny', 'Sunny', 'Sunny', 'Rainy', 'Rainy', 'Rainy', 'Rainy', 'Rainy', 'Rainy'となり、shopの日はRainyだと判定されるわけですね。

いわゆる「いかさまサイコロ推定問題」も実装してみました。といっても確率分布をいじるだけ。

class Casino(StateMachine):
    states = "FL"

    observations = "123456"

    start_probability = build_prob_table([states], (0.5, 0.5))

    transition_probability = build_prob_table(
        [states, states],
        ((0.9, 0.1),
         (0.2, 0.8)))

    emission_probability = build_prob_table(
        [states, observations],
        ((0.16, 0.16, 0.17, 0.17, 0.17, 0.17),
         (0.1,  0.1 , 0.1,  0.1,  0.1,  0.5)))

c = Casino()
rolls = ""
hidden = ""
for i in range(60):
    hidden += c.transit()
    rolls += c.emit()

print hidden
print rolls
print "".join(c.viterbi(rolls)[0])

LLFFFFFFFFLLLLFFFFFFFFLFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLL
661225134463466546136164411511233413411161362663654346666623
LLFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLL

短い範囲のいかさまはかなり取りこぼしますね。でもたぶんこれで正しい挙動のはず。