状態空間のMCMC（カルマン×ベイズ）

🎓 レベル：発展　|　重要度：A（必須）

📎 前提：カルマンフィルタ・ベイズ構造時系列　|　数理：メトロポリスヘイスティングス（ベイズ・MCMC）・ベイズ更新と逐次推論（ベイズ・逐次更新）

要点（BLUF）

線形ガウスなら、 $Q,H$ が既知のときカルマンフィルタが状態の事後を解析的に与えます（カルマンフィルタ）。問題は分散パラメータ $Q,H$ 自体が未知な場合。
このとき $Q,H$ に事前分布を置き、MCMC で事後 $p(Q,H\mid y)$ をサンプリングします。掛け算は「カルマンが返す周辺尤度 × 事前 → 事後」。
役割分担が肝：カルマンは状態（逐次更新・解析的）、MCMC はパラメータ（事後探索）。ローカルレベルとローカルトレンドモデルと同じ系列で、statsmodels の最尤の点推定が、ベイズ事後のモードに当たることを数値と図で確かめます。

1. カルマンは状態を、MCMC はパラメータを

カルマンフィルタで見たとおり、状態空間が線形ガウスで $Q,H$ が既知なら、状態の事後 $p(\alpha_t\mid y_{1:t})$ は正規分布のまま閉じて、漸化式で逐次に出ます。サンプリングは要りません。

ところが実務では $Q,H$ は未知です。ローカルレベルとローカルトレンドモデルはこれを最尤で点推定しました。その尤度こそ、カルマンフィルタが副産物として返す周辺尤度（予測誤差分解）：

p(y_{1:n}\mid Q,H) = \prod_{t=1}^{n} N\!\big(v_t;\, 0,\, F_t\big)

ここで $v_t=y_t - Z a_{t\mid t-1}$ はイノベーション、 $F_t$ はその分散——どちらも $Q,H$ を与えればカルマン漸化式から出ます（カルマンフィルタ）。状態 $\alpha_{1:n}$ はカルマンが解析的に積分消去してくれるので、尤度は $Q,H$ だけの関数になります。statsmodels の .fit() はこれを最大化して $\hat Q,\hat H$ を1点で返す。

ベイズはここに事前を掛けます：

p(Q,H\mid y_{1:n}) \;\propto\; \underbrace{p(y_{1:n}\mid Q,H)}_{\text{カルマンの周辺尤度}} \;\times\; \underbrace{p(Q)\,p(H)}_{\text{事前}}

右辺は $Q,H$ の2次元（解析的に解けない）なので、MCMC でサンプリングします（メトロポリスヘイスティングス・ハミルトニアンモンテカルロとNUTS）。役割分担が明快です。

flowchart TB
    MC["MCMC：パラメータ Q,H の事後を探索（外側）"] -->|"Q,H を提案"| KF["カルマン：状態を積分消去し周辺尤度 p(y|Q,H) を返す（内側・解析的）"]
    KF -->|"尤度 × 事前 = 事後の値"| MC
    PRI["事前 p(Q), p(H)"] --> MC
    MC --> POST["事後 p(Q,H|y)：平均・中央値・モード・信用区間"]

カルマン（内側・解析的・逐次）： $Q,H$ を1組与えるたびに、状態を逐次更新しながら尤度を返す。状態の不確実性はここで閉じて扱われます。
MCMC（外側・サンプリング）： $Q,H$ の空間を歩き回り、事後分布の形（モード・歪み・区間）を描く。

実装上の補足：本コードは PyMC で状態 $\mu_t$ とパラメータ $Q,H$ を同時にサンプルします（データ拡大法）。線形ガウスなら状態をカルマンで積分消去しても、同時サンプルから $Q,H$ の周辺だけ取り出しても、 $Q,H$ の事後は一致します。概念は「カルマンが尤度、MCMC がパラメータ」と捉えるのが見通し良いです。

2. 最尤の点 vs ベイズの事後分布

コード①：Q,H をベイズ推定し、最尤の点と重ねる

ローカルレベルとローカルトレンドモデルと同じローカルレベル系列（真値 $Q=0.5,\ H=4.0$ ）を作り、(A) statsmodels で最尤の点推定、(B) pymc で $Q,H$ の事後分布を出します。事後の平均・中央値・最頻（KDE）・95%区間を計算し、最尤の点が事後のどこに落ちるかを確かめます。水準のランダムウォークは非中心化（ベイズ構造時系列）。Windows は cores=1、数値は idata.posterior から取得（PyMC 6.0.1・ArviZ 1.2.0）。

import warnings
warnings.simplefilter("ignore")
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
import pymc as pm
import pytensor.tensor as pt
import arviz as az
from scipy.stats import gaussian_kde
from statsmodels.tsa.statespace.structural import UnobservedComponents

