Fit Hidden Markov Model

隠れマルコフモデルの適合

Baum-Welch期待値最大化アルゴリズムを使用して逐次観測データに隠れマルコフモデル（HMM）を適合させ、Viterbiアルゴリズムで最も可能性の高い隠れ状態系列をデコーディングし、情報量基準による最適な隠れ状態数の選択を行う。

使用タイミング

• 放出の系列を観測するが、基礎にある生成状態が直接観測不能な時
• データが有限数のレジーム間を切り替えるシステムから生成されていると疑う時
• 時系列を潜在フェーズにセグメント化する必要がある時（例：市場レジーム、音声の音素、生物学的配列アノテーション）
• 生成モデルの下で観測系列の確率を計算したい時
• 観測値が与えられた時の最も可能性の高い隠れ状態系列が必要な時（デコーディング）
• 最良の複雑性-適合トレードオフのために異なる隠れ状態数のモデルを比較する時

入力

必須

任意

手順

ステップ1: 隠れ状態と観測モデルの定義

1.1. 隠れ状態数K（またはステップ5のモデル選択のための候補範囲）を指定する。

1.2. データの型に基づいて放出分布族を選択する：

• 連続データ：ガウス（一変量または多変量）
• カウントデータ：ポアソンまたは負の二項
• カテゴリカルデータ：離散/多項

1.3. モデルの構成要素を定義する：

• 遷移行列 A（サイズK x K）：A[i,j] = P(z_t = j | z_{t-1} = i)
• 放出パラメータ theta_k（各状態kについて）：分布固有（例：ガウスでは平均と共分散）
• 初期状態分布 pi：pi[k] = P(z_1 = k)

1.4. 観測データが適切にフォーマットされていることを確認する：系列に欠損値がない、次元が一貫している、パラメータ数に対して十分な長さがある。

期待結果: K個の状態、選択された放出族、T >> K^2の長さのクリーンな観測データを持つ明確に指定されたHMMアーキテクチャ。

失敗時: データに欠損値がある場合、補完するか影響を受けるセグメントを除去する。TがKに対して小さすぎる場合、Kを減らすかデータを追加取得する。

ステップ2: パラメータの初期化

2.1. n_restarts回の再開のそれぞれについて初期パラメータを生成する：

• 遷移行列: ランダム確率行列（各行はディリクレ分布から抽出）またはわずかに摂動を加えた一様行列。
• 放出パラメータ: 観測値に対するK-meansクラスタリングを使用して平均を初期化；ガウス放出ではクラスタ分散を計算する。
• 初期分布: 一様またはK-meansのクラスタサイズに比例。

2.2. 最初の再開ではK-meansに基づく初期化を使用する（一般に最も強い出発点）。以降の再開ではランダム摂動を使用する。

2.3. すべての初期パラメータが有効であることを確認する：

• 遷移行列の行の和が1で、すべてのエントリが正。
• 放出パラメータが有効な定義域にある（例：共分散行列が正定値）。
• 初期分布の和が1。

期待結果: 少なくとも1つのデータ駆動初期化を含むn_restartsセットの有効な初期パラメータ。

失敗時: K-meansが収束しない場合、より多くの再開を伴う純粋にランダムな初期化を使用する。共分散行列が特異な場合、正則化定数を対角に加える。

ステップ3: パラメータ推定のためのBaum-Welch EMの実行

3.1. Eステップ（前向き-後ろ向きアルゴリズム）：

• 前向き確率alpha[t,k] = P(o_1,...,o_t, z_t=k | model)を再帰で計算する：
  • alpha[1,k] = pi[k] * b_k(o_1)
  • alpha[t,k] = sum_j(alpha[t-1,j] * A[j,k]) * b_k(o_t)
• 後ろ向き確率beta[t,k] = P(o_{t+1},...,o_T | z_t=k, model)を計算する：
  • beta[T,k] = 1
  • beta[t,k] = sum_j(A[k,j] * b_j(o_{t+1}) * beta[t+1,j])
• 状態事後確率gamma[t,k] = P(z_t=k | O, model)を計算する：
  • gamma[t,k] = alpha[t,k] * beta[t,k] / P(O | model)
• 遷移事後確率xi[t,i,j] = P(z_t=i, z_{t+1}=j | O, model)を計算する。

