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

隠れマルコフモデルの例のつづき。

さて、たとえば、友人が初日に「散歩」二日目に「買い物」三日目に「掃除」という順で行動したら、その観測結果が得られる確率はいくらでしょうか、そして、このような観測結果が得られたとき三日間の天気はどのようであったでしょうか。前の疑問に対しては前向きアルゴリズム、後の疑問についてはビタビアルゴリズムを利用できます。（観測された事象系列を結果として生じる隠された状態の最も尤もらしい並びをビタビ経路と呼びます。）これらは、同じコードで実装することができます。以前のエントリーと同様、状態遷移図は次のように書けるとします。

まずアルゴリズムを理解するために初日に「散歩」したという情報のみがある場合を考えます。
初日の状態確率は（添字は何日目かを表すとします）

Pr(雨₀)=0.6, Pr(晴れ₀)=0.4

初日「雨」で「散歩」する確率は

Pr(雨₀,散歩)=Pr(雨₀)*Pr(散歩|雨₀)=0.6*0.1=0.06

このとき次の日「雨」の確率は

Pr(雨₀,散歩,雨₁)=Pr(雨₀,散歩)*Pr(雨₁|雨₀)=0.06*0.7=0.042

同様に初日「晴れ」で「散歩」して次の日が「雨」の確率は

Pr(晴れ₀,散歩,雨₁)=0.4*0.6*0.4=0.096

つまり初日に「散歩」して次の日が「雨」であるとき

Pr(散歩,雨₁)=Pr(雨₀,散歩,雨₁)+Pr(晴れ₀,散歩,雨₁)=0.042+0.096=0.138

となります。Pr(雨₀,散歩,雨₁)<Pr(晴れ₀,散歩,雨₁) であることから、「晴れ→雨」であった可能性が高い、ということがわかります。

二日目が「晴れ」である場合についても計算すると「雨」「散歩」「晴れ」は

Pr(雨₀,散歩,晴れ₁)=0.6*0.1*0.3=0.018

「晴れ」「散歩」「晴れ」は

Pr(晴れ₀,散歩,晴れ₁)=0.4*0.6*0.6=0.144

つまり「散歩」して次の日が「晴れ」になるとき

Pr(散歩,晴れ₁)=0.018+0.144=0.162

となり、また「晴れ→晴れ」でありそうだということがわかります。

これらを併せて考えれば、「散歩」をする全体確率は

Pr(散歩)=0.138+0.162=0.3

このときビタビ経路は「晴れ→晴れ」となります。

では初日に「散歩」、二日目に「買い物」した場合はどうなるでしょう。全体確率とビタビ経路を求めます。
初日に「散歩」したとき、二日目の「雨」である確率は0.138、「晴れ」である確率は0.162であるとわかっています。
また、「雨」に至るまでの尤もな経路は「晴れ→雨」、「晴れ」に至るまでの尤もな経路は「晴れ→晴れ」だとわかっています。


（つづく？）

→ Viterbi　algorithm わかりやすくなってる！！！
→ 実装（コピペ）

	#!/usr/bin/env python
	# -*- coding:utf-8 -*-
	# Viterbi algorithm
	# http://en.wikipedia.org/wiki/Viterbi_algorithm
	#
	# > python viterbi.py
	# 0 1 2
	# Rainy: 0.06000 0.03840 0.01344
	# Sunny: 0.24000 0.04320 0.00259
	# (0.01344, ['Sunny', 'Rainy', 'Rainy'])
	# HMM
	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},
	}
	# Helps visualize the steps of Viterbi.
	def print_dptable(V):
	print " ",
	for i in range(len(V)): print "%7d" % i,
	print
	for y in V[0].keys():
	print "%.5s: " % y,
	for t in range(len(V)):
	print "%.7s" % ("%f" % V[t][y]),
	print
	def viterbi(obs, states, start_p, trans_p, emit_p):
	V = [{}]
	path = {}
	# Initialize base cases (t == 0)
	for y in states:
	V[0][y] = start_p[y] * emit_p[y][obs[0]]
	path[y] = [y]
	# Run Viterbi for t > 0
	for t in range(1,len(obs)):
	V.append({})
	newpath = {}
	for y in states:
	(prob, state) = max([(V[t-1][y0] * trans_p[y0][y] * emit_p[y][obs[t]], y0) for y0 in states])
	V[t][y] = prob
	newpath[y] = path[state] + [y]
	# Don't need to remember the old paths
	path = newpath
	print_dptable(V)
	(prob, state) = max([(V[len(obs) - 1][y], y) for y in states])
	return (prob, path[state])
	print viterbi(observations,
	states,
	start_probability,
	transition_probability,
	emission_probability)

view raw viterbi.py hosted with ❤ by GitHub