如何进行线性回归,将误码率考虑在内？

Question

我正在为一些有限大小的物理系统进行计算机模拟,之后我正在对无穷大进行外推(热力学极限).一些理论认为数据应该随着系统规模线性扩展,所以我做的是线性回归.

我的数据有噪音,但对于每个数据点,我可以估算出错误.因此,例如数据点看起来像:

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中执行此操作.

1. 我知道的第一种方式是:

   m, c, r_value, p_value, std_err = scipy.stats.linregress(x_list, y_list)
   

   我理解这给了我结果的错误栏,但这没有考虑初始数据的错误栏.

2. 我知道的第二种方式是:

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

这里我们使用每个点的误差条的倒数作为在最小二乘近似中使用的权重.因此,如果一个点不是那么可靠,那么它不会对结果造成太大影响,这是合理的.

但我无法弄清楚如何获得结合这两种方法的东西.

我真正想要的是第二种方法的作用,意思是当每个点都影响不同权重的结果时使用回归.但与此同时,我想知道我的结果有多准确,这意味着,我想知道结果系数的误码是什么.

我怎样才能做到这一点？

Answer 1

不完全确定这是不是你的意思,但是......使用pandas,statsmodels和patsy,我们可以比较普通的最小二乘拟合和加权最小二乘拟合,它使用你提供的噪声的倒数作为权重矩阵(顺便说一句,statsmodels会抱怨样本量<20.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
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, axisbg='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.savefig("so.png", dpi=100, alpha=True)
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).