3.2. Mステップ（パラメータ再推定）：

• 遷移行列の更新：A[i,j] = sum_t(xi[t,i,j]) / sum_t(gamma[t,i])
• 重み付き十分統計量を使用して放出パラメータを更新する：
  • ガウス平均：mu_k = sum_t(gamma[t,k] * o_t) / sum_t(gamma[t,k])
  • ガウス共分散：重み付き散布行列 + 正則化
  • 離散：b_k(v) = sum_t(gamma[t,k] * I(o_t=v)) / sum_t(gamma[t,k])
• 初期分布の更新：pi[k] = gamma[1,k]

3.3. 対数尤度を計算する：log P(O | model) = log sum_k(alpha[T,k])。アンダーフローを防ぐためlog-sum-expトリックを使用する。

3.4. スケーリング: 長い系列での数値的アンダーフローを防ぐためスケーリングされた前向き-後ろ向き変数を使用する。各タイムステップでalphaを正規化し、対数スケーリング係数を蓄積する。

3.5. 対数尤度の変化がconvergence_tol未満になるかmax_iterationsに達するまでEステップとMステップを繰り返す。

3.6. すべての再開にわたって、最高の最終対数尤度を持つパラメータセットを保持する。

期待結果: 各再開での反復を通じて単調非減少の対数尤度、max_iterations以内の収束。最終パラメータが有効（確率行列、正定値共分散）。

失敗時: 対数尤度が減少する場合、EステップまたはMステップにバグがある -- 数式を検証する。収束が非常に遅い場合、より良い初期化を試すかmax_iterationsを増加する。共分散が特異になる場合、正則化を増加する。

ステップ4: 最も可能性の高い状態系列のViterbiデコーディング

4.1. Viterbi変数を初期化する：

• delta[1,k] = log(pi[k]) + log(b_k(o_1))
• psi[1,k] = 0（先行なし）

4.2. t = 2,...,Tについて前向きに再帰する：

• delta[t,k] = max_j(delta[t-1,j] + log(A[j,k])) + log(b_k(o_t))
• psi[t,k] = argmax_j(delta[t-1,j] + log(A[j,k]))

4.3. 終端処理：

• z*_T = argmax_k(delta[T,k])
• 最良パス対数確率：max_k(delta[T,k])

4.4. t = T-1,...,1についてバックトレースする：

• z*_t = psi[t+1, z*_{t+1}]

4.5. デコードされた状態系列z* = (z*_1, ..., z*_T)とその対数確率を出力する。

4.6. Viterbiパス確率を前向きアルゴリズムからの総系列確率と比較して、最良パスがどの程度支配的かを評価する。

期待結果: 各エントリが{1,...,K}にある長さTの単一の最も可能性の高い状態系列。Viterbi対数確率は総対数尤度以下であるべき。

失敗時: Viterbiパスの対数確率が負の無限大の場合、ゼロであるべきでないところで遷移または放出確率がゼロになっている。log(0)を防ぐためフロア値を追加する。

ステップ5: モデル選択の実行（モデル次数間のBIC/AIC）

5.1. state_rangeの各候補隠れ状態数Kについて、完全なHMM（ステップ2-4）を適合する。

5.2. 自由パラメータ数pを計算する：

• 遷移行列：K * (K - 1)（各行は単体）
• 放出パラメータ：族に依存する（例：d次元の完全共分散ガウス：K * (d + d*(d+1)/2)）
• 初期分布：K - 1

5.3. 情報量基準を計算する：

• BIC = -2 * log_likelihood + p * log(T)
• AIC = -2 * log_likelihood + 2 * p
• AICc = AIC + 2*p*(p+1) / (T - p - 1)（小標本補正）

5.4. 最小のBIC（一貫性を優先）またはAIC（予測を優先）を持つモデルを選択する。両方を報告する。

5.5. 結果を表にまとめる：各Kについて対数尤度、パラメータ数、BIC、AIC、収束状態を示す。

5.6. 最適なKがstate_rangeの境界にある場合、範囲を拡張して再適合する。

期待結果: BIC/AICに明確な最小値があり、最適な隠れ状態数を特定する。選択されたモデルが収束済みで解釈可能な状態の意味を持つ。