# 04-03 と同じローカルレベル系列（水準=RW, 観測=水準+ノイズ）
np.random.seed(2)
n = 120
Q_true, H_true = 0.5, 4.0
alpha = np.zeros(n); alpha[0] = 30.0
for t in range(1, n):
    alpha[t] = alpha[t-1] + np.random.normal(0, np.sqrt(Q_true))
y = alpha + np.random.normal(0, np.sqrt(H_true), n)

# (A) statsmodels の最尤推定（カルマン尤度を最大化＝点推定）
res = UnobservedComponents(y, level="local level").fit(disp=False)
p = dict(zip(res.param_names, res.params))
Q_mle, H_mle = p["sigma2.level"], p["sigma2.irregular"]
print(f"最尤(MLE)  Q = {Q_mle:.3f}   H = {H_mle:.3f}")

# (B) ベイズ：Q,H に事前を置き MCMC で事後分布（状態も同時にサンプル＝状態空間MCMC）
with pm.Model() as model:
    sigma_level = pm.HalfNormal("sigma_level", sigma=2.0)   # sqrt(Q) の事前
    sigma_obs = pm.HalfNormal("sigma_obs", sigma=5.0)       # sqrt(H) の事前
    mu0 = pm.Normal("mu0", mu=30.0, sigma=10.0)
    z = pm.Normal("z", 0.0, 1.0, shape=n-1)                 # 水準RWの非中心化
    mu = pm.Deterministic("mu", mu0 + pt.concatenate([[0.0], pt.cumsum(sigma_level*z)]))
    pm.Normal("y_obs", mu=mu, sigma=sigma_obs, observed=y)
    Q = pm.Deterministic("Q", sigma_level**2)
    H = pm.Deterministic("H", sigma_obs**2)
    idata = pm.sample(draws=1000, tune=1000, chains=2, cores=1,
                      random_seed=42, target_accept=0.97, progressbar=False)

print(f"発散数 = {int(idata.sample_stats['diverging'].values.sum())}")
print(f"R-hat  Q = {float(az.rhat(idata)['Q'].values):.3f}  H = {float(az.rhat(idata)['H'].values):.3f}")

Qs = idata.posterior["Q"].values.ravel()
Hs = idata.posterior["H"].values.ravel()
def mode_kde(s):
    g = np.linspace(s.min(), s.max(), 500)
    return g[gaussian_kde(s)(g).argmax()]
print("ベイズ事後（真値 Q=0.5, H=4.0）:")
print(f"  Q  平均={Qs.mean():.3f}  中央値={np.median(Qs):.3f}  最頻(KDE)={mode_kde(Qs):.3f}  "
      f"95%区間=[{np.percentile(Qs,2.5):.3f}, {np.percentile(Qs,97.5):.3f}]")
print(f"  H  平均={Hs.mean():.3f}  中央値={np.median(Hs):.3f}  最頻(KDE)={mode_kde(Hs):.3f}  "
      f"95%区間=[{np.percentile(Hs,2.5):.3f}, {np.percentile(Hs,97.5):.3f}]")
print(f"  MLE は事後のどこか: Q_mle={Q_mle:.3f}（最頻 {mode_kde(Qs):.3f} の近く）  "
      f"H_mle={H_mle:.3f}（最頻 {mode_kde(Hs):.3f} の近く）")

# 図：Q,H の事後分布に MLE(点推定)と真値を重ねる
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
for a, s, mle, tv, name in [(ax[0], Qs, Q_mle, Q_true, "Q（状態ノイズ分散）"),
                            (ax[1], Hs, H_mle, H_true, "H（観測ノイズ分散）")]:
    a.hist(s, bins=40, density=True, color="C0", alpha=0.5, label="ベイズ事後")
    g = np.linspace(s.min(), s.max(), 300)
    a.plot(g, gaussian_kde(s)(g), color="C0", lw=1.5)
    a.axvline(mle, color="C3", lw=2, ls="-", label=f"MLE 点推定={mle:.2f}")
    a.axvline(tv, color="k", lw=1.5, ls="--", label=f"真値={tv:.2f}")
    a.axvline(np.median(s), color="C1", lw=1.5, ls=":", label=f"事後中央値={np.median(s):.2f}")
    a.set_xlabel(name); a.set_ylabel("事後密度"); a.legend(fontsize=8)
ax[0].set_title("分散パラメータ Q の事後分布")
ax[1].set_title("分散パラメータ H の事後分布")
plt.tight_layout(); plt.show()

出力：

最尤(MLE)  Q = 0.438   H = 5.215
発散数 = 0
R-hat  Q = 1.003  H = 1.003
ベイズ事後（真値 Q=0.5, H=4.0）:
  Q  平均=0.580  中央値=0.501  最頻(KDE)=0.381  95%区間=[0.161, 1.472]
  H  平均=5.288  中央値=5.225  最頻(KDE)=5.197  95%区間=[3.711, 7.188]
  MLE は事後のどこか: Q_mle=0.438（最頻 0.381 の近く）  H_mle=5.215（最頻 5.197 の近く）

