模型数据最小二乘拟合¶

这是一篇针对物理学家 (例如物理学家、天文学家) 或工程师的 statsmodels 快速入门。

为什么需要它？

因为大多数 statsmodels 是由统计学家编写的，他们使用不同的术语，有时还使用不同的方法，因此很难知道哪些类和函数是相关的，以及它们的输入和输出含义。

[1]:

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

线性模型¶

假设您有数据点，测量值为 y，位置为 x，测量误差为 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 月 3 日