** INFORMATION ON THE LEAST-SQUARES FITTING PROGRAM GENLS ** February 1992 A. BACKGROUND GENLS V.1.1 is an "errors-in-variables" generalized least squares program written by J. Ross Macdonald which allows fitting of functions linear or nonlinear in some or all of their parameters, a(j), (j < 11). It applies when there are errors in one or both of the data variables, x and y. In its current version, it handles up to N = 151 x-data values and applies only when the fitting function (constraint), F(x,y,a) = 0, can be written as either y = f(x,a) or x = g(y,a), where a is the vector of parameters. Any (up to 9 total) of the 10 possible parameters of the model may be taken fixed. Six fitting functions are provided, but others can be added if the function subprogram, GENLSF.FOR, is recompiled. Because one or more of the fitting parameters may be taken fixed, the six available fitting functions actually provide hundreds of fitting possibilities. The program provides the option to change the weightings of the x and y residuals during operation. A great many weighting choices and possibilities are provided. GENLS is operated through a series of interlocking menus and is designed to be transparent and user-friendly. Thus, possibilities for use and questions which might arise are primarily incorporated and answered in the menus. Therefore, to use the program (see below) one needs only to provide some data and start GENLS. GENLS first carries out an approximate fit based on the Deming solution of the EOV problem (1,2); this is used as input to an iterative algorithm for solving the problem accurately (3,4). Because of its use of unconventional partial derivatives, this algorithm converges to an exact solution very rapidly, and it requires the inversion of much smaller matrices than do most other EOV programs. The convergence-time advantage of GENLS thus improves rapidly as N, the number of data points, becomes large. Recompilation of the GENLS source code, which is provided, allows one to increase the number of fitting parameters and/or the number of data points beyond the above default values. B. MAIN GENLS FILES The actual GENLS source code (Fortran 77) is contained in the files GENLSA.FOR, GENLSB.FOR, GENLSP.FOR, AND GENLSF.FOR. The object versions, xx.OBJ, are also included. The executable file is GENLS.EXE. It was compiled with Microsoft Fortran V.5.1, requires the MS-DOS operating system, and thus runs on a PC/AT machine or clone. A numerical co-processor chip is not necessary but is recommended to decrease run time for large data sets and/or many free parameters. The program can be recompiled to operate on other types of computers. See the note above the first "OPEN" in GENLSB.FOR, however, before recompiling. The program requires about 165KB of RAM to run under MS-DOS. The files MKGLS.BAT and MGENLS are provided to help with recompilation, should it be needed. The files GENPYORK, GENREED, GEXP1PC, and GEXP15PC are input data files for testing the program. The last two, for exponential decay data, involve y data with 1% and 15% proportional errors, respectively, and no errors in the x values. The file GENPYORK.TST is the result of running GENPYORK in simplest fashion and is provided for comparison as a test of the program. C. RUNNING AND TESTING GENLS Before running GENLS, have available a file, say MYDAT, with input data. (See above data input files for some possible formats.) The simplest input file consists of a line of identifying information followed by two columns of data, x and y, free format. Three other input possibilities are described in Menu 0, which always appears when GENLS is started. To run GENLS for input file MYDAT, enter the command GENLS MYDAT. You will then be requested to enter the data structure identification type (1-4) in Menu 0. It is a good idea to test GENLS on your machine before starting your own fitting runs. To do so, enter GENLS GENPYORK. Then, starting with Menu 0, enter M0-4, M1-1, and M3-2. Here, M1-1 means make choice 1 when Menu 1 appears. The program will run, showing some output on the screen, and will end at Menu 5. Pick M5-7 (to make the output file GENOUTF), followed by M5-8, to exit from the program. If GENLS is working properly on your machine, GENOUTF and GENPYORK.TST should be virtually identical. In general, if the screen output rolls by too fast, you can either pause the display or make sure to save the results in a file for later examination (make choice M5-7). The final sum-of-squares of the weighted residuals is denoted by S. In the convergence-choice menu, M1A, there is a convergence factor G. Its default value is 1. When it is set to less than 1, the approach to convergence is slower but is likely to be more certain. The output denoted (A/S)*(DS/DA) involves normalized discrete derivatives for each parameter. For a full converged solution, its values should all be less than 10**(-5) and preferably less than 10**(-7). Also present in the output are parameter increments; they should be very much smaller in magnitude than the parameter values when convergence is achieved. D. OUTPUT FILES When Menu 5 is reached, you may exit or make one or several output files. After making files you can, without leaving the program, change the fitting function or weighting and make another fit. Its results can be saved in one or all of the output files. Thus, sequential results can be saved, representing the results of several fits with the same data. The file GENOUTC is always made. It is a summary of the main (menu) commands entered during the fitting(s), and thus provides a history of the run. The other three output files are produced only if selected as specific choices in Menu 5 after a fit. The other output files are GENOUTF, GENOUTP, and GENOUTT. The first saves the important results of a fit; the second (invoked by entering negative values of a choice number - see menu) saves any rough plots of residuals made after fitting (seven choices in general); and the third provides tables of input data and (final) weighting, adjusted variable values, and various residual quantities. Of particular importance are the DI(i), the composite weighted residuals (see Ref. 6 for a discussion of such residuals and Ref. 7 for the determination of their signs when there are errors in both variables). E. WEIGHTING Definitions and discussion of error and weighting models are discussed in Refs. 4, 8, 9, and 10. In GENLS, many weighting possibilities (some of which are unavailable in any other EOS program) are provided in Menu 3. See especially Ref. 4 for more weighting-model details. For example, the Independent Diagonal Weighting Model (IDWM) listed in Menu 3 involves a weighting matrix containing only diagonal elements which are independent input constants, while all the x-weights are the same and all the y- weights are the same in the Standard Weighting Model (SWM). Unless the uncertainties of the x and y data points are rather well established, they often need not be entered with the data but plausible values may be set from within the program. The Menu-3 choices set values for the estimated standard deviations of the x and y data values, SIG(X) and SIG(Y). These values can be all equal for the x data and all equal to another value for the y data (i.e the SWM). See Refs. 4, 8, and 11 for a discussion of the simple PW and GM straight-line fitting/weighting models (M3, M4 choices 1 or 3) incorporated in GENLS. In addition, GENLS includes a rather general important weighting model (M3, choice 5) which involves the possibility of both additive and power-law components; it provides, after convergence, a self-consistent set of weights (generally different for each point) which involve the final values of the x and y: X and Y, their adjusted values. In the usual case where one can express Y in terms of a function of X and a, the fitting model is Y = f(X,a), where at convergence the a-vector involves the final estimated parameter values. The above possibilities allow one to carry out ordinary least squares fits of the data if desired. For fits with x taken exact (=X), one would choose M3-3 and set SIG(X) = 0 and SIG(Y) = 1. For fits with y taken exact, one would just reverse these values. In the physical sciences, where data to be fitted often have a large range, it is frequently found that the uncertainties in the data are each approximately proportional to the absolute value of the datum itself - constant percentage error. To take proper account of such a situation, one needs to use standard deviations to form weights which are proportional to the absolute values of the converged, adjusted X and Y results (see 4, 9, and 10). For example, if it is believed that the x-errors are 1% and the y ones are 10%, one would go to Menu 3 and pick choice 5. Then set UX = UY = 0, PX =0.01, and PY = 0.1, followed by XIX = 1, XIY = 1. Sometimes, one might expect only additive errors in the x data and proportional ones for the y data. In Menu 3 one would then pick choice 5 and set UX = UY = 0, PX = a constant, PY = another constant, and take XIX = 0, XIY = 1. To include additive errors in both the x and y uncertainties, set UX AND UY non-zero. To examine such fitting in practice, make a fit of the GEXP15PC data, where it is known that the x data are exact and the y data involve 15% proportional random errors. One would pick choice 5 at Menu 3. The rest of the commands are listed in the file GENOUTC.EXP which is included with the other files. An important necessary, but not sufficient, test of proper fitting and weighting is to see how well SF, the overall standard deviation of the fit, approximates unity. The square of this quantity is the statistician's chi-square. For the GEXP15PC run discussed above, a value very close to 1.00 is found when the proper weighting is used. Note that when a fit is particularly poor, perhaps because of failure of convergence within the allowed number of iterations, one will often find that one or more of the estimated correlations between parameters is greater than 1! Then all the results must be rejected. F. GENERAL COMMENTS AND SUGGESTIONS No initial parameter choices need be made for linear-model fitting: i.e. a model linear in all parameters. GENLS is particularly easy to use for straight-line fitting. If no initial weighting changes are needed, after setting the data type in Menu 0, choose 1 from Menu 1 and 3 from Menu 3. GENLS will then carry out the fit. It is possible to transform Y data to the natural log of that data by selecting function choice 6 in Menu 2. If the input data is for exponential response, one would then choose 1 (polynomial model) with NP = 2 in Menu 2 in order to fit the result to a straight line. Under some circumstances, the exponential parameters can be estimated with very low bias by this procedure (10,12). In any event, it is often useful in obtaining initial parameter values for a direct fit to an exponential model. The transformed Y data can be restored to their original untransformed values by choice 0 in Menu 3. When the fitting model is nonlinear in some or all of its parameters, it is necessary to provide initial values for all the parameters (free and fixed). Especially for data with a wide range of variation, reasonably good estimates of the parameters are needed to ensure appropriate convergence. Various procedures are suggested by the program when convergence fails. One possibility is to turn on the Levenberg-Marquardt convergence-enhancing method by making choice 5 in Menu 1A. It only applies to the Deming part of the fitting procedure and has not yet been fully optimized. Nevertheless, it is worth trying for difficult situations. It can happen that when a complicated nonlinear model is being fitted with not-highly-accurate data (mediocre data), such final converged properties of the fit as estimated parameter standard deviations and parameter correlations may depend on the choice of convergence values (Menu 1A). Remember that these estimated quantities are calculated from the inverse of the design matrix of the fitting by linearization of the model when convergence requirements have been met. The correlation estimates, and to a lesser degree, the estimated standard deviations of the parameters, can then be sensitive functions of these convergence conditions. To reduce this sensitivity and to produce as good estimates as possible, it has been found useful, when the problem arises, to progressively make some of the requirements more stringent and to examine fit results after each such change. In particular, a value as small as 10**(-6) or 10**(-7) for H, the relative increment for numerical derivatives, may be necessary. But if H is set to too small a value (use H>=10**(-8), it may happen that the general-method part of the program does not converge but cycles until the maximum number of iterations is reached. A larger value of H should then be used, possibly along with a value of G as small as 0.3. Even though the parameter estimates and sum of squares may not change within the number of decimal places listed on the screen, as H is decreased below 10**(-4), obvious changes in the correlations and estimated standard deviations may sometimes appear. To ensure convergence, the Parameter Convergence Tolerance (General) factor, EPS1, should not usually be set to less than about 10**(-10). If it is set too small, cycling can also appear. Note that the residual plots (choice 4 in Menu 5) are scaled so that a square plot region is bounded by Ymax, Ymin, Xmax, and Xmin. These values, and the ratio of the Y-to-X range, are listed at the top of the plot. When choice 4 is made, selected plots both appear on the screen and are stored in GENOUTP. The composite residual values (choice 7 in the residuals menu) are plotted against the adjusted, not the input, values of X. Finally, the output contains a quantity SSE, a negative entropy objective function (7,13,14). It is not minimized directly as part of the fit but is nevertheless a useful measure of goodness- of-fit. It is primarily a measure of the differences of the residuals, and is zero when their absolute values are equal. Thus, it is measure of how well weighting has converted a heteroscedastic situation to a homoscedastic (equal variances) one. When SSE is large ( > 0.4) one can generally conclude that the fitting function and/or weighting are not appropriate. When it is small ( < 0.2), much better fitting is present. A value of the order of 0.1 is about as low as one ever finds. If the fitting function itself is appropriately chosen, then adjustment of weighting so that SF becomes close to 1 and SSE reaches a minimum, with a value preferably less than 0.2, is the recommended approach to picking the most reasonable weighting when other weighting information is absent. G. FOR FURTHER INFORMATION Please feel free to get in touch with Ross Macdonald if you have problems with or questions about GENLS, and especially if you have suggestions for future improvement. You may write to him at Dr. J. Ross Macdonald, Department of Physics and Astronomy, University of North Carolina, Chapel Hill, NC 27599-3255, call him at 919/967-5005, or reach him via E-mail at macd@gibbs.oit.unc.edu. Good fitting! REFERENCES 1. W. E. Deming, W. E. Deming, Statistical Adjustment of Data (Wiley, New York, 1943; Dover (reprint), New York,1964). 2. J. R. Wolberg, Prediction Analysis (Van Nostrand, New York, 1967). 3. D. R. Powell and J. R. Macdonald, "A rapidly convergent iterative method for the solution of the generalized nonlinear least squares problem," The Computer Journal 15, 148-155 (1972); erratum 16, 51 (1973). See also D. R. Powell, "A Fortran Program for Generalized Nonlinear Least Squares," Manual for the program GNLS, Texas Instruments Technical Report, 1975. 4. J. R. Macdonald and W. J. Thompson, "Least squares fitting when both variables contain errors: pitfalls and possibilities," Am. J. Phys. 69, 66-73 (1992). 5. B. C. Reed, "Linear least-squares fits with errors in both coordinates," Am. J. Phys. 57, 642-646 (1989); erratum 58, 189 (1990). 6. J. R. Macdonald, "Review of some experimental and analytical equations of state," Rev. Mod. Phys. 41, 316-349 (1969). 7. A. Celmins, "Analysis of residuals from multidimensional model fitting," Computers & Chemistry 8, 81-89 (1984). 8. D. York, "Least-squares fitting of a straight line," Can. J. Phys. 44, 1079-1086 (1966). 9. J. R. Macdonald, "Impedance spectroscopy: old problems and new developments," Electrochim. Acta 35, 1483-1492 (1990). 10. J. R. Macdonald and W. J. Thompson, "Strongly heteroscedastic nonlinear regression, "Commun. in Statistics - Simulation and Computation 20, 843-886 (1991). 11. D. S. Riggs, J. A. Guarnieri, and S. Addelman, "Fitting straight lines when both variables are subject to error," Life Sciences 22, 1305-1360 (1978). 12. W. J. Thompson and J. R. Macdonald, "Correcting parameter bias caused by taking logs of exponential data", Am. J. Phys. 59, 854-856 (1991). 13. H. Akaike, "A new look at the statistical model identification," IEEE Trans. Automation Control, AC-19, 716-723 (1974). 14. Y. Hiragi, H. Urakawa, and K. Tanabe, "Statistical procedure for deconvolving experimental data," J. Appl. Phys. 58, 5-11 (1985).