MATLAB で Viterbi アルゴリズムを実装する
先日、第8回「続・わかりやすいパターン認識」読書会 に参加しました。第 8 章の隠れマルコフモデルについて解説していただきました。そこで今回は、教科書に載っている Viterbi アルゴリズムを実装して試してみます。MATLAB で実装してみました。
遷移確率と出力確率は、教科書の (8.116) 式にならって以下のように定義します。A が遷移確率、B が出力確率で、行列の大きさから分かるように 3 つの状態と 2 つの出力記号を持つモデルになります。初期状態 rho も教科書の (8.117) 式のとおりに定めます。
>> A = [ 0.1 0.7 0.2 0.2 0.1 0.7 0.7 0.2 0.1 ]; >> B = [ 0.9 0.1 0.6 0.4 0.1 0.9 ]; >> rho = [1 1 1] / 3;
まず、これらのパラメータにしたがう系列を生成します。以下のプログラムで生成できます。引数の A, B, rho は先に定めた各パラメータを表し、n で系列長を指定します。戻り値は、s が状態系列、x が観測記号系列です。
function [s, x] = GenerateSample(A, B, rho, n) s(1) = sample(rho); x(1) = sample(B(s(1), :)); for t = 2:n s(t) = sample(A(s(t-1), :)); x(t) = sample(B(s(t), :)); end end function i = sample(p) i = min(find(rand <= cumsum(p) / sum(p))); end
この関数を次のように実行すれば、s と x を生成できます。
>> [s, x] = GenerateSample(A, B, rho, 100);
生成したデータの様子を眺めてみます。まず、s と x のそれぞれについて、histcounts 関数で度数分布を調べます。
>> histcounts(s) ans = 36 28 36 >> histcounts(x) ans = 57 43
Statistics Toolbox の crosstab 関数を用いると 2 変量のクロス集計を行えます。以下では s(t-1) と s(t) の関係、s(t) と x(t) の関係を見てみました。
>> crosstab(s(1:end-1), s(2:end)) ans = 5 22 9 3 1 24 28 5 2 >> crosstab(s(1:end-1), s(2:end)) ans = 5 22 9 3 1 24 28 5 2
Viterbi アルゴリズムの実装は以下のようになりました。行末のコメントは教科書の数式番号に対応します。確率は log にしてアンダーフローを避けています。また、argmax 関数は MATLAB に用意されていないので自作しました。
function s = Viterbi(A, B, rho, x) n = length(x); % Step 1 初期化 P(:, 1) = log(rho') + log(B(:, x(1))); % (8.28) I(:, 1) = zeros(length(rho), 1); % (8.29) % Step 2 再帰的計算 for t = 2:n Z = bsxfun(@plus, P(:, t-1), log(A)); I(:, t) = argmax(Z); % (8.31) P(:, t) = diag(Z(I(:, t), :)) + log(B(:, x(t))); % (8.30) end % Step 3 終了 s(n) = min(find(P(:, n) == max(P(:, n)))); % (8.33) % Step 4 系列の復元 for t = n-1:-1:1 s(t) = I(s(t+1), t+1); % (8.34) end end function I = argmax(Z) [R, C] = find(bsxfun(@eq, Z, max(Z))); I = R(find(diff([0; C]))); end
実装上の細かな点として、確率が最大になる状態が複数個存在する場合に、それらの中から一つを選ぶ処理が必要でした。上記のプログラムでは、(8.33) 式の行で min を取っている部分、argmax 関数で diff を取っている部分です。この実装は、単純にインデックスが最小の状態を選択するものです。
このプログラムを用いて、次のように状態系列を推定できます。推定された系列と真の系列を比較すると、正解率は 77% になりました。
>> estimated = Viterbi(A, B, rho, x); >> mean(estimated == s) ans = 0.7700
このようにして Viterbi アルゴリズムを実装できましたが、MATLAB の Statistics Toolbox には隠れマルコフモデルを扱う関数群が用意されています。Viterbi アルゴリズムによる状態系列の推定は hmmviterbi 関数で実行できます。
隠れマルコフ モデルの最も可能性の高い状態パス - MATLAB hmmviterbi - MathWorks 日本
これらの関数では、初期状態は 1 番目の状態に固定されており、rho のような指定はできません。これは、下記のウェブページの末尾に「初期状態分布の変更」として説明されているように、仮想的な初期状態を考えることで対応するようです。
隠れマルコフ モデル (HMM) - MATLAB & Simulink - MathWorks 日本
この方法にしたがって、hmmviterbi を実行してみます。まず、元の A, B, rho から A_hat, B_hat を作成します。得られる行列は以下のとおりです。
>> A_hat = [0 rho; zeros(size(A, 1), 1) A] A_hat = 0 0.3333 0.3333 0.3333 0 0.1000 0.7000 0.2000 0 0.2000 0.1000 0.7000 0 0.7000 0.2000 0.1000 >> B_hat = [zeros(1, size(B, 2)); B] B_hat = 0 0 0.9000 0.1000 0.6000 0.4000 0.1000 0.9000
次のように hmmviterbi 関数を実行します。得られる状態系列は B_hat の行番号に対応するため、関数の結果から 1 を引いています。これを自作の関数による結果と比較したところ、ちゃんと結果が一致しました。
>> estimated2 = hmmviterbi(x, A_hat, B_hat) - 1; >> all(estimated == estimated2) ans = 1
