我在做一个有限大小的物理系统的计算机模拟,然后我在做外推到无穷大(热力学极限)。一些理论认为数据应该与系统规模成线性关系,所以我做的是线性回归。
我得到的数据是有噪声的,但对于每个数据点,我可以估计误差线。例如,数据点看起来像:
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方法所做的,这意味着当每个点以不同的权重影响结果时使用回归。但同时我想知道我的结果有多精确,也就是说,我想知道得到的系数的误差条是多少。
我该怎么做?
不完全确定这是否是您的意思,但是…使用pandas、statsmodels和patsy,我们可以比较普通的最小二乘拟合和使用您提供的噪声逆作为权重矩阵的加权最小二乘拟合(顺便说一下,statsmodels会抱怨样本大小<;20)。
WLS残差:
加权拟合(
wls_fit.mse_resid
或wls_fit.scale
)残差的均方误差为0.22964824982287,拟合的r平方值为0.754。如果您需要每个可用属性和方法的列表,可以通过调用fits的
summary()
方法和/或执行dir(wls_fit)
来获得关于fits的大量数据。我编写了一个简洁的函数来执行数据集的加权线性回归,这是GSL's "gsl_fit_wlinear" function的直接转换。如果您想确切地知道函数在执行fit时正在执行的操作,则这非常有用
为了锻炼身体,你应该
它将返回线性回归系数
a
(截距)和b
(斜率)的最佳估计值,以及协方差矩阵cov_00
、cov_01
和cov_11
的元素。对a
上的误差的最佳估计是cov_00
的平方根,对b
的最佳估计是cov_11
的平方根。残差的加权和在chi2
变量中返回。重要信息:此函数接受逆方差,而不是逆标准差作为数据点的权重。
我发现this文档有助于理解和设置我自己的加权最小二乘程序(适用于任何编程语言)。
通常,学习和使用优化的程序是最好的方法,但有时了解程序的勇气是重要的。
相关问题 更多 >
编程相关推荐