考虑到误差线，如何进行线性回归？答案

【问题标题】：How to do linear regression, taking errorbars into account?考虑到误差线，如何进行线性回归？
【发布时间】：2014-02-23 12:44:30
【问题描述】：

我正在对一些有限大小的物理系统进行计算机模拟，然后我将外推到无穷大（热力学极限）。一些理论说数据应该随着系统大小线性扩展，所以我在做线性回归。

我拥有的数据很嘈杂，但对于每个数据点，我都可以估计误差线。因此，例如数据点如下所示：

x_list = [0.3333333333333333, 0.2886751345948129, 0.25, 0.23570226039551587, 0.22360679774997896, 0.20412414523193154, 0.2, 0.16666666666666666]
y_list = [0.13250359351851854, 0.12098339583333334, 0.12398501145833334, 0.09152715, 0.11167239583333334, 0.10876248333333333, 0.09814170444444444, 0.08560799305555555]
y_err = [0.003306749165349316, 0.003818446389148108, 0.0056036878203831785, 0.0036635292592592595, 0.0037034897788415424, 0.007576672222222223, 0.002981084130692832, 0.0034913019065973983]

假设我正在尝试在 Python 中执行此操作。

我知道的第一种方法是：
```
m, c, r_value, p_value, std_err = scipy.stats.linregress(x_list, y_list)
```
我知道这给了我结果的错误栏，但这没有考虑到初始数据的错误栏。

我知道的第二种方法是：

m, c = numpy.polynomial.polynomial.polyfit(x_list, y_list, 1, w = [1.0 / ty for ty in y_err], full=False)

在这里，我们使用每个点的误差条的倒数作为权重，用于最小二乘近似。所以如果一个点不是真的那么可靠，它不会对结果有很大的影响，这是合理的。

但我不知道如何获得结合这两种方法的东西。

我真正想要的是 second 方法所做的，这意味着当每个点以不同的权重影响结果时使用回归。但同时我想知道我的结果有多准确，也就是说，我想知道结果系数的误差线是什么。

我该怎么做？

【问题讨论】：

是我误解了你，还是你想使用y_err 系列作为权重矩阵？

标签： python numpy linear-regression least-squares extrapolation

【解决方案1】：

不完全确定这是否是您的意思，但是……使用 pandas、statsmodels 和 patsy，我们可以比较普通的最小二乘拟合和加权最小二乘拟合，后者使用您提供的噪声的倒数作为权重矩阵（顺便说一下，statsmodels 会抱怨样本量

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 300

import statsmodels.formula.api as sm

x_list = [0.3333333333333333, 0.2886751345948129, 0.25, 0.23570226039551587, 0.22360679774997896, 0.20412414523193154, 0.2, 0.16666666666666666]
y_list = [0.13250359351851854, 0.12098339583333334, 0.12398501145833334, 0.09152715, 0.11167239583333334, 0.10876248333333333, 0.09814170444444444, 0.08560799305555555]
y_err = [0.003306749165349316, 0.003818446389148108, 0.0056036878203831785, 0.0036635292592592595, 0.0037034897788415424, 0.007576672222222223, 0.002981084130692832, 0.0034913019065973983]

# put x and y into a pandas DataFrame, and the weights into a Series
ws = pd.DataFrame({
    'x': x_list,
    'y': y_list
})
weights = pd.Series(y_err)

wls_fit = sm.wls('x ~ y', data=ws, weights=1 / weights).fit()
ols_fit = sm.ols('x ~ y', data=ws).fit()

# show the fit summary by calling wls_fit.summary()
# wls fit r-squared is 0.754
# ols fit r-squared is 0.701

# let's plot our data
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, facecolor='w')
ws.plot(
    kind='scatter',
    x='x',
    y='y',
    style='o',
    alpha=1.,
    ax=ax,
    title='x vs y scatter',
    edgecolor='#ff8300',
    s=40
)

# weighted prediction
wp, = ax.plot(
    wls_fit.predict(),
    ws['y'],
    color='#e55ea2',
    lw=1.,
    alpha=1.0,
)
# unweighted prediction
op, = ax.plot(  
    ols_fit.predict(),
    ws['y'],
    color='k',
    ls='solid',
    lw=1,
    alpha=1.0,
)
leg = plt.legend(
    (op, wp),
    ('Ordinary Least Squares', 'Weighted Least Squares'),
    loc='upper left',
    fontsize=8)

