🔖 キーワード索引

📍 文脈ボックス — あなたが今見ているもの

あなたは 「推定法 / 状態空間モデル / 逐次推定」 の交点に位置する手法ページを見ています。 データ同化は 仮説検定 のような「真偽判定」ではなく、 真の状態に対する確率分布を逐次更新する 操作です。

同化は、 観測のみに頼る統計手法（平均・t検定 など）でも、 モデルのみに頼る古典力学的シミュレーションでもなく、 「両者をベイズの枠で結合する第 3 の道」です。

🎨 直感で掴む

あなたが GPS で現在地を知ろうとしている場面を想像してください。 ① 車の運動モデル「直前は 60 km/h で北上していたから、 1 秒後はおよそ 16.7 m 北」と 予測 できますが、 タイヤの滑りや風で誤差が貯まります。 ② 一方 GPS 観測値 は瞬時に位置を返しますが、 ±5 m 程度のノイズがあります。

データ同化は 「モデル予測の確信度」と「観測の確信度」を分散で測り、 確信度の重みで両者を平均」する処理です。 もしモデル予測の分散が小さければ予測寄りに、 観測の分散が小さければ観測寄りに最終推定が引き寄せられます。

統計学的に言えば 「事前分布（モデル）× 尤度（観測）= 事後分布（同化解）」。 これを 1 ステップごとに繰り返すのが 逐次データ同化、 期間全体の観測を一括して最適化するのが 変分データ同化（4D-Var） です。

比喩：占い師と気象観測員の合議

明日の気温を当てたい。 占い師（モデル予測）は「24℃ ± 3℃」と言い、 朝の気温計（観測）は前夜 22℃ ± 0.5℃を示しました。 観測員は「気温の自然な変動を考えて翌朝は 22.5℃ ± 0.5℃」と判定。 占い師より観測員の分散が小さいので最終予測は 22.5℃ 寄り。 もし観測機が壊れていて ±10℃ だと反対に占い師寄りになります。

📐 数式または定義

線形ガウスの状態空間モデル：

$$x_{t} = F_{t} x_{t-1} + w_{t}, \quad w_{t} \sim \mathcal{N}(0, Q_{t})$$

$$y_{t} = H_{t} x_{t} + v_{t}, \quad v_{t} \sim \mathcal{N}(0, R_{t})$$

カルマンフィルタの予測ステップ：

$$x^{f}_{t} = F_{t} x^{a}_{t-1}, \quad P^{f}_{t} = F_{t} P^{a}_{t-1} F_{t}^{\top} + Q_{t}$$

カルマンフィルタの更新ステップ：

$$K_{t} = P^{f}_{t} H_{t}^{\top} (H_{t} P^{f}_{t} H_{t}^{\top} + R_{t})^{-1}$$

$$x^{a}_{t} = x^{f}_{t} + K_{t}(y_{t} - H_{t} x^{f}_{t})$$

$$P^{a}_{t} = (I - K_{t} H_{t}) P^{f}_{t}$$

変分同化（3D-Var）のコスト関数：

$$J(x) = \tfrac{1}{2}(x - x^{f})^{\top} (P^{f})^{-1}(x - x^{f}) + \tfrac{1}{2}(y - Hx)^{\top} R^{-1}(y - Hx)$$

4D-Var はモデル制約を時間方向に含めた拡張：

$$J(x_0) = \tfrac{1}{2}(x_0 - x_0^f)^{\top}(P_0^f)^{-1}(x_0 - x_0^f) + \tfrac{1}{2}\sum_{t=0}^{T}(y_t - H_t M_t(x_0))^{\top} R_t^{-1}(y_t - H_t M_t(x_0))$$

粒子フィルタ（重要度サンプリング）の重み更新：

$$w_t^{(i)} \propto w_{t-1}^{(i)} \cdot p(y_t | x_t^{(i)})$$

🔬 数式を言葉で読み解く

第 2 式右辺の $y_t - H_t x^f_t$ は イノベーション（観測でわかった「予測のズレ」）。 これに $K_t$ を掛けて予測値に足す＝ 観測で予測を補正するのが核心です。 $R$ が大きい（観測が信用できない）と $K$ が小さくなり、 補正は弱まります。

共分散の更新 $P^a = (I - KH)P^f$ は、 「観測を取り込むほど不確実性が縮まる」ことを表現します。 $K=0$ なら何も変わらず、 $K=1$（観測完全信頼）なら $P^a = (1-H)P^f$ となり残った不確実性は $H$ が捉えられない次元のみ。

スカラー版で直感を再確認

