「隠れマルコフモデル」を少し勉強してみたので、まとめてみる。
モデルの考え方について一番参考になったのは、ここ（pdf）と、ここ（pdf）と、ここだ。
観測事象が図上部の折れ線グラフに示されている。ここで観測事象は何でもよい。株価の騰落率の時間変化でも良いし、ある店舗の売り上げの時間変化でも良い。ここで例にあげた図の観測事象は、時間に関して初期・中期・後期のそれぞれで振る舞いが異なるように見える。
これをモデルに落とし込むため隠れマルコフモデルは「隠れた定常状態（上の例では３つ）が存在し、各定常状態の確率分布（図下部の青・赤・緑で示された確率分布）に従って観測事象が表れ、かつ定常状態間は時間経過と共にマルコフ性を持つ確率過程で遷移する」という考え方をする。
マルコフ性とは
定常状態（上の例ではA・B・C）が時刻\(t = 1\)～\(n\)に、時間的に連続した系列\( Q:=\{q_1,q_2, \cdot \cdot \cdot, q_n\}\)となる場合に、観測事象が系列\( Y:=\{y_1,y_2, \cdot \cdot \cdot, y_n\}\)となる確率 \(P(Y|Q)\)は
\[P(Y|Q)=\pi_{q_0}\prod_{t=1}^{n} a_{q_t q_{t+1}}b_{q_t}(y_t)\]
ここで、
1つは、Baum-Welchアルゴリズムを用いて、上記定式化の各パラメータ\(\pi_{q_0}\)、\(a_{q_t q_{t+1}}\)、\(b_{q_t}(y_t)\)を推定すること。
もう１つは、Viterviアルゴリズムを用いて、観測事象の系列\( Y:=\{y_1,y_2, \cdot \cdot \cdot, y_n\}\)が観測される場合に、尤度関数を最大化する、定常状態系列\( Q:=\{q_1,q_2, \cdot \cdot \cdot, q_n\}\)を求めることだ。
このテストデータは、
の２つの定常状態があり、時間1～200と401～600は定常状態A、時間201～400と601～800は定常状態Bとなっているデータだ。
データ作成用のRコードは以下のとおり。
ここで、結果の内容を見ていく。
まずは、Distribution parametersの項。ここは先の\(b_{q_t}(y_t)\)の推定値、すなわち隠れた定常状態の確率分布について示されている。
・「State 1」がテストデータの定常状態B
・「State 2」がテストデータの定常状態A
を表わしている。
これらの結果に表示された平均と分散は、テストデータを用意した際のパラメータと良く一致している。
次にInitial probabilitiesの項。ここは先の\(\pi_{q_0}\)の推定値を示している。「Pi 2（=State2が初期状態として現れる確率）」がほぼ１であり、定常状態Aが最初に現れているテストデータと合致している。
最後にTransition matrixの項。これは先の\(a_{q_t q_{t+1}}\)の推定値を示している。
モデルの考え方について一番参考になったのは、ここ（pdf）と、ここ（pdf）と、ここだ。
▼隠れマルコフモデルの考え方
隠れマルコフモデルを一言でいうと、時間的に非定常な観測事象を、「（隠れた）複数の定常状態が、マルコフ性を持つ確率過程で遷移し、それぞれの定常状態の確率分布に従って観測事象が生起される」という考え方で記述しようとするモデル。といえる。文章で言っても伝わりづらいので、考え方のイメージを下の図に示す。
これをモデルに落とし込むため隠れマルコフモデルは「隠れた定常状態（上の例では３つ）が存在し、各定常状態の確率分布（図下部の青・赤・緑で示された確率分布）に従って観測事象が表れ、かつ定常状態間は時間経過と共にマルコフ性を持つ確率過程で遷移する」という考え方をする。
マルコフ性とは
「状態間の遷移確率は、過去の経緯とか関係なく、「現在がどの状態か？」のみで決定される」という性質である。例えば現在が定常状態Aであれば（過去にどういう状態遷移をしたに関わらず）状態Aから、Aに自己遷移する確率は何％、Bに遷移する確率は何％、Cに遷移する確率は何％というように決定されるという性質である。
▼隠れマルコフモデルを定式化する。
上記のような考え方を定式化すると、以下のようになる。定常状態（上の例ではA・B・C）が時刻\(t = 1\)～\(n\)に、時間的に連続した系列\( Q:=\{q_1,q_2, \cdot \cdot \cdot, q_n\}\)となる場合に、観測事象が系列\( Y:=\{y_1,y_2, \cdot \cdot \cdot, y_n\}\)となる確率 \(P(Y|Q)\)は
\[P(Y|Q)=\pi_{q_0}\prod_{t=1}^{n} a_{q_t q_{t+1}}b_{q_t}(y_t)\]
ここで、
- \(b_{q_t}(y_t)\)は、時刻\(t\)の隠れた定常状態\(q_t\)から観測事象\(y_t\)が現れる確率分布。（つまり上図での青や赤や緑で表わされた確率分布）
- \(a_{q_t q_{t+1}}\)は、時刻が\(t\)から\(t+1\)に進む際に、隠れた定常状態\(q_t\)から\(q_{t+1}\)に遷移する確率。
- \(\pi_{q_0}\)は、状態\(q_0\)が初期定常情報源となる確率。
▼隠れマルコフモデルで何が推定できるのか？
隠れマルコフモデルには主に２つのアルゴリズムを用いて以下の推定が可能である。1つは、Baum-Welchアルゴリズムを用いて、上記定式化の各パラメータ\(\pi_{q_0}\)、\(a_{q_t q_{t+1}}\)、\(b_{q_t}(y_t)\)を推定すること。
もう１つは、Viterviアルゴリズムを用いて、観測事象の系列\( Y:=\{y_1,y_2, \cdot \cdot \cdot, y_n\}\)が観測される場合に、尤度関数を最大化する、定常状態系列\( Q:=\{q_1,q_2, \cdot \cdot \cdot, q_n\}\)を求めることだ。
▼解析用のテストデータ作成
では、実際にRで隠れマルコフモデルの解析を行っていく。まずはじめに下図のような解析用のテストデータを作成しておく。このテストデータは、
- 定常状態A： 平均10、分散6の正規分布で観測事象が出現する状態
- 定常状態B： 平均12、分散7の正規分布で観測事象が出現する状態
の２つの定常状態があり、時間1～200と401～600は定常状態A、時間201～400と601～800は定常状態Bとなっているデータだ。
set.seed(1) #乱数シード作成
size <- 200 #データセットの大きさ
mean.a <- 10 #定常状態Aの平均
var.a <- 6 #定常状態Aの分散
mean.b <- 12 #定常状態Bの平均
var.b <- 7 #定常状態Bの分散
x.1 <- rnorm(size, mean.a, sqrt(var.a)) #データ作成
x.2 <- rnorm(size, mean.b, sqrt(var.b)) #データ作成
x.3 <- rnorm(size, mean.a, sqrt(var.a)) #データ作成
x.4 <- rnorm(size, mean.b, sqrt(var.b)) #データ作成
x <- c(x.1, x.2, x.3, x.4) #データ連結
plot(x, xlab = "time", ylab = "observed quantity")
size <- 200 #データセットの大きさ
mean.a <- 10 #定常状態Aの平均
var.a <- 6 #定常状態Aの分散
mean.b <- 12 #定常状態Bの平均
var.b <- 7 #定常状態Bの分散
x.1 <- rnorm(size, mean.a, sqrt(var.a)) #データ作成
x.2 <- rnorm(size, mean.b, sqrt(var.b)) #データ作成
x.3 <- rnorm(size, mean.a, sqrt(var.a)) #データ作成
x.4 <- rnorm(size, mean.b, sqrt(var.b)) #データ作成
x <- c(x.1, x.2, x.3, x.4) #データ連結
plot(x, xlab = "time", ylab = "observed quantity")
▼Baum-Welchアルゴリズム
実際にBaum-Welchアルゴリズムを用いて、上記定式化の各パラメータ\(\pi_{q_0}\)、\(a_{q_t q_{t+1}}\)、\(b_{q_t}(y_t)\)を推定する。
Rには隠れマルコフ解析用のパッケージ「RHmm」があるので、それを使えば簡単に推定できる（パッケージは前もってインストールしとく必要あり）。実際の解析用のコードは以下のとおり。
> library("RHmm") #隠れマルコフモデルのライブラリ読み込み
> hmm.fitted <- HMMFit(x, nStates = 2) #フィッティング（ステート数は2を設定）
> print(hmm.fitted$HMM) #Baum-Welchアルゴリズムでの推定量を表示
> hmm.fitted <- HMMFit(x, nStates = 2) #フィッティング（ステート数は2を設定）
> print(hmm.fitted$HMM) #Baum-Welchアルゴリズムでの推定量を表示
まずRHmmパッケージを読み込み、先のテストデータ「x」、定常状態数は２（今回はAとBの２つだから）を指定してHMMFit関数を実行する。この関数の戻り値はBaum-Welchアルゴリズムの解析結果も含むので、それをプリントする・・・だけ。
実行結果は以下のとおり。
Initial probabilities:
Pi 1 Pi 2
3.50879e-35 1
Transition matrix:
State 1 State 2
State 1 0.997176507 0.002823493
State 2 0.005380641 0.994619359
Conditionnal distribution parameters:
Distribution parameters:
mean var
State 1 11.89274 7.75078
State 2 9.99124 6.01621
Pi 1 Pi 2
3.50879e-35 1
Transition matrix:
State 1 State 2
State 1 0.997176507 0.002823493
State 2 0.005380641 0.994619359
Conditionnal distribution parameters:
Distribution parameters:
mean var
State 1 11.89274 7.75078
State 2 9.99124 6.01621
ここで、結果の内容を見ていく。
まずは、Distribution parametersの項。ここは先の\(b_{q_t}(y_t)\)の推定値、すなわち隠れた定常状態の確率分布について示されている。
・「State 1」がテストデータの定常状態B
・「State 2」がテストデータの定常状態A
を表わしている。
これらの結果に表示された平均と分散は、テストデータを用意した際のパラメータと良く一致している。
次にInitial probabilitiesの項。ここは先の\(\pi_{q_0}\)の推定値を示している。「Pi 2（=State2が初期状態として現れる確率）」がほぼ１であり、定常状態Aが最初に現れているテストデータと合致している。
最後にTransition matrixの項。これは先の\(a_{q_t q_{t+1}}\)の推定値を示している。
- 定常状態Aの時、定常状態Aに（自己）遷移する確率は0.997程度。定常状態Bに遷移する確率は0.003程度。
- 定常状態Bの時、定常状態Bに（自己）遷移する確率は0.995程度。定常状態Aに遷移する確率は0.005程度。
▼Viterviアルゴリズム
次にViterviアルゴリズムを用いて、テストデータのような観測事象が現れたとき、最も尤（もっと）もらしい（尤度関数を最大化する）、定常状態の系列（ここでは状態Aと状態Bの時間系列）を求めてみる。
解析用のRコードは以下のとおり。
> #先ほどのHMMFitの戻り値とテストデータXからViterbiアルゴリズムフィッティング
> vitervi.fitted = viterbi(hmm.fitted,x)
> #Viterviアルゴリズムで求めた、最も尤もらしい状態の系列をプロットする。
> plot(vitervi.fitted$states, xlab = "time", ylab = "Most Likely Status ID")
> vitervi.fitted = viterbi(hmm.fitted,x)
> #Viterviアルゴリズムで求めた、最も尤もらしい状態の系列をプロットする。
> plot(vitervi.fitted$states, xlab = "time", ylab = "Most Likely Status ID")
結果は、次のグラフで、テストデータ作成時に設定した状態遷移の時系列とよく合致している結果が得られている。
以上。
