カルマンフィルタ｜時系列分析

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：状態空間モデルの枠組み　|　接続：ベイズ更新と逐次推論（ベイズ・逐次更新）・正規正規モデル（ベイズ・正規共役）

要点（BLUF）

カルマンフィルタ＝状態空間モデルの隠れ状態を、観測が来るたびに逐次更新して推定するアルゴリズム。各時点で「予測（事前）→ 観測で更新（事後）」の2ステップを回します。
予測ステップは状態方程式で事前分布 $(a_{t|t-1},P_{t|t-1})$ を作り、更新ステップは観測 $y_t$ とカルマンゲイン $K_t$ で事後分布 $(a_t,P_t)$ に補正します。
これは ベイズの逐次更新（ベイズ更新と逐次推論）のガウス版——事前×尤度→事後の正規共役（正規正規モデル）を時間方向に繰り返すのと同型です。自前実装が真の隠れ水準を相関 0.914 で追え、statsmodels とも一致することを確かめます。

1. 予測ステップと更新ステップ

カルマンフィルタは、状態の事後分布を正規分布 $N(a_t, P_t)$ で持ち回ります（線形＋ガウスなら事後も常に正規）。 $a_t$ が状態の推定値（平均）、 $P_t$ がその不確実性（分散）。各時点 $t$ で2段階：

予測ステップ（事前分布）——前時点の事後 $(a_{t-1},P_{t-1})$ を状態方程式で1歩進めます。

a_{t\mid t-1} = T\,a_{t-1}, \qquad P_{t\mid t-1} = T\,P_{t-1}\,T' + R\,Q\,R'

平均は遷移行列 $T$ で進め、分散は $T$ で変換したうえ状態ノイズ $Q$ を足して増やす（時間が進むと不確実になる）。

更新ステップ（事後分布）——観測 $y_t$ が届いたら、予測とのズレで補正します。

v_t = y_t - Z\,a_{t\mid t-1}, \qquad F_t = Z\,P_{t\mid t-1}\,Z' + H

K_t = P_{t\mid t-1}\,Z'\,F_t^{-1}

a_t = a_{t\mid t-1} + K_t\,v_t, \qquad P_t = (I - K_t Z)\,P_{t\mid t-1}

ここで $v_t$ はイノベーション（予測誤差）、 $F_t$ はその分散、 $K_t$ がカルマンゲイン。「予測 $a_{t|t-1}$ に、観測のズレ $v_t$ をゲイン $K_t$ の重みで足し込む」のが心臓部です。更新後は分散 $P_t$ が予測時より小さくなる（観測で情報が増えた）。

ローカルレベル（ $Z=T=R=1$ ）ならスカラーになり、手で追えます：

a_{t\mid t-1}=a_{t-1},\quad P_{t\mid t-1}=P_{t-1}+Q,\quad F_t=P_{t\mid t-1}+H,\quad K_t=\frac{P_{t\mid t-1}}{P_{t\mid t-1}+H}

flowchart LR
    P0["前時点の事後 a_{t-1}, P_{t-1}"] --> PR["予測：a_pred=T a_{t-1}, P_pred=T P T' + R Q R'"]
    PR --> UP["更新：イノベーション v_t と ゲイン K_t で補正"]
    UP --> P1["今時点の事後 a_t, P_t"]
    P1 --> PR

2. カルマン＝ベイズ更新の連続版

カルマンフィルタの1ステップは、ベイズの逐次更新そのものです（ベイズ更新と逐次推論）。

事前分布（予測ステップ）： $\alpha_t \mid y_{1:t-1} \sim N(a_{t\mid t-1},\, P_{t\mid t-1})$
尤度（観測モデル）： $y_t \mid \alpha_t \sim N(Z\alpha_t,\, H)$
事後分布（更新ステップ）： $\alpha_t \mid y_{1:t} \sim N(a_t,\, P_t)$

正規の事前×正規の尤度＝正規の事後、という正規共役（正規正規モデル）を、時間方向に毎ステップ繰り返しているだけ。カルマンゲイン $K_t$ は共役更新の「事前と観測をどの比率で混ぜるか」の重みに相当します。前時点の事後が次時点の（予測を挟んだ）事前になる——この逐次性こそ、カルマンが大量データでも一定メモリで回る理由です。

\underbrace{N(a_{t\mid t-1},P_{t\mid t-1})}_{\text{事前}} \;\times\; \underbrace{N(Z\alpha_t,H)}_{\text{尤度}} \;\propto\; \underbrace{N(a_t,P_t)}_{\text{事後}}

コード①：1次元ローカルレベルにカルマンフィルタを自前実装

真の隠れ水準 $\mu_t$ （ランダムウォーク）を仕込み、観測 $y_t=\mu_t+\varepsilon_t$ だけから水準を推定します。predict/update を for ループで素朴に実装し、推定が真の隠れ水準を追えるか（相関・RMSE）を確かめます。あわせて、観測ノイズ $H$ を変えると定常カルマンゲインがどう動くかを見ます。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用