1 次元では $K = \sigma_f^2 / (\sigma_f^2 + \sigma_r^2)$、 これは「分散の逆数を重みとする加重平均」と等価です。 $\sigma_f^2 \gg \sigma_r^2$ → $K \approx 1$ → 観測寄り、 $\sigma_f^2 \ll \sigma_r^2$ → $K \approx 0$ → 予測寄り。

🧮 実値で計算してみる（SSDSE-B-2026 都道府県人口）

北海道の 総人口 は SSDSE-B-2026 によれば、 2021 年 5,183,000 人 → 2022 年 5,140,000 人 → 2023 年 5,092,000 人。 単純な線形外挿では 2024 年予測 $x^f = 5{,}046{,}000$ 人。 仮にモデル分散 $P^f = (40{,}000)^2$、 観測誤差 $R = (10{,}000)^2$ とし、 観測値 $y = 5{,}050{,}000$ 人だったとします。

観測の精度がモデル予測より高い（$R \ll P^f$）ため、 解析値は観測寄りの 5,049,765 人。 解析分散は $9.4\times 10^7$ と 事前分散の 1/17 に縮小し、 「観測を取り込んだ分だけ確信度が増した」ことが定量化されます。

続けて 2 ステップ目を計算

2025 年の予測は $x^f_{25} = x^a_{24} = 5{,}049{,}765$、 $P^f_{25} = P^a_{24} + Q = 9.4\times 10^7 + 1.6\times 10^9 = 1.69\times 10^9$ となります。 もし 2025 年に観測 $y_{25}=5{,}010{,}000$ が得られれば $K_{25}=1.69/1.70=0.994$、 $x^a_{25}=5{,}010{,}238$ と観測ほぼそのまま採用される計算になります。

🐍 Python 実装

SSDSE-B-2026 の都道府県人口を題材に、 ローカルレベルモデル＋カルマンフィルタを実装します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 実 SSDSE-B-2026 を読み込み（cp932 で エンコード）
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', skiprows=1, encoding='cp932')
print(df.columns.tolist()[:6])

# 北海道（地域コード R01000）の総人口時系列を抽出
hokkaido = df[df['都道府県']=='北海道'].sort_values('年度')
y = hokkaido['総人口'].values.astype(float)
years = hokkaido['年度'].values
print(years, y)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

def kalman_local_level(y, Q=40000**2, R=10000**2, x0=None, P0=1e10):
    """ ローカルレベルモデル x_t = x_{t-1} + w, y = x + v """
    n = len(y)
    x_f = np.zeros(n); x_a = np.zeros(n)
    P_f = np.zeros(n); P_a = np.zeros(n)
    x_a[0] = y[0] if x0 is None else x0
    P_a[0] = P0
    for t in range(1, n):
        # 予測ステップ
        x_f[t] = x_a[t-1]
        P_f[t] = P_a[t-1] + Q
        # 更新ステップ
        K = P_f[t] / (P_f[t] + R)
        x_a[t] = x_f[t] + K * (y[t] - x_f[t])
        P_a[t] = (1 - K) * P_f[t]
    return x_a, P_a, x_f

x_a, P_a, x_f = kalman_local_level(y)
for t, yr in enumerate(years):
    print(f"{yr}: 観測={y[t]:.0f}  予測={x_f[t]:.0f}  解析={x_a[t]:.0f}")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 47 都道府県すべてに同化を適用し 2024 年予測の信頼区間を計算
pred_2024 = {}
for pref in df['都道府県'].unique():
    sub = df[df['都道府県']==pref].sort_values('年度')
    yi = sub['総人口'].values.astype(float)
    xa, Pa, xf = kalman_local_level(yi)
    next_pred = xa[-1]
    se = np.sqrt(Pa[-1] + 40000**2)
    pred_2024[pref] = (next_pred, next_pred-1.96*se, next_pred+1.96*se)

result = pd.DataFrame(pred_2024, index=['予測','下限95%','上限95%']).T
print(result.sort_values('予測', ascending=False).head(10))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# pykalman でマルチ変数系（人口＋出生数の同化）
from pykalman import KalmanFilter
sub = df[df['都道府県']=='東京都'].sort_values('年度')
Y = sub[['総人口','出生数']].values.astype(float)
kf = KalmanFilter(transition_matrices=np.eye(2),
                  observation_matrices=np.eye(2),
                  transition_covariance=np.diag([1e8, 1e4]),
                  observation_covariance=np.diag([1e6, 1e3]))
