python - 使用python进行线性拟合的误差传播
问题描述
假设我对一些因变量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
,并且我想确定y_umn
一些未测量的 x 值的 y 值x_unm
。如果我不考虑错误,我可以很容易地在 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
?
解决方案
这是一个可以通过分析完成的问题,但这可能更适合作为数学/统计讨论。例如参见(在许多来源中):
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
Polyfit:m=1.91 b=2.15
这是完整的代码:
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()
推荐阅读
- javascript - TypeScript 和 Redux:为什么我需要添加 ` | undefined` 到我的 Reducer 状态类型?
- npm - npm install 在 node_modules 的父文件夹中下载额外的依赖项
- c# - 视频转换:从 UDP 套接字到 Mat 对象的 h.264
- vb6 - “模板”字符串中的嵌套替换
- apache-kafka - Kafka:每个 poll() 调用是否一次只使用一个主题?
- c++ - 尝试编译简单的测试程序时 CMake 失败
- django - 如何在表单中检索 id 作为隐藏值,以便将该值作为外键存储在 django 的数据库中
- sql - 需要在查询中查找和替换文本
- powershell - 如何获取通讯组列表中联系人的手机号码?
- r - 在 Rstudio 中运行简单的 xml 代码后出现错误