# 真のローカルレベル：隠れ水準=ランダムウォーク, 観測=水準+ノイズ
np.random.seed(1)
n = 150
Q_true, H_true = 0.3, 5.0
alpha = np.zeros(n); alpha[0] = 0.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)

# カルマンフィルタ（スカラー版）を自前実装。T=1, R=1 のローカルレベル
Q, H = Q_true, H_true            # ここでは真の分散を既知として与える
a, P = y[0], 1e6                 # 初期状態の平均と分散（分散大＝ほぼ無情報事前）
a_filt = np.zeros(n)             # 事後（フィルタ）状態の保存先
gain = np.zeros(n)               # カルマンゲインの記録
for t in range(n):
    # --- 予測（事前）：T=1 なので平均はそのまま、分散に Q を足す ---
    a_pred = a                  # a_{t|t-1} = T a_{t-1}
    P_pred = P + Q              # P_{t|t-1} = T P T' + R Q R'
    # --- 更新（事後）：観測 y_t で補正 ---
    v = y[t] - a_pred           # イノベーション（予測誤差）
    F = P_pred + H              # イノベーション分散 = Z P Z' + H
    K = P_pred / F              # カルマンゲイン
    a = a_pred + K * v          # a_{t|t} = a_{t|t-1} + K v
    P = (1 - K) * P_pred        # P_{t|t} = (1 - K Z) P_{t|t-1}
    a_filt[t] = a; gain[t] = K

# フィルタ推定が隠れ水準を追えたか（初期の過渡を除いて評価）
warm = 10
corr = np.corrcoef(a_filt[warm:], alpha[warm:])[0, 1]
rmse_filt = np.sqrt(np.mean((a_filt[warm:] - alpha[warm:])**2))
rmse_obs = np.sqrt(np.mean((y[warm:] - alpha[warm:])**2))
print(f"カルマンゲインの収束値 K = {gain[-1]:.3f}")
print(f"フィルタ推定と真の水準の相関 = {corr:.3f}")
print(f"RMSE  フィルタ vs 真水準 = {rmse_filt:.3f}")
print(f"RMSE  生観測 vs 真水準   = {rmse_obs:.3f}（フィルタはこれより小さい）")

# ゲインが観測ノイズで変わることを確認：H を大きくするとゲインは小さくなる
def steady_gain(Q, H):
    a, P = 0.0, 1e6
    for _ in range(500):
        P_pred = P + Q
        K = P_pred / (P_pred + H)
        P = (1 - K) * P_pred
    return K
print(f"定常ゲイン  H=1   -> K={steady_gain(0.3, 1.0):.3f}（観測を信じる）")
print(f"定常ゲイン  H=5   -> K={steady_gain(0.3, 5.0):.3f}")
print(f"定常ゲイン  H=50  -> K={steady_gain(0.3, 50.0):.3f}（観測を信じない）")

# 自前カルマンフィルタが隠れ水準を追えたかを図で確認
fig, ax = plt.subplots(figsize=(9, 4.5))
ax.plot(y, color="0.7", lw=0.8, marker="o", ms=2.5, label="観測 y_t（ノイジー）")
ax.plot(alpha, color="k", ls="--", lw=1.5, label="真の隠れ水準")
ax.plot(a_filt, color="C0", lw=1.8, label="自前カルマンフィルタ推定")
ax.set_xlabel("時点 t"); ax.set_ylabel("値"); ax.legend()
ax.set_title("自前カルマンフィルタが観測のノイズを均し隠れ水準を追う")
plt.tight_layout(); plt.show()

出力：

カルマンゲインの収束値 K = 0.217
フィルタ推定と真の水準の相関 = 0.914
RMSE  フィルタ vs 真水準 = 0.950
RMSE  生観測 vs 真水準   = 2.231（フィルタはこれより小さい）
定常ゲイン  H=1   -> K=0.418（観測を信じる）
定常ゲイン  H=5   -> K=0.217
定常ゲイン  H=50  -> K=0.075（観測を信じない）

出力の意味：フィルタ推定は真の隠れ水準と相関 0.914、RMSE は $0.950$ 。生の観測（RMSE $2.231$ ）より2倍以上正確に真の水準を当てています——カルマンが観測のノイズを均して信号を取り出した証拠です。後半のゲイン実験が核心：観測ノイズ $H$ を $1\to5\to50$ と上げると、ゲインは $0.418\to0.217\to0.075$ と小さくなります。 $K_t=P_{t|t-1}/(P_{t|t-1}+H)$ なので、 $H$ が大きい（観測が信用できない）ほど $K_t$ は小さく、新しい観測 $v_t$ を控えめにしか取り込まない＝予測（過去の蓄積）を重視する。逆に $H$ が小さければ $K_t\to1$ で観測をほぼそのまま信じる。「測定が荒いほど観測を割り引く」という、人の直観どおりの重み付けを式が自動でやってくれます。