means, covs = kf.filter(Y)
print('東京都・人口と出生数の同化解析値:')
print(pd.DataFrame(means, columns=['人口','出生数'], index=sub['年度']))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# statsmodels の状態空間モデル (UnobservedComponents) で同化
import statsmodels.api as sm
osaka = df[df['都道府県']=='大阪府'].sort_values('年度')
yi = osaka['総人口'].values.astype(float)
model = sm.tsa.UnobservedComponents(yi, level='local linear trend')
res = model.fit(disp=False)
print(res.summary())
print('カルマンスムーザ解析値:', res.smoothed_state[0])
forecast = res.forecast(steps=3)
print('3 年先までの同化予測:', forecast)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# アンサンブルカルマンフィルタを 47 都道府県人口に適用（多変量同化）
def enkf_step(X_ens, y, H, R):
    """ X_ens: (n_state, n_ens), y: 観測, H: 観測演算子, R: 観測共分散 """
    n_ens = X_ens.shape[1]
    x_mean = X_ens.mean(axis=1, keepdims=True)
    X_anom = X_ens - x_mean
    P = (X_anom @ X_anom.T) / (n_ens - 1)
    HPH = H @ P @ H.T
    K = P @ H.T @ np.linalg.inv(HPH + R)
    # 観測を摂動して観測アンサンブル生成
    Y_pert = y[:, None] + np.linalg.cholesky(R) @ np.random.randn(len(y), n_ens)
    X_new = X_ens + K @ (Y_pert - H @ X_ens)
    return X_new

# 47都道府県人口を多変量状態としてアンサンブル同化（実SSDSEデータから初期化）
prefs = df['都道府県'].unique()
n_ens = 50
latest = df.sort_values('年度').groupby('都道府県').last()['総人口'].loc[prefs].values
X_ens = latest[:, None] + 50000*np.ones((len(prefs), n_ens))  # 初期アンサンブル
H = np.eye(len(prefs)); R = (10000.0**2) * np.eye(len(prefs))
y_obs = latest.astype(float)  # 観測を最新値に置く
X_ens_new = enkf_step(X_ens.astype(float), y_obs, H, R)
print('EnKF後の解析平均（上位5県）:', dict(zip(prefs[:5], X_ens_new.mean(axis=1)[:5])))

1
2
3
4
5
6
7
8
9

# プロット：観測・予測・解析を可視化（北海道）
plt.figure(figsize=(8,4))
plt.plot(years, y, 'ko-', label='観測 (SSDSE-B)')
plt.plot(years, x_f, 'r--', label='予測 (forecast)')
plt.plot(years, x_a, 'b-', label='解析 (analysis)')
plt.fill_between(years, x_a - 1.96*np.sqrt(P_a), x_a + 1.96*np.sqrt(P_a), alpha=0.2)
plt.title('北海道 総人口のデータ同化（ローカルレベル）')
plt.xlabel('年度'); plt.ylabel('人口')
plt.legend(); plt.tight_layout(); plt.savefig('hokkaido_assimilation.png', dpi=150)

📊 多変量同化の具体例

SSDSE-B-2026 で「総人口」「出生数」「死亡数」を同時に状態とした 2 次元同化を考えます。 状態 $x_t = (人口_t, 出生_t)^\top$、 観測 $y_t$ が両者を含むケース。 関連性として「出生数が多い県は若年人口維持で総人口の減少が緩い」という共分散構造を取り込めます。

同化解析値は単に観測平均ではなく、 人口・出生数・死亡数の 動的相関 をモデル経由で取り込んでいる点が重要です。 「人口減少県では死亡数が高い」「都市部では出生は維持されるが人口は微増のみ」など、 単変量同化では見えない構造が浮かぶ。

🛠 実務ワークフロー

1. 状態の定義：何を $x$ と置くか（人口、 気温、 流速…）。 次元数が決まる。
2. 力学モデル $F, Q$ の設計：自己回帰・物理方程式・成長率モデルなど。
3. 観測モデル $H, R$ の設計：センサーごとに $H$ を組み、 $R$ は校正データから推定。
4. 初期分布 $x_0, P_0$ の設定：分かっている情報があれば事前分散を絞る。
5. 逐次更新ループ：予測 → 更新を毎ステップ実行。
6. 診断：イノベーション系列の白色性、 $\chi^2$ 残差検定で $Q, R$ を補正。
7. スムージング：オフラインで過去全期間を再推定（RTS スムーザ）。
8. 不確実性報告：解析分散・信頼区間を必ず出力。

⚠️ 落とし穴

🌐 関連手法・派生

🔬 応用事例

• 気象予報：気象庁の数値予報モデルに 4D-Var を組み込み、 約 1 億変数を持つ大気場を毎 6 時間更新。
• 津波警報：地震直後の沖合津波計と数値モデルを同化し沿岸到達高を逐次推定。
• パンデミック予測：SEIR モデルに新規感染者数を同化、 $R_t$ を推定（COVID-19 では複数研究グループが粒子フィルタを使用）。
• 交通流推定：高速道路の検出器データと交通シミュレーションを EnKF で結合し渋滞長を予測。
• ナウキャスティング：GDP を四半期遅れ統計と月次指標の動的因子モデルで毎月更新。
• 地殻変動モニタリング：GNSS 位置時系列とプレート運動モデルの同化で歪み速度推定。
• 都道府県人口推計：本ページで扱う SSDSE-B のように、 国勢調査のあいだの中間年人口をモデル＋観測で補間。

