スマートフォンに「おはよう」と話しかけたとき、マイクが拾うのは空気の振動だけです。それでも音声認識システムは、波形の背後にある「お」「は」「よ」「う」という音素の系列を当てに行き、最終的に文字列に変換します。生物学者がDNAの塩基配列 ATGC を眺めるとき、彼らはその裏側に「遺伝子コード領域」「非コード領域」「CpGアイランド」といった見えない区分があると考え、観測された塩基から状態を推定しようとします。金融アナリストが日々の株価リターンを見るとき、彼らは「市場は今、強気レジームなのか、それとも弱気レジームなのか」を読み取ろうとします。
これらに共通する構図は驚くほど単純です。観測できる時系列の背後に、観測できない離散的な状態が存在し、その状態が時間とともに遷移しながら観測を生成している。この状況をもっとも素直に確率モデル化したものが、本記事で扱う 隠れマルコフモデル (Hidden Markov Model, HMM) です。
HMMは1970年代に音声認識の文脈で本格的に発展して以来、半世紀にわたって時系列モデリングの中核を担ってきました。近年ではTransformerや状態空間モデル(Mamba)など深層学習が台頭していますが、それでもHMMは「観測の背後にある離散状態」というシンプルで強力な発想を提供し続けています。応用は今も幅広く息づいており、音素デコーダ、CpGアイランド検出、品詞タグ付け、金融レジーム検出、ウェアラブルセンサーからの活動認識など、枚挙にいとまがありません。
本記事の内容
- マルコフ性・マルコフ連鎖の復習（遷移行列、定常分布）
- 隠れマルコフモデルの直感的なイメージと、$\lambda = (\pi, A, B)$ による定式化
- HMMの3つの基本問題（評価・復号・学習）の概要
- 状態空間モデルファミリーにおけるHMMの位置づけ（カルマンフィルタとの関係）
- Pythonで小さなHMMを実装し、状態列のサンプリングと挙動を確認
- 音声認識・遺伝子解析・株価レジーム検出など実応用への接続
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
直感 — 観測の裏に潜む状態を当てに行く
部屋に閉じこもった友人と毎日電話で会話している場面を想像してください。友人の今日の気分（「機嫌が良い」「機嫌が悪い」）は、扉が閉じているので直接見ることはできません。しかし、電話越しの話題の内容（「冗談を言う」「愚痴をこぼす」「無口で短い返事しか返さない」）から、なんとなく察しがつきます。「3日連続で愚痴ばかりだから、いまは『機嫌が悪い』モードが続いているな」とか、「今日は冗談が混じってきたから、そろそろ『機嫌が良い』に切り替わってきたかな」と推測するわけです。
この状況には、HMMの本質がすべて詰まっています。直接観測できる量（話の内容）の背後に、観測できない隠れた状態（友人の気分）が存在し、その状態がゆっくり遷移しながら観測を生成している。私たちは観測の系列から、隠れ状態の系列を確率的に推定したい。
ここに、HMMが置く2つの仮定を加えると、モデルが完成します。
1つ目は マルコフ性 です。今日の気分は昨日の気分にだけ依存し、3日前の気分は（昨日の気分が決まれば）もう関係ない。気分は「直前の気分→今日の気分」という1次の遷移確率で動くと仮定します。
2つ目は 観測の条件付き独立性 です。今日の話の内容は、今日の気分にだけ依存して決まる。気分が同じなら、過去の発言や明日の発言とは独立に分布する。
この2つの仮定だけで、驚くほど強力なモデルが構築できます。気分の遷移確率、気分から発言が出る確率、初期気分の確率さえ与えてやれば、観測された発言系列から「機嫌の良し悪し」の確率分布が計算できるのです。
天気予報の例でもまったく同じ構図が成り立ちます。今日の天気（晴れ・曇り・雨）の裏に「高気圧支配下」「低気圧通過中」のような大気状態があり、状態がゆっくり遷移しながら毎日の天気を生成している。観測した天気の系列から、現在の大気状態を推定する——というのもHMMの典型的な使い方です。
直感はここまでにして、まずはHMMの土台である「マルコフ連鎖」を復習し、その上に「観測される変数」を載せていく形でHMMを定式化しましょう。
マルコフ連鎖の復習
隠れマルコフモデルの「マルコフ」は、ロシアの数学者アンドレイ・マルコフに由来します。彼が20世紀初頭に研究した マルコフ連鎖 (Markov chain) は、HMMの土台となる確率過程です。
マルコフ性
時系列の確率変数 $X_1, X_2, X_3, \dots$ がマルコフ性をもつとは、
$$ p(X_n \mid X_1, X_2, \dots, X_{n-1}) = p(X_n \mid X_{n-1}) $$
が成り立つことをいいます。「現在の状態が分かれば、未来は過去と独立になる」と言い換えてもよいでしょう。日本語では「マルコフ性」「マルコフ的」と呼ばれます。
このマルコフ性のおかげで、同時分布がきれいに因子分解できます。確率の乗法定理をそのまま適用すると、
$$ p(X_1, X_2, \dots, X_N) = p(X_1) \prod_{n=2}^{N} p(X_n \mid X_1, \dots, X_{n-1}) $$
ですが、マルコフ性を使えば右辺の条件部分が1ステップだけに縮みます。
$$ p(X_1, X_2, \dots, X_N) = p(X_1) \prod_{n=2}^{N} p(X_n \mid X_{n-1}) $$
これが 1次マルコフモデル の同時分布です。$p(X_1)$ を 初期分布、$p(X_n \mid X_{n-1})$ を 遷移確率 と呼びます。データ点が高々1ステップしか相互作用しないので、計算もパラメータも劇的に簡単になります。
遷移行列
状態空間が有限 $\{1, 2, \dots, N\}$ のときは、遷移確率を $N \times N$ の行列 $A = (A_{ij})$ にまとめられます。
$$ A_{ij} = p(X_{n+1} = j \mid X_n = i) $$
$A_{ij}$ は「いま状態 $i$ にいて、次のステップで状態 $j$ に移る確率」です。確率なので $A_{ij} \geq 0$、各行で $\sum_j A_{ij} = 1$（どこかには必ず遷移する）が要請されます。このような非負成分・行和1の行列を 確率行列 (stochastic matrix) と呼びます。
3状態の例（晴れ・曇り・雨）で書くと、
$$ A = \begin{pmatrix} 0.7 & 0.2 & 0.1 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.3 & 0.5 \end{pmatrix} $$
のような形になります。1行目は「晴れの翌日は、70%晴れ・20%曇り・10%雨」のような意味です。
時刻 $n$ における状態の確率分布をベクトル $\bm{p}_n = (p_n(1), p_n(2), \dots, p_n(N))$ と書くと、1ステップ後の分布は単に行列積で計算できます。
$$ \bm{p}_{n+1} = \bm{p}_n A $$
これを反復すれば $\bm{p}_{n+k} = \bm{p}_n A^k$ となります。マルコフ連鎖の挙動が「行列のべき乗」に集約されるのが、有限状態マルコフ連鎖の美しいところです。
定常分布
$A$ がエルゴード的（既約・非周期的）であれば、$A^k$ を反復すると行ベクトル $\bm{\pi}^*$ に収束します。これを 定常分布 と呼び、$\bm{\pi}^* = \bm{\pi}^* A$ を満たします。「十分時間が経つと、初期分布を忘れて固有の分布に落ち着く」というのが定常分布の意味です。実は $\bm{\pi}^*$ は $A^\top$ の固有値 1 に対応する左固有ベクトル（行ベクトル）として特徴づけられ、線形代数の問題に帰着します。
定常分布はHMMでも陰に陽に登場します。たとえば長期的に「Loaded のコインが選ばれる確率」は遷移行列の定常分布で決まりますし、Baum-Welch学習で初期分布を $\pi_i = \gamma_1(i)$ で更新する代わりに「定常分布で固定する」という変種もあります。また MCMC（マルコフ連鎖モンテカルロ）のサンプラーを設計するときも、定常分布が目的分布になるよう遷移核を組むのが基本戦略です。
高次マルコフ性と1次への帰着
「1次マルコフ性は強すぎる仮定ではないか？」と感じるかもしれません。たとえば文章生成で「3単語前まで」見るような2次・3次マルコフモデルは自然に思いつきます。しかし、ここには一つトリックがあります。状態を「直前の $k$ ステップをまとめた組」に拡張すれば、$k$ 次マルコフ連鎖は 1次マルコフ連鎖に帰着 できるのです。状態空間が指数的に膨らむ代償と引き換えに、理論はすべて1次の枠組みで議論できる——これがマルコフ連鎖論の標準的な扱いです。HMMでも同様に、「1次」という仮定は本質的な制約ではなく、状態の定義しだいで任意の次数を表現できます。
ここまでで、状態がそのまま観測される マルコフ連鎖の話が一通り揃いました。次に、ここに「観測層」をかぶせて、状態を直接見ない設定 — それが隠れマルコフモデル — に進みます。
隠れマルコフモデルの数学的定義
マルコフ連鎖の世界を一段拡張しましょう。マルコフ連鎖では各時刻の状態 $X_n$ が直接観測できる前提でしたが、現実では話の内容から気分を、塩基配列から遺伝子領域を、リターン系列から市場レジームを推測したい——つまり 状態は隠れていて、別の確率分布を通して観測が生成される という設定が多いのです。これを定式化したのが隠れマルコフモデルです。
モデルの構成要素
時刻を $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)$。「状態 $j$ にいるときに、観測 $x$ が出る確率」です。これは状態ごとに違うので、状態数だけ分布があります。
これらをまとめて $\lambda = (\pi, A, B)$ と書き、これがHMMのパラメータです。観測が離散値（語彙サイズ $K$）なら $B$ は $N \times K$ の行列、観測が連続値ならガウス分布などのパラメトリック分布を用い、特にガウスを使うものを Gaussian 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) $$
パラメータ $\lambda$ で書き直せば、
$$ p(z_{1:T}, x_{1:T} \mid \lambda) = \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$ どうしの間には直接のエッジがありません。観測同士の相関はすべて隠れ状態を介してのみ生まれます。条件付き独立性が「$x_{t+1:T}$ は $z_t$ で条件付ければ $x_{1:t}$ と独立」のようにきれいに成り立ち、これが Forward-Backward 計算の正当性を担保します。
1対K符号化
実装上は、隠れ状態 $z_t$ をスカラーで持つこともあれば、$N$ 次元の one-hot ベクトルとして持つこともあります。one-hot 表現では、
$$ \bm{z}_t = (z_{t,1}, z_{t,2}, \dots, z_{t,N}), \quad z_{t,k} \in \{0, 1\}, \quad \sum_k z_{t,k} = 1 $$
として、「いまの状態」を $z_{t,k} = 1$ で表します。この表現を使うと、同時分布の指数部分が
$$ p(\bm{z}_t \mid \bm{z}_{t-1}, A) = \prod_{i=1}^{N} \prod_{j=1}^{N} A_{ij}^{z_{t-1,i} z_{t,j}} $$
のような形でコンパクトに書けるため、教科書（特にPRMLや須山ベイズ）ではこちらの記法もよく使われます。指数の $z_{t-1,i} z_{t,j}$ は「両方の状態に同時にいる」場合だけ 1 になる indicator で、$A_{ij}$ 自身が選ばれるという意味です。スカラー表記とone-hot表記は本質的に同じものを書いているだけなので、書籍を読み比べるときはどちらの流儀かをまず確認するとよいでしょう。
パラメータ数
HMMの自由なパラメータ数を数えておきましょう。初期分布 $\pi$ は $N$ 個の値で和 1 の制約があるので $N-1$ 個、遷移行列 $A$ は各行で和 1 の制約があるので $N(N-1)$ 個、離散観測の出力分布 $B$ は同じく $N(K-1)$ 個。合計で $N – 1 + N(N-1) + N(K-1) = N^2 + NK – 1$ 個です。$N=10, K=100$ なら 1099 個。深層学習モデルと比べると驚くほど少なく、これがHMMが少ない教師信号でも学習しやすい理由です。連続観測（Gaussian HMM、$D$ 次元）なら、各状態に平均 $D$ 次元と共分散行列 $D(D+1)/2$ 個を持つので、$N(D + D(D+1)/2)$ がここに加わります。
定義が揃ったところで、HMMで実際にやりたい計算が何なのかを整理しましょう。
HMMの3つの基本問題
HMMが応用される場面は、Rabinerの古典的なチュートリアル（1989年）以来、次の3つの問題に集約されてきました。評価 (Evaluation)、復号 (Decoding)、学習 (Learning) の3つです。それぞれ「与えられるもの」と「求めたいもの」が違います。
問題1: 評価 — モデルと観測の整合度を測る
モデル $\lambda = (\pi, A, B)$ と観測系列 $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$ のうち、観測をもっともよく説明するモデルを選ぶというシーンが代表的です。音声認識で「単語ごとに学習したHMMを並べて、いま聞こえた音声系列がどの単語のモデルから出たかを尤度比較で当てる」のがまさにこれです。
素朴に状態列の和を取ると、項数は $N^T$ 通り。$T = 100, N = 10$ で $10^{100}$ 通りに膨らみ、宇宙の年齢でも終わりません。これを多項式時間 $O(N^2 T)$ に落とすのが Forwardアルゴリズム です（詳細は後述のリンク先記事）。
問題2: 復号 — 隠れ状態の系列を推定する
観測 $x_{1:T}$ から、もっともらしい隠れ状態系列 $\hat z_{1:T}$ を推定する問題です。「もっともらしい」には2通りの解釈があり、それぞれ異なるアルゴリズムに対応します。
- 個別最大事後確率推定 (Posterior Marginals / MPM): 各時刻 $t$ について $\arg\max_i p(z_t = i \mid x_{1:T})$ を独立に選ぶ。各時刻の周辺事後を最大化する素直な戦略です。Forward-Backwardアルゴリズムで $\gamma_t(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})$ で系列全体として最尤の経路を求める。Viterbiアルゴリズムで動的計画法的に計算します。
両者は一般に 異なる解 を返します。MPMは時刻ごとに独立に最大化するため、遷移確率 $A_{ij} = 0$ の禁則を破る系列を返す可能性があります（例：$A_{12} = 0$ なのに $\hat z_t = 1, \hat z_{t+1} = 2$）。Viterbiは経路として整合する系列を保証しますが、各時刻で見ると最尤ではないこともあります。応用に応じて使い分けが必要です。
問題3: 学習 — パラメータを推定する
観測 $x_{1:T}$（または複数の系列）だけから、パラメータ $\lambda = (\pi, A, B)$ を推定する問題です。隠れ状態が観測できないため、最尤推定はEMアルゴリズムの形を取り、HMMに特化したものが Baum-Welchアルゴリズム として知られています。
Baum-Welchの各反復は、現在のパラメータで $\gamma_t(i)$ と $\xi_t(i, j) = p(z_t=i, z_{t+1}=j \mid x_{1:T})$ を計算するEステップと、これらを「重み付き頻度」とみなしてパラメータを更新するMステップに分かれます。観測対数尤度が反復ごとに単調非減少することが保証される、EMの典型的な構造です。
これら3つの問題は独立ではなく、内側で同じ Forward / Backward を共有します。本記事ではHMMの 全体像と数学的定義 までを扱い、Forward-Backward / Viterbi / Baum-Welch の具体的な再帰式・収束証明・実装は、記事末尾でリンクする深掘り版にゆずります。
問題の整理ができたところで、HMMを「状態空間モデル」というより広い枠組みの中に位置づけてみましょう。
状態空間モデルの中でのHMM
HMMは、より広い 状態空間モデル (State Space Model, SSM) ファミリーの一員として理解するのが見通しのよい捉え方です。状態空間モデルとは、
- 状態方程式: $z_t = f(z_{t-1}) + \text{noise}$
- 観測方程式: $x_t = g(z_t) + \text{noise}$
という形で、隠れ状態の時間発展と、状態から観測への生成過程を分けてモデル化する一般的な枠組みです。HMMは「状態空間が 離散 で、遷移と観測が 確率行列／確率分布 で与えられる」特別な場合に対応します。
逆に、状態空間が 連続 で、状態方程式と観測方程式がともに 線形ガウス であれば、それは カルマンフィルタ が扱う線形ガウス状態空間モデルになります。
| 離散状態 | 連続状態 | |
|---|---|---|
| 線形・ガウス | HMM（離散） | カルマンフィルタ |
| 一般 | HMM（離散） | EKF / UKF / 粒子フィルタ |
Forward-Backwardアルゴリズムは、連続状態の世界では カルマンフィルタ + RTSスムーザー に対応します。Viterbiは最尤軌道推定（線形ガウスでは平均軌道に一致）に、Baum-Welchは線形ガウスでの EMによるシステム同定 に対応します。実は「離散状態空間の動的計画法」と「連続状態空間の解析的Bayes更新」は、本質的に同じ枠組みの2バージョンなのです。
非線形まで広げると、拡張カルマンフィルタ (EKF) や アンセンテッドカルマンフィルタ (UKF)、粒子フィルタ が登場し、さらに深層化すると Mamba などの最新の状態空間モデル（Selective SSM）に繋がっていきます。HMMはその「離散・小規模」端で、いまも生き続けている古典です。
この大きな見取り図を頭に入れた上で、実際にPythonで小さなHMMを動かしてみましょう。
Python実装 — HMMを動かしてみる
理論を「動くもの」として体感するために、ここではHMMから状態列と観測列をサンプリング（順問題）し、続いて Forward 計算で観測尤度を求める部分までを実装します。詳細な Forward-Backward / Viterbi / Baum-Welch の log-domain 実装と数値安定化、株価レジーム検出の実例は、記事末尾でリンクする深掘り記事にまとめています。
HMMからのサンプリング
まずは典型例として、Rabiner のチュートリアルでもおなじみの「公平／不正なコイン (Fair / Loaded)」設定を使います。観測は「表 (H=0) / 裏 (T=1)」、隠れ状態は「Fair (0) / Loaded (1)」。Loaded のときは表が出やすく、状態自身は粘性をもって遷移する設定です。
import numpy as np
import matplotlib.pyplot as plt
# 真のパラメータ
pi_true = np.array([0.5, 0.5]) # 初期分布: Fair/Loaded を 50/50 で開始
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):
"""HMMから状態列zと観測列xを長さTぶんサンプリング"""
N, K = B.shape
z = np.zeros(T, dtype=int)
x = np.zeros(T, dtype=int)
z[0] = rng.choice(N, p=pi) # 初期状態を 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=300, rng=rng)
print("最初の20ステップの隠れ状態:", z_true[:20])
print("最初の20ステップの観測 :", x_obs[:20])
print(f"Loaded だった時刻の割合: {(z_true == 1).mean():.3f}")
このコードを走らせると、まず初期状態が $\pi$ に従って引かれ、続く各時刻で「遷移行列 $A$ で次状態を引き、その状態の出力分布 $B$ で観測を引く」という2段階のサンプリングが繰り返されます。生成された状態列 z_true は Fair と Loaded の間を粘性をもって行き来し、観測列 x_obs はその裏側を反映して「H が連続する区間」と「H/T が拮抗する区間」が交互に現れるはずです。(z_true == 1).mean() の値は遷移行列の定常分布の Loaded 成分（およそ $0.05 / (0.05 + 0.10) = 0.33$）に近づきます。
可視化 — 隠れ状態と観測
サンプリングした系列を可視化して、状態と観測の関係を視覚的に確認します。
fig, axes = plt.subplots(2, 1, figsize=(11, 5), sharex=True)
# 上段: 隠れ状態 (真の状態列)
axes[0].step(range(len(z_true)), z_true, where='post', color='black')
axes[0].set_ylabel('Hidden state z')
axes[0].set_yticks([0, 1])
axes[0].set_yticklabels(['Fair', 'Loaded'])
axes[0].set_title('Hidden state sequence (ground truth)')
# 下段: 観測列 (コインの表裏)
axes[1].step(range(len(x_obs)), x_obs, where='post', color='steelblue')
axes[1].set_ylabel('Observation x')
axes[1].set_yticks([0, 1])
axes[1].set_yticklabels(['H', 'T'])
axes[1].set_xlabel('time t')
axes[1].set_title('Observed coin flips')
plt.tight_layout()
plt.savefig('hmm_sampling.png', dpi=140, bbox_inches='tight')
plt.show()
このグラフから2つの重要な観察が得られます。第1に、上段の隠れ状態は階段状にゆっくり変化しています。これは $A_{11} = 0.95, A_{22} = 0.90$ という強い自己ループ確率の結果で、状態が一度入ると平均的に $1/(1-0.95) = 20$ ステップ程度滞在することを意味します。第2に、下段の観測列はそれよりずっと細かく上下していて、ぱっと見では Fair / Loaded の区別がつきません。観測だけ見せられて隠れ状態を推定する難しさが、この図から直観的に伝わるはずです。
このように、HMMが解こうとしている問題の本質は「下段の観測しか手に入らない状況で、上段の隠れ状態を当てに行く」ことなのです。
Forward アルゴリズムによる尤度計算
サンプリングできたところで、観測列 $x_{1:T}$ の尤度 $p(x_{1:T} \mid \lambda)$ を計算してみましょう。Forward 再帰の本体は次のたった数行で書けます。$\alpha_t(j) = p(x_{1:t}, z_t = j)$ を確率の積で素朴に持つと underflow するので、ここでは小さい $T$ で動作確認するだけにとどめ、本格的な log-domain 実装は深掘り記事で扱います。
def forward(pi, A, B, obs):
"""前向き再帰で alpha_t(j) を計算し、最終的に観測尤度を返す
パラメータ
----------
pi : (N,) 初期分布
A : (N, N) 遷移行列
B : (N, K) 離散観測の出力分布
obs: (T,) int 観測系列(整数ID列)
"""
T = len(obs)
N = len(pi)
alpha = np.zeros((T, N))
# 初期化: alpha_1(j) = pi_j * B_j(x_1)
alpha[0] = pi * B[:, obs[0]]
# 再帰: alpha_{t+1}(j) = [sum_i alpha_t(i) * A_ij] * B_j(x_{t+1})
for t in range(1, T):
alpha[t] = (alpha[t - 1] @ A) * B[:, obs[t]]
# 終端: 観測尤度は alpha_T(j) を全状態で和を取ったもの
likelihood = alpha[-1].sum()
return alpha, likelihood
# 短い系列で動かしてみる
alpha, lik = forward(pi_true, A_true, B_true, x_obs[:30])
print(f"観測尤度 p(x_1:30 | lambda) = {lik:.6e}")
print(f"対数尤度 log p(x_1:30 | lambda) = {np.log(lik):.4f}")
# 全状態を列挙する素朴な総和と一致するかチェック (N^T 通りを走査)
from itertools import product
def brute_force_likelihood(pi, A, B, obs):
"""すべての状態系列を列挙して尤度を計算 (T<=10 程度の検算用)"""
N = len(pi)
T = len(obs)
total = 0.0
for z in product(range(N), repeat=T):
p = pi[z[0]] * B[z[0], obs[0]]
for t in range(1, T):
p *= A[z[t - 1], z[t]] * B[z[t], obs[t]]
total += p
return total
lik_short = forward(pi_true, A_true, B_true, x_obs[:8])[1]
lik_bf = brute_force_likelihood(pi_true, A_true, B_true, x_obs[:8])
print(f"Forward : {lik_short:.10e}")
print(f"Brute force : {lik_bf:.10e}")
print(f"相対誤差 : {abs(lik_short - lik_bf) / lik_bf:.3e}")
ここでのポイントは、再帰式 $\alpha_{t+1}(j) = [\sum_i \alpha_t(i) A_{ij}] B_j(x_{t+1})$ をベクトル化して alpha[t] = (alpha[t-1] @ A) * B[:, obs[t]] の1行で書いている点です。Numpyのブロードキャストにより、$N$ 次元ベクトルに $N \times N$ 行列を掛けて再び $N$ 次元、そこに $N$ 次元の観測確率を要素ごとに掛けて、1ステップが完了します。各時刻で $O(N^2)$、全体で $O(N^2 T)$ の計算量です。さらに後半では「Forward再帰の答え」と「全 $N^T$ 状態を素朴に列挙した総和」が（数値誤差の範囲で）一致することを $T=8$ で検算しています。素朴な総和は $T=8, N=2$ で 256 通りに過ぎませんが、$T=100$ にすればもう $10^{30}$ 通りなので、現実には Forward 再帰なしには計算不可能です。
短い系列であれば確率の積をそのまま持っても問題ありませんが、$T$ が数百を超えると lik の値が容易に $10^{-100}$ を下回り、float64の表現範囲を出てしまいます。この underflow を防ぐために対数領域で計算する logsumexp ベースの実装が、深掘り記事の本題になります。
計算量の実測
Forward 再帰の計算量が $O(N^2 T)$ であることを、状態数 $N$ を振って実測してみます。
import time
def measure_forward_time(N, T=2000, n_trials=5, seed=0):
"""状態数 N, 長さ T の Forward を n_trials 回平均"""
rng = np.random.default_rng(seed)
pi = rng.dirichlet(np.ones(N))
A = rng.dirichlet(np.ones(N), size=N)
K = 4 # 観測語彙サイズは固定
B = rng.dirichlet(np.ones(K), size=N)
obs = rng.integers(0, K, size=T)
times = []
for _ in range(n_trials):
t0 = time.perf_counter()
forward(pi, A, B, obs)
times.append(time.perf_counter() - t0)
return np.median(times)
Ns = [2, 4, 8, 16, 32, 64]
elapsed = [measure_forward_time(N) for N in Ns]
for N, e in zip(Ns, elapsed):
print(f"N={N:3d}: {e*1000:7.2f} ms (per-step {e*1e6/2000:.2f} us)")
# log-log でフィッティング
log_N = np.log(Ns)
log_e = np.log(elapsed)
slope = np.polyfit(log_N, log_e, 1)[0]
print(f"\n実測スロープ (log-log): {slope:.2f} (理論値: 2.00)")
実行すると、状態数 $N$ を倍にするごとに実行時間がおよそ4倍に増え、log-log プロットの傾きが理論値の 2.0 に近づいているはずです（小さい $N$ では Numpy のオーバーヘッドが支配的で傾きがやや小さくなることがあります）。これは「各時刻で $N \times N$ の行列・ベクトル積をする」という Forward 再帰の構造をそのまま反映した結果で、HMMが扱える状態数の上限を見積もるのにも使えます。たとえば $N = 1000, T = 10^4$ なら $N^2 T = 10^{10}$ 演算で、現代のCPUで数十秒〜数分のオーダー。大規模な音声認識では、ビーム探索などで実効的な $N$ を絞り込んで、この限界を回避するテクニックが発達しました。
計算量グラフの可視化
実測した計算時間と理論曲線 $O(N^2)$ を重ねてプロットし、両者がきれいに揃うことを視覚的に確認します。
fig, ax = plt.subplots(figsize=(7, 5))
ax.loglog(Ns, elapsed, 'o-', color='steelblue', lw=2, markersize=8,
label='measured (Forward)')
# 理論曲線 N^2 をフィッティング (定数倍だけ揃える)
Ns_arr = np.array(Ns)
theory = elapsed[0] * (Ns_arr / Ns_arr[0]) ** 2
ax.loglog(Ns, theory, '--', color='crimson', lw=2, label=r'theoretical $O(N^2)$')
ax.set_xlabel('Number of states N')
ax.set_ylabel('Forward time [s] (T=2000)')
ax.set_title(f'Forward complexity: measured slope = {slope:.2f}')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('hmm_complexity.png', dpi=140, bbox_inches='tight')
plt.show()
このプロットでは、実測点（青、丸）と理論曲線 $\propto N^2$（赤、破線）が log-log 軸上でほぼ平行になり、Forward 再帰の理論的計算量がそのまま現実の実行時間に表れていることが確認できます。線がぴったり重ならないのは、Numpy のBLAS呼び出しオーバーヘッドや L1/L2 キャッシュ効率など、漸近解析では捨象される定数因子の効果です。それでも $N$ が大きい領域では理論どおりの 2 乗スケーリングが見えるのは、HMM の動的計画法構造が極めて素直に実装に反映されている証拠と言えます。
隠れ状態と出力の関係を確認
HMMが「状態ごとに違う観測分布をもつ」という直感を、サンプリング結果から検証します。
# 状態ごとの観測ヒストグラム
fig, axes = plt.subplots(1, 2, figsize=(9, 4), sharey=True)
for state in [0, 1]:
mask = (z_true == state)
counts = np.bincount(x_obs[mask], minlength=2)
freq = counts / counts.sum()
axes[state].bar(['H', 'T'], freq, color=['steelblue', 'crimson'])
axes[state].set_title(f'State {state} ({"Fair" if state==0 else "Loaded"})')
axes[state].set_ylim(0, 1.0)
axes[state].set_ylabel('Empirical frequency')
for i, f in enumerate(freq):
axes[state].text(i, f + 0.02, f"{f:.2f}", ha='center')
plt.tight_layout()
plt.savefig('hmm_emission.png', dpi=140, bbox_inches='tight')
plt.show()
このヒストグラムから、第1に、Fair 状態（左）では H/T がほぼ 50/50 に近い経験頻度を示し、第2に、Loaded 状態（右）では H が圧倒的に多く、Loaded の出力分布 $B_2 = (0.8, 0.2)$ を反映していることが確認できます。サンプル数が十分なら、経験頻度は真の出力分布に収束していくはずです。実装したHMMが、設定した遷移行列・出力分布どおりに動いていることが、この検証から確かめられました。
遷移行列のべき乗と定常分布
最後に、マルコフ連鎖パートで触れた定常分布が、実際に遷移行列のべき乗で得られることを数値的に確認します。
# 遷移行列を何度も自身に掛けると、行が一定のベクトルに収束する
A_pow = A_true.copy()
print("A^1 =")
print(np.round(A_pow, 4))
for k in [2, 5, 20, 100]:
A_pow = np.linalg.matrix_power(A_true, k)
print(f"\nA^{k} =")
print(np.round(A_pow, 4))
# 定常分布: A^T の固有値1に対応する左固有ベクトル
eigvals, eigvecs = np.linalg.eig(A_true.T)
idx = np.argmin(np.abs(eigvals - 1.0)) # 固有値1に最も近いもの
stationary = np.real(eigvecs[:, idx])
stationary /= stationary.sum() # 確率として正規化
print(f"\n定常分布 pi* = {np.round(stationary, 4)}")
print(f"理論値 (Fair確率) = {0.10 / (0.05 + 0.10):.4f}")
このコードの出力を見ると、$A^k$ の行は $k$ が大きくなるにつれて区別がつかなくなり、すべての行が定常分布 $\bm{\pi}^* \approx (0.667, 0.333)$ にそろってきます。これは「どこから出発しても、十分時間が経てば同じ分布に落ち着く」というエルゴード性の現れです。理論値 $\pi^*_{\text{Loaded}} = A_{12} / (A_{12} + A_{21}) = 0.05 / 0.15 = 0.333$ とも一致しており、HMMの基盤としてマルコフ連鎖が確かに正しく実装されていることが確認できます。サンプリングで「Loaded だった時刻の割合」がおよそ 0.33 に近かったのも、この定常分布の現れにほかなりません。
理論と実装の感触が掴めたところで、HMMが実際にどんな分野で活躍しているかを概観しましょう。
HMMの応用 — 半世紀の蓄積
HMMは1970年代以降、現実の様々なドメインで使われてきました。それぞれで「何が観測で、何が隠れ状態か」のマッピングを押さえると、応用への展開が見えてきます。
音声認識 (GMM-HMM時代)
音声波形を10ms程度のフレームに区切り、各フレームから MFCC (メル周波数ケプストラム係数) を抽出して観測ベクトル $x_t$ とします。隠れ状態は音素（あるいはトライフォン）に対応し、各状態の観測分布をガウス混合 (GMM) で表現する GMM-HMM が、2010年代前半まで音声認識の主役でした。Forward 再帰で単語仮説のスコアリングを行い、Viterbiでビーム探索的にデコードする——という基本骨格は、現在の DNN-HMM や End-to-End モデル（Transformer / RNN-Transducer）にも形を変えて受け継がれています。
生物情報学 (CpGアイランド検出・遺伝子予測)
DNA配列はATGCの4記号で表せます。CpGアイランド（CG含有量の高い領域、遺伝子のプロモーターに多い）を「CpGアイランド内 / 外」の2状態HMMで検出するのが古典的な応用です。離散観測 $\{A, T, G, C\}$ なので、本記事で扱った離散HMMがそのまま使えます。これを拡張すると、エクソン／イントロン／非翻訳領域を識別する 遺伝子予測HMM になり、GENSCAN や HMMER といったツールに発展しました。タンパク質配列のドメイン解析でも プロファイルHMM が広く使われています。
自然言語処理 (品詞タグ付け・固有表現抽出)
単語列 $w_1, \dots, w_T$ を観測、品詞列 $t_1, \dots, t_T$ を隠れ状態とした HMM は、品詞タグ付け (POS tagging) の標準アプローチでした。$A_{ij}$ は「名詞のあとに動詞」のような構文遷移、$B_j(w)$ は「動詞は『go, run, eat』のような単語を出す」という語彙確率を捉えます。現代ではBERT等のニューラルモデルが主流ですが、低資源言語や軽量タグ付けでは今も実用的です。同じ枠組みで 固有表現抽出 (NER) にも応用され、CRFへと発展しました。
金融 (レジームスイッチ検出)
日々のリターン系列 $r_t = \log(p_t / p_{t-1})$ を観測、市場レジーム（「強気・弱気」「低ボラ・高ボラ」）を隠れ状態とする Gaussian HMM で、市場の状態変化を抽出できます。Hamilton の Markov Switching Model として計量経済学でも長く使われており、リスク管理やポートフォリオの動的調整、レジーム依存のオプション価格付けに応用されます。学習されたボラティリティの大小から状態を「高ボラ／低ボラ」に解釈し、ポジションサイズを調整する戦略が典型例です。
人間活動認識
スマートフォンやウェアラブルの加速度・ジャイロセンサーの時系列から、「歩行・走行・座位・立位・階段昇降」などの活動を推定するのに HMM がよく使われます。観測がベクトル値なので Gaussian HMM、状態遷移が物理的に滑らか（突然「走行」から「座位」にはなりにくい）なので、HMM の自己ループ確率の強さがちょうど機能します。ヘルスケア・スポーツ・高齢者見守りなど、応用ドメインも多彩です。
手書き文字認識・ジェスチャー認識
ペンタブレットの座標時系列 $(x_t, y_t)$ や、Kinect 等の骨格座標時系列から、書かれた文字や行ったジェスチャーを推定するのにもHMMが古くから使われました。文字の各ストロークを状態とみなし、ストローク遷移と座標分布のモデル化を行います。
ロボット工学・SLAM
ロボットの環境地図と自己位置の同時推定 (SLAM) では、離散グリッドベースで状態を持つ場合に HMM の Forward-Backward が使えます。連続版はカルマンフィルタや粒子フィルタ (Monte Carlo Localization) になり、両者は理論的に同じ家系の枠組みです。屋内ロボットの「いまどの部屋にいるか」の推定なども、部屋を離散状態とすればHMMが自然に使えます。
HMMを使う際の実務的な注意
HMMを実応用するときに気をつけたいポイントもいくつかあります。第1に、状態数 $N$ の決定 は最尤推定だけでは決まらず、AIC・BIC・クロスバリデーションのようなモデル選択指標を使う必要があります。状態数を増やしすぎると過学習で「データの細かいノイズに対応する状態」が生まれ、解釈性が落ちます。第2に、Baum-Welchは局所最適に陥る ため、複数の初期値で学習して最尤のものを選ぶのが定石です。Gaussian HMM ではK-means初期化、離散HMMでは Dirichlet からの乱択初期化がよく使われます。第3に、ラベルスイッチング（学習後に「状態0が Loaded、状態1が Fair」になることもその逆もある）はEMの本質的な不定性なので、結果の解釈時に出力分布の特徴から人間がラベルを当て直すワークフローを組みます。第4に、観測が非定常（時間とともに分布自体がドリフトする）な場合は、固定パラメータのHMMでは追従できないので、オンライン学習や変分推論版のHMMを検討します。
このように、HMMの守備範囲は驚くほど広く、現代でも実用ツールとして生き続けています。深層学習が台頭しても、「観測の背後の離散状態を効率的に推論する枠組み」という骨格は今なお魅力的なのです。
まとめ
本記事では、隠れマルコフモデル (HMM) の全体像を直感・定式化・応用の三層で解説しました。
- 直感: 観測の背後に「観測できない離散状態」が存在し、状態がマルコフ連鎖的に遷移しながら観測を生成する確率モデル。気分と発言、市場レジームとリターン、音素と音声特徴量、などに自然に対応する。
- マルコフ連鎖の復習: 状態が直接観測できる場合は、初期分布と遷移行列 $A$ で完全に決まる。エルゴード的なら長時間平均が定常分布に収束。
- HMMの数学的定義: パラメータは $\lambda = (\pi, A, B)$。同時分布は $p(z_{1:T}, x_{1:T}) = \pi_{z_1} B_{z_1}(x_1) \prod_t A_{z_{t-1}z_t} B_{z_t}(x_t)$ という単純な経路積で書ける。
- 3つの基本問題: 評価（観測尤度の計算 — Forward）、復号（隠れ状態列の推定 — Forward-Backward / Viterbi）、学習（パラメータ推定 — Baum-Welch / EM）。すべて動的計画法で多項式時間に解ける。
- 状態空間モデルファミリー: 離散状態×離散時間がHMM、連続状態×線形ガウスがカルマンフィルタ、非線形に拡張するとEKF/UKF/粒子フィルタ、深層化するとMambaなどのSSMに繋がる。
- Python実装: HMMからの順サンプリングと、Forward再帰による短い系列の尤度計算を確認。状態ごとの観測ヒストグラムから「Fair vs Loaded」の出力分布の違いも確認できた。
- 応用: 音声認識・遺伝子予測・品詞タグ付け・金融レジーム検出・活動認識など、半世紀にわたって実用されてきた古典かつ現役のモデル。
ここまでで、HMMの「全体像と数学的な骨格」が掴めたはずです。次のステップとして、Forward-Backwardアルゴリズム・Viterbiアルゴリズム・Baum-Welch(EM)の 再帰式の導出と log-domain 実装、そして hmmlearn による株価レジーム検出までを一気通貫で扱った深掘り記事を用意しています。実装で手を動かして本当に「動く」HMMを書きたい方は、こちらをどうぞ。
また、関連する以下の記事も参考にしてください。