データへの最小二乗法によるモデルのフィッティング¶
これは、物理学者、天文学者、またはエンジニアなどの物理科学者向けに、statsmodelsの簡単な紹介です。
なぜこれが必要なのか？
それは、statsmodelsのほとんどが統計学者によって書かれており、彼らは異なる用語や時には異なる手法を使用しているため、どのクラスや関数が関連しているのか、またその入力や出力が何を意味するのかを理解するのが難しいからです。
[1]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
線形モデル¶
次のようなデータポイントがあり、位置 x での測定値 y と測定誤差 y_err があるとします。
このデータに直線モデルを適合させるために、statsmodels をどのように使用できるのでしょうか？
詳細な議論については、Hogg et al. (2010), 「Data analysis recipes: Fitting a model to data」 を参照してください… ここでは、彼らが表 1 に示した例のデータを使用します。
したがって、モデルは f(x) = a * x + b であり、図 1では、再現したい結果が印刷されています… このデータに対する「標準加重最小二乗法フィット」の最適パラメータとパラメータ誤差は以下の通りです：
a = 2.24 ± 0.11b = 34 ± 18
[2]:
data = """
x y y_err
201 592 61
244 401 25
47 583 38
287 402 15
203 495 21
58 173 15
210 479 27
202 504 14
198 510 30
158 416 16
165 393 14
201 442 25
157 317 52
131 311 16
166 400 34
160 337 31
186 423 42
125 334 26
218 533 16
146 344 22
"""
try:
from StringIO import StringIO
except ImportError:
from io import StringIO
data = pd.read_csv(StringIO(data), delim_whitespace=True).astype(float)
# Note: for the results we compare with the paper here, they drop the first four points
data.head()
[2]:
| x | y | y_err | |
|---|---|---|---|
| 0 | 201.0 | 592.0 | 61.0 |
| 1 | 244.0 | 401.0 | 25.0 |
| 2 | 47.0 | 583.0 | 38.0 |
| 3 | 287.0 | 402.0 | 15.0 |
| 4 | 203.0 | 495.0 | 21.0 |
直線をフィットさせるには、加重最小二乗法クラス WLS を使用します。パラメータは次のように設定されます：
exog=sm.add_constant(x)endog=yweights=1 / sqrt(y_err)
exog は x を列として、もう一つの列に 1 を加えた 2 次元の配列である必要があることに注意してください。この 1 の列を加えることで、モデル y = a * x + b をフィットさせることになります。これを省略すると、モデル y = a * x をフィットさせることになります。
また、cov_type='fixed scale' オプションを使用して、statsmodels に絶対スケールでの測定誤差があることを伝える必要があります。これを指定しないと、statsmodels は重みをデータポイント間の相対的な重みとして扱い、内部でスケーリングを再調整して最適なモデルの chi**2 / ndf = 1 となるようにします。
[3]:
exog = sm.add_constant(data["x"])
endog = data["y"]
weights = 1.0 / (data["y_err"] ** 2)
wls = sm.WLS(endog, exog, weights)
results = wls.fit(cov_type="fixed scale")
print(results.summary())
WLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.400
Model: WLS Adj. R-squared: 0.367
Method: Least Squares F-statistic: 193.5
Date: Wed, 28 Aug 2024 Prob (F-statistic): 4.52e-11
Time: 10:25:16 Log-Likelihood: -119.06
No. Observations: 20 AIC: 242.1
Df Residuals: 18 BIC: 244.1
Df Model: 1
Covariance Type: fixed scale
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 213.2735 14.394 14.817 0.000 185.062 241.485
x 1.0767 0.077 13.910 0.000 0.925 1.228
==============================================================================
Omnibus: 0.943 Durbin-Watson: 2.901
Prob(Omnibus): 0.624 Jarque-Bera (JB): 0.181
Skew: -0.205 Prob(JB): 0.914
Kurtosis: 3.220 Cond. No. 575.
==============================================================================
Notes:
[1] Standard Errors are based on fixed scale
scipy.optimize.curve_fitとの比較¶
[4]:
# You can use `scipy.optimize.curve_fit` to get the best-fit parameters and parameter errors.
from scipy.optimize import curve_fit
def f(x, a, b):
return a * x + b
xdata = data["x"]
ydata = data["y"]
p0 = [0, 0] # initial parameter estimate
sigma = data["y_err"]
popt, pcov = curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))
print("a = {0:10.3f} +- {1:10.3f}".format(popt[0], perr[0]))
print("b = {0:10.3f} +- {1:10.3f}".format(popt[1], perr[1]))
a = 1.077 +- 0.077
b = 213.273 +- 14.394
自作のコスト関数との照合¶
[5]:
# You can also use `scipy.optimize.minimize` and write your own cost function.
# This does not give you the parameter errors though ... you'd have
# to estimate the HESSE matrix separately ...
from scipy.optimize import minimize
def chi2(pars):
"""Cost function."""
y_model = pars[0] * data["x"] + pars[1]
chi = (data["y"] - y_model) / data["y_err"]
return np.sum(chi ** 2)
result = minimize(fun=chi2, x0=[0, 0])
popt = result.x
print("a = {0:10.3f}".format(popt[0]))
print("b = {0:10.3f}".format(popt[1]))
a = 1.077
b = 213.274
非線形モデル¶
[6]:
# TODO: we could use the examples from here:
# http://probfit.readthedocs.org/en/latest/api.html#probfit.costfunc.Chi2Regression