モデルのデータへの最小二乗フィッティング¶
これは、物理科学者（例えば、物理学者、天文学者）やエンジニア向けの 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.11 * b = 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()
/tmp/ipykernel_4678/2751177292.py:28: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
data = pd.read_csv(StringIO(data), delim_whitespace=True).astype(float)
[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 = y * weights = 1 / sqrt(y_err)
exog は、x を列として持ち、1の追加列を持つ2次元配列でなければならないことに注意してください。この1の列を追加することは、モデル y = a * x + b をフィッティングしたいことを意味し、これを省略することは、モデル y = a * x をフィッティングしたいことを意味します。
また、statsmodels に、絶対スケールを持つ測定誤差が実際にあることを伝えるために、cov_type='fixed scale' オプションを使用する必要があります。そうしない場合、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: Thu, 03 Oct 2024 Prob (F-statistic): 4.52e-11
Time: 15:50:45 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