📏 同化の評価指標

同化器が「自分の不確実性を正しく主張しているか」は スプレッド/RMSE 比で判断します。 1 未満なら過信（共分散が小さすぎ→インフレ）、 1 超なら過分散（共分散が大きすぎ→デフレ）です。

📜 歴史と発展

データ同化の起源は 1960 年代のアポロ計画とカルマンによる線形最小分散推定の提案にさかのぼります。 当初は航空宇宙の軌道決定に応用され、 1970 年代以降は気象学（Kalman, Bucy）、 海洋学（Evensen の EnKF）、 数値天気予報（ECMWF の 4D-Var 業務化 1997 年）と裾野を広げました。 2000 年代以降は 地球システム同化（大気・海洋・陸面・氷雪を同時に同化）へと進化しています。

• 1960 — Kalman, Bucy: 線形ガウス系の最適フィルタ理論
• 1986 — Lorenc: 気象データ同化を変分問題に定式化
• 1994 — Evensen: アンサンブルカルマンフィルタを海洋同化に適用
• 1997 — ECMWF: 4D-Var を業務天気予報で運用開始
• 2000s — 粒子フィルタの高次元への拡張（SIR, SMC²）
• 2010s — 機械学習との融合（同化＋ニューラルネットでの emulator 学習）
• 2020s — ハイブリッド EnVar、 アンサンブル変分、 同化＋深層学習（DA4ML, ML4DA）

⚙ 共分散インフレーションと局所化

EnKF では有限アンサンブルゆえに「共分散の縮小（covariance collapse）」が起こり、 同化が観測を受け付けなくなります。 対策として：

• 乗法インフレーション：$P^f \leftarrow (1+\alpha) P^f$（典型 $\alpha=0.05\sim0.20$）
• 加法インフレーション：$P^f \leftarrow P^f + \Sigma_\text{add}$
• 局所化：距離関数 $L(r)$（Gaspari-Cohn 関数）でハダマード積を取り遠方の偽相関を抑える

$$P^f_\text{loc} = L \circ P^f$$

🔗 関連用語（前提・並列・発展）

前提

並列

発展

📚 関連グループ教材

• 📚 標本抽出と中心極限定理 — 推定の前提となる分布論
• 📚 ロジスティック回帰と GLM — 観測モデルの設計
• 📚 仮説検定 — 同化残差の検定
• 📚 機械学習 — モデル予測の代替

❓ FAQ

Q1. データ同化は機械学習とどう違う？

機械学習は大量データから関数形を学習しますが、 同化は 明示的な物理/力学モデルを保持 し、 観測でその状態のみを更新します。 ハイブリッドアプローチでは ML が代理モデルを担い、 同化が状態を逐次更新する分業が研究されています。

Q2. データ同化と回帰の違いは？

回帰は静的な「説明変数→目的変数」の関係を推定しますが、 同化は 時間発展する状態 を観測で逐次補正します。 状態空間モデルでは回帰係数自体を状態として同化することも可能です（時変パラメータ）。

Q3. 線形カルマンが「最適」なのはどんな意味で？

線形ガウス系において、 KF はベイズ事後分布の 正確な平均と共分散を返します。 非線形・非ガウスでは「線形最小分散」という弱い意味でしか最適ではありません。

Q4. 観測が一切無いステップではどうなる？

予測ステップのみが実行され、 $P^f$ が単調増加します。 観測が再開した瞬間に大きな更新（カルマンゲイン高）が起き、 解析誤差は再び縮小します。

Q5. 50 状態×100 観測のシステムで何を計算する？

$P^f$ は 50×50、 $HPH^\top + R$ は 100×100。 後者の逆行列が支配的な計算量。 EnKF を使えば $P^f$ を陽に持たず、 アンサンブルとの行列積で計算量を削減できます。

Q6. SSDSE-B のような年次データでも同化は有効？

はい。 年次の観測でも「年次から年次へのモデル」と「年次観測」を同化すれば、 国勢調査年と年次推計の差を自然に統合できます。 観測欠損期もモデルだけで予測し、 観測が出た年に補正が入ります。

🗺 概念マップ

データ同化は 状態空間モデル の解析手法、 ベイズ推定 の逐次版、 最適化問題 としても定式化される三つの顔を持ちます。 用語の地図全体 も参照してください。