出力の意味：核心は最後の3行です。最尤の点推定 $Q_{\text{mle}}=0.438$ は、ベイズ事後の最頻（モード） $0.381$ の近くに落ち、 $H_{\text{mle}}=5.215$ も事後モード $5.197$ にほぼ重なります。これは偶然でなく必然——弱情報事前のもとで事後 $\propto$ 尤度 $\times$ 事前 $\approx$ 尤度なので、事後を最大化する点（モード＝MAP）が尤度を最大化する点（MLE）にほぼ一致するからです。「最尤はベイズ事後のモードの近似」が数値で見えました。

それ以上に重要なのは、ベイズが点でなく分布を返すこと。分散の事後は右に歪みます（ $Q$ で平均 $0.580$ > 中央値 $0.501$ > 最頻 $0.381$ ）。最尤の1点だけ見ると「左の裾の方の値」を握っていることになり、 $Q$ の不確実性（95%区間 $[0.161,\,1.472]$ は7倍以上の幅）を見落とします。真値はどうか—— $Q=0.5$ は事後中央値 $0.501$ にほぼ命中、 $H=4.0$ は最尤・ベイズとも $5.2$ 付近と高めに出ましたが（このデータ実現が偶然 $H$ を高く支持した）、それでも真値 $4.0$ は95%区間 $[3.711,\,7.188]$ の内側に収まっています。収束も良好（発散 0、R-hat $1.003$ ）。

3. なぜ事後分布を持つと嬉しいか

最尤の点推定 $\hat Q,\hat H$ で予測区間を作ると、「 $Q,H$ がぴったり $\hat Q,\hat H$ だ」と決めつけた区間になります。だが上で見たとおり $Q$ の事後は7倍幅の不確実性を持つ。ベイズはこのパラメータ不確実性を予測に伝播させます（ベイズ構造時系列の事後予測）。短い系列・分散が識別しにくいモデルほど、この差は無視できません。最尤区間は楽観的に、ベイズ区間は正直に広がります。

また、 $Q,H$ に事前を入れられるのも実利です。「観測ノイズはこの程度のはず」という装置仕様の知識を $H$ の事前に入れれば、データが少なくても破綻しにくい。最尤はこの注入がやりにくい。

4. 数式の直観

カルマン＝尤度マシン、MCMC＝パラメータ探索機：カルマンは「 $Q,H$ をくれれば、状態をうまく消して尤度を返す」装置。MCMC は「その尤度×事前を頼りに $Q,H$ 空間を歩く」装置。役割が直交しているので、組み合わせて状態空間のベイズ推定になります。
最尤はほぼ事後モード（MAP）：事前が平坦に近ければ $\arg\max p(\theta\mid y)=\arg\max p(y\mid\theta)$ 。だから最尤はベイズの特別な場合（MAP に無情報事前）と読めます（点推定と損失関数）。違いは「点を返すか、分布を返すか」。
分散の事後は右に歪む：分散は正の量で、上に裾を引きます。だから平均 > 中央値 > 最頻。「代表値をどれで報告するか」で値が変わるので、ベイズでは区間ごと示すのが誠実です（事後分布の要約）。
逐次（カルマン）と探索（MCMC）の時間軸が違う：カルマンは時間方向に1パスで状態を更新（ベイズ更新と逐次推論）。MCMC はパラメータ空間を反復探索。前者は $O(n)$ 、後者はサンプル数×尤度評価。線形ガウスでは両者を入れ子にするのが最効率です。

⚠️ よくある誤解・落とし穴

「カルマンとMCMCはどちらか一方」ではない：線形ガウスの状態空間ベイズでは両方使います——カルマンで状態を消して尤度を作り、MCMC でパラメータを回す。対立でなく分業です。
「最尤＝事後モードは常に成立」ではない：事前が強い、または尤度が歪む・多峰だと、MLE と MAP はずれます。本例で一致したのは弱情報事前＋単峰の尤度だから。事前感度は要確認（事前の選び方と感度分析）。
「点推定が真値に近ければ十分」ではない： $H$ は最尤・ベイズとも真値 $4.0$ から $5.2$ へずれました。単一実現では珍しくない。区間が真値を含むかで評価すべきで、点だけ見ると過信します。
「線形ガウスでなくてもこの手が効く」ではない：状態が非ガウス・遷移が非線形だと、カルマンは厳密尤度を返せません。**粒子フィルタ（逐次モンテカルロ）**で尤度を近似し、その上で MCMC（粒子MCMC）や、状態を直接サンプルする方法が要ります（本章の範囲外）。確率的ボラティリティ（SV）モデルが典型例（第6章ボラティリティ目次）。
「状態を同時サンプルすれば速い」ではない：データ拡大法は状態 $n$ 本ぶん次元が増え、 $\sigma$ と増分の funnel で詰まりやすい。非中心化や、状態をカルマンで積分消去して低次元の $Q,H$ だけ MCMC する設計が安定・高速なことが多い（階層モデルの実例と再パラメータ化）。