NSM-024
最尤推定（MLE）と尤度関数
作成: 2026-05-26 / NSM-024 / プロジェクト: Mathematics in Neuroscience プレゼン準備
本ノートの読み方
本ノートは二部構成です。
前半（A1〜A6）: やさしい入門 — 数式を極力使わず、コイン投げの例などで尤度の直感をつかむ。数学の予備知識なしで読める。後半（B1〜B7）: 数学者向け厳密版 — CTMC の観測データに対する尤度関数の形式的定義・Colquhoun-Hawkes 形式・最適化の数値的側面・漸近統計理論を扱う。
前半だけ読んで「雰囲気をつかむ」でも、後半から直接入っても構いません。
【本ノートでの用語規約】
確率 (probability) : パラメータ Q を固定し、データ（dwell time）の実現値を変数として見た値。$P(\text{data} \mid Q)$尤度 (likelihood) : データを観測した後で、Q を変数として同じ式を見た値。$L(Q) := P(\text{data} \mid Q)$（数値は同じだが「何を変数と見るか」が逆）最尤推定 (MLE) : 尤度 $L(Q)$ を最大にする $\hat{Q}$ を選ぶ推定法対数尤度 : $\ell(Q) := \log L(Q)$。最大化の問題として等価だが数値的に扱いやすい最適化条件式 : $\partial \log L / \partial \theta = 0$。A6 ではこの式の記号（∂, θ, = 0）を 1 つずつ読み解き、なぜ数値最適化が必要かを解説する
目次
前半 — やさしい入門
A1. 実験で何が見えて、何が見えないのか
単一イオンチャネルの実験では、パッチクランプ記録から「チャネルが開いている時間」と「閉じている時間」の系列が得られます。
観測できるもの : dwell time の系列 $t_1, t_2, \ldots, t_n$（開閉の滞在時間）直接観測できないもの : チャネルが内部でどの状態にいるか、遷移速度行列 Q の値
ところが私たちが本当に知りたいのは Q です。チャネルがどの速さで状態を遷移しているかを知れば、薬の効果や変異の影響を定量化できます。dwell time から Q を逆推定するのが「逆問題」であり、その標準的な手法が最尤推定（MLE）です。
実験で見えるもの・見えないもの
内部（直接観測不可）
状態 C1
状態 C2
状態 O
Q 行列（遷移速度）
← 推定したい！
確率分布として
dwell time を生成
観測データ（パッチクランプ）
O: t₁
C: t₂
O: t₃
C: t₄
O: t₅
dwell time 系列
← これだけが分かる
t₁, t₂, ..., tₙ
図 1. 実験で観測できるもの（dwell time 系列）と観測できないもの（内部状態・Q 行列）の関係。MLE は右から左への逆向きの推定を行う。
A2. 尤度とは何か — 視点の反転
「尤度（ゆうど、likelihood）」という概念の核心は視点の反転 にあります。
確率の視点 : Q が決まっている。そのとき「dwell time $t$ が観測される確率（密度）は何か？」→ 変数はデータ $t$、パラメータ Q は固定。尤度の視点 : データ $t$ はすでに観測されている。そのとき「どの Q がこのデータを生み出す最もらしいパラメータか？」→ 変数は Q、データは固定。
核心的な一文 : 尤度は「Q を変えたとき、観測されたデータが得られる確率（密度）がどう変わるか」を測る関数です。数値上は確率密度と同じ式ですが、「何を変数と見るか」が違います。
コイン投げの例
コインを 10 回投げたら 8 回表が出たとします。コインが歪んでいるかどうかを判定したい。
表の確率を $p$ とすると、10 回中 8 回表の確率は $\binom{10}{8} p^8 (1-p)^2$
$p = 0.5$（公平なコイン）のとき: 約 0.044
$p = 0.8$ のとき: 約 0.302
データは変わらない（8 回表）。$p$ を変えたときの「確率の大きさ」が尤度 $L(p)$ です。$p = 0.8$ の尤度がずっと高いため、「8 回表を説明するには $p = 0.8$ がもっともらしい」という結論になります。
コイン投げの尤度 L(p) — 10 回投げて 8 回表
表の確率 p
0
0.2
0.4
0.6
0.8
1.0
尤度 L(p)
p̂ = 0.8（MLE）
L(0.8) ≈ 0.302
L(0.5) ≈ 0.044
図 2. 10 回投げて 8 回表の場合の尤度 L(p)。最大値は p = 0.8 のとき（赤点）。p = 0.5 の尤度は相対的に低い（紫点）。これが最尤推定の直感。
A3. なぜ積の形になるのか
イオンチャネルの実験では、複数の dwell time $t_1, t_2, \ldots, t_n$ が独立に観測されます。
マルコフ性（無記憶性）から「次の遷移は過去の履歴に依存しない」ので、各 $t_k$ は独立です。独立な事象が「同時に起こる」確率は積になります:
$$L(Q) = f(t_1 \mid Q) \cdot f(t_2 \mid Q) \cdots f(t_n \mid Q) = \prod_{k=1}^{n} f(t_k \mid Q)$$
ここで $f(t_k \mid Q)$ は「パラメータ Q のモデルのもとで dwell time $t_k$ が観測される確率密度」です。
dwell time 系列と尤度の積
t₁ = 2.1ms
t₂=0.8ms
t₃ = 3.5ms
t₄=1.2ms
· · ·
tₙ=4.0ms
観測された dwell time 系列（独立）
f(t₁|Q)
f(t₂|Q)
f(t₃|Q)
f(t₄|Q)
· · ·
f(tₙ|Q)
×
×
×
×
L(Q) = ∏ f(tₖ|Q) ← 尤度関数
図 3. 観測された各 dwell time $t_k$ に対して密度 $f(t_k \mid Q)$ を計算し、すべての積をとったものが尤度 $L(Q)$。マルコフ性（独立性）が積の形を正当化する。
A4. 最尤推定の考え方
「尤度 $L(Q)$ を最大にする Q を選ぶ」というのが最尤推定（Maximum Likelihood Estimation）の考え方です。
$$\hat{Q} = \arg\max_Q L(Q)$$
言い換えると: 「観測されたデータが、もっとも高い確率（密度）で出てくるような Q を採用する」。これが自然で直感的な推定原理です。
A5. なぜ log を取るのか
実際の計算では $L(Q)$ を直接最大化するのではなく、$\log L(Q)$（対数尤度）を最大化します。理由は 3 つあります。
数値安定性 : $n = 1000$ の観測データがあると、$L(Q) = \prod f(t_k|Q)$ は非常に小さい数の積になり浮動小数点のアンダーフローが起きる。積 → 和 : $\log L(Q) = \sum_{k=1}^{n} \log f(t_k|Q)$ となり和の形で計算できる。等価性 : $\log$ は単調増加関数なので、$L(Q)$ と $\log L(Q)$ の最大化は同じ $\hat{Q}$ を与える（$\arg\max$ は変わらない）。
A6. 最適化条件式 ∂log L/∂θ = 0 を読み解く
教科書確認済み
最尤推定の核心は次の式です。
$$\frac{\partial \log L(Q)}{\partial \theta} = 0$$
「対数尤度 $\log L(Q)$ をパラメータ $\theta$ について偏微分してゼロに等しくする」—— これは 尤度の山の頂上を探す方程式 です。一つひとつの記号を読み解いていきましょう。
θ
log L
頂上
傾き = 0
傾き > 0
傾き < 0
θ̂
（最尤推定値）
∂ log L / ∂θ = 0
この式を解く = 頂上を探す
図 6. log L(θ) の山型グラフ。最大点（θ̂）では接線の傾き = 0、すなわち ∂log L / ∂θ = 0 が成立する。
A6-1. 記号を 1 つずつ読む
記号 名前 意味
∂ （パーシャル）偏微分記号 変数が複数あるとき、1 個だけ動かしたときの傾き 。通常の微分 d と区別するために使う。 log L(Q) 対数尤度 尤度 L(Q) の自然対数。積 → 和に変換した、扱いやすい形（A5 参照）。 θ （theta、シータ）パラメータ Q を構成する遷移速度定数（α, β, γ, ... など）の総称。複数あるとき全体を θ と書く。 = 0 ゼロ条件 「傾きがゼロ」= 「山の頂上 or 谷底」。尤度は山型なので最大点。
log L(α, β) の等高線図（2 変数の場合）
最大点（∂/∂α=0 かつ ∂/∂β=0）
α
β
α 方向だけ動かす
= ∂ log L / ∂α
β 方向だけ動かす
= ∂ log L / ∂β
図 7. 2 変数関数 log L(α, β) の等高線図。偏微分 ∂/∂α は「α だけ動かしたときの傾き」、∂/∂β は「β だけ動かしたときの傾き」。最大点ではどちらも 0 になる。
A6-2. なぜ「= 0」が最大化を意味するか
1 変数関数 $y = f(x)$ の最大点・最小点では 接線の傾き = 0 （微分の基本）
多変数関数では すべての方向への偏微分 = 0 が必要条件
$\log L(Q)$ は $L(Q)$ と同じ最大点をもつ（$\log$ は単調増加、A5 参照）
だから「$\log L$ の偏微分 = 0」を解けば 最尤推定値 $\hat{Q}$ が見つかる
A6-3. 「数値最適化」とは — 手では解けないから
$\partial \log L / \partial \theta = 0$ は通常 手で解けません 。Q が複雑な行列指数（$e^{Qt}$）を含むため、連立方程式の解析解が存在しないからです。そこでコンピュータで反復計算して近似解を求めます。
手法 仕組み
勾配降下法 （gradient descent）現在地での傾き（勾配）の方向に少しずつパラメータを動かす。シンプルだが収束が遅い場合も。 Newton-Raphson 法 1 階微分に加え 2 階微分（曲率）も使い、より大きなステップで速く収束する。 EM 法（Baum-Welch） 隠れマルコフモデル専用。隠れ状態の期待値と Q を交互に更新（E ステップ・M ステップ）。
勾配降下法による数値最適化のイメージ
θ
log L
初期値
θ̂（収束）
各ステップ: 傾きを計算 → 上り方向へ移動 → 繰り返す → 頂上で ∂ log L/∂θ ≈ 0
図 8. 勾配降下法による数値最適化。初期値から出発し、各ステップで ∂log L/∂θ を計算してパラメータを更新、頂上（θ̂）に収束する。
A6-4. イオンチャネルの言葉で
「dwell time の系列が観測された。どの遷移速度定数 α, β なら、このデータがもっとも『ありえそう』に出るか？」をコンピュータの反復計算で探す。
具体例: 2 状態モデル（Closed ⇌ Open）
Q 行列のパラメータ: $\theta = (\alpha, \beta)$（$\alpha$: 開口速度、$\beta$: 閉口速度）
観測された平均開口時間が 5 ms なら → MLE は $\alpha \approx 200\,\text{s}^{-1}$ 付近
実際には多次元の $\theta$ に対して $\partial \log L / \partial \alpha = 0$ かつ $\partial \log L / \partial \beta = 0$ を同時に満たす $(\hat{\alpha}, \hat{\beta})$ を数値探索する
後半 — 数学者向け厳密版
B1. 形式的設定
教科書確認済み
連続時間マルコフ連鎖（CTMC）$X(t)$ を仮定する（NSM-001 参照）。状態空間 $\mathcal{S}$ は有限で、開状態集合 $\mathcal{O}$ と閉状態集合 $\mathcal{C}$ に分割される: $\mathcal{S} = \mathcal{O} \cup \mathcal{C}$。
記号 意味
$Q$ 遷移速度行列（$n \times n$、行和ゼロ、非対角非負） $Q_{OO}$ 開状態間の遷移速度の部分行列（$|\mathcal{O}| \times |\mathcal{O}|$） $Q_{OC}$ 開状態 → 閉状態への遷移速度（$|\mathcal{O}| \times |\mathcal{C}|$） $Q_{CC}, Q_{CO}$ 同様に閉状態ブロック $t_O^{(k)}, t_C^{(k)}$ $k$ 番目の開・閉 dwell time $\boldsymbol{\phi}_O$ 開状態突入時の条件付き初期占有分布（行ベクトル）
観測データは dwell time の交互系列: $t_O^{(1)}, t_C^{(1)}, t_O^{(2)}, t_C^{(2)}, \ldots, t_O^{(N)}, t_C^{(N)}$ を想定する（単純化として burst 内観測を考える）。
B2. dwell time の確率密度
教科書確認済み — Colquhoun & Hawkes (1981) / (1982)
単純 2 状態モデルの場合
状態 O（1 つ）から状態 C（1 つ）へ速度 $q_{OC}$ で遷移する場合、開状態 dwell time は指数分布に従う:
$$f_O(t) = q_{OC} \cdot e^{-q_{OC} t}, \quad t \ge 0$$
一般の多状態モデル — Colquhoun-Hawkes 形式
標準的だが手元で要確認
開状態集合内に複数の状態がある場合、開状態 dwell time の密度は行列指数を使った形になる:
開状態 dwell time の密度（Colquhoun-Hawkes）
$$f_O(t) = \boldsymbol{\phi}_O \, e^{Q_{OO} t} \, (-Q_{OO}) \, \mathbf{1}$$
ここで:
$\boldsymbol{\phi}_O$: 開状態に突入したときの（条件付き）初期占有分布（行ベクトル）
$e^{Q_{OO} t}$: 開状態内に「まだ滞在している」確率行列（$Q_{OO}$ は負定値の部分行列）
$(-Q_{OO}) \mathbf{1}$: 開状態から閉状態への脱出速度ベクトル（$Q_{OC} \mathbf{1}$ と等価）
$Q_{OO}$ の固有値はすべて負の実部を持つ（CTMC の既約性・有限状態空間から）ため、$e^{Q_{OO} t}$ は $t \to \infty$ で $0$ に収束し、$f_O(t)$ は正しく確率密度になる。固有値分解すると指数の和（mixture）の形が現れる（NSM-003 参照）。
B3. 尤度関数の構造
教科書確認済み — Colquhoun & Hawkes (1982), Qin et al. (1996)
独立観測の積
burst 内の観測が独立と仮定すると、$N$ 個の開閉サイクルに対する尤度は:
$$L(Q) = \prod_{k=1}^{N} f_O\!\left(t_O^{(k)}\right) \cdot f_C\!\left(t_C^{(k)}\right)$$
連続観測列の尤度（行列積形式）
開閉状態の切り替わりを通じた条件付き確率の連鎖を明示的に書くと:
$$L(Q) = \boldsymbol{\phi}_O \; e^{Q_{OO} t_O^{(1)}} \, Q_{OC} \; e^{Q_{CC} t_C^{(1)}} \, Q_{CO} \; e^{Q_{OO} t_O^{(2)}} \, Q_{OC} \; \cdots \; \mathbf{1}$$
この形式はブロック行列の積として書けるため、「forward algorithm」（Baum-Welch の特殊ケース）で効率的に計算できる。
尤度 L(Q) の Q パラメータ依存性（模式図）
Q のスカラーパラメータ θ（例: q₁₂）
尤度 L(θ)
Q̂ = arg max L(Q)
（最尤推定値）
θ̂
接線の傾き = 0
（∂L/∂θ = 0）
図 4. 尤度 $L(Q)$ のパラメータ $\theta$（$Q$ の各エントリ）依存性の模式図。最大点 $\hat{Q}$ が最尤推定値。実際は多次元パラメータ空間での最大化になる。
B4. 対数尤度
教科書確認済み
対数尤度 $\ell(Q) := \log L(Q)$ は積を和に変換する:
$$\ell(Q) = \sum_{k=1}^{N} \left[ \log f_O\!\left(t_O^{(k)}\right) + \log f_C\!\left(t_C^{(k)}\right) \right]$$
対数尤度が有用な理由:
計算安定性 : $10^{-300}$ 以下の積が生じないよう log 変換する（アンダーフロー防止）スコア関数 : $\nabla_\theta \ell(Q)$ が「スコア」として最適化に使われるフィッシャー情報量 : $I(\theta) = -\mathbb{E}[\nabla^2 \ell(Q)]$ が漸近分散の逆数を与える
対数尤度 ℓ(θ) = log L(θ) と最尤推定の最適化条件
対数尤度 ℓ(θ)
θ̂
∂ℓ/∂θ = 0
スコア関数 ∂ℓ/∂θ（= 微分）
0
θ̂（スコア = 0）
一階条件: ∂ℓ/∂θ = 0
θ
θ
図 5. 左: 対数尤度 $\ell(\theta) = \log L(\theta)$ の概形（山型）。右: そのスコア関数 $\partial \ell / \partial \theta$。スコアが 0 を横切る点が一階条件 = 最尤推定値 $\hat{\theta}$。
B5. 最尤推定の最適化
標準的だが手元で要確認（数値手法の詳細）
一階最適化条件（スコア方程式）
Q のパラメータ $\theta = (q_{12}, q_{21}, \ldots)$ に対し:
$$\frac{\partial \ell(Q)}{\partial \theta_j} = 0 \quad \text{for all } j$$
解析的な解は一般に存在しないため、数値最適化が必要。
数値最適化手法
手法 特徴 適用場面
勾配上昇法（gradient ascent） 実装簡単・収束は遅い 凸問題・初期探索 Newton-Raphson 法 2 次収束・ヘッセ行列が必要 小〜中規模問題 EM 法（Baum-Welch） 隠れ状態 HMM の特殊版・単調増加を保証 隠れ状態が多い場合 BFGS（準 Newton 法） ヘッセ行列を近似・実用的に高速 一般的な MLE 実装
制約条件
Q 行列の条件（$q_{ij} \ge 0$（$i \ne j$）、$q_{ii} = -\sum_{j \ne i} q_{ij}$）を守りながら最適化する必要がある。独立パラメータは非対角の非負エントリ $\{q_{ij}\}_{i \ne j}$ のみ（行和ゼロ条件で $q_{ii}$ は決まる）。
B6. 統計的性質
教科書確認済み（漸近論の標準結果）
漸近正規性
観測数 $N \to \infty$ のとき（規則的な識別可能なモデルで）:
$$\sqrt{N}\left(\hat{\theta} - \theta^*\right) \xrightarrow{d} \mathcal{N}\!\left(\mathbf{0},\; I(\theta^*)^{-1}\right)$$
ここで $\theta^*$ は真のパラメータ、$I(\theta)$ はフィッシャー情報行列:
$$I(\theta) = -\mathbb{E}\!\left[\frac{\partial^2 \ell(Q)}{\partial \theta_i \partial \theta_j}\right]$$
クラメール・ラオ下界
任意の不偏推定量 $\tilde{\theta}$ に対して: $\operatorname{Var}(\tilde{\theta}) \ge I(\theta)^{-1}$（クラメール・ラオ下界）。MLE はこの下界を漸近的に達成する（漸近有効性）。
一致性
$N \to \infty$ で $\hat{\theta} \xrightarrow{p} \theta^*$（確率収束）。ただし識別可能性（同定可能性）が前提。
識別可能性（Identifiability）について
複数の異なる Q 行列が全く同じ尤度を生成する（=識別不可能）場合、MLE は一意でない。CTMC の同定可能性は非自明な問題であり、状態数・制約構造に依存する。
NSM-005 （逆問題の定式化）を参照。
B7. 実装の数値安定性
標準的だが実装の詳細は手元確認を要する
行列指数の計算
尤度計算の核心は $e^{Q_{OO} t}$ の計算。効率的な手法:
対角化法 : $Q_{OO} = V \Lambda V^{-1}$ のとき $e^{Q_{OO}t} = V e^{\Lambda t} V^{-1}$（NSM-022 参照）Scaling-and-Squaring 法 : $e^{Qt} = (e^{Qt/2^k})^{2^k}$ として計算する一般的アルゴリズム（MATLAB の expm が使用）Padé 近似 : 行列指数の有理関数近似（scaling-and-squaring と組み合わせる）
対数尤度のスケール問題
長い観測系列では $L(Q)$ が極めて小さくなる。解決策:
対数尤度の逐次加算: $\ell = \sum_k \log f(t_k|Q)$
行列積の対数変換（log-sum-exp トリック）
スケーリング因子の明示的管理（Rabiner の HMM アルゴリズムに類似）
勾配計算
数値微分（有限差分）か自動微分（JAX 等）を使う。解析的な $\partial e^{Qt} / \partial q_{ij}$ の計算は行列微分の理論を要する（NSM-022 参照）。
参考文献
Colquhoun, D. & Hawkes, A. G. (1981). On the stochastic properties of single ion channels. Proceedings of the Royal Society of London B , 211, 205–235. → 単一チャネル dwell time の確率密度（行列指数形式）の原著定式化。 教科書確認済み Colquhoun, D. & Hawkes, A. G. (1982). On the stochastic properties of bursts of single ion channel openings and of clusters of bursts. Philosophical Transactions of the Royal Society of London B , 300, 1–59. → burst 内の尤度計算と Kolmogorov 方程式の逆問題への応用。 教科書確認済み Colquhoun, D., Hawkes, A. G. & Srodzinski, K. (1996). Joint distributions of apparent open and shut times of single-ion channels and maximum likelihood fitting of mechanisms. Philosophical Transactions of the Royal Society of London A , 354, 2555–2590. → dead time 補正を含む尤度計算（HJC theory）の最重要論文。 教科書確認済み Qin, F., Auerbach, A. & Sachs, F. (1996). Estimating single-channel kinetic parameters from idealized patch-clamp data containing missed events. Biophysical Journal , 70, 264–280. → 見落とし事象を含む MLE の定式化。 標準的だが手元未確認 Cox, D. R. & Hinkley, D. V. (1974). Theoretical Statistics. Chapman & Hall. → MLE の漸近理論（一致性・漸近正規性・クラメール・ラオ下界）の標準的教科書。 教科書確認済み（定理は標準的）
作成: 2026-05-26 / NSM-024 / プロジェクト: Mathematics in Neuroscience プレゼン準備