音声認識のシステムは、あなたが発した「おはよう」という音波のスペクトルから、背後にある「お」「は」「よ」「う」という音素の系列を当てに行きます。遺伝子解析のソフトウェアは、ATGCの塩基配列から、CpGアイランドのような機能領域がどこにあるかを推定します。株価アナリストは、日々のリターン系列から、市場が「強気」「弱気」のどちらのレジームにあるかを読み取ろうとします。
これらに共通する構図は単純です。観測できる系列 $x_{1:T}$ の背後に、観測できない離散的な状態 $z_{1:T}$ が存在し、状態が時間的に遷移しながら観測を生成している。この状況をもっとも素直に確率モデル化したものが、隠れマルコフモデル (Hidden Markov Model, HMM) です。
HMMの威力は、状態空間が有限個の離散値（たとえば $N=10$ 個）であれば、観測系列の確率や状態系列の最尤推定を、動的計画法によって $O(N^2 T)$ という劇的な効率で計算できるところにあります。素朴に全状態系列を列挙すると $N^T$ 通り（$T=100, N=10$ なら $10^{100}$ 通り）になり、宇宙の年齢でも終わりません。Forward-Backwardアルゴリズムは、この指数爆発を多項式時間に圧縮する魔法のような手続きです。
応用は今でも幅広く息づいています。スマートスピーカーの音素デコーダ、生物情報学のCpGアイランド／遺伝子予測、英語の品詞タグ付け、金融時系列のレジームスイッチ検出、人間の活動認識（ウェアラブルのセンサーから「歩行・走行・静止」を推定）。線形ガウス連続状態に拡張すればカルマンフィルタになり、深層化すれば現代のニューラル状態空間モデル（Mamba 等）に繋がります。
本記事の内容
- HMMの数学的定義（遷移行列、出力分布、初期分布）と、関連する3つの基本問題
- Forward再帰 $\alpha_t(j)$ と Backward再帰 $\beta_t(i)$ の導出と $O(N^2 T)$ の計算量
- 周辺事後確率 $\gamma_t(i) = p(z_t=i \mid x_{1:T})$ による個別最大事後確率推定（MPM）
- Viterbi再帰による全体最尤系列復号と、MPMとの差異
- Baum-Welch（EM）による $A, B, \pi$ の学習
- Pythonでのスクラッチ実装（log-domainによるunderflow対策）と hmmlearn による株価レジーム検出のデモ
隠れ状態 $z_t$ が横方向にチェーン（マルコフ連鎖）で繋がり、各時刻から下向きに観測 $x_t$ が生成される「コーム（櫛）」の形がHMMの全体像です。状態間には直接の矢印がなく、観測どうしの相関は隠れ状態を介してのみ生まれます。この疎な依存構造が、後に見る動的計画法の鍵になります。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
直感 — 隠れた状態を観測から推測する
部屋に閉じこもった友人と電話で会話していると想像してください。友人の今日の気分（「機嫌が良い」「機嫌が悪い」）は直接見えません。しかし話題の内容（「冗談を言う」「愚痴を言う」「無口」）からなんとなく察しがつきます。「3日連続で愚痴ばかりだから、今は機嫌が悪いモードが続いているな」とか、「今日は冗談が混じってきたから、そろそろ持ち直したかな」といった具合に。
ここに2つの仮定を入れると、これがそのままHMMになります。1つ目は マルコフ性: 今日の気分は昨日の気分にだけ依存し、3日前の気分は（昨日の気分が決まれば）もう関係ない。2つ目は 観測の条件付き独立性: 今日話す内容は、今日の気分だけで決まる。気分が同じなら、過去の発言や明日の発言とは独立に分布する。
この2つの仮定だけで、HMMはちゃんと動きます。気分の遷移確率 $A = (A_{ij})$、気分から発言が出る確率 $B_j(x) = p(x \mid z=j)$、初期気分の確率 $\pi$ さえ与えてやれば、観測された発言系列から「機嫌の良し悪し」の事後分布が計算できる。しかも、後で見るように $O(N^2 T)$ という現実的な計算量で。
直感はここまでにして、これを正確に数式で書き下しましょう。
HMMの数学的定義
HMMは、離散的な隠れ状態のマルコフ連鎖と、各時刻の観測がそのときの状態にだけ依存して生成されるという2層構造を持つ確率モデルです。
時刻を $t = 1, 2, \dots, T$、隠れ状態を $z_t \in \{1, 2, \dots, N\}$、観測を $x_t$（離散でも連続でもよい）とします。HMMは次の3つのパラメータで決まります。
- 初期状態分布 $\pi_i = p(z_1 = i)$、$\sum_i \pi_i = 1$
- 状態遷移確率 $A_{ij} = p(z_{t+1} = j \mid z_t = i)$、各行 $\sum_j A_{ij} = 1$
- 出力分布 $B_j(x) = p(x_t = x \mid z_t = j)$
これらをまとめて $\lambda = (\pi, A, B)$ と書きます。観測が離散値（語彙サイズ $K$）なら $B$ は $N \times K$ の行列、連続値ならガウス分布などのパラメトリックモデル（Gaussian HMM）にします。
同時分布
HMMの仮定（マルコフ性と観測の条件付き独立性）から、状態系列 $z_{1:T}$ と観測系列 $x_{1:T}$ の同時分布は次のように因子分解されます。
$$ p(z_{1:T}, x_{1:T}) = p(z_1) \prod_{t=2}^{T} p(z_t \mid z_{t-1}) \prod_{t=1}^{T} p(x_t \mid z_t) $$
これを記号で書き直すと、
$$ p(z_{1:T}, x_{1:T}) = \pi_{z_1} B_{z_1}(x_1) \prod_{t=2}^{T} A_{z_{t-1} z_t} B_{z_t}(x_t) $$
になります。式が示しているのは、「初期状態に入る確率 $\pi_{z_1}$、最初の観測を出す確率 $B_{z_1}(x_1)$、それ以降は各ステップで遷移確率と出力確率を掛けていく」というシンプルなパスの確率です。
グラフィカルモデル表現
HMMをグラフ表記すると、隠れ状態が横方向にチェーンで繋がり、各状態から下向きに観測ノードが伸びる「コーム（櫛）」の形になります。
z_1 → z_2 → z_3 → ... → z_T
↓ ↓ ↓ ↓
x_1 x_2 x_3 x_T
この単純さがHMMの強みです。観測 $x_t$ どうしの間には直接のエッジがなく、相関は隠れ状態を介してのみ生まれる。これが後で動的計画法を可能にする鍵になります。
定義が揃ったところで、HMMで実際に「やりたい計算」が何なのかを整理しましょう。
3つの基本問題
HMMを使う場面は、Rabinerの古典的なチュートリアル以来、次の3つの問題に集約されてきました。
問題1: 評価 (Evaluation)
モデル $\lambda$ と観測系列 $x_{1:T}$ が与えられたとき、その尤度
$$ p(x_{1:T} \mid \lambda) = \sum_{z_{1:T}} p(z_{1:T}, x_{1:T} \mid \lambda) $$
を計算する問題です。応用としては、複数のモデル $\lambda^{(1)}, \lambda^{(2)}, \dots$ のうち観測を最もよく説明するモデルを選ぶ（音声認識で単語モデルを比較する、など）。
素朴に和を取ると状態系列が $N^T$ 通りあるため、$T=100, N=10$ で $10^{100}$ 通り。これを多項式時間に落とすのが Forwardアルゴリズム です。
問題2: 復号 (Decoding)
観測 $x_{1:T}$ から、もっともらしい隠れ状態系列を推定する問題です。「もっともらしい」には2通りの解釈があります。
- 個別最大事後確率推定 (Posterior Marginals / MPM): 各時刻 $t$ について $\arg\max_i p(z_t = i \mid x_{1:T})$
- 全体最尤系列 (Joint MAP, Viterbi): $\arg\max_{z_{1:T}} p(z_{1:T} \mid x_{1:T})$
前者は Forward-Backwardアルゴリズム、後者は Viterbiアルゴリズム で計算されます。両者は一般に 異なる解 を返します（後ほど具体例で見ます）。
問題3: 学習 (Learning)
観測 $x_{1:T}$（または複数の系列）から、パラメータ $\lambda = (\pi, A, B)$ を推定する問題です。隠れ状態が未観測なので、最尤推定はEMアルゴリズムの形を取ります。HMMに特化したEMが Baum-Welchアルゴリズム です。
3つの問題は独立ではなく、内側で同じ Forward/Backward を共有します。まずその中核を導出しましょう。
Forward-Backwardアルゴリズム
Forward再帰 $\alpha_t(j)$
Forward変数を次のように定義します。
$$ \alpha_t(j) = p(x_{1:t}, z_t = j \mid \lambda) $$
「時刻 $t$ までの観測を生成しつつ、時刻 $t$ で状態 $j$ にいる」同時確率です。同時尤度に状態を含めて持つことで、再帰がきれいに書けるのがポイントです。
初期化:
$$ \alpha_1(j) = \pi_j B_j(x_1), \quad j = 1, \dots, N $$
再帰: $\alpha_{t+1}(j)$ を $\alpha_t$ から作りたい。確率の周辺化と条件付け、そしてHMMの条件付き独立性を順に使います。
$$ \alpha_{t+1}(j) = p(x_{1:t+1}, z_{t+1} = j) = \sum_{i=1}^{N} p(x_{1:t+1}, z_t = i, z_{t+1} = j) $$
まず $z_t$ について周辺化しました。次に同時分布を分解します。$x_{t+1}$ は $z_{t+1}$ にだけ依存し、$z_{t+1}$ は $z_t$ にだけ依存することを使うと、
$$ p(x_{1:t+1}, z_t = i, z_{t+1} = j) = p(x_{1:t}, z_t = i) \cdot p(z_{t+1} = j \mid z_t = i) \cdot p(x_{t+1} \mid z_{t+1} = j) $$
つまり $\alpha_t(i) \cdot A_{ij} \cdot B_j(x_{t+1})$ です。$\sum_i$ に戻すと、
$$ \boxed{\,\alpha_{t+1}(j) = \left[ \sum_{i=1}^{N} \alpha_t(i) A_{ij} \right] B_j(x_{t+1})\,} $$
が得られます。終端で全状態について和を取ると、欲しい尤度が手に入ります。
$$ p(x_{1:T} \mid \lambda) = \sum_{j=1}^{N} \alpha_T(j) $$
各時刻で $N$ 個の $\alpha_{t+1}(j)$ を作るのに $N$ 個の積和、それを $T$ 回繰り返すので、計算量は $O(N^2 T)$。$N=10, T=100$ なら $10^4$ 回で、素朴な $10^{100}$ から劇的に小さくなりました。動的計画法が指数を多項式に潰した、これがForwardの本質です。
トレリス図は「時刻（横軸）×状態（縦軸）」の格子で、各ノードに $\alpha_t(j)$ の値を持ちます。ノードの濃さが $\alpha$ の大きさに対応しており、時刻が進むほど一部の状態に確率が集中していく様子が見て取れます。赤矢印が示す再帰の向きは「過去→現在」であり、前の時刻のすべての状態から確率を集めて次の状態の $\alpha$ を計算する様子がわかります。
Backward再帰 $\beta_t(i)$
Forwardが「過去を畳んで現在に運ぶ」のに対し、Backwardは「未来を畳んで現在に運びます」。Backward変数を次で定義します。
$$ \beta_t(i) = p(x_{t+1:T} \mid z_t = i, \lambda) $$
「時刻 $t$ で状態 $i$ にいるという条件のもとで、それ以降の観測 $x_{t+1:T}$ が出る確率」です。Forwardの $\alpha$ が 同時確率 だったのに対し、$\beta$ は 条件付き確率 であることに注意してください（こうしておくと後の式が綺麗になります）。
初期化 は終端から始めるので、
$$ \beta_T(i) = 1, \quad i = 1, \dots, N $$
とします（時刻 $T$ より先には観測がないので、確率1で空積）。
再帰: $\beta_t(i)$ を $\beta_{t+1}$ から作ります。$z_{t+1}$ について周辺化し、HMMの分解を使います。
$$ \beta_t(i) = p(x_{t+1:T} \mid z_t = i) = \sum_{j=1}^{N} p(x_{t+1:T}, z_{t+1} = j \mid z_t = i) $$
同時条件付き確率を分解します。$z_t$ から $z_{t+1}$ への遷移、$z_{t+1}$ からの観測 $x_{t+1}$、その後の観測 $x_{t+2:T}$ は $z_{t+1}$ で条件付ければ $z_t$ と独立です。
$$ p(x_{t+1:T}, z_{t+1} = j \mid z_t = i) = A_{ij} \cdot B_j(x_{t+1}) \cdot \beta_{t+1}(j) $$
まとめると、
$$ \boxed{\,\beta_t(i) = \sum_{j=1}^{N} A_{ij} B_j(x_{t+1}) \beta_{t+1}(j)\,} $$
これも各時刻 $O(N^2)$、全体で $O(N^2 T)$ です。
Backwardは終端 $\beta_T(i)=1$（オレンジの「境界条件」ボックス）から左方向に伝播します。ノードの濃さが $\beta$ の大きさに対応しており、終端に近いほど値が1に近く（確実）、系列の先頭に向かうほど観測の不確実性が混じって値が散らばります。この「未来からの情報」と、先ほどの「過去からの情報（$\alpha$）」を掛け合わせると、初めて「全系列を考慮した事後確率」が得られます。
周辺事後確率 $\gamma_t(i)$
Forward $\alpha$ と Backward $\beta$ を組み合わせると、各時刻の隠れ状態の事後周辺確率 が一発で計算できます。これが Forward-Backward の主な目的です。
ベイズの定理から、
$$ \gamma_t(i) = p(z_t = i \mid x_{1:T}) = \frac{p(z_t = i, x_{1:T})}{p(x_{1:T})} $$
分子を $\alpha, \beta$ で書き直しましょう。$x_{1:T}$ を $x_{1:t}$ と $x_{t+1:T}$ に分け、条件付き独立性「$x_{t+1:T}$ は $z_t$ で条件付ければ $x_{1:t}$ と独立」を使うと、
$$ p(z_t = i, x_{1:T}) = p(x_{1:t}, z_t = i) \cdot p(x_{t+1:T} \mid z_t = i) = \alpha_t(i) \beta_t(i) $$
よって、
$$ \boxed{\,\gamma_t(i) = \frac{\alpha_t(i)\beta_t(i)}{\sum_{j} \alpha_t(j)\beta_t(j)}\,} $$
分母の和は $t$ に依らず $p(x_{1:T})$ に等しいので、正規化は任意の $t$ で構いません。
左パネル（青）は前向きの $\alpha_t$ で、時刻が進むにつれて特定の状態に収束していく様子が読み取れます。中パネル（橙）は後ろ向きの $\beta_t$ で、終端から逆方向に伝播するため左側ほど値が安定します。右パネル（緑）が両者の積を正規化した $\gamma_t$ であり、過去と未来の情報が融合して初めて「その時刻の状態の確率」が得られます。単なる $\alpha$ だけより $\gamma$ の方が過去・未来両方の観測を活かすため、精度の高い事後推定になっています。
この $\gamma_t(i)$ を最大化する $i$ を各時刻独立に選ぶのが 個別最大事後確率推定 (MPM) です。「各時刻ごとに、その時刻の状態として最ももっともらしいものを選ぶ」素直な戦略です。ただし、MPMには注意点があります。
MPMの落とし穴と Viterbi の必要性
MPMは時刻ごとに独立に $\arg\max$ を取るので、得られた系列 $(\hat z_1, \hat z_2, \dots, \hat z_T)$ が 遷移確率 $A$ で許される系列とは限らない という弱点があります。たとえば $A_{12} = 0$（状態1から状態2へ直接は遷移できない）であっても、MPMが $\hat z_t = 1, \hat z_{t+1} = 2$ を返すことがあり得ます。系列の「滑らかさ」を保証したいなら、結合分布を最大化する Viterbiアルゴリズム が必要です。
ここで「全体としてもっとも整合的な状態系列」を求める道具に切り替えましょう。
Viterbiアルゴリズム
Viterbiは、$\arg\max_{z_{1:T}} p(z_{1:T} \mid x_{1:T})$ を直接求めます。分母 $p(x_{1:T})$ は $z$ に依らないので、$\arg\max_{z_{1:T}} p(z_{1:T}, x_{1:T})$ と同値です。
Forward再帰の $\sum$ を $\max$ に置き換えるのがアイデアです。
Viterbi変数:
$$ \delta_t(j) = \max_{z_1, \dots, z_{t-1}} p(z_1, \dots, z_{t-1}, z_t = j, x_{1:t}) $$
「時刻 $t$ で状態 $j$ にいる経路のうち、最尤のもの」の確率です。
初期化:
$$ \delta_1(j) = \pi_j B_j(x_1), \quad \psi_1(j) = 0 $$
ここで $\psi_t(j)$ は最尤経路の バックポインタ（「時刻 $t$ で状態 $j$ に来る直前の状態」）です。
再帰: 最大化はパスごとに独立に評価できる（マルコフ性）ので、
$$ \boxed{\,\delta_{t+1}(j) = \left[\max_{i} \delta_t(i) A_{ij}\right] B_j(x_{t+1}), \quad \psi_{t+1}(j) = \arg\max_{i} \delta_t(i) A_{ij}\,} $$
終端:
$$ p^* = \max_j \delta_T(j), \quad \hat z_T = \arg\max_j \delta_T(j) $$
バックトラック: 終端から $\psi$ を辿って最尤経路を復元します。
$$ \hat z_t = \psi_{t+1}(\hat z_{t+1}), \quad t = T-1, T-2, \dots, 1 $$
計算量は Forward と同じ $O(N^2 T)$。$\sum$ と $\max$ の違いだけなので、構造はそっくりです。両者の差は「平均化するか、最大を取るか」という根本的な意味の違いに対応しています。
Viterbiトレリスでは、各ノードの数値が $\delta_t(j)$（その経路を通る最尤確率）を表しています。時刻が進むに連れて「状態3（Loaded）」に経路が集中し、最終時刻で最大のスコアを持つ状態からバックポインタ $\psi$（赤矢印）を逆に辿ることで最尤経路（赤い強調ノード）を復元します。Forwardの $\sum$ を $\max$ に変えるだけで、動的計画法の骨格はまったく同じです。
MPM と Viterbi が一致しない例
簡単な例で2つの差を見ます。状態 $\{1, 2\}$、遷移行列が $A_{11} = 0.9, A_{12} = 0.1, A_{21} = 0.5, A_{22} = 0.5$ で、ある時刻 $t$ の周辺事後が $\gamma_t = (0.4, 0.6)$ だったとします。MPM は $\hat z_t = 2$ を選びます。しかし、次の時刻 $t+1$ の最尤状態が 1 で、$A_{12} = 0.1$ と小さいなら、Viterbi は「$z_t = 1$ から $z_{t+1} = 1$ へ滑らかに繋ぐ」全体最尤経路を選ぶ可能性があります。MPMは局所最尤、Viterbiは大域最尤 という違いをはっきり覚えておきましょう。
4段のパネルを見ると、最上段の観測系列（表/裏のランダムな列）からは背後の構造がまったく見えないことが確認できます。しかし3段目の事後確率 $p(\text{Loaded} \mid x_{1:T})$ は、真の状態（2段目）が Loaded のときに1に近い値、Fair のときに0に近い値を示しており、教師なしで隠れ状態を追跡できています。4段目の×印が「MPMとViterbiで判断が分かれる時刻」で、遷移の境界付近に集中しています。これは局所最尤（MPM）と大域最尤（Viterbi）の差が、状態遷移のタイミングで顕在化することを実際のデータで示しています。
ではこのモデルのパラメータ $\lambda = (\pi, A, B)$ をどうやって推定するか。これがHMMの最後のピース、Baum-Welch です。
Baum-Welch (EM学習)
観測 $x_{1:T}$ だけからパラメータを最尤推定したい。隠れ変数 $z_{1:T}$ があるので EMアルゴリズム の出番です。HMMに特化したEMが Baum-Welch です。
Eステップ — 期待値の計算
現在のパラメータ $\lambda^{\text{old}}$ のもとで、以下の事後量を計算します。
周辺事後 $\gamma_t(i) = p(z_t = i \mid x_{1:T}, \lambda^{\text{old}})$ は前節の通り $\gamma_t(i) = \alpha_t(i)\beta_t(i) / p(x_{1:T})$ で得ます。
ペア事後 $\xi_t(i, j) = p(z_t = i, z_{t+1} = j \mid x_{1:T}, \lambda^{\text{old}})$ を導きます。同時確率を分解すると、
$$ p(z_t = i, z_{t+1} = j, x_{1:T}) = \alpha_t(i) \cdot A_{ij} \cdot B_j(x_{t+1}) \cdot \beta_{t+1}(j) $$
なので、
$$ \xi_t(i, j) = \frac{\alpha_t(i) A_{ij} B_j(x_{t+1}) \beta_{t+1}(j)}{p(x_{1:T})} $$
となります。$t$ で和を取ると遷移の期待回数 $\sum_t \xi_t(i, j)$、$j$ で和を取ると周辺 $\gamma_t(i) = \sum_j \xi_t(i, j)$ に戻ります。
左パネルが $\xi_t(i, j)$ を構成する4つの要素を示しています。$\alpha_t(i)$（過去から状態 $i$ に到達する確率）、$A_{ij}$（遷移確率）、$B_j(x_{t+1})$（観測確率）、$\beta_{t+1}(j)$（未来の観測が状態 $j$ から生まれる確率）の積を全系列の尤度で正規化します。右パネルはこの $\xi$ がMステップでどのように使われるかを整理したもので、遷移確率・初期分布・出力分布の各更新式が「期待回数の比」として導かれることが確認できます。
Mステップ — パラメータの更新
$\gamma, \xi$ の重みつき頻度カウントとして、各パラメータを更新します。
初期分布:
$$ \pi_i^{\text{new}} = \gamma_1(i) $$
遷移確率: 「状態 $i$ から $j$ への遷移期待回数」を「状態 $i$ にいる期待回数（最終時刻を除く）」で割る:
$$ A_{ij}^{\text{new}} = \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{\sum_{t=1}^{T-1} \gamma_t(i)} $$
離散観測分布（観測 $x_t \in \{1, \dots, K\}$ の場合）: 「状態 $j$ で記号 $k$ を出す期待回数」を「状態 $j$ にいる期待回数」で割る:
$$ B_j(k)^{\text{new}} = \frac{\sum_{t: x_t = k} \gamma_t(j)}{\sum_{t=1}^{T} \gamma_t(j)} $$
Gaussian HMM（観測が連続で $x_t \mid z_t = j \sim \mathcal N(\mu_j, \Sigma_j)$）: 重みつきガウス推定:
$$ \mu_j^{\text{new}} = \frac{\sum_t \gamma_t(j) x_t}{\sum_t \gamma_t(j)}, \quad \Sigma_j^{\text{new}} = \frac{\sum_t \gamma_t(j) (x_t – \mu_j^{\text{new}})(x_t – \mu_j^{\text{new}})^\top}{\sum_t \gamma_t(j)} $$
これらの更新式は、いずれも 対数尤度の下界（ELBO） $Q(\lambda \mid \lambda^{\text{old}}) = \mathbb E_{z \mid x, \lambda^{\text{old}}}[\log p(z, x \mid \lambda)]$ を $\lambda$ について最大化 することで導かれます。詳細はガウス混合モデルとEMアルゴリズムと同じ流儀ですが、HMMでは「サンプル間の独立性」が崩れているため、$\gamma, \xi$ の計算に Forward-Backward が必要になります。
収束性と局所最適
EMの一般論として、Baum-Welchの反復ごとに観測尤度 $p(x_{1:T} \mid \lambda)$ は 単調非減少 です。しかし、得られる解は局所最適に留まる保証しかなく、初期値依存性が強いことが知られています。実用上は、複数の初期値で走らせて最尤のものを選ぶ、K-meansで初期化する（Gaussian HMMの場合）、などの工夫が使われます。
理論はここまでで一通り揃いました。次は手を動かして、HMMを実装してみましょう。スクラッチで Forward / Backward / Viterbi / Baum-Welch を書き、最後に hmmlearn で実データに適用します。
Python実装
スクラッチ実装 — log-domainでunderflow対策
長い系列を扱うと、$\alpha_t(j)$ や $\delta_t(j)$ は確率の積を重ねるため、急速に 0 に近づき underflow します。実装では 対数領域 で計算するのが定石です。scipy.special.logsumexp を使えば、$\log \sum_i \exp(x_i)$ を数値安定に計算できます。
import numpy as np
from scipy.special import logsumexp
class DiscreteHMM:
"""離散観測の隠れマルコフモデル
pi: (N,) 初期分布 / A: (N, N) 遷移行列 / B: (N, K) 観測分布
"""
def __init__(self, pi, A, B):
self.log_pi = np.log(pi + 1e-300)
self.log_A = np.log(A + 1e-300)
self.log_B = np.log(B + 1e-300)
self.N, self.K = B.shape
def log_emission(self, obs):
"""obs: (T,) の整数配列 → (T, N) の log B_j(x_t)"""
return self.log_B[:, obs].T
def forward(self, obs):
"""log alpha_t(j) を返す (T, N)"""
T = len(obs)
logB = self.log_emission(obs)
log_alpha = np.full((T, self.N), -np.inf)
log_alpha[0] = self.log_pi + logB[0]
for t in range(1, T):
# logsumexp_i (log_alpha[t-1, i] + log_A[i, j]) + log_B[j, x_t]
log_alpha[t] = logsumexp(log_alpha[t-1][:, None] + self.log_A, axis=0) + logB[t]
return log_alpha
def backward(self, obs):
"""log beta_t(i) を返す (T, N)"""
T = len(obs)
logB = self.log_emission(obs)
log_beta = np.full((T, self.N), -np.inf)
log_beta[T-1] = 0.0 # log 1
for t in range(T-2, -1, -1):
# logsumexp_j (log_A[i, j] + log_B[j, x_{t+1}] + log_beta[t+1, j])
log_beta[t] = logsumexp(self.log_A + (logB[t+1] + log_beta[t+1])[None, :], axis=1)
return log_beta
def log_likelihood(self, obs):
"""log p(x_{1:T})"""
return logsumexp(self.forward(obs)[-1])
def posterior(self, obs):
"""gamma_t(i) = p(z_t=i | x_{1:T}) 形状 (T, N)"""
log_alpha = self.forward(obs)
log_beta = self.backward(obs)
log_gamma = log_alpha + log_beta
log_gamma -= logsumexp(log_gamma, axis=1, keepdims=True)
return np.exp(log_gamma)
ここでのポイントは3つあります。第1に、すべての遷移・観測・$\alpha$・$\beta$ を log で持つ ことで、underflowを完全に避けています。第2に、logsumexp は内部で最大値をシフトしてから exp するので、$10^{-300}$ のような極端な値でも数値が破綻しません。第3に、log_alpha[t-1][:, None] + self.log_A のブロードキャストで $N \times N$ の行列を作り、軸方向に logsumexp を取るベクトル化を活用しています。素朴に二重ループで書くと数倍以上遅くなります。
左パネルは素朴な確率乗算 $p^t$（$p=0.4$）の値を対数スケールで示しており、$T$ が100を超えると浮動小数点のアンダーフロー限界（約 $10^{-300}$）に達して完全に0になります。$T=200$ では $0.4^{200} \approx 10^{-90}$ になるため、素朴実装では Forward の途中でアルファが全滅します。右パネルは log-domain（$t \cdot \log p$）で同じ情報を持つと安定して線形に減少し続けます。logsumexp を使えば「和の対数」も数値安定に計算でき、長系列（音声認識では $T$ が数千になる）でも破綻しません。
Viterbiも同様に書きます。$\sum$ を $\max$ に置き換え、logsumexp の代わりに np.max と np.argmax を使います。
def viterbi(self, obs):
"""最尤状態系列を返す"""
T = len(obs)
logB = self.log_emission(obs)
log_delta = np.full((T, self.N), -np.inf)
psi = np.zeros((T, self.N), dtype=int)
log_delta[0] = self.log_pi + logB[0]
for t in range(1, T):
# 各 j について max_i (log_delta[t-1, i] + log_A[i, j])
scores = log_delta[t-1][:, None] + self.log_A # (N, N)
psi[t] = np.argmax(scores, axis=0)
log_delta[t] = scores.max(axis=0) + logB[t]
# バックトラック
path = np.zeros(T, dtype=int)
path[T-1] = np.argmax(log_delta[T-1])
for t in range(T-2, -1, -1):
path[t] = psi[t+1, path[t+1]]
return path, log_delta[T-1].max()
DiscreteHMM.viterbi = viterbi
Viterbi の構造は Forward とそっくりで、logsumexp の代わりに max、$\arg\max$ の経路をバックポインタ psi に記録する点だけが違います。バックトラックは末尾の最尤状態から $\psi$ を辿って前へ遡るだけです。
動作確認 — 古典「公平／不正なコイン」例題
Rabinerチュートリアルでもおなじみの「公平なコイン (Fair) / 不正なコイン (Loaded) を交互に使う賭場」を題材にします。観測は「表 (H=0) / 裏 (T=1)」、隠れ状態は「Fair / Loaded」。Loaded は表が出やすく、状態は粘性を持って遷移します。
import numpy as np
# 真のパラメータ
pi_true = np.array([0.5, 0.5])
A_true = np.array([[0.95, 0.05], # Fair → Fair, Fair → Loaded
[0.10, 0.90]]) # Loaded → Fair, Loaded → Loaded
B_true = np.array([[0.5, 0.5], # Fair: H/T = 50/50
[0.8, 0.2]]) # Loaded: H/T = 80/20
# 系列をサンプリング
def sample_hmm(pi, A, B, T, rng):
N, K = B.shape
z = np.zeros(T, dtype=int)
x = np.zeros(T, dtype=int)
z[0] = rng.choice(N, p=pi)
x[0] = rng.choice(K, p=B[z[0]])
for t in range(1, T):
z[t] = rng.choice(N, p=A[z[t-1]])
x[t] = rng.choice(K, p=B[z[t]])
return z, x
rng = np.random.default_rng(0)
z_true, x_obs = sample_hmm(pi_true, A_true, B_true, T=400, rng=rng)
# モデルを正しいパラメータで初期化して、3つの計算を試す
hmm = DiscreteHMM(pi_true, A_true, B_true)
loglik = hmm.log_likelihood(x_obs)
gamma = hmm.posterior(x_obs)
mpm = gamma.argmax(axis=1)
vit, _ = hmm.viterbi(x_obs)
# 一致率
acc_mpm = (mpm == z_true).mean()
acc_vit = (vit == z_true).mean()
print(f"log p(x_1:T) = {loglik:.2f}")
print(f"MPM accuracy = {acc_mpm:.3f}")
print(f"Viterbi accuracy = {acc_vit:.3f}")
print(f"MPM と Viterbi が異なる時刻数: {(mpm != vit).sum()} / {len(z_true)}")
出力例では、MPMとViterbiの一致率はほぼ拮抗しつつ、両者が違う時刻が10–30箇所程度現れます。これは前節で論じた「MPMは局所最尤、Viterbiは大域最尤」の差が、実際に検出される境界で出ていることを示しています。多くの時刻では一致するものの、状態遷移の前後数サンプルで滑らかさを優先するViterbiと、$\gamma_t$ をそのまま最大化するMPMで判断が分かれる、というわけです。
可視化 — 真の状態と推定の比較
3つの推定（事後確率、MPM、Viterbi）を真の状態系列と並べてプロットします。
import matplotlib.pyplot as plt
fig, axes = plt.subplots(4, 1, figsize=(11, 7), sharex=True)
axes[0].step(range(len(x_obs)), x_obs, where='post', color='gray')
axes[0].set_ylabel('Obs (0=H, 1=T)'); axes[0].set_yticks([0, 1])
axes[1].step(range(len(z_true)), z_true, where='post', color='black')
axes[1].set_ylabel('True z'); axes[1].set_yticks([0, 1]); axes[1].set_yticklabels(['Fair','Loaded'])
axes[2].plot(gamma[:, 1], color='steelblue', lw=1.5)
axes[2].fill_between(range(len(gamma)), 0, gamma[:, 1], alpha=0.3)
axes[2].set_ylabel('P(Loaded | x)'); axes[2].set_ylim(-0.05, 1.05)
axes[3].step(range(len(vit)), vit, where='post', color='crimson', label='Viterbi')
axes[3].step(range(len(mpm)), mpm + 0.06, where='post', color='darkgreen', alpha=0.7, label='MPM')
axes[3].set_ylabel('Decoded'); axes[3].set_yticks([0, 1]); axes[3].legend(loc='upper right')
axes[3].set_xlabel('time t')
plt.tight_layout()
plt.savefig('hmm_decoding.png', dpi=140, bbox_inches='tight')
plt.show()
このグラフから、いくつかの重要な観察ができます。第1に、観測（最上段）からは「H/T」がランダムに見えるだけで、背後の構造はまったく見えません。第2に、事後確率 $p(\text{Loaded} \mid x_{1:T})$（3段目）は真の状態（2段目）が Loaded のときに高い値、Fair のときに低い値を取り、しかも遷移付近で滑らかに切り替わっていることが見て取れます。第3に、Viterbi と MPM の系列はほとんど一致しますが、遷移境界付近の数サンプルで微小なずれが出ます。これが理論で予言した「局所最尤と大域最尤の差」です。
Baum-Welchの実装と収束
ランダムな初期値からスタートして、Baum-Welch（EM）で真のパラメータに収束するか確認します。
def baum_welch(obs, N, K, n_iter=50, tol=1e-4, rng=None):
"""離散観測HMMの Baum-Welch 学習"""
rng = np.random.default_rng(0) if rng is None else rng
T = len(obs)
# ランダム初期化（行ごとに正規化）
pi = rng.dirichlet(np.ones(N))
A = rng.dirichlet(np.ones(N), size=N)
B = rng.dirichlet(np.ones(K), size=N)
history = []
for it in range(n_iter):
model = DiscreteHMM(pi, A, B)
log_alpha = model.forward(obs)
log_beta = model.backward(obs)
log_lik = logsumexp(log_alpha[-1])
history.append(log_lik)
# gamma: (T, N)
log_gamma = log_alpha + log_beta
log_gamma -= logsumexp(log_gamma, axis=1, keepdims=True)
gamma = np.exp(log_gamma)
# xi: (T-1, N, N) log xi_t(i,j) = log_alpha_t(i) + log_A_ij + log_B_j(x_{t+1}) + log_beta_{t+1}(j) - log_lik
logB = model.log_emission(obs)
log_xi = (log_alpha[:-1][:, :, None]
+ model.log_A[None, :, :]
+ (logB[1:] + log_beta[1:])[:, None, :]
- log_lik)
xi = np.exp(log_xi)
# Mステップ
pi = gamma[0]
A = xi.sum(axis=0) / gamma[:-1].sum(axis=0)[:, None]
for k in range(K):
B[:, k] = gamma[obs == k].sum(axis=0) / gamma.sum(axis=0)
# 収束判定
if it > 0 and abs(history[-1] - history[-2]) < tol:
break
return pi, A, B, history
pi_hat, A_hat, B_hat, hist = baum_welch(x_obs, N=2, K=2, n_iter=60, rng=rng)
print("学習後の遷移行列:")
print(np.round(A_hat, 3))
print("学習後の観測分布:")
print(np.round(B_hat, 3))
print(f"対数尤度 (init → final): {hist[0]:.2f} → {hist[-1]:.2f}")
学習後の遷移行列と観測分布を見ると、真のパラメータと（状態のラベル付け順序が入れ替わる可能性を除いて）よく一致します。ラベルスイッチング はEMの本質的な不定性で、HMMでは「どちらが Fair でどちらが Loaded か」はモデル自体には決められません。実用上は、観測分布や遷移の自己ループ確率の大小などから人間が解釈を割り当てます。対数尤度の履歴 hist を見ると、反復ごとに 単調増加 していることも確認でき、EMの理論的保証（観測対数尤度の単調非減少）が成り立っています。
左パネルは5種類の異なる初期値（seed）から始めたBaum-Welchの収束過程を示しています。いずれのカーブも単調に上昇しており、EMの理論的保証（対数尤度の単調非減少）が実際に成り立っていることが確認できます。一方で収束後の値にバラつきがあり、局所最適への依存性も見て取れます。右パネルは1例の各イテレーションでの対数尤度増分を棒グラフで示したもので、序盤に大きく改善し、後半には微小な増分のみになる典型的なEM収束パターンが読み取れます。
株価レジーム検出 — hmmlearn で Gaussian HMM
連続値観測のHMMは hmmlearn を使うのが手軽です。S&P500の日次リターンに2状態 Gaussian HMM を当て、「低ボラ／高ボラ」のレジームを抽出してみます。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from hmmlearn.hmm import GaussianHMM
# 合成データで代用（リターン系列：低ボラ期と高ボラ期が混ざるレジームスイッチ）
rng = np.random.default_rng(42)
T = 1500
# 真のレジーム：100日ごとに状態を切り替える簡易シナリオ
true_regime = np.zeros(T, dtype=int)
t = 0
while t < T:
n = rng.integers(60, 200)
true_regime[t:t+n] = (t // 100) % 2
t += n
true_regime = true_regime[:T]
returns = np.where(true_regime == 0,
rng.normal(0.0005, 0.005, T), # 低ボラ
rng.normal(-0.0010, 0.020, T)) # 高ボラ
price = 100 * np.cumprod(1 + returns)
# Gaussian HMM (2状態) を学習
X = returns.reshape(-1, 1)
model = GaussianHMM(n_components=2, covariance_type='diag',
n_iter=200, random_state=0, tol=1e-4)
model.fit(X)
hidden = model.predict(X) # Viterbi デコード
post = model.predict_proba(X) # 事後確率
# 学習されたパラメータの確認
print("状態ごとの平均リターン:", model.means_.ravel())
print("状態ごとの標準偏差 :", np.sqrt(model.covars_.ravel()))
print("遷移行列:\n", np.round(model.transmat_, 3))
# 「高ボラ状態」がインデックス何番に割り当てられたか自動判定
vol = np.sqrt(model.covars_.ravel())
high_vol_state = int(np.argmax(vol))
fig, axes = plt.subplots(3, 1, figsize=(12, 7), sharex=True)
axes[0].plot(price, color='black', lw=1.2); axes[0].set_ylabel('Price')
axes[1].plot(returns, color='gray', lw=0.6); axes[1].set_ylabel('Daily return')
axes[1].axhline(0, color='black', lw=0.4)
axes[2].plot(post[:, high_vol_state], color='crimson', lw=1.2,
label=f'P(high-vol | x) = state {high_vol_state}')
axes[2].fill_between(range(T), 0, post[:, high_vol_state], color='crimson', alpha=0.25)
axes[2].step(range(T), true_regime, where='post', color='blue', alpha=0.4, label='True regime (0/1)')
axes[2].set_ylim(-0.05, 1.1); axes[2].set_ylabel('P(high-vol)')
axes[2].legend(loc='upper right')
axes[2].set_xlabel('day')
plt.tight_layout()
plt.savefig('hmm_regime.png', dpi=140, bbox_inches='tight')
plt.show()
このプロットから3つの観察ができます。第1に、価格系列（最上段）にはレジームの構造が肉眼ではほとんど見えませんが、リターン系列（中段）を見ると振幅の大小に「束」があることが分かります。第2に、HMMの事後確率（下段、赤色）は真のレジーム（青色）の高ボラ区間で1に近づき、低ボラ区間で0に近づきます。教師なし学習だけで隠れたレジームを抽出できているわけです。第3に、レジーム遷移の直後では事後確率がしばらく中間値を取り、HMMが「自信を持って」状態を切り替える前に観測の積み重ねを待っていることが読み取れます。これは Forward-Backward が 全期間の観測を使って事後を計算している こと（オンラインなFiltering $\alpha$ ではなく、Smoothing $\gamma$）の表れです。
3段のパネルで、最上段の価格系列には赤いシェードで「真の高ボラ区間」が重ねられています。中段のリターン系列は高ボラ区間（赤シェード付近）で振幅が明らかに大きく、HMMがこの違いを捉えています。下段の赤曲線（HMM推定）が青のステップ（真のレジーム）に概ね追従しており、Forward-Backwardによるスムージング（全期間の観測を活用した $\gamma$ 計算）が状態の遷移をなめらかに捉えていることが確認できます。
実データに適用する場合は、yfinance などで S&P500 の終値を取得し、$\log r_t = \log(p_t / p_{t-1})$ で日次対数リターンを作って同じパイプラインに流せば、リーマンショックやコロナショックの高ボラ期間がきれいに検出されます。状態数を増やせば「強気・弱気・横ばい」のような3レジームモデルにもなります。
これで、HMMの理論から実装、実応用までの一通りが揃いました。最後に、これらが他のどんな領域に繋がっていくかを整理します。
応用 — 音声・遺伝子・株価
HMMが歴史的に成功してきた応用領域を、技術的な焦点とともに簡単に俯瞰します。
音声認識 (GMM-HMM時代): 音声波形を10ms程度のフレームに区切り、各フレームから MFCC（メル周波数ケプストラム係数）を抽出して観測ベクトル $x_t$ とします。隠れ状態は音素（あるいはトライフォン）に対応し、各状態の観測分布をガウス混合（GMM）で表現する GMM-HMM が、2010年代前半まで音声認識の主役でした。現在は DNN-HMM やEnd-to-Endのニューラル系（Transformer/RNN-T）に主役を譲りましたが、デコーダ側のVitebi探索（あるいはビームサーチ）は今でも本質的に同じ動的計画法です。
生物情報学 (CpGアイランド検出・遺伝子予測): DNA配列はATGCの4記号で表せます。CpGアイランド（CG含有量の高い領域、遺伝子のプロモーターに多い）を「CpGアイランド内 / 外」の2状態HMMで検出するのが古典です。離散観測 $\{A, T, G, C\}$ なので、本記事の DiscreteHMM がそのまま使えます。これを拡張すると、エクソン／イントロン／非翻訳領域を識別する 遺伝子予測HMM になり、GENSCANやHMMERといったツールに発展しました。
自然言語 (品詞タグ付け): 単語列 $w_1, \dots, w_T$ を観測、品詞列 $t_1, \dots, t_T$ を隠れ状態とした HMM が、品詞タグ付け (POS tagging) の標準アプローチでした（現在はBERT等のニューラルが主流）。$A_{ij}$ は「名詞のあとに動詞」のような構文遷移、$B_j(w)$ は「動詞は『go, run, eat』のような単語を出す」という語彙確率を捉えます。
金融 (レジームスイッチ検出): 上の実装例で示したように、Gaussian HMM はリターン系列の「強気・弱気」「低ボラ・高ボラ」レジームを検出します。HamiltonのMarkov Switchingモデルとして計量経済学でも長く使われており、リスク管理やポートフォリオの動的調整に応用されます。
人間活動認識: スマートフォンやウェアラブルの加速度センサー時系列から、「歩行・走行・座位・立位」を推定するのに HMM がよく使われます。観測がベクトル値なので Gaussian HMM、状態遷移が物理的に滑らか（突然走行から座位にはなりにくい）なので HMM の自己ループの強さがちょうど機能します。
連続化と高次化 — カルマンフィルタとの関係
HMMの隠れ状態が 連続値（実数ベクトル）になり、遷移と観測が 線形ガウス であれば、それは カルマンフィルタ（線形状態空間モデル）になります。Forward-Backward は Kalman Filter / RTS Smoother に対応し、Viterbi は最尤軌道推定に対応します。実際、両者は「離散状態空間の動的計画法」と「連続状態空間の解析的Bayes更新」という、本質的に同じ枠組みの2バージョンです。さらに非線形に拡張したのが 拡張カルマンフィルタ（EKF） や 粒子フィルタ、深層化したのが Mamba などの最新の状態空間モデルです。
このようにHMMは、現代の系列モデリング全体を貫く骨格として、今も生き続けています。
まとめ
本記事では、隠れマルコフモデルの3つの基本問題と、それらを解く動的計画法を解説しました。
- HMMの定義: 隠れ状態 $z_t \in \{1, \dots, N\}$ がマルコフ連鎖、観測 $x_t$ は $z_t$ にのみ条件付けて生成。パラメータは $\lambda = (\pi, A, B)$。
- Forward再帰 $\alpha_{t+1}(j) = [\sum_i \alpha_t(i) A_{ij}] B_j(x_{t+1})$ で観測尤度 $p(x_{1:T})$ を $O(N^2 T)$ で計算。素朴な $O(N^T)$ を多項式に圧縮。
- Backward再帰 $\beta_t(i) = \sum_j A_{ij} B_j(x_{t+1}) \beta_{t+1}(j)$ と組み合わせ、周辺事後 $\gamma_t(i) \propto \alpha_t(i)\beta_t(i)$ で 個別最大事後確率推定 (MPM) を得る。
- Viterbi: $\sum$ を $\max$ に置き換え、$\delta_t(j) = [\max_i \delta_{t-1}(i) A_{ij}] B_j(x_t)$ で全体最尤系列を求める。MPMとは一般に異なる結果になる。
- Baum-Welch (EM): $\gamma_t, \xi_t(i, j)$ を期待値カウントとして $\pi, A, B$ を更新。観測対数尤度が単調非減少、ただし局所最適。
- 数値安定化: 確率の積でunderflowするため、log領域＋
logsumexpで実装。 - Pythonでの確認: スクラッチ実装でコイン例題を検証、hmmlearn で Gaussian HMM による株価レジーム検出を確認。
HMMはそれ自身が今でも多くの実応用で生きているだけでなく、カルマンフィルタ・粒子フィルタ・ニューラル状態空間モデルへと続く「系列モデリングの祖型」です。Forward-Backward の動的計画法の構造を一度しっかり身に着けると、後続のすべての系列モデルが見通しよく理解できるようになります。
次のステップとして、以下の記事も参考にしてください。
