The attached patch adds the options yerrorbars and xyerrorbars to the set fit command.
If the xyerrorbars is provided the fit command expects a using spec of the type
x:y:sx:sy to fit y=f(x), where sy is the standard deviation of the y-values and sx
the standard deviation of the independend variable x. For the fit the
"effective variance method" (Jay Orear, Am. J. Phys., Vol. 50, No. 10, October 1982)
is then used to incorporate the standard deviation of the x-values in the estimation
of the weights s used to calculate the Chi^2. By default (option yerrorbars) the
weights are simply s = sy (as befor this patch), but with the xyerrorbars option
the weights are calculated using the formula s**2 = (f'(x)*sx)**2 + sy**2, where
f'(x) is the derivative of the function f at x. The results of this new method were
sucessfully compared to results using ROOT, which uses a similar approach.
Bastian Märkisch
2014-02-17
Thanks for the quick response. Below are my findings for Orear's data set, which more or less agree with your analysis. I guess the conclusion is that we just get stuck in a local minimum and the code can thus go to CVS right now. We really should do a chi² plot of the region...
Summary of the fit results:
a1 a2 a1_err a2_err ------------------------------------------------ initial values 1.0000e-003 5.900e+005 dy only 9.8993e-004 5.893e+005 3.8731e-005 2.289e+004 (*)dy and dx 1.0162e-003 5.937e+005 1.9933e-004 1.116e+005 dy and dx, errorscaling 1.0162e-003 5.937e+005 1.7022e-004 9.527e+004 dy and dx, no prescale 1.0096e-003 5.900e+005 1.9840e-004 1.110e+005 ------------------------------------------------ published 1.0163e-003 5.937e+005 exact 1.0731e-003 6.250e+005 ------------------------------------------------ Minuit 1.0727e-003 6.250e+005 2.5254e-004 1.399e+005
Note that our fit result (*) and the result published in J. Orear, Am. J. Phys. 52, 278 (1984) agree.
This is also an example where--due to the large difference in the order of magnitude of the parameters--prescaling is helpful to obtain convergence!
Ethan Merritt
2014-02-17
Is there a version of this patch that applies against current CVS?
Bastian Märkisch
2014-02-18
This patch should apply to current CVS. It includes a demo fitxyerr.dem
that uses Orear's data and Pearson's data with York's weights. The documentation still needs to be reworked.
Ethan Merritt
2014-02-19
Thanks. That works nicely other some trivial compiler warnings.
Help me to understand the output of the demo.
Demo 1)
How can it be that both the RMS and CHI2 values are better for the xyerr fit, but the p value is worse?
Demo 2)
"This is also an example where--due to the large difference in
the order of magnitude of the parameters--prescaling is required
to obtain convergence." I have lost track - is patch #507 still relevant or does it implement the same thing as the current "set fit prescale"?
I gather that the fit is not exact even with prescaling, but in the demo the residual error is on the same order of magnitude as the convergence limit. However, if you tighten the convergence limit to "set fit limit 1e-09" then it reaches the same end point and claims an epsilon of e-10. Doesn't this mean there is something wrong with the function evaluation rather than with the convergence?
Bastian Märkisch
2014-02-19
Demo 1)
How can it be that both the RMS and CHI2 values are better for the xyerr fit, but the p value is worse?
For the y-errors fit I get p = 3.5e-5, whereas for the xy-error fit it is 0.15, which is much better indeed, i.e. closer to 0.5, and in the 'acceptable' range of lets say [0.06:0.94].
Demo 2)
"This is also an example where--due to the large difference in
the order of magnitude of the parameters--prescaling is required
to obtain convergence." I have lost track - is patch #507 still relevant or does it implement the same thing as the current "set fit prescale"?
Patch #507 is included in CVS since 2013-03-12. Question is, if it shouldn't be included in 4.6.5, too.
I gather that the fit is not exact even with prescaling, but in the demo the residual error is on the same order of magnitude as the convergence limit. However, if you tighten the convergence limit to "set fit limit 1e-09" then it reaches the same end point and claims an epsilon of e-10. Doesn't this mean there is something wrong with the function evaluation rather than with the convergence?
The fact that it does not reach the exact solution is why I started the discussion again. Without further investigation (e.g. a real Chi² plot of the region), I can only assume that the algorithm gets stuck in a local minimum. As Alexander pointed out, the calculation of the WSSR seems to be correct. Also, the endpoint is the same as in Orear's paper, which probably tells us that the starting point is just not good enough for both algorithms.
Although the problem looks easy at first glance, it is really very difficult: the correlation between the parameters is very high ( gnuplot says >0.994 - as does a colleagues' Mathematica script ), and -without prescaling- the region of convergence must be very elongated due to the very different scale of the parameters.
Alexander Täschner
2014-02-19
The fact that it does not reach the exact solution is why I started the discussion >again. Without further investigation (e.g. a real Chi² plot of the region), I can >only assume that the algorithm gets stuck in a local minimum.
Attached is a Chi² plot of the region. The contour lines equal the minimum Chi²
+1, +0.1 and +0.01. The points show the exact solution and Orears one. From the
plot it is clear, that finding the absolut minimum is not really easy.
Bastian Märkisch
2014-02-22
Bastian Märkisch
2014-02-28
Alexander, thanks for the nice Chi² plot. It would be great to be able to do plots like this after a fit with using a chisq()
function in gnuplot. Any ideas how to implement this? Of course we would have to limit it to one or two parameters.
From the plot it seems that it is really the fit algorithm that gets stuck and that it is not a problem of the effective variance calculation per se. An modified version of the patch is now included in CVS. The help text for gnuplot.doc still needs work and is not included yet.
Bastian Märkisch
2014-02-28
Bastian Märkisch
2014-03-04