Program: FIT
Purpose: Fit a set of data points (Xi,Yi) to a model, that depends
nonlinearly on its parameters, using a Chi-squared
minimization.
Category: FITTING, UTILITY
File: fit.c
Author: M.G.R. Vogelaar
Keywords:
X= Enter X values:
Examples:
X=file(fitdata.dat,1,1:) , read first column of ASCII file
X=1:3:0.05 , generate x=1 to 3, step size 0.05
Y= Enter Y values:
Examples:
Y=image(spectr1, -127 0 128 0)
If the Y-data contains blank values then the
corresponding weights will be set to 0.
GRDEVICE= Plot device: [List of devices]
Destination of plot, Screen or Hardcopy.
Cursor interaction is possible in X11 and GIDS
============= INSIDE THE FITTING LOOP: ===========
WEIGHTS= Enter weights (W_i=1/VAR_i): [1.0 for all]
The weights are the inverse of the variances in Y_i.
The maximum number of weights is the same as the number
of data points and is initialized to 1.0 (default). If
there are less weights specified than the number of
data points, then the missing weights are set to 1.0.
Pressing the left mouse button while in cursor mode
(CURSOR=YES) toggles between the original weight and
weight 0.0 for a marked data point. Weights that are
set to 0 with the mouse can only get a new value
if they are unmarked. Data X values that can not be
evaluated by the expression evaluation are excluded
by setting the weight to 0.0.
FX= Enter expression:
Describe your model as a function of X.
See description for allowed functions and constants.
At this point keywords are prompted that
represent a model parameter (best fit parameter) in
your expression. The maximum number of parameters is 31.
The length of a parameter may not exceed 16 characters.
The expression is case insensitive.
E.g. FX=A*EXP(-(X-B)*(X-B)/(2.0*SIG*SIG)) generates
A=, B= and SIG= as best-fit parameters and keywords.
This function is also an example of a model that depends
NONlinearly on its parameters.
WARNING: Do not create variable names equal to one
of the listed keywords because otherwise you can get
wrong initial estimates.
....=
....=
....=
CURSOR= Exclude data with cursor? Y/[N]
After the data is displayed in a window, it is
possible to mark data using the cursor. A marked data
point gets a green colour and weight 0. If you mark it
again, it will become yellow again and gets
its original weight (as in WEIGHTS=).
QUIT= Want to quit? Y/[N]
End loop over keywords and fits and end program.
FILENAME= Name of output ASCII file: [No output to file]
If a name is specified, an ASCII file is created
where fit parameters are listed in a row. If you
press carriage return, there will be no output to
an ASCII file.
APPEND= File exists, ok to append? [Y]/N
The file specified in FILENAME= already exists. You
can append to this file with APPEND=Y. If APPEND=N
you will be prompted for another name.
==== Next keywords are asked in a loop, to be ====
==== aborted with INSET= ====
=========== HIDDEN KEYWORDS: ===========
Plotting:
---------
** PGMOSAIC= View surface sub divisions in x,y: [1,1]
View surface can contain a number of plots in
in X and Y direction (mosaic). Default is 1 plot in
both X- and Y direction.
** PGBOX= Corners of box Xl,Yl,Xh,Yh: [default by application]
It is possible to overrule the calculated
PGPLOT box size with PGBOX= The coordinates (x,y) of
the lower point are given first.
** PGWIDTH= Give line width 1..21: [2]
** PGHEIGHT= Give character height: [1.0]
** PGFONT= Give font 1..4: [2]
1 single stroke 'normal' font
2 roman font
3 italic font
4 script font
** YERRRANGE= Enter Y range for residues: [calculated Ylo, Yhi]
The difference between model and original data is
displayed in a second window. The range in Y can be
adjusted with this keyword. If you change the default
values after a fit, use CLEAR=YES to erase the plot
first.
** CLEAR= Clear plot page? Y/[N]
Fitting:
--------
** TOLERANCE= Fractional tolerance for Chi-squared: [0.001]
Fitting is stopped when successive iterations fail
to produce a decrement in reduced Chi-squared less
than TOLERANCE= The value cannot be less than a
certain minimum set by the system. This means that
maximum accuracy can be obtained with TOLERANCE=0.0
** LAB= Mixing parameter: [0.001]
The Mixing parameter, LAB= determines the initial
weight of steepest descent method relative to the
Taylor method. LAB= should be a small value (i.e. 0.001).
LAB can only be zero when the partial derivatives are
independent of the parameters. In fact in this case LAB=
should be exactly equal to zero.
** MAXITS= Maximum number of iterations: [100]
Stop the fitting process if it needs more than
this number of iteration.
Description: Program FIT fits model parameters to data that you
enter with keywords X= and Y= The 'best' fit
parameters are those that minimizes Chi-squared.
The model is entered as an expression in which you
define the variable names yourself. The model may
depend non-linearly on its parameters. If for
instance you entered FX=A*EXP(-(X-B)*(X-B)/(2.0*SIG*SIG))
then you are prompted to enter initial estimates for
A=, B= and SIG=
The following functions are recognized:
SQRT, SIN, ASIN, COS, ACOS, TAN, ATAN, EXP, LN, LOG,
SINH, COSH, TANH, RAD, DEG, SINC.
The following constants are recognized:
PI 3.14159....
C speed of light (SI)
H Planck (SI)
K Boltzmann (SI)
G gravitation (SI)
S Stefan-Boltzman (SI)
M mass of sun (SI)
P parsec (SI)
BLANK Universal undefined value
To prevent that these constants conflict with any variables
you must enclose the constants by braces, e.g. FX={C}*X
The maximum number of parameters is 31 and the length of
a parameter name may not exceed 16 characters.
The (initial) function is plotted (using your initial
estimates) so that you have some control over how close the
initial estimates are to the data. In most cases the fit will
succeed, even with very rough initial estimates, but
sometimes the fit is more sensitive and you have to think
a bit harder about the initial values.
The method used for non-linear models is the Marquardt
method (Marquardt D.W. 1963, Journal of the Society for
Industrial and applied Mathematics, vol 11 pp. 431-441).
The method is a maximum neighbourhood method which
performs an optimum interpolation between the Taylor
series fitting method and the gradient method.
The Mixing parameter, LAB= determines the initial
weight of the gradient method relative to the Taylor method.
The program modifies this value so only a start value
is required. A value LAB=0.001 works satisfactorily.
The fitting process needs a condition for stopping.
When successive iterations fail to produce a decrement
in reduced(!) Chi-squared less than TOLERANCE=
the fitting process is stopped. Iterating to machine
accuracy is possible if you set TOLERANCE=0, but this
will not be statistically very meaningful. Therefore
the default value for the tolerance is 0.001.
Usually the Marquardt method needs symbolic expressions
for the partial derivatives. This program is able to
create these expressions after the model is entered.
If you want to check the symbolic partial derivatives
generated by the program, set Hermes in 'test' mode and
watch the log file.
If a fit is ended successfully, the fitted parameters are
displayed in the log file together with their uncertainties
If N is the number of independent data points and n is the
number of fitted parameters, then a reduced Chi-squared_r
can be calculated from the Chi-squared:
Chi-squared_r = Chi-squared/(N-n-1)
The least squares fit maximizes the probability
distribution function for Chi-squared by minimizing
the function:
( y_i-y0[x_i] )^2
Chi-squared = SIGMA {-----------------}
sd_i^2
where sd_i is the standard deviation of data point i.
Note that 1/sd_i^2 is the definition of the weights
in WEIGHTS=
The maximum of the probability distribution function is
also calculated. It describes the probability that a
RANDOM set of N data points would yield a value of
Chi-squared as large or larger when compared with the
parent function (Bevington, data reduction and Error Analysis
for the Physical Sciences).
Notes: .......
Example:
fit
(start fit) 06/02/97 09:59:54
FIT Version 1.0 (Feb 6 1997)
FIT X=file(fitdata.dat,1,1:)
FIT Y=file(fitdata.dat,2,1:)
FIT GRDEVICE=x11
FIT FX=a*exp(-(x-b)*(x-b)/(2.0*sig*sig))
FIT SIG=1
FIT B=1.5
FIT A=3.3
==================================RESULTS=======================
Fitted to the data was FX = A*EXP(-(X-B)*(X-B)/(2.0*SIG*SIG))
The fitted parameters (with errors):
SIG = 1.975 (+/- 0.0960817)
B = 1.56574 (+/- 0.0469892)
A = 3.41704 (+/- 0.0168025)
least squares fit converged to TOL=0 in 5 iterations,
The least squares fit converged to TOL=0 using a mixing parameter LAB=0.001.
You can rerun this with new values for FX=, TOLERANCE=, MAXITS=, LAB=
and new initial estimates.
=================================================================
FIT QUIT=y
(end fit) 06/02/97 10:02:38
FIT +++ FINISHED +++
Updates: Jan 28, 1997: VOG, Document created.
Nov 21, 2006: VOG, Included output to Ascii file