書記が数学やるだけ#939 非ガウス型状態空間モデル
非線形の状態空間モデルについて。
問題
説明
非ガウス型状態空間モデルにおける予測・フィルター・平滑化は以下のように一般化される。
粒子フィルター（逐次モンテカルロ法）は、ベイズ統計学における事後分布の再帰的推定を、数式による厳密解から数値シミュレーションによる近似解で表す手法である。
解答
1. 一期先予測 (One-step Ahead Prediction)
定義: 時刻 t−1 までの観測値に基づき、時刻 t の状態を推定するプロセス。
導出: チャップマン・コルモゴロフ方程式（Chapman-Kolmogorov equation）を用いる。
意味: システムモデル p(xt ∣xt−1 ) を介して、前時刻の事後分布を現時刻へ遷移（伝播）させる。
2. フィルタリング (Filtering)
定義: 最新の観測値 yt を取り込み、一期先予測分布を修正して現在の事後分布を得るプロセス。
導出: ベイズの定理（Bayes' theorem）に基づく。
意味: 観測モデル（尤度） p(yt ∣xt ) によって予測の不確実性を縮小させ、現在の状態に関する最新の知見を構成する。
3. 平滑化 (Smoothing)
定義: 全期間（1 から T まで）の観測値 y1:T を用いて、過去の任意の時刻 t (t<T) の状態を再推定するプロセス。
導出: 固定区間平滑化（Fixed-interval smoothing）では、逆向きの再帰計算を行う。
意味: 未来の観測情報 yt+1:T を過去へ逆伝播させることで、フィルタリング時に生じていた推定のラグや不確実性を低減し、より精緻な状態軌跡を復元する。
非ガウス型状態空間モデルの事例について，まずはコーシー分布ノイズを含むデータについて「線形ガウスを仮定するKF」と「非ガウスを許容するPF」の両方で推定する。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import cauchy, norm
# --- 1. 地震波風データの生成 (Non-Gaussian Data) ---
np.random.seed(7)
T = 300
x_true = np.zeros(T)
y = np.zeros(T)
for t in range(1, T):
# システムノイズ: 95%は微小変動、5%でコーシー分布によるジャンプ
if np.random.rand() > 0.95:
v = cauchy.rvs(loc=0, scale=2.0)
else:
v = np.random.normal(0, 0.02)
x_true[t] = x_true[t-1] + v
y[t] = x_true[t] + np.random.normal(0, 0.5)
# --- 2. 粒子フィルタ (PF) の推定 ---
N = 2000
particles = np.zeros(N)
weights = np.ones(N) / N
x_est_pf = np.zeros(T)
for t in range(T):
# 予測: システムモデルにCauchy分布を仮定
particles += cauchy.rvs(loc=0, scale=0.1, size=N)
# 更新: 観測尤度 (Gaussian)
likelihood = norm.pdf(y[t], loc=particles, scale=0.5)
weights *= likelihood
weights += 1e-300
weights /= np.sum(weights)
# 推定
x_est_pf[t] = np.sum(particles * weights)
# リサンプリング
if 1.0 / np.sum(weights**2) < N/2:
indices = np.random.choice(range(N), size=N, p=weights)
particles = particles[indices]
weights = np.ones(N) / N
# --- 3. 線形カルマンフィルタ (KF) の推定 ---
x_kf = 0.0
P_kf = 1.0
Q_kf = 0.1 # ガウス分布を仮定したプロセスノイズ
R_kf = 0.25 # 観測ノイズ
x_est_kf = np.zeros(T)
for t in range(T):
# 予測
x_pred = x_kf
P_pred = P_kf + Q_kf
# 更新
K = P_pred / (P_pred + R_kf)
x_kf = x_pred + K * (y[t] - x_pred)
P_kf = (1 - K) * P_pred
x_est_kf[t] = x_kf
# --- 4. 可視化 (2つのパネルに分割) ---
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True)
# 上段: 線形カルマンフィルタの限界
axes[0].plot(x_true, color='black', lw=2, label='True Signal')
axes[0].plot(x_est_kf, color='green', lw=1.5, label='Linear KF Estimate')
axes[0].scatter(range(T), y, color='gray', s=5, alpha=0.3, label='Observations')
axes[0].set_title("Failure of Linear Model (Kalman Filter) on Abrupt Changes")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 下段: 粒子フィルタによる高精度推定
axes[1].plot(x_true, color='black', lw=2, label='True Signal')
axes[1].plot(x_est_pf, color='red', lw=1.5, label='Particle Filter Estimate')
axes[1].scatter(range(T), y, color='gray', s=5, alpha=0.3, label='Observations')
axes[1].set_title("Success of Non-Gaussian Model (Particle Filter)")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()線形モデル（カルマンフィルタ）では，地震波のような急激なステップ状の変化（ジャンプ）に対し、推定線がなだらかに追従してしまい、大きなラグが生じる。一方で非ガウス型状態空間モデル（粒子フィルタ）では，コーシー分布をシステムモデルに組み込んでいるため，大きな変化を「あり得る事象」として許容し、垂直に近い形でジャンプに即応する。
次に，観測方程式が2乗であるために正負の判断が困難になる状況で，両者の挙動がいかに異なるかを示す。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# 再現性のためのシード設定
np.random.seed(42)
# --- 1. データの生成 (Highly Non-linear & Periodic Model) ---
T = 100
sigma_v = 1.0
sigma_w = 1.0
x_true = np.zeros(T)
y = np.zeros(T)
x_now = 0.1
for t in range(T):
v = np.random.normal(0, sigma_v)
w = np.random.normal(0, sigma_w)
# 非線形な状態遷移
x_now = x_now/2 + 25*x_now/(1 + x_now**2) + 8*np.cos(1.2*t) + v
x_true[t] = x_now
# 二乗を含む観測方程式
y[t] = (x_now**2)/20 + w
# --- 2. 粒子フィルター (PF) の推定 ---
N = 1000
particles = np.random.normal(0, 5, N)
weights = np.ones(N) / N
x_est_pf = np.zeros(T)
for t in range(T):
# 予測
v = np.random.normal(0, sigma_v, N)
particles = particles/2 + 25*particles/(1 + particles**2) + 8*np.cos(1.2*t) + v
# 更新 (尤度計算)
likelihood = norm.pdf(y[t], loc=(particles**2)/20, scale=sigma_w)
weights *= likelihood
weights += 1e-300
weights /= np.sum(weights)
# 推定 (加重平均)
x_est_pf[t] = np.sum(particles * weights)
# リサンプリング
if 1.0 / np.sum(weights**2) < N/2:
indices = np.random.choice(N, size=N, p=weights)
particles = particles[indices]
weights = np.ones(N) / N
# --- 3. 線形カルマンフィルター (KF) の推定 ---
x_kf = 0.0
P_kf = 10.0
Q_kf = 5.0 # 無理やり追従させるための大きめのシステムノイズ
R_kf = 1.0
x_est_kf = np.zeros(T)
for t in range(T):
# 予測 (線形近似 F=1)
x_pred = x_kf
P_pred = P_kf + Q_kf
# 更新 (線形近似 H=1, y=x+w と仮定)
K = P_pred / (P_pred + R_kf)
x_kf = x_pred + K * (y[t] - x_pred)
P_kf = (1 - K) * P_pred
x_est_kf[t] = x_kf
# --- 4. 結果の可視化 (2つのパネル) ---
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True)
# 上段: 線形KF (失敗例)
axes[0].plot(x_true, color='black', lw=2, linestyle='--', label='True State (Non-linear)')
axes[0].plot(x_est_kf, color='green', lw=1.5, label='Linear KF Estimate (Misfit)')
axes[0].scatter(range(T), y, color='blue', s=10, alpha=0.3, label='Observed (y)')
axes[0].set_title("Failure of Linear Model (Kalman Filter)")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 下段: 粒子フィルター (成功例)
axes[1].plot(x_true, color='black', lw=2, linestyle='--', label='True State (Non-linear)')
axes[1].plot(x_est_pf, color='red', lw=1.5, label='Particle Filter Estimate')
axes[1].scatter(range(T), y, color='blue', s=10, alpha=0.3, label='Observed (y)')
axes[1].set_title("Success of Particle Filter on Highly Non-linear System")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()線形カルマンフィルター (KF) は、事後分布を常に単峰性のガウス分布 N(x^t ,Pt ) と仮定する。観測モデルがのような非単調関数である場合、尤度関数は x の正負にピークを持つ多峰性を示す。しかしKFは分布を一つの平均値に集約しようとするため、多峰的な尤度を適切に処理できず、平均値が一方の局所解（プラス側）にトラップされてしまう。
一方で粒子フィルター (PF) は、事後分布を経験分布として表現する。尤度 p(yt ∣xt ) が正負の領域で等しく高い値を持つとき、粒子群は自然に二峰性の分布を形成する。この「解の分岐」を保持したまま次時刻の遷移モデル p(xt+1 ∣xt ) を適用することで、時間発展の過程でシステムモデルと矛盾する側のピークが自然に淘汰され、非線形性が強い状況下でも正しい状態へと収束が可能となる。
本記事のもくじはこちら：
