我有一个函数f(x) = a/x
,我有一组包含f(x) +- df(x)
和值的数据x +- dx
。我如何告诉 gnuplot 进行加权拟合a
?
我知道fit
接受这个using
术语,这适用于df(x)
,但它不适用于dx
。似乎 gnuplot 将我所遇到的错误视为我的函数 ( )x
的整个 RHS 的错误。a/x +- dx
如何进行加权拟合f(x) +- df(x) = a/(x +- dx)
以找到最佳值a
?
我有一个函数f(x) = a/x
,我有一组包含f(x) +- df(x)
和值的数据x +- dx
。我如何告诉 gnuplot 进行加权拟合a
?
我知道fit
接受这个using
术语,这适用于df(x)
,但它不适用于dx
。似乎 gnuplot 将我所遇到的错误视为我的函数 ( )x
的整个 RHS 的错误。a/x +- dx
如何进行加权拟合f(x) +- df(x) = a/(x +- dx)
以找到最佳值a
?
从 5.0 版开始,gnuplot 明确规定要考虑自变量的不确定性
fit f(x) datafile using 1:2:3:4 xyerror
使用“Orear 的有效方差法”。
(上面的命令需要格式为 的数据x y dx dy
。)
你正在拟合一个方程,如:
z = a/(x +- dx)
这可以等效地写为:
z = a/x +- dz
一个合适的dz。
我认为(如果我的微积分和统计数据正确),您可以通过以下方式从 x 和 dx 计算 dz:
dz = partial_z/partial_x*dx
假设 dx 很小。
对于这种情况,这会产生:
dz = -a/x**2*dx
所以现在您有一个包含 2 个变量 ( z = a/x - (a/x**2)*dx
) 的函数,您希望它适合 1 个常数。当然,我可能是错的......我已经有一段时间没有玩过这些东西了。
这里一个简单的例子就足以证明 gnuplot 正在做你想做的事:
使用以下玩具模型数据构造一个平面文本文件 data.dat:
#f df x dx
1 0.1 1 0.1
2 0.1 2 0.1
3 0.1 3 0.1
4 0.1 4 0.1
10 1.0 5 0.1
仅仅通过观察数据,我们预计一个好的模型将是直接比例 f = x,在 x = 5,f = 10 处有一个明显的异常值。我们可以告诉 gnuplot 以两种非常不同的方式拟合这些数据。以下脚本 weightedFit.gp 演示了差异:
# This file is called weightedFit.gp
#
# Gnuplot script file for demonstrating the difference between a
# weighted least-squares fit and an unweighted fit, using mock data in "data.dat"
#
# columns in the .dat are
# 1:f, 2:d_f, 3: x, 4: d_x
# x is the independent variable and f is the dependent variable
# you have to play with the terminal settings based on your system
# set term wxt
#set term xterm
set autoscale # scale axes automatically
unset log # remove any log-scaling
unset label # remove any previous labels
set xtic auto # set xtics automatically
set ytic auto # set ytics automatically
set key top left
# change plot labels!
set title "Weighted and Un-weighted fits"
set xlabel "x"
set ylabel "f(x)"
#set key 0.01,100
# start with these commented for auto-ranges, then zoom where you want!
set xr [-0.5:5.5]
#set yr [-50:550]
#this allows you to access ASE values of var using var_err
set fit errorvariables
## fit syntax is x:y:Delta_y column numbers from data.dat
#Fit data as linear, allowing intercept to float
f(x)=m*x+b
fW(x)=mW*x+bW
# Here's the important difference. First fit with no uncertainty weights:
fit f(x) 'data.dat' using 3:1 via m, b
chi = sprintf("chiSq = %.3f", FIT_WSSR/FIT_NDF)
t = sprintf("f = %.5f x + %.5f", m, b)
errors = sprintf("Delta_m = %.5f, Delta_b = %.5f", m_err, b_err)
# Now, weighted fit by properly accounting for uncertainty on each data point:
fit fW(x) 'data.dat' using 3:1:2 via mW, bW
chiW = sprintf("chiSqW = %.3f", FIT_WSSR/FIT_NDF)
tW = sprintf("fW = %.5f x + %.5f", mW, bW)
errorsW = sprintf("Delta_mW = %.5f, Delta_bW = %.5f", mW_err, bW_err)
# Pretty up the plot
set label 1 errors at 0,8
set label 2 chi at 0,7
set label 3 errorsW at 0,5
set label 4 chiW at 0,4
# Save fit results to disk
save var 'fit_params'
## plot using x:y:Delta_x:Delta_y column numbers from data.dat
plot "data.dat" using 3:1:4:2 with xyerrorbars title 'Measured f vs. x', \
f(x) title t, \
fW(x) title tW
set term jpeg
set output 'weightedFit.jpg'
replot
set term wxt
生成的绘图 weightedFit.jpg 讲述了这个故事:绿色拟合没有考虑数据点的不确定性,并且对于数据来说是一个糟糕的模型。蓝色拟合说明了异常值中的较大不确定性,并正确拟合了比例模型,获得斜率 1.02 +/- 0.13 和截距 -0.05 +/- 0.35。
因为我今天刚加入,所以我缺乏发布图片所需的“10 名声望”,所以你只需要自己运行脚本来查看是否合适。在工作目录中拥有脚本和数据文件后,请执行以下操作:
gnuplot> 加载“weightedFit.gp”
您的 fit.log 应如下所示:
*******************************************************************************
Thu Aug 20 14:09:57 2015
FIT: data read from 'data.dat' using 3:1
format = x:z
x range restricted to [-0.500000 : 5.50000]
#datapoints = 5
residuals are weighted equally (unit weight)
function used for fitting: f(x)
f(x)=m*x+b
fitted parameters initialized with current variable values
iter chisq delta/lim lambda m b
0 1.0000000000e+01 0.00e+00 4.90e+00 2.000000e+00 -2.000000e+00
1 1.0000000000e+01 0.00e+00 4.90e+02 2.000000e+00 -2.000000e+00
After 1 iterations the fit converged.
final sum of squares of residuals : 10
rel. change during last iteration : 0
degrees of freedom (FIT_NDF) : 3
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 1.82574
variance of residuals (reduced chisquare) = WSSR/ndf : 3.33333
Final set of parameters Asymptotic Standard Error
======================= ==========================
m = 2 +/- 0.5774 (28.87%)
b = -2 +/- 1.915 (95.74%)
correlation matrix of the fit parameters:
m b
m 1.000
b -0.905 1.000
*******************************************************************************
Thu Aug 20 14:09:57 2015
FIT: data read from 'data.dat' using 3:1:2
format = x:z:s
x range restricted to [-0.500000 : 5.50000]
#datapoints = 5
function used for fitting: fW(x)
fW(x)=mW*x+bW
fitted parameters initialized with current variable values
iter chisq delta/lim lambda mW bW
0 2.4630541872e+01 0.00e+00 1.78e+01 1.024631e+00 -4.926108e-02
1 2.4630541872e+01 0.00e+00 1.78e+02 1.024631e+00 -4.926108e-02
After 1 iterations the fit converged.
final sum of squares of residuals : 24.6305
rel. change during last iteration : 0
degrees of freedom (FIT_NDF) : 3
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 2.86534
variance of residuals (reduced chisquare) = WSSR/ndf : 8.21018
p-value of the Chisq distribution (FIT_P) : 1.84454e-005
Final set of parameters Asymptotic Standard Error
======================= ==========================
mW = 1.02463 +/- 0.1274 (12.43%)
bW = -0.0492611 +/- 0.3498 (710%)
correlation matrix of the fit parameters:
mW bW
mW 1.000
bW -0.912 1.000
有关文档,请参阅http://gnuplot.info/。干杯!