使用python进行线性拟合的误差传播答案

【问题标题】：Error propagation in a linear fit using python使用python进行线性拟合的误差传播
【发布时间】：2018-10-13 04:41:38
【问题描述】：

假设我对一些因变量 y 相对于一些自变量 x 进行了多次测量。我还记录了每次测量中的不确定性dy。例如，这可能看起来像

import numpy as np
x = np.array([1, 2, 3, 4])
y = np.array([4.1, 5.8, 8.1, 9.7])
dy = np.array([0.2, 0.3, 0.2, 0.4])

现在假设我希望测量值遵循线性关系 y = mx + b，并且我想确定一些未测量的 x 值 x_unm 的 y 值 y_umn。如果不考虑错误，我可以很容易地在 Python 中执行线性拟合：

fit_params, residuals, rank, s_values, rcond = np.polyfit(x, y, 1, full=True)
poly_func = np.poly1d(fit_params)

x_unm   # The unmeasured x value
y_unm = poly_func(x_unm)  # The unmeasured x value

我对这种方法有两个问题。首先是np.polyfit 没有包含每个点的错误。其次，我不知道y_unm 的不确定性是什么。

有谁知道如何用不确定性拟合数据以使我能够确定y_unm 中的不确定性？

【问题讨论】：

标签： python numpy statistics linear-regression

【解决方案1】：

这是一个可以通过分析解决的问题，但这可能更适合作为数学/统计讨论。例如参见（在许多来源中）：

https://www.che.udel.edu/pdf/FittingData.pdf

拟合误差可以通过解析计算。重要的是要注意，在考虑测量误差时，拟合本身是不同的。

在 python 中，我不确定处理错误的内置函数，但这里有一个使用 scipy.optimize.fmin 进行卡方最小化的示例

#Calculate Chi^2 function to minimize
def chi_2(params,x,y,sigy):
    m,c=params
    return sum(((y-m*x-c)/sigy)**2)

data_in=(x,y,dy)
params0=[1,0]

q=fmin(chi_2,params0,args=data_in)

为了比较，我使用了这个、你的 polyfit 解决方案和解析解决方案，并为你提供的数据绘制了图表。

给定技术的参数结果：

带 fmin 的加权卡方： m=1.94609996 b=2.1312239

分析： m=1.94609929078014 b=2.1312056737588647

多合体：米=1.91 b=2.15

Linear fits to given data

这里是完整的代码：

import numpy as np
from scipy.optimize import fmin
import matplotlib.pyplot as plt

x = np.array([1, 2, 3, 4])
y = np.array([4.1, 5.8, 8.1, 9.7])
dy = np.array([0.2, 0.3, 0.2, 0.4]) 

#Calculate Chi^2 function to minimize
def chi_2(params,x,y,sigy):
    m,c=params
    return sum(((y-m*x-c)/sigy)**2)

data_in=(x,y,dy)
params0=[1,0]

q=fmin(chi_2,params0,args=data_in)

#Unweighted fit to compare

a=np.polyfit(x,y,deg=1)

#Analytic solution
sx=sum(x/dy**2)
sx2=sum(x**2/dy**2)
s1=sum(1./dy**2)
sy=sum(y/dy**2)
sxy=sum(x*y/dy**2)

ma=(s1*sxy-sx*sy)/(s1*sx2-sx**2)
ba=(sx2*sy-sx*sxy)/(sx2*s1-sx**2)

xplt=np.linspace(0,5,100)
yplt1=xplt*q[0]+q[1]


yplt2=xplt*a[0]+a[1]

yplt3=xplt*ma+ba

plt.figure()
plt.plot(xplt,yplt1,label='Error Weighted',color='black')
plt.plot(xplt,yplt2,label='Non-Error Weighted',color='blue')
plt.plot(xplt,yplt3,label='Error Weighted Analytic',linestyle='--',color='red')
plt.errorbar(x,y,yerr=dy,fmt='ko')
plt.legend()
plt.show()

【讨论】：