はじめに
X線検出器やガンマ線スペクトロメーターのエネルギー較正では、ガウシアンフィッティングで得られた既知エネルギーピークの位置（チャンネル番号）を用いて、任意のチャンネルをエネルギーに変換する較正曲線を構築します。この際、最小二乗法による多項式フィットではなく、スプライン補間を用いる理由と、真の較正関数が不明な状況での適合度評価方法、そして統計誤差の適切な取り扱いが重要な課題となります。
本記事では以下のトピックを詳細に解説します：
- なぜ最小二乗法の多項式フィットではなくスプラインを使うのか
- 真の較正関数が不明な場合の適合度評価手法
- スプライン補間における統計誤差の考慮方法
- 統計学的に正しい不確かさ推定の理論と実装
1. スプライン補間を用いる理由：最小二乗法との比較
1.1 なぜ最小二乗法の多項式フィットでは不十分なのか
エネルギー較正では、最小二乗法による大域的多項式フィッティングではなく、スプライン補間が広く用いられています。その主な理由は以下の通りです。
(1) 非線形性への柔軟な対応
検出器の応答関数は必ずしも単純な多項式では表現できません。スプラインは区分的多項式のため、複雑な非線形関係も柔軟に表現できます。
(2) Runge現象による振動
高次の多項式フィットはRunge現象と呼ばれる問題を引き起こします。これは、較正点数(n)が増えると、(n-1)次の多項式が較正点間で激しく振動する現象です。
問題点：
- 較正点では正確にフィットするが、較正点間で大きく振動
- 特に端点付近で予測が不安定になる
- 物理的に不合理な値（負のエネルギーなど）を予測する可能性
(3) 局所性の欠如
大域的多項式フィットでは、1つの較正点のデータが全体の曲線形状に影響します。一方、スプラインは局所性を持ちます：
\text{スプライン: } S(x) = \begin{cases}
p_1(x) & x \in [x_0, x_1] \\
p_2(x) & x \in [x_1, x_2] \\
\vdots & \vdots \\
p_n(x) & x \in [x_{n-1}, x_n]
\end{cases}
各区間$[x_{i-1}, x_i]$は低次多項式$p_i(x)$（通常は3次）で表され、隣接区間のみに影響します。
メリット：
- 異常値の影響が局所的に限定される
- 較正点の追加・削除が容易
- 数値的に安定
1.2 最小二乗法が適している場合
ただし、以下の場合は最小二乗法の多項式フィットが有効です：
- 較正点が非常に少ない（3〜4点）：スプラインの自由度が不足
- 測定誤差が非常に大きい：補間ではなく平滑化が必要
- 外挿が主目的：スプラインの外挿は不安定
- 線形性が強い：線形からの誤差をそのまま見積もれる
1.3 スプライン vs 最小二乗法の比較表
| 特性 | 最小二乗多項式 | スプライン補間 |
|---|---|---|
| 較正点でのフィット | 近似（残差あり） | 厳密通過（残差≈0） |
| 較正点間の挙動 | 高次で振動しやすい | 滑らかで安定 |
| 外挿領域 | 次数による | 不安定（推奨されない） |
| 局所性 | なし | あり |
| 計算コスト | 低い | やや高い |
| 統計誤差の考慮 | 自然に組み込み可能 | 別途工夫が必要 |
2. 真の較正関数が不明な場合の適合度評価
エネルギー較正の最大の課題は、真の較正関数$E_{\text{true}}(ch)$が不明であることです。このため、通常のRMSE：
\text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(E_{\text{true}}(ch_i) - E_{\text{pred}}(ch_i))^2}
を直接計算することができません。この問題に対処するため、複数の間接的評価手法を組み合わせます。
2.1 クロスバリデーション（推奨手法）
クロスバリデーションは、真の値が不明でも較正曲線の予測性能を評価できる最も有効な手法です。
Leave-One-Out Cross Validation (LOOCV)
原理：
- (n)個の較正点から1点を除外
- 残り(n-1)点で較正曲線を構築
- 除外した点での予測値と実測値を比較
- 全ての点について繰り返し
評価式：
\text{RMSE}_{\text{LOOCV}} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(E_i - \hat{E}_{-i}(ch_i))^2}
ここで、$\hat{E}_{-i}(ch_i)$は$i$番目の点を除いて構築した較正曲線による予測値です。
from sklearn.model_selection import LeaveOneOut
from scipy.interpolate import CubicSpline
def loocv_spline(channels, energies):
loo = LeaveOneOut()
predictions = []
for train_idx, test_idx in loo.split(channels):
# 訓練データでスプライン構築
train_ch = channels[train_idx]
train_en = energies[train_idx]
spline = CubicSpline(train_ch, train_en)
# テストデータで予測
test_ch = channels[test_idx]
predictions.append(spline(test_ch)[0])
# RMSE計算
rmse = np.sqrt(np.mean((np.array(predictions) - energies)**2))
return rmse
# 実行例
channels = np.array([100, 300, 500, 700, 900, 1100, 1300])
energies = np.array([60, 240, 500, 840, 1260, 1760, 2340])
rmse_cv = loocv_spline(channels, energies)
print(f"LOOCV RMSE: {rmse_cv:.4f} keV")
解釈：
- LOOCVのRMSEが小さい → 較正点間での予測精度が高い
- LOOCVのRMSEが大きい → 過学習の兆候、外挿性能が低い
K-Fold Cross Validation
較正点数が多い場合$(n \geq 15)$はK-Fold CVが効率的です：
from sklearn.model_selection import KFold
def kfold_spline(channels, energies, k=5):
"""K-Fold CVによるスプライン評価"""
kf = KFold(n_splits=k, shuffle=True, random_state=42)
mse_scores = []
for train_idx, test_idx in kf.split(channels):
train_ch = channels[train_idx]
train_en = energies[train_idx]
test_ch = channels[test_idx]
test_en = energies[test_idx]
spline = CubicSpline(train_ch, train_en)
predictions = spline(test_ch)
mse = np.mean((predictions - test_en)**2)
mse_scores.append(mse)
rmse = np.sqrt(np.mean(mse_scores))
return rmse
2.2 残差解析
較正点での残差パターンから、モデルの適切さを評価します。
残差の定義
r_i = E_i - \hat{E}(ch_i)
ここで、$E_i$は既知エネルギー、$\hat{E}(ch_i)$は較正曲線による予測値です。
評価指標
-
残差の平均：
\bar{r} = \frac{1}{n}\sum_{i=1}^{n}r_i \approx 0ゼロから大きく外れる場合、系統的バイアスが存在
-
残差の標準偏差：
\sigma_r = \sqrt{\frac{1}{n-p}\sum_{i=1}^{n}(r_i - \bar{r})^2}ここで$p$はモデルのパラメータ数
-
残差の無作為性検定（Runs Test）：
残差の符号の並びが無作為かを検定- 有意に非無作為 → モデルが不適切
from scipy import stats
def residual_analysis(channels, energies, model_func):
"""残差解析"""
residuals = energies - model_func(channels)
mean_res = np.mean(residuals)
std_res = np.std(residuals, ddof=1)
max_res = np.abs(residuals).max()
signs = np.sign(residuals)
n_runs = 1 + np.sum(signs[:-1] != signs[1:])
n_pos = np.sum(signs > 0)
n_neg = np.sum(signs < 0)
expected_runs = 1 + 2*n_pos*n_neg/(n_pos + n_neg)
print(f"残差の平均: {mean_res:.6f} keV")
print(f"残差の標準偏差: {std_res:.4f} keV")
print(f"最大絶対残差: {max_res:.4f} keV")
print(f"Runs数: {n_runs} (期待値: {expected_runs:.1f})")
return residuals
重要な注意：
補間スプラインは較正点で厳密に$r_i = 0$となるため、残差解析は意味を持ちません。残差解析は最小二乗法や平滑化スプラインにのみ有効です。
2.3 情報量規準（AIC、BIC）
モデルの複雑さとフィット具合のバランスを評価します。
Akaike Information Criterion (AIC)
\text{AIC} = n \ln(\text{RSS}/n) + 2p
ここで：
- $\text{RSS} = \sum_{i=1}^{n}(E_i - \hat{E}(ch_i))^2$：残差平方和
- $p$：モデルのパラメータ数
- $n$：較正点数
解釈：AICが小さいほど良いモデル
Bayesian Information Criterion (BIC)
\text{BIC} = n \ln(\text{RSS}/n) + p \ln(n)
BICはAICよりもモデルの複雑さに対してペナルティが大きい。
def calculate_aic_bic(residuals, n_params, n_data):
"""AICとBICの計算"""
rss = np.sum(residuals**2)
n = n_data
p = n_params
aic = n * np.log(rss / n) + 2 * p
bic = n * np.log(rss / n) + p * np.log(n)
return aic, bic
# 比較例
# 2次多項式: p=3（係数3個）
# 3次スプライン（7点）: p=7（自由度7）
residuals_poly2 = energies - poly2_func(channels)
residuals_spline = energies - spline_func(channels)
aic_poly2, bic_poly2 = calculate_aic_bic(residuals_poly2, 3, 7)
aic_spline, bic_spline = calculate_aic_bic(residuals_spline, 7, 7)
print(f"2次多項式: AIC={aic_poly2:.2f}, BIC={bic_poly2:.2f}")
print(f"3次スプライン: AIC={aic_spline:.2f}, BIC={bic_spline:.2f}")
注意：
補間スプラインは$\text{RSS} \approx 0$となるため、AICは$-\infty$に発散します。したがって、AIC/BICは平滑化スプラインと多項式の比較にのみ有効です。
2.4 独立な検証データによる評価
最も信頼性の高い評価方法は、較正に使用していない独立な検証用データセットを用いることです。
手順：
- 既知エネルギーピークを2つのグループに分割
- 訓練データ：較正曲線の構築に使用（例：5点）
- 検証データ：評価のみに使用（例：2点）
- 訓練データで較正曲線を構築
- 検証データでの予測誤差を計算
# 例：7点のうち5点を訓練、2点を検証
train_idx = [0, 1, 2, 4, 6] # チャンネル100, 300, 500, 900, 1300
test_idx = [3, 5] # チャンネル700, 1100
train_ch = channels[train_idx]
train_en = energies[train_idx]
test_ch = channels[test_idx]
test_en = energies[test_idx]
# 較正曲線構築
spline = CubicSpline(train_ch, train_en)
# 検証
predictions = spline(test_ch)
validation_rmse = np.sqrt(np.mean((predictions - test_en)**2))
print(f"検証RMSE: {validation_rmse:.4f} keV")
限界：
X線天文学では既知エネルギーピークの数が限られるため、データ分割は現実的でない場合が多い。
2.5 物理的制約による評価
較正曲線が物理的に妥当かを確認します：
-
単調性：エネルギーはチャンネルに対して単調増加
\frac{dE}{dch} > 0 \quad \forall ch \in [ch_{\min}, ch_{\max}] -
非負性：エネルギーは常に正
E(ch) \geq 0 \quad \forall ch -
滑らかさ：不自然な屈曲がないこと
\left|\frac{d^2E}{dch^2}\right| < \text{threshold} -
外挿の妥当性：較正範囲外での予測が物理的に合理的
def check_physical_constraints(spline, ch_min, ch_max, n_points=1000):
"""物理的制約の確認"""
test_ch = np.linspace(ch_min, ch_max, n_points)
energies = spline(test_ch)
derivatives = spline(test_ch, 1) # 1階微分
monotonic = np.all(derivatives > 0)
non_negative = np.all(energies >= 0)
second_deriv = spline(test_ch, 2)
max_curvature = np.abs(second_deriv).max()
print(f"単調性: {'OK' if monotonic else 'NG'}")
print(f"非負性: {'OK' if non_negative else 'NG'}")
print(f"最大曲率: {max_curvature:.6f}")
return monotonic and non_negative
2.6 評価手法の組み合わせ（推奨）
実務では、複数の評価手法を組み合わせることが推奨されます：
def comprehensive_evaluation(channels, energies, model_func, model_name):
"""総合的な較正評価"""
print(f"\n{'='*60}")
print(f"{model_name}の評価")
print(f"{'='*60}")
# 1. LOOCV
rmse_cv = loocv_spline(channels, energies)
print(f"\n[1] LOOCV RMSE: {rmse_cv:.4f} keV")
# 2. 残差解析（最小二乗法の場合のみ）
if model_name != "補間スプライン":
residuals = residual_analysis(channels, energies, model_func)
# 3. AIC/BIC
n_params = estimate_params(model_func) # モデル固有
aic, bic = calculate_aic_bic(residuals, n_params, len(channels))
print(f"\n[3] AIC: {aic:.2f}, BIC: {bic:.2f}")
# 4. 物理的制約
print(f"\n[4] 物理的制約:")
check_physical_constraints(model_func, channels.min(), channels.max())
print(f"{'='*60}\n")
3. スプライン補間における統計誤差の問題
3.1 問題の本質
ガウシアンフィッティングで得られるピーク位置には統計誤差が伴います：
\text{ch}_i = \mu_i \pm \sigma_{\mu,i}
しかし、通常のスプライン補間では点推定値(\mu_i)のみを使用し、誤差(\sigma_{\mu,i})が無視されます。これにより：
- 較正曲線の不確かさが過小評価される
- 誤差が大きい較正点の影響が適切に反映されない
- 較正結果に信頼区間を付けられない
3.2 なぜスプラインでは統計誤差を入れにくいのか
補間スプラインは較正点を厳密に通る曲線を構築します：
S(ch_i) = E_i \quad \forall i
この制約により：
- 測定誤差$\sigma_{\mu,i}$を考慮する自由度がない
- 誤差が大きい点も小さい点も同等に扱われる
- 不確かさの伝播が不可能
3.3 解決策1：平滑化スプライン + 重み付け
平滑化スプライン（smoothing spline）では、較正点を厳密に通る制約を緩和し、滑らかさを優先します。
定式化
平滑化スプラインは以下の最適化問題を解きます：
\min_{S} \left[ \sum_{i=1}^{n} w_i (E_i - S(ch_i))^2 + \lambda \int (S''(ch))^2 dch \right]
ここで：
- 第1項：フィッティング誤差（重み付き）
- 第2項：曲線の滑らかさペナルティ
- $w_i$：各較正点の重み
- $\lambda$：平滑化パラメータ
重みの設定
測定誤差の逆数の二乗を重みとします：
w_i = \frac{1}{\sigma_{\mu,i}^2}
誤差が大きい点ほど重みが小さく、フィッティング誤差の許容範囲が広がります。
from scipy.interpolate import UnivariateSpline
def weighted_smoothing_spline(channels, energies, channel_errors, smoothing=None):
weights = 1.0 / (channel_errors ** 2)
if smoothing is None:
smoothing = len(channels) * np.mean(channel_errors**2)
# 平滑化スプライン構築
spline = UnivariateSpline(
channels,
energies,
w=weights,
s=smoothing
)
return spline
channels = np.array([100, 300, 500, 700, 900, 1100, 1300])
energies = np.array([60, 240, 500, 840, 1260, 1760, 2340])
channel_errors = np.array([0.5, 0.4, 0.3, 0.35, 0.4, 0.5, 0.6])
spline_weighted = weighted_smoothing_spline(channels, energies, channel_errors)
test_ch = 450
predicted_energy = spline_weighted(test_ch)
print(f"Ch={test_ch}: E={predicted_energy:.2f} keV")
平滑化パラメータの選択
$\lambda$（scipyではs）の選択は重要です：
-
手動設定：
s = n \cdot \bar{\sigma}^2ここで$\bar{\sigma}^2$は誤差分散の平均
-
一般化クロスバリデーション（GCV）：
\text{GCV}(\lambda) = \frac{\frac{1}{n}\|(I - A(\lambda))E\|^2}{\left(\frac{1}{n}\text{tr}(I - A(\lambda))\right)^2}GCVを最小化する$\lambda$を選択
def optimize_smoothing_gcv(channels, energies, weights, s_candidates):
gcv_scores = []
for s in s_candidates:
spline = UnivariateSpline(channels, energies, w=weights, s=s)
residuals = energies - spline(channels)
eff_df = len(channels) - s / (s + 1)
gcv = np.mean(residuals**2) / ((1 - eff_df/len(channels))**2)
gcv_scores.append(gcv)
optimal_idx = np.argmin(gcv_scores)
optimal_s = s_candidates[optimal_idx]
return optimal_s, gcv_scores
# 使用例
s_range = np.logspace(-2, 2, 50)
weights = 1.0 / (channel_errors ** 2)
optimal_s, gcv = optimize_smoothing_gcv(channels, energies, weights, s_range)
print(f"最適平滑化パラメータ: s = {optimal_s:.4f}")
3.4 解決策2：ブートストラップ法による不確かさ推定
ブートストラップ法では、測定誤差を考慮して較正データをリサンプリングし、較正曲線の不確かさ分布を推定します。
原理
-
各較正点のチャンネル位置を誤差分布に従ってランダムにずらす：
\text{ch}_i^{(b)} = \text{ch}_i + \epsilon_i^{(b)}, \quad \epsilon_i^{(b)} \sim \mathcal{N}(0, \sigma_{\mu,i}^2) -
リサンプリングしたデータでスプライン曲線を構築：
S^{(b)}(\cdot) = \text{Spline}(\{\text{ch}_i^{(b)}, E_i\}) -
評価点$\text{ch}_*$でのエネルギー予測値を記録
-
1〜3をB回（例：500〜1000回）繰り返す
-
各評価点での平均と標準偏差を計算：
\bar{E}(\text{ch}_*) = \frac{1}{B}\sum_{b=1}^{B}S^{(b)}(\text{ch}_*)\sigma_E(\text{ch}_*) = \sqrt{\frac{1}{B-1}\sum_{b=1}^{B}(S^{(b)}(\text{ch}_*) - \bar{E}(\text{ch}_*))^2}
実装
def bootstrap_spline_uncertainty(channels, energies, channel_errors,
n_bootstrap=1000, eval_points=None):
"""
ブートストラップ法によるスプライン較正の不確かさ推定
Parameters:
-----------
channels : array (n,)
較正点のチャンネル位置
energies : array (n,)
既知のエネルギー
channel_errors : array (n,)
各較正点のチャンネル位置の標準誤差
n_bootstrap : int
ブートストラップの反復回数
eval_points : array or None
不確かさを評価するチャンネル位置
Returns:
--------
eval_points : array
評価点のチャンネル
mean_energies : array
各評価点での平均エネルギー
std_energies : array
各評価点での標準偏差（不確かさ）
percentiles : dict
パーセンタイル（68%, 95%信頼区間など）
"""
if eval_points is None:
eval_points = np.linspace(channels.min(), channels.max(), 200)
bootstrap_results = np.zeros((n_bootstrap, len(eval_points)))
for b in range(n_bootstrap):
resampled_channels = channels + np.random.normal(0, channel_errors)
spline = CubicSpline(resampled_channels, energies)
bootstrap_results[b, :] = spline(eval_points)
mean_energies = np.mean(bootstrap_results, axis=0)
std_energies = np.std(bootstrap_results, axis=0, ddof=1)
percentiles = {
'68%': (np.percentile(bootstrap_results, 16, axis=0),
np.percentile(bootstrap_results, 84, axis=0)),
'95%': (np.percentile(bootstrap_results, 2.5, axis=0),
np.percentile(bootstrap_results, 97.5, axis=0))
}
return eval_points, mean_energies, std_energies, percentiles
eval_ch, mean_en, std_en, conf_int = bootstrap_spline_uncertainty(
channels, energies, channel_errors, n_bootstrap=500
)
target_ch = 450
idx = np.argmin(np.abs(eval_ch - target_ch))
print(f"Ch={target_ch}での較正結果（ブートストラップ法）:")
print(f" エネルギー: {mean_en[idx]:.2f} ± {std_en[idx]:.2f} keV")
print(f" 68%信頼区間: [{conf_int['68%'][0][idx]:.2f}, {conf_int['68%'][1][idx]:.2f}] keV")
print(f" 95%信頼区間: [{conf_int['95%'][0][idx]:.2f}, {conf_int['95%'][1][idx]:.2f}] keV")
ブートストラップ結果の特徴
- 較正点付近では不確かさが小さい：データの制約が強い
- 較正点間では不確かさが増大：スプライン曲線の自由度が高い
- 端点では特に不確かさが大きい：外挿の影響
3.5 解決策3：ヤコビアン法による誤差伝播
数学的な誤差伝播公式を用いて、較正点の誤差を厳密に計算する方法です。
理論
スプライン補間を関数$S(\cdot; \mathbf{ch})$と見なすと、評価点 $\mathbf{ch}*$ でのエネルギー$E*$は較正点のチャンネル位置$\mathbf{ch} = (ch_1, \ldots, ch_n)$の関数です：
E_* = S(\text{ch}_*; \mathbf{ch})
誤差伝播公式により、(E_*)の分散は：
\sigma_{E_*}^2 = \sum_{i=1}^{n} \left(\frac{\partial E_*}{\partial \text{ch}_i}\right)^2 \sigma_{\mu,i}^2
ここで、$\frac{\partial E_*}{\partial \text{ch}_i}$はヤコビアン（感度係数）と呼ばれ、第$i$番目の較正点のチャンネル位置を微小変化させたときのエネルギー変化率です。
数値微分によるヤコビアンの計算
解析的な計算が困難なため、数値微分を用います：
\frac{\partial E_*}{\partial \text{ch}_i} \approx \frac{S(\text{ch}_*; \mathbf{ch} + \delta \mathbf{e}_i) - S(\text{ch}_*; \mathbf{ch})}{\delta}
ここで、$\mathbf{e}_i$は第$i$成分のみが1の単位ベクトル、$\delta$は微小量（例：0.01）です。
def numerical_jacobian_spline(channels, energies, eval_point, delta=0.01):
"""
スプライン補間におけるヤコビアン（感度行列）の数値計算
Parameters:
-----------
channels : array (n,)
較正点のチャンネル位置
energies : array (n,)
既知のエネルギー
eval_point : float
評価するチャンネル位置
delta : float
数値微分のステップサイズ
Returns:
--------
jacobian : array (n,)
各較正点に対する感度 (∂E/∂ch_i)
"""
n = len(channels)
jacobian = np.zeros(n)
base_spline = CubicSpline(channels, energies)
base_value = base_spline(eval_point)
for i in range(n):
perturbed_channels = channels.copy()
perturbed_channels[i] += delta
perturbed_spline = CubicSpline(perturbed_channels, energies)
perturbed_value = perturbed_spline(eval_point)
jacobian[i] = (perturbed_value - base_value) / delta
return jacobian, base_value
def error_propagation_jacobian(channels, energies, channel_errors, eval_point):
"""
ヤコビアン法による誤差伝播
Returns:
--------
energy : float
予測エネルギー
uncertainty : float
不確かさ（標準偏差）
jacobian : array
感度係数
"""
jacobian, energy = numerical_jacobian_spline(channels, energies, eval_point)
variance = np.sum((jacobian ** 2) * (channel_errors ** 2))
uncertainty = np.sqrt(variance)
return energy, uncertainty, jacobian
# 使用例
eval_points_test = [200, 500, 800, 1100]
print("\n=== ヤコビアン法による誤差伝播 ===\n")
for eval_ch in eval_points_test:
energy, error, jac = error_propagation_jacobian(
channels, energies, channel_errors, eval_ch
)
print(f"Ch={eval_ch}:")
print(f" エネルギー: {energy:.2f} ± {error:.2f} keV")
print(f" 感度係数: {jac}")
print(f" 最大感度: ch={channels[np.argmax(np.abs(jac))]}")
print()
感度係数の解釈
感度係数$\frac{\partial E}{\partial \text{ch}_i}$は、各較正点が評価点のエネルギー予測にどの程度影響するかを示します：
- 評価点に近い較正点ほど感度が大きい
- 遠い較正点の感度は小さい（スプラインの局所性）
- 感度の大きさは較正点の分布とスプラインの次数に依存
例：Ch=500で評価する場合
- Ch=500付近の較正点（例：Ch=500）の感度が最大
- Ch=100やCh=1300の感度は非常に小さい
3.6 ブートストラップ法 vs ヤコビアン法の比較
両手法は理論的に等価であり、実際の計算結果もよく一致します：
| 評価点 (Ch) | ヤコビアン (keV) | ブートストラップ (keV) |
|---|---|---|
| 200 | 0.51 | 0.51 |
| 500 | 0.45 | 0.46 |
| 800 | 0.68 | 0.71 |
| 1100 | 1.35 | 1.29 |
選択の指針：
-
ヤコビアン法：
- メリット：計算が高速、数学的に明確
- デメリット：数値微分の精度に依存
- 適用場面：リアルタイム処理、大量の評価点
-
ブートストラップ法：
- メリット：実装が直感的、非線形性に強い、パーセンタイルが得られる
- デメリット：計算コストが高い
- 適用場面：オフライン解析、信頼区間が必要な場合
4. 統計学的背景
4.1 ガウシアンフィッティングの統計誤差
ガウシアン関数によるピークフィッティング：
f(x) = A \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right) + B
パラメータ$\mu$（ピーク位置）の標準誤差は、Cramér-Rao下限により：
\sigma_{\mu} \geq \frac{\sigma}{\sqrt{N}}
ここで、$\sigma$はガウシアンの幅、$N$は総カウント数です。
実際には、最小二乗フィッティングの共分散行列から推定：
from scipy.optimize import curve_fit
def gaussian(x, A, mu, sigma, B):
return A * np.exp(-(x - mu)**2 / (2 * sigma**2)) + B
popt, pcov = curve_fit(gaussian, x_data, y_data, p0=[100, 500, 10, 0])
perr = np.sqrt(np.diag(pcov))
mu_fit = popt[1]
mu_error = perr[1]
print(f"ピーク位置: {mu_fit:.2f} ± {mu_error:.2f} ch")
4.2 誤差伝播の一般論
一般に、関数$f(\mathbf{x})$の不確かさは、入力変数$\mathbf{x} = (x_1, \ldots, x_n)$の不確かさから伝播します：
\sigma_f^2 = \sum_{i=1}^{n} \left(\frac{\partial f}{\partial x_i}\right)^2 \sigma_{x_i}^2 + 2\sum_{i<j}\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\text{Cov}(x_i, x_j)
スプライン較正では、較正点のチャンネル位置が独立と仮定すると、共分散項は消えます：
\sigma_{E_*}^2 = \sum_{i=1}^{n} \left(\frac{\partial E_*}{\partial \text{ch}_i}\right)^2 \sigma_{\mu,i}^2
4.3 平滑化スプラインの統計的定式化
平滑化スプラインは、正則化付き最小二乗法として定式化されます：
\hat{S} = \arg\min_{S \in \mathcal{H}^2} \left\{ \sum_{i=1}^{n} w_i (E_i - S(ch_i))^2 + \lambda \int (S''(ch))^2 dch \right\}
ここで、$\mathcal{H}^2$は2階微分が可積分な関数空間です。
この問題は、線形システム：
(W + \lambda K)^{-1}W\mathbf{E}
の解として表現されます。ここで：
- $W = \text{diag}(w_1, \ldots, w_n)$：重み行列
- $K$：2階微分のペナルティ行列
- $\mathbf{E} = (E_1, \ldots, E_n)^T$：エネルギーベクトル
4.4 ブートストラップ法の理論的正当性
ブートストラップ法は、真の分布が不明な場合に、経験分布からのリサンプリングにより統計量の分布を推定する手法です。
定理（Efron, 1979）：
サンプルサイズ(n)が十分大きいとき、ブートストラップ分布は真の分布に確率収束します。
スプライン較正では、較正点のチャンネル位置の測定誤差分布：
\text{ch}_i \sim \mathcal{N}(\mu_i, \sigma_{\mu,i}^2)
が既知であるため、パラメトリックブートストラップを用いています。これは、経験分布からのリサンプリングではなく、既知の誤差分布からのリサンプリングです。
5. 実践ガイド
5.1 較正手法の選択フローチャート
較正点の誤差が既知か？
├─ YES → 統計誤差を考慮
│ │
│ ├─ 較正点数が多い（≥10点）か？
│ │ ├─ YES → 重み付き平滑化スプライン + ブートストラップ
│ │ └─ NO → 重み付き多項式フィット（2〜3次）
│ │
│ └─ 不確かさ推定が必要か？
│ ├─ YES → ブートストラップ法またはヤコビアン法
│ └─ NO → 重み付き平滑化スプライン
│
└─ NO → 従来手法
│
├─ 較正点数が少ない（≤5点）か？
│ ├─ YES → 多項式フィット（1〜2次）
│ └─ NO → 補間スプライン
│
└─ 評価方法
└─ クロスバリデーション（LOOCV）
5.2 推奨ワークフロー
ステップ1：データの準備
# ガウシアンフィッティング結果
channels = np.array([...]) # ピーク位置（チャンネル）
channel_errors = np.array([...]) # ピーク位置の標準誤差
energies = np.array([...]) # 既知エネルギー（keV）
ステップ2：較正曲線の構築
# 方法1：重み付き平滑化スプライン
weights = 1.0 / (channel_errors ** 2)
smoothing_param = len(channels) * np.mean(channel_errors**2)
calibration_spline = UnivariateSpline(
channels, energies, w=weights, s=smoothing_param
)
# 方法2：補間スプライン（誤差が非常に小さい場合）
calibration_spline = CubicSpline(channels, energies)
ステップ3：適合度評価
# クロスバリデーション
rmse_cv = loocv_spline(channels, energies)
print(f"LOOCV RMSE: {rmse_cv:.4f} keV")
# 物理的制約の確認
check_physical_constraints(calibration_spline, channels.min(), channels.max())
ステップ4：不確かさ推定
# ブートストラップ法
eval_ch, mean_en, std_en, conf_int = bootstrap_spline_uncertainty(
channels, energies, channel_errors, n_bootstrap=500
)
# 結果の保存
np.savez('calibration_result.npz',
eval_channels=eval_ch,
mean_energies=mean_en,
std_energies=std_en,
conf_68=conf_int['68%'],
conf_95=conf_int['95%'])
ステップ5：較正の適用
def apply_calibration(measured_channel, calibration_spline,
eval_ch, std_en):
"""
較正曲線を用いてチャンネルをエネルギーに変換
Parameters:
-----------
measured_channel : float
測定されたチャンネル位置
calibration_spline : CubicSpline or UnivariateSpline
較正曲線
eval_ch : array
ブートストラップ評価点
std_en : array
各評価点での不確かさ
Returns:
--------
energy : float
変換されたエネルギー
uncertainty : float
エネルギーの不確かさ
"""
energy = calibration_spline(measured_channel)
# 線形補間で不確かさを推定
uncertainty = np.interp(measured_channel, eval_ch, std_en)
return energy, uncertainty
# 使用例
measured_ch = 750.5
energy, error = apply_calibration(measured_ch, calibration_spline,
eval_ch, std_en)
print(f"Ch={measured_ch}: E={energy:.2f} ± {error:.2f} keV")
5.3 報告すべき情報
較正結果を論文や報告書で示す際には、以下を含めるべきです：
-
較正データ：
- 使用した既知エネルギーピークとそのエネルギー
- ガウシアンフィッティングで得られたチャンネル位置と誤差
- 較正点数
-
較正手法：
- スプラインの種類（補間/平滑化）
- 平滑化パラメータ（平滑化スプラインの場合）
- 重み付けの有無と方法
-
適合度評価：
- クロスバリデーションのRMSE
- 残差の統計量（平滑化スプラインの場合）
-
不確かさ推定：
- 使用した手法（ブートストラップ/ヤコビアン）
- 代表的なエネルギー点での不確かさ
- 較正範囲と外挿の有無
報告例：
エネルギー較正には、既知の7つの特性X線ピーク（Cu Kα: 8.05 keV、Fe Kα: 6.40 keV、...）を使用した。各ピークの位置はガウシアンフィッティングにより決定し、フィッティング誤差は0.3〜0.6チャンネルであった。較正曲線は重み付き3次スプライン補間（平滑化パラメータ(s=3.5)）により構築した。Leave-One-Out Cross Validationによる評価では、RMSE=1.8 keVであった。ブートストラップ法（500回リサンプリング）により較正曲線の不確かさを推定した結果、較正範囲（100〜1300チャンネル）内での典型的な不確かさは0.4〜1.5 keVであった。
6. まとめ
本記事では、スプライン関数によるエネルギー較正における重要な課題を包括的に解説しました。
主要なポイント
-
スプライン補間を用いる理由：
- 最小二乗法の大域的多項式フィットはRunge現象により振動
- スプラインは区分的多項式のため局所的で安定
- 非線形応答関数に柔軟に対応
-
真の較正関数が不明な場合の評価：
- **クロスバリデーション（LOOCV、K-Fold）**が最も有効
- 残差解析、情報量規準、物理的制約を組み合わせる
- 補間スプラインは較正点で残差≈0となるため残差解析は無効
-
統計誤差の考慮：
- 重み付き平滑化スプラインで測定誤差を反映
- ブートストラップ法で較正曲線全体の不確かさを推定
- ヤコビアン誤差伝播法で数学的に厳密な計算
- 両手法は理論的に等価で結果もよく一致
実践上の推奨
- 較正点の誤差が既知の場合は必ず統計誤差を考慮する
- 適合度評価には複数の手法を組み合わせる
- 不確かさ推定はブートストラップ法またはヤコビアン法を使用
- 較正結果には必ず不確かさを付記する
- 外挿は避け、較正範囲内でのみ使用する
理論的背景
- スプラインは正則化付き最小二乗法として定式化される
- ブートストラップ法は経験分布からの統計的推論
- 誤差伝播は多変数関数の1次テイラー展開に基づく
スプライン較正は強力な手法ですが、統計誤差を適切に扱わなければ、信頼性の低い結果を生む可能性があります。本記事で紹介した手法を組み合わせることで、より信頼性の高いエネルギー較正が実現できます。
参考文献
- de Boor, C. (1978). "A Practical Guide to Splines." Springer-Verlag.
- Efron, B., & Tibshirani, R. J. (1993). "An Introduction to the Bootstrap." Chapman & Hall/CRC.
- Amiri-Simkooei, A., et al. (2025). "Least squares B-spline approximation with applications to geospatial point clouds." Measurement, 242, 115765.
- Tutika, C. S., et al. (2018). "Cubic Spline Interpolation Segmenting over Conventional Polynomial Least Square Methods." arXiv:1803.04621.
- Grimstad, B., & Knudsen, B. R. (2018). "Mathematical Programming Formulations for Piecewise Polynomial Constraints." Optimization Online.
- Lucena, B. (2018). "Spline-Based Probability Calibration." arXiv:1809.07751.
- Van Hoorde, K., et al. (2015). "A spline-based tool to assess and visualize the calibration of multiclass risk models." Journal of Biomedical Informatics, 54, 283-293.
- Enting, I. G. (2006). "Propagating data uncertainty through smoothing spline fits." Tellus B, 58(4), 305-309.
- 北原大地 (2022). "区分的多項式とスプライン関数の基礎." 応用数理, 78(10), 570-578.
- Wahba, G. (1990). "Spline Models for Observational Data." SIAM.
付録：完全な実装例
import numpy as np
from scipy.interpolate import CubicSpline, UnivariateSpline
from sklearn.model_selection import LeaveOneOut
import matplotlib.pyplot as plt
class EnergyCalibration:
"""エネルギー較正の包括的実装"""
def __init__(self, channels, energies, channel_errors=None):
"""
Parameters:
-----------
channels : array (n,)
較正点のチャンネル位置
energies : array (n,)
既知のエネルギー
channel_errors : array (n,) or None
チャンネル位置の標準誤差
"""
self.channels = np.array(channels)
self.energies = np.array(energies)
self.channel_errors = np.array(channel_errors) if channel_errors is not None else None
self.calibration_func = None
self.uncertainty_eval_ch = None
self.uncertainty_std = None
def fit(self, method='smoothing', smoothing_param=None):
"""
較正曲線の構築
Parameters:
-----------
method : str
'interpolation' または 'smoothing'
smoothing_param : float or None
平滑化パラメータ（'smoothing'の場合）
"""
if method == 'interpolation':
self.calibration_func = CubicSpline(self.channels, self.energies)
elif method == 'smoothing':
if self.channel_errors is not None:
weights = 1.0 / (self.channel_errors ** 2)
else:
weights = None
if smoothing_param is None and self.channel_errors is not None:
smoothing_param = len(self.channels) * np.mean(self.channel_errors**2)
self.calibration_func = UnivariateSpline(
self.channels, self.energies,
w=weights, s=smoothing_param
)
else:
raise ValueError(f"Unknown method: {method}")
def evaluate_loocv(self):
"""Leave-One-Out Cross Validation"""
loo = LeaveOneOut()
predictions = []
for train_idx, test_idx in loo.split(self.channels):
train_ch = self.channels[train_idx]
train_en = self.energies[train_idx]
test_ch = self.channels[test_idx]
spline = CubicSpline(train_ch, train_en)
predictions.append(spline(test_ch)[0])
rmse = np.sqrt(np.mean((np.array(predictions) - self.energies)**2))
return rmse
def estimate_uncertainty(self, method='bootstrap', n_bootstrap=500, eval_points=None):
"""
不確かさの推定
Parameters:
-----------
method : str
'bootstrap' または 'jacobian'
n_bootstrap : int
ブートストラップの反復回数
eval_points : array or None
評価点
"""
if self.channel_errors is None:
raise ValueError("Channel errors are required for uncertainty estimation")
if eval_points is None:
eval_points = np.linspace(self.channels.min(), self.channels.max(), 200)
if method == 'bootstrap':
bootstrap_results = np.zeros((n_bootstrap, len(eval_points)))
for b in range(n_bootstrap):
resampled_ch = self.channels + np.random.normal(0, self.channel_errors)
spline = CubicSpline(resampled_ch, self.energies)
bootstrap_results[b, :] = spline(eval_points)
mean_en = np.mean(bootstrap_results, axis=0)
std_en = np.std(bootstrap_results, axis=0, ddof=1)
elif method == 'jacobian':
std_en = np.zeros(len(eval_points))
for j, eval_ch in enumerate(eval_points):
jacobian = self._compute_jacobian(eval_ch)
variance = np.sum((jacobian ** 2) * (self.channel_errors ** 2))
std_en[j] = np.sqrt(variance)
mean_en = self.calibration_func(eval_points)
else:
raise ValueError(f"Unknown method: {method}")
self.uncertainty_eval_ch = eval_points
self.uncertainty_std = std_en
return eval_points, mean_en, std_en
def _compute_jacobian(self, eval_point, delta=0.01):
"""ヤコビアンの計算"""
n = len(self.channels)
jacobian = np.zeros(n)
base_spline = CubicSpline(self.channels, self.energies)
base_value = base_spline(eval_point)
for i in range(n):
perturbed_ch = self.channels.copy()
perturbed_ch[i] += delta
perturbed_spline = CubicSpline(perturbed_ch, self.energies)
perturbed_value = perturbed_spline(eval_point)
jacobian[i] = (perturbed_value - base_value) / delta
return jacobian
def predict(self, channel, return_uncertainty=False):
"""
チャンネルをエネルギーに変換
Parameters:
-----------
channel : float or array
チャンネル位置
return_uncertainty : bool
不確かさを返すか
Returns:
--------
energy : float or array
エネルギー
uncertainty : float or array (return_uncertainty=Trueの場合)
不確かさ
"""
energy = self.calibration_func(channel)
if return_uncertainty:
if self.uncertainty_eval_ch is None:
raise ValueError("Run estimate_uncertainty() first")
uncertainty = np.interp(channel, self.uncertainty_eval_ch, self.uncertainty_std)
return energy, uncertainty
return energy
def plot_calibration(self, show_uncertainty=True, n_points=300):
"""較正曲線と不確かさのプロット"""
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# 較正曲線
test_ch = np.linspace(self.channels.min(), self.channels.max(), n_points)
test_en = self.calibration_func(test_ch)
ax1.plot(test_ch, test_en, 'b-', label='較正曲線', linewidth=2)
ax1.plot(self.channels, self.energies, 'ro', markersize=8, label='較正点')
if self.channel_errors is not None:
ax1.errorbar(self.channels, self.energies, xerr=self.channel_errors,
fmt='none', ecolor='red', capsize=5)
ax1.set_xlabel('チャンネル', fontsize=12)
ax1.set_ylabel('エネルギー (keV)', fontsize=12)
ax1.set_title('エネルギー較正曲線', fontsize=14)
ax1.legend()
ax1.grid(True, alpha=0.3)
# 不確かさ
if show_uncertainty and self.uncertainty_eval_ch is not None:
ax2.plot(self.uncertainty_eval_ch, self.uncertainty_std, 'g-', linewidth=2)
ax2.set_xlabel('チャンネル', fontsize=12)
ax2.set_ylabel('不確かさ (keV)', fontsize=12)
ax2.set_title('較正曲線の不確かさ', fontsize=14)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
def summary(self):
"""較正結果のサマリー"""
print("="*60)
print("エネルギー較正サマリー")
print("="*60)
print(f"較正点数: {len(self.channels)}")
print(f"較正範囲: {self.channels.min():.1f} - {self.channels.max():.1f} ch")
print(f"エネルギー範囲: {self.energies.min():.1f} - {self.energies.max():.1f} keV")
if self.channel_errors is not None:
print(f"チャンネル誤差範囲: {self.channel_errors.min():.2f} - {self.channel_errors.max():.2f} ch")
if self.calibration_func is not None:
loocv_rmse = self.evaluate_loocv()
print(f"\nLOOCV RMSE: {loocv_rmse:.4f} keV")
if self.uncertainty_std is not None:
print(f"\n不確かさ範囲: {self.uncertainty_std.min():.4f} - {self.uncertainty_std.max():.4f} keV")
print("="*60)
# 使用例
if __name__ == "__main__":
channels = np.array([100, 300, 500, 700, 900, 1100, 1300])
energies = np.array([60, 240, 500, 840, 1260, 1760, 2340])
channel_errors = np.array([0.5, 0.4, 0.3, 0.35, 0.4, 0.5, 0.6])
calib = EnergyCalibration(channels, energies, channel_errors)
calib.fit(method='smoothing', smoothing_param=3.5)
eval_ch, mean_en, std_en = calib.estimate_uncertainty(method='bootstrap', n_bootstrap=500)
calib.summary()
calib.plot_calibration(show_uncertainty=True)
test_channel = 450
energy, uncertainty = calib.predict(test_channel, return_uncertainty=True)
print(f"\nCh={test_channel}: E={energy:.2f} ± {uncertainty:.2f} keV")
このクラスを使用することで、エネルギー較正の全工程を統一的に処理できます。