非線形フィットでは、success=Trueやパラメータの値だけで採用を決めない。least_squaresへ残差ベクトルを渡し、範囲制約、パラメータ尺度、終了理由、残差の形を確認する。残差が小さくても、モデルが物理的に正しいことや、パラメータが一意に決まることまで示したわけではない。
最小例で確かめる
以下は説明用に作った小さなデータである。実測データや実行速度の測定結果ではない。コード全体をexample.pyとして保存すれば、入力ファイルを別途用意せずに実行できる。assertは、この例で成り立つべき形や値を確認するために入れてある。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
t = np.linspace(0., 4., 21)
def model(p, t):
a, k, c = p
return a * np.exp(-k * t) + c
y = model([2.5, 0.8, 0.3], t) + 0.04 * np.sin(5 * t) # 決定的な模擬値
def residual(p):
return model(p, t) - y
options = dict(bounds=([0., 0., 0.], [10., 3., 2.]),
x_scale=[2., 1., 0.5], max_nfev=1000)
fit = least_squares(residual, [1., 0.2, 0.1], **options)
budget = least_squares(residual, [1., 0.2, 0.1], max_nfev=1,
bounds=options["bounds"], x_scale=options["x_scale"])
def no_offset(p):
return p[0] * np.exp(-p[1] * t) - y
wrong_model = least_squares(no_offset, [1., 0.2], bounds=([0.,0.],[10.,3.]))
rms = np.sqrt(np.mean(fit.fun**2))
wrong_rms = np.sqrt(np.mean(wrong_model.fun**2))
print("success / status:", fit.success, fit.status)
print("parameters:", fit.x.round(6).tolist())
print("RMS:", round(float(rms), 6))
print("without offset RMS:", round(float(wrong_rms), 6))
print("budget-limited success / status:", budget.success, budget.status)
assert fit.success and np.isfinite(fit.x).all()
assert fit.fun.shape == y.shape
np.testing.assert_allclose(2 * fit.cost, np.sum(fit.fun**2))
assert rms < 0.04 and wrong_rms > rms
assert not budget.success and budget.status == 0
fig, axes = plt.subplots(2, 1, figsize=(7, 6), sharex=True, layout="constrained")
axes[0].plot(t, y, "o", label="synthetic observations", color="#26778c")
axes[0].plot(t, model(fit.x, t), label="with offset", color="#a83262")
axes[0].plot(t, y + wrong_model.fun, "--", label="without offset", color="#e38619")
axes[0].set(ylabel="Signal [arb. unit]", title="Fit and residual diagnostics")
axes[0].legend()
axes[1].axhline(0, color="0.5", linewidth=1)
axes[1].plot(t, fit.fun, "o-", color="#a83262", label="with offset")
axes[1].plot(t, wrong_model.fun, "s--", color="#e38619", label="without offset")
axes[1].set(xlabel="Time [s]", ylabel="Model - observation")
axes[1].legend()
fig.savefig("PY084-residuals.png", dpi=160)
plt.close(fig)
実行結果
success / status: True 1
parameters: [2.512365, 0.811382, 0.306866]
RMS: 0.027453
without offset RMS: 0.081391
budget-limited success / status: False 0
残差をスカラーへ先にまとめない
least_squaresには、各観測点についてモデル値から観測値を引いた1次元配列を返す関数を渡す。自分で二乗和を計算して一つの値にすると、ソルバーへ渡す問題の形が変わってしまう。標準のlinear損失ではcostが残差二乗和の半分になるので、例ではこの関係も検査した。
模擬値は指数減衰に定数項と決定的な正弦のずれを加えたもので、実測値でも乱数による試験結果でもない。初期値は生成に使った値と変え、推定後に残差を確認している。推定されたパラメータが生成時の値と完全一致しないのは、意図的なずれが入っているためである。
制約と尺度を問題に合わせる
boundsは各パラメータの下限・上限であり、ここでは振幅、減衰率、定数項を非負に制限した。物理的に許される範囲が分かるときに役立つが、都合のよい形へ固定するための道具ではない。解が境界へ張り付くなら、モデル、単位、初期値、境界の妥当性を点検する。
x_scaleはパラメータの代表的な尺度を与える。値の桁が大きく違うパラメータを同じように動かすと、探索が進みにくい場合がある。ただし尺度を指定すれば必ず大域的な最適解へ到達するわけではない。初期値を複数変えて結果が安定するか、識別に必要な観測範囲があるかも調べたい。
終了理由と残差の形を確認する
予算をmax_nfev=1に制限したbudgetは、有限なパラメータを持っていてもsuccess=False、status=0となる。数字が返ったことを成功判定にしてはいけない。fit.messageには詳細な終了理由があるので、処理記録に残すと後から診断しやすい。許容誤差による停止も、残差が十分小さいという意味と同一ではない。
図の下段は残差を時刻に対して描いたものだ。定数項を省いたモデルではRMSが大きくなり、時刻に沿った偏りも見える。この例は省略の影響を示すために作ったデータなので、実測データで同じ形が出れば原因は必ず定数項、と断定できるわけではない。系統的な形は仮定を見直す手掛かりである。
小さい残差だけでモデルを選ばない
パラメータを増やせば残差が減る場合があるが、過剰適合や強いパラメータ相関が生じることもある。未使用データでの予測、残差の分散、測定誤差、モデルの根拠を併せて考える。least_squaresの結果から、そのまま厳密な信頼区間が得られるわけでもない。
観測ごとに精度が違うなら、残差を標準偏差などで割る重み付けを検討する。ただし重みの根拠を持つことが必要である。外れ値用の損失を採用すると目的関数の意味も変わるので、通常の二乗和と区別して報告する。まず入力の有限性、残差のshape、終了状態、残差の図を確認し、その先でモデルの妥当性を判断しよう。
動作確認環境と参考資料
Linux・CPython 3.12.14、NumPy 2.3.5、pandas 2.2.3、SciPy 1.17.0、Matplotlib 3.10.8の環境で掲載コードを実行した。使用するライブラリはコード冒頭のimportを参照してほしい。公式資料の最新版と、この実行確認版は区別している。数値の末尾や表の表示幅は環境によって変わることがある。