[Gnuplot-info] Outputs of fit: details of how they're defined

[Daniel H. <vi5...@ya...>]

Dear All,

I find myself using Gnuplot (5.0 patchlevel 3) in teaching
least-squares fitting in a Bayesian framework.  As such, I need to
know some of the details of how the various outputs from the fit
command are defined.

>From a mixture of reading the manual and experimenting with running
various commands, I've come to the following tentative view of what's
happening:

Let's say that I have run

set fit noerrorscaling
set fit errorvariables
set fit covariancevariables

fit some_equation(L,M,x) 'some_data_set.dat' using ($1):($3):($2):($4) errors x,z via L,M

where some_equation(L,M,x) is a theoretical formula, with adjustable
parameters L and M, which is supposed to predict the quantity in the
third column of the data file, when the value of the quantity in the
first column of the data file is x; the second column of the data file
contains the standard uncertainties in the measured values in the
first column of the data file 'some_data_set.dat'; and the fourth
column of the data file contains the standard uncertainties in the
measured values in the third column of the data file
'some_data_set.dat'.

Then I _think_ I expect to find (with apologies if my Unicode
equations don't display correctly):

- in the variables L and M, the values l and m of the L and M
parameters respectively, which, among all the (L,M) pairs that the
Marquardt-Levenberg algorithm has tried, produce the lowest value of
the χ⁲ statistic;
- in the variable FIT_WSSR, the value Χ⁲ of the χ⁲ statistic that is
achieved by setting L = l and M = m;
- in the variable L_err, the value (presumably calculated in a
finite-difference way) of
√(2/((∂⁲(χ⁲)/∂L⁲)ₗ¸ₘ-((∂⁲(χ⁲)/(∂L∂M))ₗ¸ₘ)⁲/(∂⁲(χ⁲)/∂M⁲)ₗ¸ₘ));
- in the variable M_err, the value (presumably calculated in a
finite-difference way) of
√(2/((∂⁲(χ⁲)/∂M⁲)ₗ¸ₘ-((∂⁲(χ⁲)/(∂L∂M))ₗ¸ₘ)⁲/(∂⁲(χ⁲)/∂L⁲)ₗ¸ₘ)); and
- in the variable FIT_COV_M_L, the value (presumably calculated in a
finite-difference way) of
2(∂⁲(χ⁲)/(∂L∂M))ₗ¸ₘ/((∂⁲(χ⁲)/∂L⁲)ₗ¸ₘ(∂⁲(χ⁲)/∂M⁲)ₗ¸ₘ-((∂⁲(χ⁲)/(∂L∂M))ₗ¸ₘ)⁲).

Please can anyone confirm whether I've correctly assessed what the fit
routine is putting in those output variables?

As an aside, I believe the following two statements to be true when
this least-squares fitting process is done in a Bayesian framework:

- As long as the prior probability density over parameter space,
evaluated at the point (L,M) = (l,m), is non-zero, the leading-order
Laplace's method approximation to the standard deviation of the
posterior probability distribution for L, marginalized over M, is
√(2/((∂⁲(χ⁲)/∂L⁲)ₗ¸ₘ-((∂⁲(χ⁲)/(∂L∂M))ₗ¸ₘ)⁲/(∂⁲(χ⁲)/∂M⁲)ₗ¸ₘ)).

- As long as the prior probability density over parameter space,
evaluated at the point (L,M) = (l,m), is non-zero, the leading-order
Laplace's method approximation to the standard deviation of the
posterior probability distribution for M, marginalized over L, is
√(2/((∂⁲(χ⁲)/∂M⁲)ₗ¸ₘ-((∂⁲(χ⁲)/(∂L∂M))ₗ¸ₘ)⁲/(∂⁲(χ⁲)/∂L⁲)ₗ¸ₘ)).

I had to (re)invent those two statements myself, by means of an
enormously long-winded derivation: I haven't seen them in any
published source.  However, I can't believe that I'm the first to make
those statements, and I have a sneaking suspicion that the person who
wrote the phrase "asymptotic standard error", both in the Gnuplot
manual and in the comments of file <src/fit.c>, in the Gnuplot source
code tree, _has_ seen those two statements or something very like them
in a published source.  If that person is on this mailing list, it'd
be great if s/he could provide a reference, please.

Thanks very much.

-- 

Kind regards,

Dan