plt.tight_layout()
fig.set_size_inches(6.40, 5.12)
plt.show()

WLS 残差：

[0.025624005084707302,
 0.013611438189866154,
 -0.033569595462217161,
 0.044110895217014695,
 -0.025071632845910546,
 -0.036308252199571928,
 -0.010335514810672464,
 -0.0081511479431851663]

加权拟合（wls_fit.mse_resid或wls_fit.scale）的残差均方误差为0.22964802498892287，拟合的r平方值为0.754.

如果您需要每个可用属性和方法的列表，您可以通过调用他们的summary() 方法和/或执行dir(wls_fit) 来获取有关拟合的大量数据。

【讨论】：

你确定参数weights应该设置为1/y_err吗？ statsmodel 页面上的示例使用 weights=1/(w**2)。他们写道：“在这个例子中，w 是误差的标准差。WLS 要求权重与误差方差的倒数成正比。” statsmodels.org/stable/examples/notebooks/generated/wls.html
@Jim 我不确定！我不能确定 statsmodels 是否在 7 年前就 WLS 中的误差方差提出了建议，所以我只是使用我在 intro stats 中学到的知识来指定模型。您是否尝试过使用 weights = 1 / (w ** 2) 运行示例？很高兴以任何一种方式纠正答案。
我还没有尝试以 1/(w**2) 的重量运行它。我觉得应该差不多吧。而且我不确定 1(w**2) 是否正确。我想在过去的某个时候，我决定使用 1/(w**2) 进行分析，但那是很久以前的事了。

【解决方案2】：

我写了一个简洁的函数来执行数据集的加权线性回归，它是GSL's "gsl_fit_wlinear" function的直接翻译。如果您想确切地知道函数在执行拟合时正在做什么，这很有用

def wlinear_fit (x,y,w) :
    """
    Fit (x,y,w) to a linear function, using exact formulae for weighted linear
    regression. This code was translated from the GNU Scientific Library (GSL),
    it is an exact copy of the function gsl_fit_wlinear.
    """
    # compute the weighted means and weighted deviations from the means
    # wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i)
    W = np.sum(w)
    wm_x = np.average(x,weights=w)
    wm_y = np.average(y,weights=w)
    dx = x-wm_x
    dy = y-wm_y
    wm_dx2 = np.average(dx**2,weights=w)
    wm_dxdy = np.average(dx*dy,weights=w)
    # In terms of y = a + b x
    b = wm_dxdy / wm_dx2
    a = wm_y - wm_x*b
    cov_00 = (1.0/W) * (1.0 + wm_x**2/wm_dx2)
    cov_11 = 1.0 / (W*wm_dx2)
    cov_01 = -wm_x / (W*wm_dx2)
    # Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2
    chi2 = np.sum (w * (y-(a+b*x))**2)
    return a,b,cov_00,cov_11,cov_01,chi2

为了锻炼身体，你会这样做

a,b,cov_00,cov_11,cov_01,chi2 = wlinear_fit(x_list,y_list,1.0/y_err**2)

这将返回线性回归的系数a（截距）和b（斜率）的最佳估计值，以及协方差矩阵cov_00、cov_01 和@987654328 的元素@。 a 上误差的最佳估计是 cov_00 的平方根，b 上的误差是 cov_11 的平方根。残差的加权和在chi2 变量中返回。

重要提示：此函数接受逆方差，而不是逆标准差作为数据点的权重。

【讨论】：

【解决方案3】：

sklearn.linear_model.LinearRegression 支持在fit 期间指定权重：

x_data = np.array(x_list).reshape(-1, 1)  # The model expects shape (n_samples, n_features).
y_data = np.array(y_list)
y_err  = np.array(y_err)

model = LinearRegression()
model.fit(x_data, y_data, sample_weight=1/y_err)

这里的样本权重指定为1 / y_err。可能有不同的版本，通常最好将这些样本权重裁剪为最大值，以防y_err 变化很大或有小的异常值：

sample_weight = 1 / y_err
sample_weight = np.minimum(sample_weight, MAX_WEIGHT)

应根据您的数据确定MAX_WEIGHT 的位置（通过查看y_err 或1 / y_err 分布，例如，如果它们有异常值，则可以对其进行裁剪）。

【讨论】：

【解决方案4】：

我发现this 文档有助于理解和设置我自己的加权最小二乘例程（适用于任何编程语言）。

通常学习和使用优化的例程是最好的方法，但有时了解例程的内容很重要。

【讨论】：