NSM-030

Q行列の最尤推定 — 逆問題の数値的解法

作成日: 2026-06-03 / §6 逆問題パート / Mathematics in Neuroscience プレゼン準備

このノートの中心命題:
単一チャネル記録から観測された dwell time データ $\{t_O^{(k)}, t_C^{(k)}\}$ に対し、 尤度関数 $L(Q)$ を最大化することで Q 行列を推定するのが最尤推定（MLE）である。 これは「順問題（Q → $P(t)$）」の逆である「逆問題（データ → Q）」の数値的解法であり、 Baum-Welch（EM アルゴリズム）や数値最適化によって実装される。

1. 問いの設定 — 逆問題としての MLE

Q 行列を既知とすれば、Kolmogorov 方程式の解 $P(t) = e^{Qt}$ によって 任意の時刻における遷移確率が計算できる（順問題）。 しかし実験では Q 行列は未知であり、観測されるのは dwell time データである。

問いは次のように立式される:

「どのような Q 行列であれば、この実験データが得られる確率が最も高いか」
= 最尤推定（Maximum Likelihood Estimation, MLE）の問い
問題の種類既知未知手法
順問題 Q 行列 $P(t)$、dwell time 分布 $P(t) = e^{Qt}$、固有値分解
逆問題（MLE） 観測 dwell time データ Q 行列（速度定数 $\theta$） 尤度最大化、数値最適化

2. 尤度関数 $L(Q)$ の構成

2.1 尤度関数の定義

単一チャネル記録では、チャネルの開口時間（open time）$t_O^{(k)}$ と 閉口時間（shut time）$t_C^{(k)}$ が交互に観測される。 $N$ イベントの独立な観測に対する尤度関数は次式で定義される:

$$L(Q) = \prod_{k=1}^{N} f_O\!\left(t_O^{(k)}\right) \cdot f_C\!\left(t_C^{(k)}\right)$$
記号意味
$N$観測されたイベント（開閉サイクル）の総数
$t_O^{(k)}$$k$ 番目の open time（開口持続時間）
$t_C^{(k)}$$k$ 番目の shut time（閉口持続時間）
$f_O(t)$open time の確率密度関数（Q 行列の $Q_{OO}$ ブロックから導出）
$f_C(t)$shut time の確率密度関数（Q 行列の $Q_{CC}$ ブロックから導出）
積（∏）の意味:
各イベントはマルコフ性（無記憶性）により独立であるため、全イベントの同時確率は 個別の確率密度の積で表される。$N = 1000$ イベントであれば 1000 個の密度の積となり、 数値的に非常に小さな値になる。このため実際の計算では対数尤度（§4）を用いる。

2.2 確率密度関数 $f_O, f_C$ の形

open time 分布と shut time 分布はそれぞれ「指数の和（mixture）」の形をとる（NSM-003 参照）:

$$f_O(t) = \boldsymbol{\pi}_O\, e^{Q_{OO} t}\,(-Q_{OO})\mathbf{1}_O, \qquad f_C(t) = \boldsymbol{\pi}_C\, e^{Q_{CC} t}\,(-Q_{CC})\mathbf{1}_C$$

ここで $\boldsymbol{\pi}_O, \boldsymbol{\pi}_C$ は各状態集合の開始確率ベクトル、 $\mathbf{1}$ は全成分 1 のベクトルである。これらの式は Q 行列のブロック構造に依存する（§3 参照）。

3. $Q_{AA}$ ブロック行列の構造

3.1 Q 行列のブロック分割

n 個の状態を「Open 状態の集合 $\mathcal{O}$」と「Closed 状態の集合 $\mathcal{C}$」に分割する。 このとき Q 行列は次の $2 \times 2$ ブロック構造を持つ:

