出力から状態変数を予測する
系列 $x^n$ が観測されている場合に、状態変数の系列 $s^n$ の値を予測する方法を考えます。最尤法の考え方に従って条件付き確率を最大化するのが自然な考え方です。
(s^n)^* = \mathrm{arg} \max_{s^n} \mathbb{P}(s^n | x^n)
この条件付き確率の計算は複雑になってしまうので、同時確率に直してしまったほうが扱いやすいです。
(s^n)^* = \mathrm{arg} \max_{s^n} \mathbb{P}(s^n | x^n)
= \mathrm{arg} \max_{s^n} \mathbb{P}(s^n, x^n) / \mathbb{P}(x^n) = \mathrm{arg} \max_{s^n} \mathbb{P}(s^n, x^n)
最後の等式は、$\mathbb{P}(x^n)$ による除算が、$s^n$ に依存しないことによります。
このまま、あらゆる列 $s^n$ について確率を評価し、尤度が最大となるものを選ぶとすると計算量が指数的に増大してしまうので、最大化問題を動的計画法の考え方を用いて分解します。
ビタービ・アルゴリズム
具体的には、次のような分解を考えます。任意の $t$ について、
\mathbb{P}(s^t, x^t) = \mathbb{P}(s_t, x_t | s^{t-1}, x^{t-1}) \cdot \mathbb{P}(s^{t-1}, x^{t-1}) = \mathbb{P}(s_t, x_t | s_{t-1}) \cdot \mathbb{P}(s^{t-1}, x^{t-1})
この両辺を $s^{t-1}$ について最大化すると、
\begin{align*}
\max_{s^{t-1}} \mathbb{P}(s^t, x^t) &= \max_{s^{t-1}} \left[ \mathbb{P}(s_t, x_t | s_{t-1}) \cdot \mathbb{P}(s^{t-1}, x^{t-1}) \right] \\
&= \max_{s_{t-1}} \left[ \mathbb{P}(s_t, x_t | s_{t-1}) \cdot \max_{s^{t-2}} \mathbb{P}(s^{t-1}, x^{t-1}) \right]
\end{align*}
そこで、$\psi_t(s_t) \equiv \max_{s^{t-1}} \mathbb{P}(s^t, x^t)$ と定義します。関数の引数が $s_t$ だけなのは、$t-1$ 期までの状態については最適値を取るように選ばれているためです。すると、次のような漸化式が得られます。
\psi_t(s_t)
= \max_{s_{t-1}} \mathbb{P}(s_t, x_t | s_{t-1}) \cdot \psi_{t-1}(s_{t-1})
= b(s_t, x_t) \cdot \max_{s_{t-1}} \left[ a(s_{t-1}, s_t) \cdot \psi_{t-1} (s_{t-1}) \right]
$\psi_t$ は動的計画法の枠組みでは"value function"と呼ばれます。対応する"policy function"にあたる関数を、次のように定義しておきます。
\Psi_t(s_t) \equiv \mathrm{arg} \max_{s_{t-1}} \left[ a(s_{t-1}, s_t) \cdot \psi_{t-1}(s_{t-1}) \right]
この関数が、$s_t$ を与えられた時に $s_{t-1}$ を返す関数であること、つまり未来から過去へ遡る関数である点に注意します。これは、最適「値」の計算には関係ありませんが、後に最適「解」を得るために必要になります。
仮に、$\psi_t, \Psi_t$ が任意の $t$ について求められたとします。$t=n$ の場合を考えて、
\psi_n(s_n)
= \max_{s^{n-1}} \mathbb{P}(s^n, x^n)
両辺を $s_n$ について最大化すると、
\max_{s_n} \psi_n(s_n)
= \max_{s_n} \max_{s^{n-1}} \mathbb{P}(s^n, x^n) = \max_{s^n} \mathbb{P}(s^n, x^n)
となるので、題意の最適化問題に帰着します。左辺の問題を解くことによって、最終期の最適解 $s_{n}^{*}$ が求まります。ここから、$\Psi_t$ 関数を利用して、各期の解を求めることができます。
s_{t}^{*} = \Psi_{t+1}(s_{t+1}^{*}) \;\;\;\; \text{for $t = n-1, n-2, \dots, 1$}
したがって、$\psi_t, \Psi_t$ の両関数を求めることで、題意の問題を解くことができるわけですが、これらの関数は逐次的な手順によって求めることができます。まず、初期条件は次のように定まります。
\psi_1(s_1) = \max_{s^{0}} \mathbb{P}(s^1, x^1) = \rho(s_1) b(s_1, x_1)
便宜上$s^0$ という変数を使ってはいますが特に意味はなく、第1期の状態変数の確率分布はすでに与えられています。
これを、上で求めた漸化式を逐次当てはめていけば全ての解が得られます。
\begin{align*}
\psi_t(s_t)
&= b(s_t, x_t) \cdot \max_{s_{t-1}} \left[ a(s_{t-1}, s_t) \cdot \psi_{t-1} (s_{t-1}) \right] \\
\Psi_t(s_t)
&= \mathrm{arg} \max_{s_{t-1}} \left[ a(s_{t-1}, s_t) \cdot \psi_{t-1} (s_{t-1}) \right]
\end{align*}
計算上は、$\psi_t, \Psi_t$ はいずれもベクトルで表現することができます。このうち、$\Psi_t$ については、最後に系列を復元する際に利用するので、どこかに格納しておく必要があります。
シミュレーション
次の関数は、総当り法によって最適な状態変数の系列とその尤度を見つけます。
S1 <- function(x, A, B, rho) {
N <- length(x) # size of data
J <- nrow(B) # number of possible state
maxP <- -Inf
amax <- NA
s <- rep_len(1, N)
v <- 0
while (v < J^N) {
prob_s <- 1 # P(s)
prob_xs <- 1 # P(x|s)
for (n in 1:N) {
if (n == 1) {
prob_s <- rho[s[1]]
} else {
prob_s <- prob_s * A[s[n-1], s[n]]
}
prob_xs <- prob_xs * B[s[n], x[n]]
}
p <- prob_s * prob_xs
if (p > maxP) {
maxP <- p
amax <- s
}
# next s
s[1] <- s[1] + 1
for (n in 2:N)
if (s[n-1] > J) {
s[n] <- s[n] + 1
s[n-1] <- 1
}
v <- v + 1
}
return(list(s = amax, prob = maxP))
}
次の関数はビタービ・アルゴリズムによって最適な状態系列とその尤度を計算します。
S2 <- function(x, A, B, rho) {
# Viterbi algorithm
N <- length(x) # size of data
J <- nrow(B) # number of possible state
value.fun <- rho * B[, x[1]] # the previous value function
policy.fun <- list() # keep the full history
for (n in 2:N) {
v <- diag(value.fun) %*% A
pol <- apply(v, 2, which.max)
policy.fun <- c(policy.fun, list(pol))
value.fun <- apply(v, 2, max) * B[, x[n]]
}
# compute the final value and the state
maxP <- max(value.fun)
s <- which.max(value.fun)
# recover the path
for (n in (N-1):1) {
s <- c(policy.fun[[n]][s[1]], s)
}
return(list(s = s, prob = maxP))
}
いずれも、本と同じ結果となります（171頁）。
x <- c(1, 2, 1)
S1(x, A, B, rho)
S2(x, A, B, rho)
#$s
#[1] 2 3 1
#$prob
#[1] 0.07938
ビタービ・アルゴリズムの方が断然速いです。
library(microbenchmark)
x <- floor(runif(7)*ncol(B))+1
microbenchmark(S1(x, A, B, rho), S2(x, A, B, rho))
#Unit: microseconds
# expr min lq mean median uq max neval
# S1(x, A, B, rho) 65249.613 67735.8585 71777.8841 71294.1740 73889.037 113615.508 100
# S2(x, A, B, rho) 460.986 497.1215 645.5264 573.2395 613.223 2314.339 100