失敗時: 明確な最小値がない場合（BICが単調減少）、モデルが誤指定されている可能性がある -- 異なる放出族を検討する。すべてのモデルの対数尤度が低い場合、データがHMM構造に従わない可能性がある。

ステップ6: ホールドアウトデータと事後デコーディングによる検証

6.1. データを訓練セットと検証セットに分割する（例：80/20、または複数の系列が利用可能な場合はそれを使用）。

6.2. 訓練データでモデルを適合する。前向きアルゴリズムを使用してホールドアウトデータの対数尤度を計算する（パラメータを再適合しない）。

6.3. 事後デコーディング（Viterbiの代替）：

• 各タイムステップで最大の事後確率を持つ状態を割り当てる：z^_t = argmax_k(gamma[t,k])
• これは正しくデコードされる状態の期待数を最大化する（同時パス確率を最大化するViterbiと対照的）。

6.4. Viterbiと事後デコーディングを比較する：

• 2つのデコード系列間の一致率を計算する。
• 不一致の領域は曖昧な状態割り当てを示す。

6.5. 状態の解釈可能性を評価する：

• 各状態の放出パラメータ（平均、分散、離散分布）を調べる。
• 状態がドメインの文脈で意味のあるレジームに対応することを確認する。
• 状態の滞在時間（Aの対角により示される）が妥当であることを確認する。

6.6. 観測値あたりのホールドアウト対数尤度を計算し、モデル次数間で比較して訓練セットのモデル選択を確認する。

期待結果: ホールドアウト対数尤度が訓練対数尤度に合理的に近い（重大な過適合なし）。Viterbiとデコーディングがタイムステップの90%以上で一致する。状態が異なる解釈可能な放出分布を持つ。

失敗時: ホールドアウト尤度が訓練よりもはるかに悪い場合、モデルが過適合している -- Kを減らすか正則化を増加する。状態が解釈できない場合、異なる初期化または異なる放出族を試す。

バリデーション

• 各再開でBaum-Welch反復を通じて対数尤度が単調非減少
• 遷移行列が行確率的（行の和が1、すべてのエントリが非負）
• 放出パラメータが有効な定義域にある（正定値共分散、有効な確率分布）
• Viterbiパス対数確率が総系列対数確率を超えない
• BIC/AIC曲線が選択されたモデル次数で明確な最小値を示す
• ホールドアウト対数尤度がモデルが訓練セットを超えて汎化することを確認する
• 前向きと後ろ向きの確率計算が一致する：P(O) = sum_k(alpha[T,k]) = sum_k(pi[k] * b_k(o_1) * beta[1,k])

よくある落とし穴

• EMでの局所最適: Baum-Welchアルゴリズムは局所最大に収束し、必ずしも大域最大ではない。常に複数のランダム再開を使用して最良のものを選択する。
• 数値的アンダーフロー: 前向き-後ろ向き確率は系列長に対して指数的に縮小する。ゼロへのアンダーフローを防ぐために対数空間計算またはスケーリング変数を使用する。
• 状態数過多による過適合: 追加の隠れ状態はそれぞれO(K + d^2)のパラメータを追加する。モデル選択には（尤度だけでなく）BICを使用し、ホールドアウトデータで検証する。
• ラベルスイッチング: 隠れ状態は置換を除いてのみ同定可能。再開間でモデルを比較する際、インデックスではなく放出パラメータで状態をマッチングする。
• 退化状態: 1つの状態が単一の観測値を説明するために崩壊する（分散がほぼゼロのガウス）ことがある。共分散行列への正則化がこれを防ぐ。
• Viterbiと事後デコーディングの混同: Viterbiは単一の最良同時パスを与え、事後デコーディングは各タイムステップでの最良の周辺状態を与える。異なる質問に答え、大幅に異なることがある。
• 状態滞在時間の無視: 標準HMMに内在する幾何分布の滞在時間は、長いレジーム持続時間を持つデータに適さない場合がある。滞在時間が非幾何分布の場合は隠れセミマルコフモデルを検討する。

関連スキル

• マルコフ連鎖のモデル化 -- 隠れ層の基礎にある遷移構造を理解するための前提知識
• 確率過程のシミュレーション -- テスト用の合成HMMデータの生成や適合モデルからの事後予測チェックのシミュレーションに使用可能