$$Q = \begin{pmatrix} Q_{OO} & Q_{OC} \\ Q_{CO} & Q_{CC} \end{pmatrix}$$
ブロックサイズ意味
$Q_{OO}$ $n_O \times n_O$ Open 状態内の遷移（Open → Open）
$Q_{OC}$ $n_O \times n_C$ Open → Closed への遷移速度定数
$Q_{CO}$ $n_C \times n_O$ Closed → Open への遷移速度定数
$Q_{CC}$ $n_C \times n_C$ Closed 状態内の遷移（Closed → Closed）

$f_O(t)$ の計算には $Q_{OO}$ のみ、$f_C(t)$ の計算には $Q_{CC}$ のみが使われる。 これはチャネルが「Open である間は $Q_{OO}$ で記述される部分過程に従う」という事実を反映している。

Q_OO Open 内遷移 Q_OC Open→Closed Q_CO Closed→Open Q_CC Closed 内遷移 f_O(t) の計算に使用 f_C(t) の計算に使用 状態間の遷移速度定数 Q 行列全体の行和 = 0 （確率保存）を保つ
図1. Q 行列の Open/Closed ブロック分割。$f_O$ の計算には $Q_{OO}$、$f_C$ の計算には $Q_{CC}$ が使用される。

4. 最尤推定 — 対数尤度の最大化

4.1 なぜ対数をとるか

積 $L(Q) = \prod_{k=1}^N f_O(t_O^{(k)}) \cdot f_C(t_C^{(k)})$ は $N$ が大きいとき数値的にアンダーフロー（ゼロへの丸め誤差）を起こす。 対数をとることで積が和に変換され、数値的安定性が得られる:

$$\log L(Q) = \sum_{k=1}^{N} \left[\log f_O\!\left(t_O^{(k)}\right) + \log f_C\!\left(t_C^{(k)}\right)\right]$$

対数は単調増加であるため、$L(Q)$ を最大化することと $\log L(Q)$ を最大化することは等価である。

4.2 最尤方程式

Q 行列の各速度定数を $\theta = (q_{ij})_{i \neq j}$ とまとめたとき、 対数尤度を最大化する $\theta$ は次の方程式の解である:

$$\frac{\partial \log L(Q)}{\partial \theta} = 0$$

この方程式は解析的に解けないため、数値最適化手法が必要である（§5 参照）。 推定量 $\hat{Q}$（$L(Q)$ を最大にする Q 行列）を最尤推定量という。

推定するパラメータの次元:
n 状態モデルの Q 行列は非対角成分 $n(n-1)$ 個の速度定数を持つが、 実際には 0 に固定されるものも多く（非直接遷移）、自由パラメータ数は モデルのトポロジーによって決まる。たとえば 2 状態モデルでは $\theta = (\alpha, \beta)$ の 2 パラメータのみである。

5. 数値手法

5.1 Baum-Welch アルゴリズム（EM アルゴリズム）

単一チャネル記録において、チャネルが複数の Closed 状態を持つ場合、 各時点でどの Closed 状態にいるかは直接観測されない（隠れ状態）。 このような潜在変数モデルに適した最適化手法が EM アルゴリズム（Baum-Welch）である。

ステップ名称操作
E ステップ 期待値計算 現在の Q 推定値のもとで「各時点で各状態にいる確率」を計算（前向き・後ろ向きアルゴリズム）
M ステップ 最大化 E ステップで得た期待値を使い、対数尤度を最大化する Q を更新
収束判定 — $\log L$ の変化が閾値以下になるまで E・M を繰り返す
Baum-Welch の収束保証:
EM アルゴリズムは各反復で対数尤度を単調増加させることが保証される。 ただし大域最適解ではなく局所最適解に収束する可能性があり、 初期値の選択が重要である。

5.2 Scaling-and-squaring 法（行列指数関数の数値計算）

尤度計算には $e^{Q_{OO} t}$、$e^{Q_{CC} t}$ の評価が必要であるが、 マクローリン展開を無限級数のまま評価することは数値的に非効率である。 実用的な計算には scaling-and-squaring 法が用いられる:

