3

我有一个函数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

4

3 回答 3

3

从 5.0 版开始,gnuplot 明确规定要考虑自变量的不确定性

fit f(x) datafile using 1:2:3:4 xyerror

使用“Orear 的有效方差法”。

(上面的命令需要格式为 的数据x y dx dy。)

于 2015-08-21T09:26:20.703 回答
1

你正在拟合一个方程,如:

 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 个常数。当然,我可能是错的......我已经有一段时间没有玩过这些东西了。

于 2012-10-16T17:19:33.813 回答
0

这里一个简单的例子就足以证明 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/。干杯!

于 2015-08-20T18:21:27.757 回答