Pythonで隠れマルコフモデルのFilteringの例

http://en.wikipedia.org/wiki/Viterbi_algorithm
の天気の例で使えるFiltering。

隠れマルコフモデルの例
隠れマルコフモデルの例　その２

フィルタリング

• 現在までのすべての観測があるときに，現在の状態を推定する．
• i.e. P(X_t|Y_0:t)を求める．

P(X_t|Y_0:t)をガリガリ式変形する（マルコフ性，ベイズ，周辺化でごにょごにょする）と，

が得られる．P(X_t|X_t-1)は遷移確率．P(Y_t|X_t)は出力確率．cは正規化項．で，再帰的にP(X_t|Y_0:t)を求める．以下実装．

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# 状態変数
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},
   }

# Forward Algorithm
def forward(y1, prior, states, trans_p, emit_p):
    post = {}
    for x1 in states:
        post[x1] = emit_p[x1][y1] * sum([ trans_p[x][x1] * prior[x] for x in states ])
    s = sum(post.values())
    for k,v in post.items():
        post[k] = v/s
    return post

# Filtering
def filtering(observations, states, start_p, trans_p, emit_p):
    prior = start_p
    T = len(observations)
    for t in range(T):
        post = forward( observations[t], prior, states, trans_p, emit_p )
        prior = post
    return post

print filtering(observations, states, start_probability, transition_probability, emission_probability)

上記の実行結果は

{'Rainy': 0.86342066067036916, 'Sunny': 0.1365793393296309}

つまり， walk→shop→clean という観測列があったら3日目の天気は雨が86%，晴れが14%だよっ，ということ．