はじめに
Wahba の定理に基づく曲線信頼区間の推定原理と、Integrated Wiener process と cubic smoothing spline の深い関係を物理屋向けに解説をしてみます。
Improper Priors, Spline Smoothing and the Problem of Guarding Against
Model Errors in Regression,
Grace Wahba,
Journal of the Royal Statistical Society. Series B (Methodological), Vol. 40, No. 3
(1978), 364-372.
「1回積分された Wiener 過程（Integrated Wiener process）を GP の prior に選ぶと、その posterior の期待値は必ず cubic smoothing spline になる」
— Wahba, Kimeldorf–Wahba, Rasmussen & Williams ほか
この事実は、一見すると「え、本当？」となるくらい意外です。
- 一方では 純粋に変分法の世界で定義される smoothing spline：
\sum_i (y_i - f(x_i))^2 + \lambda \int (f''(x))^2 dx - もう一方では 確率過程（Gaussian Process）の世界で定義される prior・posterior：
f \sim \mathcal{GP}(\text{IWPカーネル}),\quad y_i = f(x_i) + \epsilon_i
この二つが、数学的には同じものを定義している、というのが Wahba の定理の本質です。
この記事では、解析力学や場の理論の素養を持つ「物理屋」を想定し、
- smoothing spline を 作用最小問題として理解する
- Integrated Wiener process を「曲率エネルギーを持つ場」として理解する
- Wahba の定理：
GP の posterior 平均 = cubic smoothing spline をつなぐ - さらに GP が与える 曲線全体の信頼区間（不確かさ） の意味
- 実際の応用（例：エネルギーキャリブレーション、時系列解析など）
までを整理してみます。
(注) 専門外です。誤りやコメント大歓迎です。下記の記事の補足説明を試みました。
あたりをしっかりと読み込むべきなんでしょう。
1. 問題設定：ばらつきのある測定値から「なめらかな曲線と信頼区間」を得たい
物理実験でも天文観測でも、典型的な問題はこうです：
- 観測値：
(x_i, y_i),\quad i=1,\dots,N - 誤差（標準偏差）：
\sigma_i \quad(\text{既知とする}) - モデル：ある「なめらかな」関数 $f(x)$ があって、
y_i = f(x_i) + \epsilon_i,\quad \epsilon_i \sim \mathcal{N}(0,\sigma_i^2)
やりたいことは 2 つ：
- 「いい感じのなめらかな推定曲線」 $\hat{f}(x)$ を得る
- その曲線に対する 信頼区間（不確かさ） を、x の全域で見積もる
古典的な cubic smoothing spline は (1) だけを与える枠組みですが、Gaussian Process Regression (GPR) を使うと、(2) まで自然に手に入る、というのがポイントです。
2. smoothing spline の変分原理：曲率エネルギーを最小化する
まずは確率の話を忘れて、純粋に変分問題として smoothing spline を考えます。
2.1 smoothing spline の定義
cubic smoothing spline は、次の汎関数を最小にする関数 $f$ として定義されます：
J[f]
= \underbrace{\sum_{i=1}^N \frac{(y_i - f(x_i))^2}{\sigma_i^2}}_{\text{データとのズレ（χ²）}}
+ \underbrace{\lambda \int_a^b (f''(x))^2 dx}_{\text{曲率エネルギー（なめらかさペナルティ）}}.
- 第1項：観測点での誤差二乗和（重み付き）
- 第2項：曲率（2階微分）の 2 乗の積分 → 曲がりすぎを罰するエネルギー
- (\lambda>0)：両者のトレードオフを制御するスムージングパラメータ
$\lambda\to 0$ なら「ほぼ補間（interpolation）」、
$\lambda\to\infty$ なら「ほぼ直線」に近づきます。
2.2 まずは「曲率エネルギー」だけ最小化すると？
本質を見るために、まずデータ項を捨てて
J_0[f] = \int_a^b (f''(x))^2 dx
だけを最小化する状況を考えます。
これは、解析力学で言えば
- 「$S[q]=\int L(q,\dot{q})dt$ を最小にする軌道 q(t)」の「L が f'' まで依存する版」
だと思って良いでしょう。
3. 変分法で導く：なぜ「4階微分方程式」→ cubic になるのか
3.1 変分を入れる：f --> f + εη
変分として
f(x) \to f(x) + \epsilon \eta(x),
を考えます。$\eta(x)$ は任意の（十分滑らかな）テスト関数とします。
2階微分は
f''(x) \to f''(x) + \epsilon \eta''(x).
これを $J_0$ に代入すると：
J_0[f+\epsilon\eta]
= \int_a^b (f''+\epsilon\eta'')^2 dx
= \int_a^b \bigl(f''^2 + 2\epsilon f''\eta'' + O(\epsilon^2)\bigr) dx.
したがって、一次の変分は
\delta J_0
= \left.\frac{d}{d\epsilon} J_0[f+\epsilon\eta]\right|_{\epsilon=0}
= 2 \int_a^b f''(x)\,\eta''(x)\,dx.
3.2 部分積分を2回やる
$\eta''$ が付いていると扱いづらいので、部分積分を 2 回行います：
1回目：
\int_a^b f''\eta'' dx
= \bigl[f''\eta'\bigr]_a^b - \int_a^b f'''\eta' dx.
2回目：
\int_a^b f'''\eta' dx
= \bigl[f'''\eta\bigr]_a^b - \int_a^b f''''\eta dx.
合わせると：
\int_a^b f''\eta'' dx
= \bigl[f''\eta' - f'''\eta\bigr]_a^b + \int_a^b f''''(x)\,\eta(x)\,dx.
自然境界条件（natural boundary）（例えば $f''(a)=f''(b)=0$ など）を課すと、境界項はゼロとみなせます。
このとき、
\delta J_0 = 2 \int_a^b f''''(x)\, \eta(x)\,dx.
3.3 任意の η に対して δJ=0 ⇒ f''''(x)=0
変分法の基本定理：
$\eta(x)$ が任意のテスト関数なので、$\delta J_0=0$ が常に成り立つためには、被積分関数がゼロでなければならない
したがって、
f''''(x) = 0.
が必要条件です。
これは 4階常微分方程式 で、その一般解は
f(x) = ax^3 + bx^2 + cx + d
すなわち各区間で「三次多項式」となります。
さらに knot（接続点）で
- $f, f', f''$ が連続
になるようにすると、cubic spline になります。
✅ 結論：
曲率エネルギーを最小化すると、4階の Euler–Lagrange 方程式 $f''''=0$ が出てきて，
その解は「区間ごとの三次多項式（cubic）」になる。
4. smoothing spline：データ項を足しても本質は同じ
現実には、データとのズレも重要なので、
J[f]
= \sum_{i=1}^N \frac{(y_i - f(x_i))^2}{\sigma_i^2}
+ \lambda \int_a^b (f''(x))^2 dx
を最小化します。
ここでは
- 第1項：観測点での誤差 → 離散的な条件（デルタ関数的）
- 第2項：曲率エネルギー → 連続的ななめらかさの制約
として働きます。
厳密な導出では、Euler–Lagrange 方程式側に「デルタ関数が載った 4階微分方程式」が出てきますが、結論としては：
- 各区間の内部では $f''''(x)=0$ → 三次多項式
- 観測点で 2階導関数に jump 条件が入る
→ それを満たすように係数が決まる
という結果になり、最終的に natural cubic smoothing spline が得られます。
5. Gaussian Process の視点：Integrated Wiener process を prior にする
ここから 確率モデル の世界に入り、Wahba の定理に繋げます。
5.1 Gaussian Process の基本
Gaussian Process（GP）は、ざっくり言うと：
任意の有限集合 ((x_1,\dots,x_N)) に対し
((f(x_1),\dots,f(x_N))) が 多変量ガウス分布に従うような「関数の分布」
のことです。
- 平均関数：(m(x))
- 共分散関数（カーネル）：(k(x,x'))
を指定すれば、GP が定まります。
5.2 Integrated Wiener process（IWP）のイメージ
標準 Wiener 過程（Brown 運動）$W(t)$ に対して、
-
一回積分された Wiener 過程：
X(t) = \int_0^t W(s),ds -
もしくは「$f''(x)$ が白色雑音」のような形式的表現：
f''(x) = \xi(x), \quad \xi(x):\ \text{White noise}
を持つような GP を考えます（定義の流儀は本によってやや違うようですが、ここでは「曲率がランダムな過程」というイメージとします）。
このとき、path の事前分布は形式的に
p[f] \propto \exp\left(
-\frac{1}{2\sigma_f^2}\int_a^b (f''(x))^2 dx
\right)
という 「曲率エネルギーを持つガウス場」 の形になります。
✅ 重要ポイント： Integrated Wiener process の prior の log 密度は、smoothing spline の曲率エネルギーと同じ functional になっている。
6. 観測を入れる：posterior の負 log は smoothing spline の functional に一致する
観測モデルを
y_i = f(x_i) + \epsilon_i,\quad
\epsilon_i \sim \mathcal{N}(0, \sigma_i^2)
とします（ガウス雑音）。
Bayes の定理：
p(f\mid y) \propto p(y\mid f)\,p(f).
両辺の対数を取ってマイナスを付けると、
物理屋にはお馴染みの「作用」の形になります：
- \log p(f\mid y) = -\log p(y\mid f) - \log p(f) + \text{const}.
-
観測 likelihood は独立なガウスなので：
-\log p(y\mid f) = \frac{1}{2}\sum_i\frac{(y_i-f(x_i))^2}{\sigma_i^2} + \text{const}. -
prior は IWP：
-\log p(f) = \frac{1}{2\sigma_f^2}\int_a^b(f''(x))^2 dx + \text{const}.
両者を足すと：
- \log p(f\mid y)
= \frac{1}{2}\sum_i\frac{(y_i-f(x_i))^2}{\sigma_i^2}
+ \frac{1}{2\sigma_f^2}\int_a^b(f''(x))^2 dx
+ \text{const}.
$\lambda = 1/\sigma_f^2$ とおけば、
- \log p(f\mid y) = \frac{1}{2} J[f] + \text{const}.
ここで $J[f]$ は先ほどの smoothing spline の functional そのものです。
✅ つまり：posterior の最大化（MAP 推定）は，smoothing spline の変分問題の解と完全に一致する。
7. Wahba の定理：posterior の期待値も smoothing spline と一致する
線形ガウスモデルの世界（GP＋ガウス雑音）では、
- posterior はガウス過程のまま
- ガウス分布では 平均とモード（最尤点）は一致 する
ので、
\hat{f}_{\text{MAP}}(x)
= \arg\max_f p(f\mid y)
= \mathbb{E}[f(x)\mid y]
= \mu_{\text{post}}(x).
この $\mu_{\text{post}}(x)$ が cubic smoothing spline になっている、というのが Wahba の定理の中核です。
✅ 結論：Integrated Wiener process を prior にした Gaussian Process の posterior 平均は，cubic smoothing spline と一致する。
8. GP がもたらす「曲線全体の信頼区間」
ここまでで、
- 「最もありそうな曲線」 = smoothing spline = GP posterior mean
まで分かりました。
GP の真価は、ここから先です：
曲線に対する「不確かさ（variance）」まで自然に得られる。
8.1 posterior 分布は GP のまま
GP の枠組みでは、
f \mid y \sim \mathcal{GP}(\mu_{\text{post}}(x), k_{\text{post}}(x,x')).
となります。
したがって、
- 任意の点 $x_*$ での値は
f(x_*) \mid y \sim \mathcal{N}(\mu_{\text{post}}(x_*),, k_{\text{post}}(x_*,x_*)) - 任意の有限個の点 $(x_1,\dots,x_M)$ でのベクトルは、多変量ガウス
$
(f(x_1),\dots,f(x_M))^\top \sim \mathcal{N}\bigl(\boldsymbol{\mu},, \mathbf{K}\bigr).
$
ここから
-
点ごとの 95%信頼区間：
\mu_{\text{post}}(x) \pm 1.96~\sqrt{k_{\text{post}}(x,x)} - ある区間での平均値・積分値の分布（線形汎関数を通したガウス）
なども、そのまま計算できます。
8.2 smoothing spline との関係
smoothing spline 単体は
- 最適な一本の関数 $f_{\text{ss}}(x)$ を返すだけで、
- 「その周りにどれだけ他の関数があり得るか」は教えてくれない
一方、IWP-GP による解釈では、
- prior：曲率に対してガウスな重みがかかった「場のパス積分」
- posterior：観測で更新された「場の揺らぎ」
となるので、
smoothing spline は「曲率エネルギーを持つ場の classical solution（期待値）」であり，その周りに「量子的な揺らぎ」に相当する分散構造がある
という、物理屋的に非常に美しい描像が得られます。
9. 物理的なアナロジー：最小作用の classical 解と量子的揺らぎ
ここまでの話を、場の理論的な比喩でまとめると：
- 作用：
（前半：場の自由項、後半：観測との相互作用）
S[f] = \frac{1}{2\sigma_f^2}\int (f''(x))^2 dx + \frac{1}{2}\sum_i \frac{(y_i - f(x_i))^2}{\sigma_i^2} - classical 解：
\frac{\delta S}{\delta f} = 0 \quad\Rightarrow\quad \text{smoothing spline} - path integral：
Z[y] = \int \mathcal{D}f\ \exp\bigl(-S[f]\bigr) - 期待値：
\langle f(x)\rangle = \frac{1}{Z[y]} \int \mathcal{D}f\ f(x)~\exp\bigl(-S[f]\bigr) = \text{posterior mean} = \text{smoothing spline} - 二点関数：
\langle f(x)f(x')\rangle - \langle f(x)\rangle\langle f(x')\rangle = k_{\text{post}}(x,x')
という構造で捉えることもできる。
IWP-GP は「4階の微分作用素を持つガウス場」のパス積分問題であり、その classical solution が cubic smoothing spline と言える。
10. 応用例（実務編）：どこで役に立つのか？
10.1 キャリブレーション曲線（エネルギースケールなど）
X線マイクロカロリメータや TES などの 非線形なエネルギー応答 をキャリブレーションするとき：
-
複数の既知エネルギー線から、チャネル(PHA, filt_valueなど)→エネルギーの関係を推定したい
-
高次多項式フィットだと外挿で暴れやすい
-
スプラインなら「局所的に三次」で柔らかく補間できる
-
さらに GP の枠組みを使えば
- エネルギー軸の不確かさ（信頼区間）も同時に推定できる
- 系統誤差 vs 統計誤差の議論がしやすい
このとき、まさに
- smoothing spline = エネルギーキャリブレーション曲線の最良推定
- GP posterior covariance = キャリブレーション曲線の不確かさ
として使えます。
10.2 スペクトル・光度曲線のスムージング
-
雑音だらけの light curve や spectrum を「なめらか」に表示したいとき
-
強引に polynomial を上げていくと高周波の振動（Runge 現象）が出る
-
スプライン + GP を用いれば
- なめらかな曲線表示
- 強調したい部分の信頼区間表示（データが少ない領域では広い）
など、物理的解釈を保ちながら視覚的にわかりやすい図が得られます。
10.3 衛星軌道・姿勢の平滑化（Kalman フィルタとのつながり）
-
衛星の位置・速度・姿勢は、通常 カルマンフィルタ＋スムーザー で推定される
-
線形ガウス系においては、カルマンスムーザーと GP 回帰は数学的に等価
-
IWP-GP は「加速度が白色雑音のモデル」と考えることができ、
- smoothing spline = その軌道のなめらかな推定
- GP covariance = 軌道の不確かさ
という位置付けになります。
11. 実装イメージ
11.1 IWP カーネルの GP（概念）
IWP に対応する「spline カーネル」として、Rasmussen & Williams の式 (6.28) に相当する
k(x, x') = \sigma_f^2 \left[
\frac{1}{3} \min(x,x')^3
+ \frac{1}{2} \min(x,x')^2 |x-x'|
\right]
のようなものを使うと、
この GP の posterior 平均が上記 smoothing spline と一致します。
実際には、
- 共分散行列 $K$ をこのカーネルから構成
- 観測ノイズを加えた $K_y = K + \mathrm{diag}(\sigma_i^2)$ を作る
- GP の公式：
\mu_{\text{post}}(x_*) = k_*^T K_y^{-1} yを使う\mathrm{Var}[f(x_*)] = k(x_*,x_*) - k_*^T K_y^{-1} k_*
という手順になります。
※ NIST の interpolate.py 内の GPRSpline クラスは、まさにこの「IWP-GP に基づく smoothing spline」を実装している例です。
12. まとめ：Wahba の定理が教えてくれること
-
変分法の世界：
- 曲率エネルギー $\int (f'')^2 dx$ の最小化から 4階微分方程式 $f''''=0$ が導かれる
- データ項を加えた最適化問題の解が cubic smoothing spline
-
確率過程（GP）の世界：
- Integrated Wiener process (IWP) を prior に選ぶと、その log 密度は $\int (f'')^2 dx$ に比例
- ガウス雑音観測と組み合わせた posterior の負 log が smoothing spline の functional と一致
- 線形ガウス系では posterior の平均 = モード = smoothing spline
-
Wahba の定理：
- IWP を prior にした GP の posterior 平均は cubic smoothing spline そのもの
- posterior 共分散は「曲線全体の信頼区間」を与える
-
物理的解釈：
- smoothing spline は「曲率エネルギーを持つ1次元場の classical solution」
- GP の揺らぎは、その場の「量子的（統計的）な揺らぎ」を表す
- 解析力学の最小作用原理・場のパス積分と同じ構造を持つ
結論：
コンピュータ科学や統計の文脈で語られる「Gaussian Process Regression」や「smoothing spline」は、物理屋から見ると 「高階微分作用素を持つガウス場のパス積分問題」 そのものであり、Wahba の定理はその classical 解（spline）と揺らぎ（信頼区間）がどのように対応しているかを教えてくれる定理である。