模型對數據的最小平方法擬合¶

這是對物理科學家（例如物理學家、天文學家）或工程師快速介紹 statsmodels 的方法。

為什麼需要這個？

因為大多數 statsmodels 是由統計學家編寫的，他們使用不同的術語，有時使用不同的方法，這使得很難知道哪些類別和函數是相關的，以及它們的輸入和輸出意味著什麼。

[1]:

import numpy as np
import pandas as pd
import statsmodels.api as sm

線性模型¶

假設您在位置 x 處有測量值 y 的數據點，以及測量誤差 y_err。

您如何使用 statsmodels 將直線模型擬合到此數據？

如需詳細討論，請參閱 Hogg 等人 (2010)，“數據分析秘訣：將模型擬合到數據”… 我們將使用他們在表 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。添加此列 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:                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

上次更新時間：2024 年 10 月 03 日