$$e^{Qt} \approx \left(e^{Qt/2^s}\right)^{2^s}$$

スケーリングパラメータ $s$ を十分大きく取ると $\|Qt/2^s\|$ が小さくなり、 マクローリン展開の少数項で $e^{Qt/2^s}$ を精度良く近似できる。 その後 $2^s$ 回の行列の二乗（squaring）によって元のスケールに戻す。

Step 1. スケーリング: $s$ を選び $A = Qt/2^s$ とおく（$\|A\|$ が小さくなる）
Step 2. Padé 近似またはテイラー展開で $e^A$ を計算（少数項で十分な精度）
Step 3. スクエアリング: $e^{Qt} = (e^A)^{2^s}$ を繰り返し二乗で計算
Python での実装:
scipy.linalg.expm(Q * t) が内部的に scaling-and-squaring 法を使用している。 MATLAB では expm(Q * t) が対応する。 手動で実装する必要はなく、これらの関数を直接利用できる。

6. 推定の全体的な流れ

実験データ {t_O^(k), t_C^(k)} k=1,…,N 尤度関数 L(Q) の構成 ∏ f_O(t_O^(k)) · f_C(t_C^(k)) e^{Q_AA t} の計算 ∂logL/∂θ = 0 を数値最適化 Baum-Welch（EM）または勾配法 scaling-and-squaring 法 e^{Qt} ≈ (e^{Qt/2^s})^{2^s} 最尤推定量 Q̂ を得る argmax_Q logL(Q) 固有値・時定数 τ_k = -1/λ_k・P_open の計算
図2. Q 行列の最尤推定の全体的な流れ。

6.1 推定後の解析

最尤推定量 $\hat{Q}$ が得られた後、以下の量が計算される:

量計算方法意味
固有値 $\lambda_k$ $\hat{Q}$ の固有値分解 指数分布の速度パラメータ
時定数 $\tau_k$ $\tau_k = -1/\lambda_k$ 各成分の特徴的な時間スケール
$P_{\mathrm{open}}$ 定常分布 $\pi \hat{Q} = 0$ の Open 成分の和 チャネルが開いている時間の割合
dwell time 分布 $f_O(t) = \boldsymbol{\pi}_O e^{\hat{Q}_{OO} t}(-\hat{Q}_{OO})\mathbf{1}$ 実験データと比較して推定の妥当性を検証
まとめ — MLE による Q 行列推定の要点:
  1. 観測 dwell time データから尤度関数 $L(Q) = \prod f_O(t_O^{(k)}) \cdot f_C(t_C^{(k)})$ を構成する
  2. $f_O, f_C$ の計算には Q 行列の Open/Closed ブロック（$Q_{OO}, Q_{CC}$）が必要であり、$e^{Q_{AA} t}$ を評価する
  3. $e^{Q_{AA} t}$ の数値計算には scaling-and-squaring 法（expm）を用いる
  4. $\log L(Q)$ の最大化は Baum-Welch（EM アルゴリズム）等の数値手法で解く
  5. 推定された $\hat{Q}$ から固有値・時定数・$P_{\mathrm{open}}$ を計算し、実験との整合性を検証する

参考文献

  1. Colquhoun, D. & Hawkes, A. G. (1995). The principles of the stochastic interpretation of ion channel mechanisms. In: Single-Channel Recording, 2nd ed., pp. 397–482. → 尤度関数・ブロック行列・MLE の枠組みの主要典拠
  2. Qin, F., Auerbach, A. & Sachs, F. (1997). Maximum likelihood estimation of aggregated Markov processes. Proc. R. Soc. Lond. B, 264, 375–383. → 隠れ Markov モデルとしての単一チャネル記録の MLE
  3. Moler, C. & Van Loan, C. (2003). Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review, 45(1), 3–49. → scaling-and-squaring 法の詳細

関連項目

最終更新: 2026-06-03 / §6 逆問題パート / Mathematics in Neuroscience プレゼン準備