|
|
|
|
SATAN Help
HOW TO USE THE FIT PACKAGE SATAN provides an elaborate FIT Package for data of ANALYZERS and PSEUDO ANALYZERS. The following kinds of data can be used as input for the SATAN fit package: ANALYZERS may contain: Digital data on an
equidistant grid. The x values may start with any value (integer, real, positive
or negative), and the step size may be any positive real number. Analog data with
constant bin size. The x range may start with any value (integer, real, positive
or negative) and the bin size may be any positive real number. PSEUDOANALYZERS may contain: Digital data with
constant or variable step size. The x values can have any value (integer, real,
positive or negative). Analog data with
constant or variable bin size. The x range may start with any value (integer, real, positive
or negative) and the bin size may be any positive real number. The value of the
bin size can differ from bin to bin. The width of the bin is defined by the
range from the lower limit of the current bin to the lower limit of the next
bin. The fit handling with data of pseudo analyzers: 1. The data of pseudo analyzers are defined in the dataset by: ==> A: X Y,LN1 D The fit handling with data of analyzers: The fit to data in logarithmic representation can be performed: Errors can be defined Symmetrical vertical error bars defined in the dataset are taken into account.
Asymmetrical vertical error bars defined in the dataset are averaged and symmetrized
before entering into the fit. Two-dimensional data (array or analyzer) are also transferred to the SATAN fit
package when the command GREAD / FIT2D was entered and
when they are shown by GDISP. A function y=f(x) (e.g. a
polynomial) may be fitted to channels with z>0. The statistical weight due to the
number of counts is considered. If a polygon boundary is defined, only the data inside
this boundary are taken into account. (See also the command ARIDGE which offers more extended possibilities for
data of two-dimensional analyzers!) For fitting user-defined functions, a specific user-defined subroutine has to be
linked to GRAF (see below). The SATAN Fit Package This section provides a survey of the data fitting facilities. Following a general
description, the program executing the maximum-likelihood fit is shortly explained
("The Maximizer MLFIT")), some aspects of fitting and their realization are
described ("The Fitting Procedure") and, finally, a few applications are shown
("Examples"). In SATAN, a fit package is offered which permits a great variety and flexibility in data fitting. The programs are based on the maximizer 'MLFIT' [ V. Blobel, MLFIT: A Program to Find Maxima of Likelihood Functions, DESY Report 71/18, April 1971]. For a practically unlimited number of fit parameters this subroutine searches for local maxima of the logarithm of a likelihood function in parameter space by combining two methods, namely, the quadratic approximation developed from a Taylor's expansion and a gradient search algorithm. It allows for the fixing of parameters, giving limits that must not be exceeded by the parameters and performing a special error analysis derived from the true variation of the likelihood function near the maximum. In addition, any fit parameter may be correlated to another reference parameter by fixing the corresponding difference, sum, quotient or product. Two fit methods are offered which correspond to different definitions of the likelihood function: o The method of least squares which uses a likelihood function given by the
weighted sum of squares ('chi square'); The respective likelihood function may be plotted versus one or two fit parameters with a linear or logarithmic axis. The following fit functions are directly supported by commands: o a polynomial and the square root of a polynomial of any order, both optionally
multiplied by the exponential function of another polynomial, In addition to this set of standard functions the user may define his own fit function in a resident procedure, where the function value and first derivatives with respect to the fit parameters are calculated. The total fit function may be specified as the sum of up to ten standard or user-defined functions added initially or step by step after fitting successively parts of the spectrum. This provides a general handling of background fitting since any function may act as background for a superior one. Some fit utilities include the following. o Definition of spectrum windows to be exclusively considered by the fit The following chapter summarizes the principle of operation of the maximizer MLFIT. For detailed information refer to the original publication (flike,flsqu,fpois). Estimation of a parameter set "a = (a1 ... an)" of a theoretical function "f(x,a)" describing the experimental data "yi(xi)" is done by constructing and maximizing a likelihood function "L(a)", which represents the relative likelihood for the fit parameters to assume particular values. Two different definitions of likelihood functions are offered which correspond to fitting modes known as the least squares fit and the Poisson distribution fit. The logarithms "F(a) = ln{L(a)}" are given by the expressions F(a) = -1/2 SUM{ wi [f(xi,a) - yi]**2 } LEAST SQUARES
where the weighting factors "wi" are determined by the experimental wi = ( Delta yi )**(-2), and F(a) = SUM{ yi ln f(xi,a) } - SUM{ f(xi,a) } POISSON DISTRIBUTION
For a given parameter set the logarithm of the likelihood function is expanded into a Taylor series by considering the gradient "g" and the matrix of the second derivatives "G" neglecting higher-order terms. F(a + Delta a) = F(a) + [g]t Delta a - 1/2·[Delta a]t·G·Delta a An iteration step Delta a = K**(-1) g is calculated using a matrix K = G + eps1(2**l-1)·([g]t·G·g / [g]t·g)·I (I = unit matrix). The quantities "l" and "eps1" determine the step direction in parameter space; "l=0" corresponds to quadratic approximation, for "eps1(2**l-1) approximately 1" the step is carried out in the direction of the gradient. Usually starting with the value "l0=0" the parameter "l" is modified according to the result of the last step, which is classified by comparing the true variation of the likelihood function "Delta Ft", the estimated variation "Delta Fe" in case of quadratic behavior, and a third value Delta Fc = MAX{ eps2, eps3 |F(a + Delta a)| }
into three categories: o Divergent: Final convergence is assumed on the following conditions: o "Delta Fe" is less than "Delta Fc", A step carried out into a forbidden region is simply treated as divergent; a
condition variable is set to one of the following reason codes: Generally the number of terms processed in an iteration equals the number of data points; a smaller value indicates the running number of the last fit point before occurrence of an error condition. Symmetric errors "sigma" and correlation coefficients "rho" for the fit parameters are derived from the inverted matrix of second derivatives of the likelihood function (error matrix). They may be used for considering error propagation; the error of a quantity "f" depending on two variables "x" and "y" is defined according to the formula sigmaf**2 = sigmax**2 · (df/dx)**2 + where the expressions sigmax**2, sigmay**2 VARIANCES rhoxy, sigmax, sigmay COVARIANCES are represented by the diagonal and non-diagonal elements of the error matrix, respectively. More precise and generally asymmetric errors are obtained by variation of the likelihood function in the near of the maximum. This error analysis is time consuming, because additional iterations are performed. Fixing of a fit parameter is indicated by setting the corresponding partial
derivative to zero. If two parameters are strongly correlated, some matrix elements may
become very small during inversion. To avoid absurd results due to rounding errors, the
matrix is only partly inverted if certain matrix elements decrease by more than a factor
of "tau". The parameter corresponding to this matrix element is not modified
during the current iteration, i.e. the parameter is treated as fixed. The MLFIT parameters which may be modified interactively if required at any
critical fitting condition are listed below together with their default values: eps1 0.1 gradient search parameter A frequent change of the gradient search parameter "l" between zero and
larger values during the iterations indicates a strong non-quadratic behavior of the
likelihood function; in such a case it may be useful to start with "l>0", or
to increase "ld" and decrease "eps1" simultaneously. If the accuracy
of the likelihood function is limited by large rounding errors, it is recommended to
increase eps2 and eps3; the latter parameter is essential and may be increased by orders
of magnitude. In the following, some aspects of data fitting and their realization are described. Initialization The command FINIT resets all fit variables and modes to defined starting values. Fit Windows Generally, all data shown after a display command are considered by a fit. By FWIN up to 30 parts of the spectrum may be specified either by cursor or terminal prompting loop, or via display windows or analyzer conditions defining the actual fit regions. Fit windows may extend beyond the displayed spectrum segment; irrespectively of this fact only data within the display region are fitted. Display window names and analyzer conditions are merely used as synonyms for channel limits; a fit window defined by this means is not affected by any display command. Note that unless explicitly demanded within the commands FIT or SET, bins with zero contents are not included by the fit. Temporary windows (valid only for a single command) may be given within the commands FIT (execution of fit), FLIST (evaluation of the fit function) and FSUM (sum and moments of raw and fit data). Calibration The x-values assigned to each bin may be redefined by modifying the parameter FXCAL of an analyzer with the command AMODIFY or by applying an arithmetic operation to the x values of a pseudo analyzer. Besides the major application to get an implicit transformation of fit results, The calibration may be used to avoid an overflow in numerical computations by reducing the x-values by orders of magnitude or to define the bin center according to the mean value of the contained events. For analog data, a channel "xi" is expected to include events with an observed value lying in the interval between "xi" and "xi+1"; accordingly linear calibration is effective with coefficients set to "a = 1" and "b = 0.5" (middle of the bin). For histograms (analog data), the x values are derived from the mean of the channel numbers contained by a bin. Fit Functions The various standard fit functions are defined by FPOL (polynomial or square root of a polynomial, optionally multiplied by an exponential function of a polynomial), FLEG (sum of legendre polynomials, optionally consisting of even or odd terms only and normalized), FEXP (sum of linear exponentials), FPEAKS (Gaussian or Lorentzian peaks), FBOX (convolution of rectangular box with a peak function), FPTAIL (convolution of an exponential with a peak function), and FPVOIGT (convolution of a Lorentzian peak with a Gaussian function). The command FMY acts as an interface to a user-defined fit function, which can be loaded dynamically. Background Any desired function can be considered as a background function since the command FLAST allows the total fit function to be composed of up to ten standard or user-defined components. By FLAST/KEEP the last-defined fit function is stacked (no matter if the fit was executed or not), i.e. the function specified next is taken as an additional term; successive fits adding components step by step may be executed on different windows. Each stacked function gets an internal serial number by which it is identified for deleting or displaying singly. Automatic renumbering is done after a component i is deleted by FLAST/DEL(i). Keeping the last fit as background, the parameters of the respective function(s) are not fixed automatically. Peaks Gaussian or Lorentzian peaks are specified by the command FPEAKS providing three ways: o limits and FWHM are found by the automatic peak search routine, Whenever peaks are defined by limits, initial values for the peak parameters area, center and full width at half maximum are derived from sum, mean and variance of the net data in the window referred to. If only centers are given as prompting input, the full width at half maximum may be specified, however, a default value is computed from the data. The product of width and peak height yields the initial area. The resulting function is displayed (unless the keyword NODISP is given). Lower limits of areas and widths are set to zero; limits for the peak centers are set equal to the window limits or at least half the corresponding width above or below the center, respectively. Modified peak functions, resulting from the convolution of a Gaussian with a rectangular and an exponential function, are specified by the commands FPBOX and FPTAIL, respectively. Parameter Values FPAR is the general command to assign initial values and permitted limits to fit parameters. Each parameter is addressed by a running number (corresponding to the description of the respective fit function). Alternatively, the commands FAREA, FPOS and FWIDTH are available for peak functions, where areas, centers and full widths at half maximum are specified by peak numbers. After definition of a function, fit parameters are set to arbitrary starting values (generally zero). However, parameter values of a user-defined function may be initialized within the procedure $MYFIT which is called by FMY, and guess values for parameters of peak functions (FPEAKS) are calculated from sum and moments of corresponding spectrum regions. In case of a structured likelihood function with several local maxima it might be useful to plot this function versus one or two fit parameters to find appropriate starting values for the fit. This is supported by the command FLIKE which evaluates the likelihood function of the current fit environment and stores the result in an analyzer; a reasonable step size of parameter variation avoids an excessive consumption of computing time. At any time, the fit function may be evaluated and displayed with the actual set of parameter values. Fixing Fit Parameters A fit parameter is allowed to be held constant during iterations by fixing its absolute value or the difference, sum, quotient or product with respect to any reference parameter (which for its part must not be relatively fixed). The guess value for a relatively fixed parameter may be given in relative mode, i.e. by a difference, sum, quotient or product which internally is transposed into an absolute value. Fixing and varying are done by FPAR or in case of peak parameters by FAREA, FPOS and FWIDTH. Execution of Fit The fit is executed by the command FIT using the mode of least squares or Poisson distribution, respectively. In case of least-squares fit weighting factors for the data are represented by the reciprocal squares of the experimental errors specified by FERR; if no errors are given, all data have constant weights. The fit may be very sensitive to the choice of appropriate weighting factors. Bins with zero contents are excluded from fit by default but may be demanded to be considered. Least squares weights if required are set to 1 unless a non zero error is given; however, for spectra with small data rates the Poisson fit mode is adequate. The execution of fit is stopped by reaching convergence, by a maximum number of iterations or by exceeding a CPU-time limit; it may of course be restarted with the actual fit parameter set. The fit may also be interrupted by pressing the "BREAK" button on the SATAN interface window. In case of critical fit conditions some MLFIT parameters may be altered by means of the command FMLFIT. A detailed output after each iteration comprising true and estimated values of the likelihood function, number of terms processed, "l"-parameter, condition flag in case of divergence and current fit-parameter values is given on the keyword LIST of the FIT command. Results Values, errors and correlations of the fit parameters are listed, prepared for GSAVE or stored into global parameters by FRESULT. The command FLIST produces an output of fit-function values and experimental data. Under some circumstances it may be preferable to determine sum and moments of raw data or raw data minus background fit in an interval rather than peak fitting. FSUM calculates the sum and first to third moments of a distribution "yi(xi)" due to the following expressions: s = SUM{ yi } SUM
Graphs of the total fit function, selected components or single peaks are shown by
FDISP with an optional accuracy defined by a number of
polygon lines or cubic splines composing the line to be drawn. With FSTORE, fit data evaluated for each bin of the displayed
region (or the difference of experimental and fit data, optionally divided by Errors of Experimental and Fit Data Experimental errors (and hence weighting factors for least squares fit) may be specified by FERR to be undefined (corresponding to constant weighting factors), statistical (square root of observed counting rate), proportional (given by a percentual amount of the counting rate) or external (contained in a dataset). Initially, errors are undefined. In this case, no errors can be evaluated for the fit parameters. Symmetric errors and correlations of fit parameters are taken from the matrix of second derivatives; in addition, an extended error analysis generally yielding asymmetric errors may be demanded for particular parameters (flagged by FPAR, FAREA, FPOS, FWIDTH, respectively) by the keyword ERRANAL within the command FIT. Errors for fit data (FLIST, FSTORE) or net integrals of raw data (FSUM/FITERR) are derived from symmetric errors and correlation coefficients (which are essential for considering error propagation) and listed, provided that experimental errors have been defined. Problems In the following, a few critical fit conditions are described. In any case, the fit should be restarted using the keyword LIST of the command FIT to get a comprehensive output for each iteration which facilitates the solution of the problem. o An error condition occurs during iteration. The column 'TERMS' of the iteration output indicates the fit point processed last before occurence of the respective error denoted in the column 'COND' by a number: 2 fit parameter exceeds specified limit Depending on the error condition, different actions are recommended: - exclude the respective fit point by FWIN; o The fit converges very slowly ('further iteration recommended'). - Looping of the "l" parameter (column 'L' of the iteration
output) between 0 and larger values indicates a strong non-quadratic behavior of the
likelihood function. In case of error condition_2 check the fit parameter limits,
especially for peak positions (FPAR/LIST, FPOS/LIST); open parameter intervals which are set too
narrow (e.g. FPAR.../NOLIM). Otherwise retry the fit
with an increased MLFIT parameter "l0" (default 0), or increased "ld"
(default 2) and decreased "eps1" (default 0.1), set by FMLFIT. o The fit has reached convergence, but the result looks unreasonable. Fitting a small number of fit points (compared to the number of free fit
parameters), frequently parameters are fixed internally due to strong correlations (and
thus have unreasonable values), which is indicated by the number at column 'FIX' of the
iteration output or by missing fit-parameter errors. Decrease the parameter
"tau" (default 1E-5) by an order of magnitude with FMLFIT/TAU(...)
and restart fitting. Fit of Legendre polynomials to an angular distribution GDISP ... display command Fit of a single Gaussian to displayed data GDISP ... display command Fit of a spectrum of Gaussian peaks plus background GDISP ... display command Fit of two correlated peaks with relatively fixed parameters GDISP ... display command Fit 8 peaks, all widths are required to be equal GDISP ... Fit with the user-defined function GDISP ... display command Example for user-defined fit functions $MYFIT: PROC(C_NAME,I_NPAR,D_X,D_FPAR,I_SWIT,D_Y,D_DERIV,I_DERIV,I_RC); Link the GRAF package together with your $MYFIT function by
Independently of the above described fit
package, a
simple polynomial fit is available by the line symbol "F" in a GRAF
dataset. This fit is always performed to the representation of the data, like the fit
package does, if the option |
