James L. Dye1 ond Vincent A. Nicely Michigan Stote University E O S ~Lansing, Michigan 48823
I(
A General Purpose Curvefitting Program for Class and Research Use
The need for statistical treatment of data has been recognized for many years. It is standard procedure to require undergraduates to fit straight lines to their laboratory data by the method of leastsquares (1, 8 ) . The increasing availability of programable calculators and more sophisticated digital computers has reduced considerably the effort required to make such analyses. Although several articles in THIS JOURNAL (3-6)and in others (7-9) have described methods for solving more complex problems, the work required to program them discourages their widespread use in undergraduate laboratory courses, and often even in research. The ready availability of computational procedures and devices should insure that statistical analysis of data be a normal part of research. Such analysis is particularly important when a complex equation (or equations) containing adjustable parameters is to be fitted to the data. In such cases it is important to have error estimates on the constants and information about the mutual dependence of the constants. Unfortunately, a "good-fit" to the data does not necessarily mean that the parameters are individually well-determined. Because of the complexity of general leastsquares procedures, however, it is frequently only the linear problems which are subjected to statistical analysis. More complex problems are often "linearized" by selecting experimental conditions which yield a series of straight lines whose slopes and/or intercepts can be fitted by another straight line. A general curvefitting program which provides statistical information can provide a test of the goodness of fit of an equation without imposing special constraints upon the data collection. In addition, it can provide for the proper weighing of the individual data points, a procedure which is often not followed, even for the case of simple equations. Finally, the use of numerical integration procedures can provide a fit to the data of differential equations such as those frequently encountered in chemical kinetics. This paper describes the structure, availability, and use of a curvefitting program which was developed with these objectives in mind. Objectives of the Program
Chemists readily accept the need to use complex instruments in their research even though they might not have a detailed knowledge of the working parts of such instruments. To be sure, they understand the general nature of the equipment, but, for example, a
' T o whom correspondence should he ddressed.
lack of knowledge of the electronics behiid the front panel does not prevent the chemist from using an nmr instrument. It is important to know that the instrument is functioning properly. The usual test of this for any instrument is the output which is obtained by using a known standard sample; that is, the instrument is checked or calibrated by its response to a known input. It is our contention that complex but general computer programs can be used by chemists who are not familiar with programing techniques. Of course, just as a single-purpose instrument can be made exceedingly easy to use, a program designed for a specific computation can also be written so that the user need only punch the data in specific columns on cards. Just as an instrument designed for a specific measurement can be made very efficient, so a single-purpose computer program can be made more efficient than a multipurpose program. For repeated application to a given problem, then, specific curvefitting routines designed for the job a t hand will continue to be useful. The program described in this paper is intended for application to a wide variety of curvefitting and equation-solving problems. Because the program has been designed to handle problems of a general nature with reasonable efficiency, it is relatively complex. The user need not be concerned with this complexity, however, since only a short subroutine needs to be changed to switch from one equation to another. The program is designed to handle problems which are either linear or nonlinear in the adjustable parameters [L(a) is a linear function of a if (bL/ba) is independent of a]. In order to handle nonlinear problems it is necessary to proceed from a set of initial estimates of the parameters which must be provided by the user, and should be carefully chosen in order that convergence to the "correct" set of parameters occurs. The surest way to gain confidence in the fitting procedure is to "calibrate" the system with calculated "data" which are known to be exactly fitted by the equation being tested or which have known standard deviations from the equation. Thus, while the "black-box" approach relieves the user of the difficult task of writing a complex computer program, it simultaneously requires that he "calibrate" the program with a "known" problem which closely resembles the "unknown" problem to be handled later. Incidentally, by using calculated "data" having various known amounts of standard error, it is possible to determine with what accuracy the real data must be measured in order to yield given standard deviations of the adjustable parameters. Obviously, such a procedure could be very useful in the design of experiments. Volume 48, Number 7, July 1971
/
443
Description of the Program Design Objectives and Procedures
An experienced programer can readily adapt a Fortran program written for use on a particular computer to a form compatible with another computer. In order to avoid conversion problems, this program was written in conformity with USASI standards. It has been run successfully on CDC 3600 and 6500 computers as well as on IBM 7094 and 360-75 computers. It requires about 13,000 words of memory on the CDC computers and about 16,000 on the IBM computers. Three options are available in this program. 1 ) Fit an arbitrary function of not more than 20 parameters and four variables (dependent plus independent) to 8. dats set containing not more than 99 points (an expanded version which can handle up to 20 parameters and 500 dsta points for the case of two variables has also been written). 2) Solve arbitrary algebraic equations having not more than 20 unknowns. 3) Minimize (or maximize) an arbitrary functional having not more than 20 pmuneters.
The program also has a provision for solving sets of not more than 10 simultaneous initial-value first-order differential equations. This procedure uses the RungeKutta method (10, If) to provide values of the dependent variable to be fitted to the data. This feature is of particular value for problems in chemical kinetics since the time rate of change of a variable is the natural mathematical form but the variable itself is observed as a function of time. The program has been designed to minimize the time required by the user while making efficient use of computer time. The user must specify ten control constants, his data (including estimates of relative variance) and a short subroutine in standard form which gives the 'equation(s) to be used for the particular problem. The user can select either of two minimization methods; whether and how weights are to be used; whether auxiliary constants are to be read; the maximum number of iterations to be allowed; the convergence criterion; and, when differential equations are used, the integration accuracy required. When a functional of the form y = f(x) is fit to the data, the user can obtain a print-plot of the experimental and calculated points. Alternatively, the value of the residual can be plotted against one of the variables. Two methods can be used to carry out the options described above. Method A uses the matrix technique described by Wentworth (3) modified to include the iteration step described by Pitha and Jones (7) as their method V. Method B uses the procedure of Powell (7). Method A has been found to converge more rapidly than Method B in a number of cases involving curvefitting. However, since this may not be true in general, the user is encouraged to select the procedure most suited to his particular problem. What the Program Does
Under either option 1 or 2, the program minimizes the functional
in which n is the number of points, F, is a residual 444
/
Journol of Chemical Edumtion
defined in such a way that it would approach zero for all i as the parameters approach their "best" values if the data were completely free of errors. W' is a weight which can either be set to unity for all points (which is referred to as "not using weights") or else can be calculated by the program as Wi
=
(2)
in which
In eqn. (3), m is the number of variables, X xis the kth variable and the partial derivative is evaluated a t the ith data point. The quantity u*: is the variance of the lcth variable at the ith point. The use of weights as given by eqns. (2) and (3) has been discussed by Wentworth (3). It gives maximum likelihood estimates of the parameters when the data have erron described by independent normal distributions. When weights are calculated in this way, the fit obtained is independent of the choice of independent and dependent variables since for each variable the minimization assumes variances in the amount s~ecified (3). When the program uses method A, an estimate of the variancecovariance matrix is available at the end of the calculation. Multiple correlation coefficients (13) and estimates of standard error for each parameter are d e rived from this matrix. All of these quantities are printed at the end of a calculation. Since such information is applicable only for the curvefitting problem, it is not printed when either option 2 or option 3 is selected. When option 3 is chosen, the program minimizes 9 = F (parameters)
(4)
with respect to the parameters; in which F is any function which depends on parameters. Option 3 requires that NOPT = 1 and IMETH # 0 (see later). Use of the Progrom
Each data set to be used is preceded by cards which give the necessary control instructions, constants for the data set, and initial estimates of the parameters. The card input is diagrammed in Figure 1 to show the format and logic used. Definitions of the names used are given in the table. The structure of the read statement and the format used, together with the sample input supplied with the program, should permit the user to construct a data deck to fit his problem. The program must be given information about the equation to be fitted to the data as well as the data themselves. This is done via a short subroutine (EQN) which calculates a value of Ft (here called RESID) for use in eqn. (1) by the main program. Figure 2 gives a listing of EQN for the particular example to be discussed now. The example chosen is given in both the differential and the integrated form to illustrate the use of both procedures. (The program uses only one of these a t a time as specified by the control card.) Suppose that the total absorbance A(t), of two species is followed a t a given wavelength as a function of time, t, and, further, that the two species disappear
Definition of Terms for Figure 1
stort
NOPT READ NOPT. IYETH. ITUAX. I W T . 1RX.ISMIN
IMETH
IPLT, NCST. T E S T . KVAR
ITMAX
815, F10.0, 1 5
The number of seta of mesaurablas t o be uaed in the sum given by (1). 0 < NOPT QO. When NOPT is aero or negstive, the emoution ia terminated. The ohoice of minimination methoda: nero chooses method one; nonzero chooses method two. The msuimum number of iterations t o be permitted i f omvereenoe is not obtained.
