読者です 読者をやめる 読者になる 読者になる

y_uti のブログ

統計、機械学習、自然言語処理などに興味を持つエンジニアの技術ブログです

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