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