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

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

使用タイミング

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

入力

必須

observations

n_hidden_states

emission_type

"gaussian"

"discrete"

"poisson"

"multinomial"

任意

initial_params

n_restarts

max_iterations

convergence_tol

state_range

[n_hidden_states]

covariance_type

"full"

"full"

"diagonal"

"spherical"

regularization

手順

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

1.1. 隠れ状態数

K

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

失敗時: データに欠損値がある場合、補完するか影響を受けるセグメントを除去する。

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])

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

alpha

3.5. 対数尤度の変化が

convergence_tol

max_iterations

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

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

5.1.

state_range

K

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

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データの生成や適合モデルからの事後予測チェックのシミュレーションに使用可能