コード②：自前実装と statsmodels の一致を確認

同じ系列・同じ真の分散で、statsmodels の UnobservedComponents(y, level="local level") を分散を固定して走らせ（.smooth([H, Q])）、その filtered_state が自前実装と一致するかを最大絶対差で確かめます。一致すれば「自前のループは本物のカルマンフィルタだ」と裏が取れます。

import warnings
warnings.simplefilter("ignore")
import numpy as np
from statsmodels.tsa.statespace.structural import UnobservedComponents

# 04-02 コード1 と同じ系列・同じ真の分散を再生成
np.random.seed(1)
n = 150
Q_true, H_true = 0.3, 5.0
alpha = np.zeros(n); alpha[0] = 0.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)

# 自前実装（真の分散を既知として与える）
Q, H = Q_true, H_true
a, P = y[0], 1e6
a_self = np.zeros(n)
for t in range(n):
    a_pred, P_pred = a, P + Q
    K = P_pred / (P_pred + H)
    a = a_pred + K * (y[t] - a_pred)
    P = (1 - K) * P_pred
    a_self[t] = a

# statsmodels：分散を固定して同じカルマンフィルタを走らせる（smooth で実行）
mod = UnobservedComponents(y, level="local level")
# パラメータ順は [sigma2.irregular(=H), sigma2.level(=Q)]
res = mod.smooth([H_true, Q_true])
a_sm = res.filtered_state[0]      # フィルタ水準

max_abs_diff = np.max(np.abs(a_self - a_sm))
print(f"パラメータ名 = {res.param_names}")
print(f"自前実装と statsmodels の filtered_state 最大絶対差 = {max_abs_diff:.2e}")
print(f"自前 a[-1]      = {a_self[-1]:.4f}")
print(f"statsmodels[-1] = {a_sm[-1]:.4f}")

出力：

パラメータ名 = ['sigma2.irregular', 'sigma2.level']
自前実装と statsmodels の filtered_state 最大絶対差 = 7.63e-06
自前 a[-1]      = 4.4465
statsmodels[-1] = 4.4465

出力の意味：自前実装と statsmodels のフィルタ水準は最大絶対差 $7.63\times10^{-6}$ ——実質一致です（差は初期分散の与え方の微差のみ。statsmodels は厳密拡散初期化、こちらは $P_0=10^6$ の近似拡散）。param_names から、statsmodels のパラメータ順が [sigma2.irregular(=H),\ sigma2.level(=Q)] だと分かります。たった十数行の predict/update ループが、ライブラリの本格実装と同じ結果を出す——カルマンフィルタの中身が透けて見えたはずです。次ノートローカルレベルとローカルトレンドモデルでは、この $Q,H$ を既知とせず**データから推定（復元）**します。

3. 数式の直観

カルマンゲイン＝信頼の配分： $K_t=P_{t|t-1}/(P_{t|t-1}+H)$ は「予測の不確実 $P_{t|t-1}$ 」対「観測の不確実 $H$ 」の綱引き。予測が頼りない（ $P$ 大）か観測が正確（ $H$ 小）なら $K\to1$ （観測重視）、逆なら $K\to0$ （予測重視）。逆分散加重平均そのものです（正規正規モデル）。
イノベーションが情報の源： $v_t=y_t-Z a_{t|t-1}$ は「予測で説明できなかった分」。モデルが正しければ $v_t$ は白色（無相関・平均0）になり、これが残差診断の基礎です（モデル選択と残差診断の白色性）。
分散は増えて減る：予測で $+Q$ と増え、更新で観測ぶん減る。この増減の釣り合いが定常ゲインに落ち着き、長く回すと $K_t$ は一定値へ収束します（コード①で確認）。

⚠️ よくある誤解

「ゲインは大きいほど良い」ではない： $K_t$ が大きいと観測に振り回され、ノイズまで追ってしまう。 $K_t$ は $Q,H$ から最適に決まる量で、大小に善悪はありません。観測がノイジーなら小さいのが正しい。
「フィルタは未来も使う」ではない：カルマンフィルタは各時点で過去と現在の観測だけ（ $y_{1:t}$ ）を使うオンライン推定。未来の観測まで使って遡って改善するのは**平滑化（スムーザー）**で、別物です（平滑化と欠測補間）。
「初期値で結果が決まる」ではない：初期 $(a_0,P_0)$ は数ステップで洗い流されます（拡散初期化なら特に）。だからコード①で警戒期間 warm を除いて評価しました。
「線形・ガウスでなくても厳密」ではない：カルマンフィルタが厳密最適なのは線形＋ガウスのとき。非線形なら拡張カルマン（EKF）やUKF・粒子フィルタなどの近似が要ります（本章の範囲外）。
「 $Q,H$ は気にしなくていい」ではない：ゲインは $Q,H$ の比で決まるので、これを誤ると追従が速すぎ/遅すぎになります。実務では $Q,H$ を最尤やベイズで推定します（ローカルレベルとローカルトレンドモデル）。