Sherril D. Christian
The University of Oklahoma Norman
Graphical Least Squares Analysis
I n spite of the conceptual simplicity and utility of the method of least squares, many experimenters have only a cursory acquaintance with the technique and its applications in parameter estimation problems. This may he due in part to the somewhat cumbersome formulas which arise in fitting any but the simplest functions with the method. A geometrical or graphical approach to least squares analysis is apt to seem more palatable to the nonmathematician than a rigorous analytical treatment of the subject. In addition, graphical least squares methods can be used profitably in treating problem that would he difficult to solve in any other way. I n this paper, a graphical method of solving least squares problems will be described and illustrated with several examples. Theory
Basically, the principle of least squares is to minimize the sum of the weighted squares of the residuals (i.e., deviations between observations of a variable subject to error and values calculated from a functional relation involving other variables and adjustable parameters). Ordinarily, the values of the adjustable parameters which lead to a minimum in the sum of the weighted squares of the residuals are found by differentiating the sum with respect to each parameter in turn, and equating the derivatives to zero. The resulting set of simultaneous equations is solved for all the parameters (1,s). Recently, Wentworth has given a useful summary of procedures for fitting linear and nonlinear functions by least squares (3). The simplest problem involving least squares is that of finding the value of a single parameter to represent a collection of values of an observable, measured under conditions as nearly reproducible as possible. In this case, the least squares solution may be obtained by minimizing the function
mean square deviation, and will be utilized rather than the sum of the squares of the residuals because of the close relation between r(a) and the standard deviation. It is easy to show that the value of a that minimizes r(a) is g, the average value of y. Further, the standard deviation of the parameter a is related to that of the y values by uaz = uyZ/N,and for large N it is expected that a, will be approximately equal to the minimum r(a). In many curve-fitting problems, it is necessary to estimate the standard deviations of the parameters as well as of the ohservahles; hence, it is important that the preceding relations he well understood. One method of solving the averaging problem described in the previous paragraph is to plot r(a), calculated from equation (I), versus a, and to seek a,!., the value of a that yields a minimum in r(a). This is undeniably a difficult way to solve an easy problem; nonetheless, it is instructive to consider the graphical solution. In Figure 1, r(a) is plotted versus a for a collection of 100 numbers which were taken from a table of random numbers, normally distributed, having a standard error of unity and a mean value of zero (4). It is expected that a at the minimum will be nearly zero and that r(a) will be approximately 1. The actual values deduced from the figure are a ~ = . 0.011 and r,;, = 1.021. The question now arises, is there a simple relation between the variation of r(a) in the vicinity of the minimum and the standard deviation in a? To answer this question, the value of r at a = a,:, rt ua may be expressed as
Note that the first term under the radical is just N
r2(aen), that the second term vanishes (since where N represents the number of measured values of an observable, y, and a is the only parameter involved.' The function to be minimized, r(a), is termed the root
C y,
=
:= 1 .-A
N a ~ , ) and that uaZ=rZ(a,;,)/N, provided N is large. Therefore, the preceding equation reduces to
' If the yr values are not all of equal weight, it is necessary
to incorporate weight factors into the equation defining ?(a).
Throuehout. it will be assumed that all the st's are weighted equall;, and that no observable other then y is&bject to experimental error. It is not difficult, however, to extend the present treatment to problems in which the y,'s are not all weighted equally. See references (1-3) for discussions of the proper use of weighta.
604 / Journal of Chemical Education
where the binomial expansion has been applied in extracting the square root. According to equation (2), if r(a,;,) is incremented by the fractional part of its value 1/2(N - l), a number should be obtained which
on the parameters a, b, mum.
. . . in the vicinity of the mini-
Applications
The graphical least squares method, involving a single adjustable constant, has been useful in solving a variety of nonlinear parameter estimation problems. For example, vapor density data for mixtures of two interacting species frequently can he explained in terms of one hetero-association constant (6). I n the system trifluoroacetic acid-acetone, Ling, Christian, and Affsprung have shown that the single hetero-association reaction Figure 1.
Root mean square deviation curve for collection of random
numbers.
is equal to r(a) at a = a,i, +. re. Hence, a. can be estimated from the points of intersection of the line 1/2(N - l ) ] with the curve r(a) r(a) = ?(a,;.) [l versus a. The horizontal line in Figure 1, representing r = r(ami.) [l 1/2(N - l ) ] = 1.026, intersects the curve at points having values of a of 0.100 and -0.078, from which the estimate of aa is 0.089. However, from the relation aa2 = avZ/Nit would be predicted that a= = 0.100 for the collection of 100 numbers drawn from the random number table. Thus, reasonable est,imates of a , a, and ami. have been inferred graphically, without utilizing any information about the characteristics of the parent distribution from which the 100 numbers were selected. Fort,unately, equation (2) can be used without modification to estimate a. in any least squares problem in which the root mean square deviation has the form
+
is sufficient to account for the interaction of vapors of the two compounds (6). At low pressures, values of the measured total pressure, p, and the individual densities of the acetone and the trifluoroacetic acid are related by the equations
+
8
..
where f*(a) is a function of the single parameter a and one or more observables measured concomitantly with y,, which is assumed to he the only variable subject to error. In this case, conventional least squares analysis indicates that
If r(a & a,) is evaluated by means of equation (3), using Taylor's .. series to expand ft(a an) and substiN
*
tuting aga . /x . (djd/ba)2for aa2,equation (2) is again de.-L
duced. I n other words, both a, and ancan he estimated graphically, from a plot similar to Figure 1, for data obeying any functional relation of the form y = f(a, x, z, . . .) ; where x, z, . . . are observables (assumed free from error), y is an observable subject to error and a is the single parameter to be adjusted in minimizing r(a). When functions involving more than one parameter are to be fitted by least squares, an absolute minimum in the root mean square deviation must again be located. The values and standard errors of all the parameters may be estimated from the dependence of r
where rAand a F are the "formal pressures" of acetone and trifluoroacetic acid, defined in terms of the measured densities, p.& and p ~ and , the molecular weights, MAand MF, of the monomers of the two compounds; pA and p, are the monomer partial pressures; KF, is the known dimerization constant of trifluoroacetic acid, determined from measurements of the vapor density of pure trifluoroacetic acid; and KAFis the hetero-association constant to he determined by least squares analysis. From equations (4) it is possible to obtain a calculated value of the total pressures, p,OaLC, for each data set by choosing a trial value of KAF,solving equations (4b) and (44 simultaneously for p~ and p ~ and , substituting these values of the monomer partial pressures into equation (4a). The root mean square deviation may be written in the form
are the observed values of the where p,, aAi, and total pressure and the formal pressures for the ith data set. In defining r(KAF)it is assumed that all the error resides in the measured pressures and that all the pi values should be weighted equally. Obviously, equation (5) is formally equivalent to equation (3); therefore, a plot of r(KAF) versus KAFmay be analyzed by the method described in the previous section to estimate KAF, axAsand up. The task of calculating points on the ~ ( K A Fversus ) KAF curve would be somewhat tedious if performed with a desk calculator; however, a digital computer can do the required calculations very rapidly. Figure 2 shows the computed root mean square deviation plot for 38 data sets for the acetone-trifluoroacetic acid system a t 30'; the horizontal line represents r-. [l 1/2(N-l)]. Using equation ( 2 ) , as before, it is found that u ~ , ,= 0.0044 mm-I, KAF = 0.178 mm-' and up = 0.111 mm.
+
Volume 42, Number 7 7, November 7 965
/
605
An example of the application of the computational scheme to infrared spectral data has been reported by Worley (8). He investigated the hydration of 2pyrrolidone in CClr by measuring the formal molality of the lactam, fp, the formal molality of water, fw, and the water activity, a*. From previous measurements on dry solutions of 2-pyrrolidone in CCla he had determined the value of the dimerization constant, KP,. His hydration data at 25" were consistent with the equations 0.1101.
' 0.17
.
.
.
1
0.1 8
.
,
,
.
I 0.1 9
KAF (mm-'I Figure 2. Root mean square deviation curve for aselone-triflurrroocetis acid system.
Equations (4) are easily modified to allow the introduction of additional hetero-association constants. For example, it might be plausible in treating a given mixed vapor system to assume the presence of both the 1 :1 complex (AF) and the 1:2 complex (AF2). I n this case, the equations for p and ?rA would each include the additional term KAF,~?PF'and the equation for TF would require the add~tlonof ~KAF,PAPI~, where KAF, is the hetero-association constant for formation of the 1:2 complex from the monomers. The procedure for calculating r would be exactly the same as in the one parameter case, although in this problem, ppa'O would be a function of KAp,KAP, ?rAi, and a=,; and r would be a function of both KAp and K A ~ .Note that in a two-parameter curve-fitting problem, the factor in the denominator of the equation defining r is N - 2, not N - 1 as in the one-parameter problem ( I ) . The versatility of the parameter estimation procedure described above may be illustrated by noting some of the other problems to which it has been applied in this laboratory. While equations (4) were written in a form suitable for analyzing binary vapor density data, similar equations have been formulated for use in calculating hetero-association constants from: l ) vapor pressure data for mixtures of two solutes in a nonvolatile solvent (7); 2) infrared spectral data for two interacting solutes in solution; 3) partition and water solubility data for solutions of polar solutes in organic solvents; and 4) colligative property data.
Figure 3. Root mean square deviation contour curve= for the 2-pyrrolidone-roter system.
606
/
lournol of Chemical Educofion
where Afw is the spectrally determined increase in the solubility of water in CClr (at a given water activity) due to the presence of the lactam; mp and mw are the monomer molalities of 2-pyrrolidone and water; mw0 is the water monomer molality at unit water activity; and Kpw and K P , are ~ hetero-association constants for formation of the lactam monomer monohydrate and dimer monohydrate. Equations (5b) and (5c) readily may be solved for mp and mW, for a given choice of KPW and KP,w. If the monomer molalities are substituted into equation ( 5 4 , Afw may be calculated for each point. I n this way, the root mean square deviation of Afw may be computed forvarious combinations of KPWand Kp,w. Since two parameters are involved, it is convenient to draw contour lines of constant r(Kpw, K p 2 ~on ) a plot of Kpwversus Kp3w. Figure 3 indicates results of calculations on 60 data sets for the 2-pyrrolidone-water-CC14 system; the minimum value of r is 5.51 X molal and the lines represent the contour curves r = r,i.. [I 1/2(N - l ) ] and r = rmi,, [1 1/(N - 1) 1. From the location of the horizontal and vertical tangents to the inner contour line 3.3 molal-' and U K ~ = , it is estimated that aKpv 900 rn~lal-~. The best values of the constants are KPW 31.6 molal-I and KpZwc 14,100m ~ l a l - ~ .The outer contour line is not used in estimating the standard deviations, but is included in the figure to give some indication of the shape of the surface. An example of a two parameter, nonlinear, least squares analysis has been reported by Bennett and Franklin (9). They treated the kinetic data of Newman, Lee and Garrett (10) (on the liberation of nitrogen during the decomposition of benzazide in dioxane) in t,he form
+
-
-
+
where V is the cumulative volume of Nz produced between time zero and time t, k is the first order rate constant and V , is the total volume liberated up to time infinity. Although Newman, et al., assumed V , to be known and determined k as a single least squares parameter, Bennett and Franklin treated both V , and k as adjustable parameters in the analysis. They computed the value of k which minimizes the sum of the squares of the residuals for each of several choices of V* and concluded that the best value of k is that for which the smallest sum of squares is obtained. In this way they determined the values V- = 95.72 ml, k = 0.648 hr.-', and Q = 0.002 ml; whereas Newman,
Figure 4.
Contour curves for dot.
of Newman, Lee, and
Garrett (10).
et al., had originally reported V- = 95.20 ml, k = 0.674 hr-l, and a, = 0.005 hr-'. The present least squares scheme is easily adapted to the problem of obtaining best values of V, and k. Note that the method used by Bennett and Franklin involves the implicit assumption that values of loglo [(V- - V)/V,] should be weighted equally; however, one might alternately, and perhaps more reasonably, assume that the values of V deserve equal weights. Figure 4 shows the contour curves r = r,,,i,,. (1 1/2(N - I ) ] obtained by using each of the two weighting procedures, as well as points corresponding tor,,,. and the values of the parameters reported by Bennett and Franklin and by Newman, et al. (denoted by BF and N). Curve I results from the analysis in which log,, [(V, - V)/V,] values are weighted equally and curve I1 is inferred by assigning equal weights to all the V values. From the location of curve I it is estimated that k = 0.648 hr-'. , a,. = 0.004 hr-'. V- = 95.73 ml, and a"_ = 0.08 ml; from curve 11, the values k = 0.669 hr-l, uk = 0.008 hr-I, V, = 94.98 ml, and a~~ = 0.40 ml are calculated. As expected, the solution obtained from curve I differs only slightly from that reported by Bennett and Franklin, except that the value of rrx is approximately twice their value. This is a conseauence of the fact that the two least squares para&eters are not independently variable; -thus, the standard error i n k is somewhat greater than that calculated assuming V* to be known with certainty. The analysis based on curve I1 (which, in the opinion of the author is the preferable one for the data being treated) yields values of the parameters more nearly in agreement with those reported by Newman, et al. Note that the solutions deduced from curves I and I1 differ significantly; this example illustrates that the choice of weighting factors may considerably influence
+
the values of parameters obtained by least squares analysis. The preceding examples indicate that least squares solutions can be obtained simply and rapidly by using a digital computer to calculate root mean square deviations and plotting r versus values of the adjustable paramet,ers to locate the minimum. The technique described is equally applicable to linear and to nonlinear functions, and completely avoids problems of instability which frequently arise in matrix methods for obtaining least squares solutions. Since the method eliminates the necessity of formulating and solving the normal equations, programming a given problem for the computer is generally quite simple. If more than two variable parameters are involved, it may be somewhat difficult to locate the absolute minimum in r. However, we have successfully handled problems involving several adjustable parameters by employing the so-called "partan" (parallel tangent) method for seeking maxima or minima (11). A number of useful techniques for finding maxima or minima are discussed by Wilde (12). This work was supported by grants from the National Institutes of Health (GM-08972) and the Petroleum R e search Fund of the American Chemical Society (PRF 1481-A5). The author is indebted to Mr. Delbert Mueller and Dr. Chii Ling for performing many of the computations and to the University of Oklahoma Computer Laboratory for providing generous blocks of computer time. Literature Cited (1) DEMING,W. E., "Statistical Adjustment of Data," John Wiley and Sonns, Inc., New York, 1943. (Also available as a Dover reprint, 1964.) F. A,, "Introduction to the (2) Mooo, A. M., A N D GRAYBILL, Theory of Statistics," 2nd ed., McGraw Hill Book Company, Inc., New York, 1963. (3) WENTWORTH, W. E., J. CHEM.EDUC.,42.96 (1965). (4) DEMING,W. E., Op. cit., p. 252. (5) LING, C. "The Determination of Hetero-association Cons t ~ n t by s a. New Vapor Density Method," PhD Dissertation, University of Oklahoma, 1964. H. E., in prep(6) LING,C., CHRISTIAN,S. D., AND AFFSPRUFQ, aration. (7) TAHA, A,, "A New Method for the Study of Hydrogen Bonding in Solutions of Volatile Compounds in Nonvolatile, Nonpolar Solvents," PhD Dissertation, University of Oklahoma, 1965. (8) WOHLEY,J. D., "A Spectroscopic Investigation of the Self Association and Hydration of Lactams in Carbon Tetrachloride," PhD Dissertation, University of Oklahoma, 1964. .. ( 9 ) BENNETT, C. A., AND F R A N ~ L ~ NN. , L., t i ~ t a t i s t i eAndy& s~ in Chemistry and the Chemical Industry," John Wiley and Sons, Inc., New York, 1954, p. 235. A. B., J . Am. (10) NEWMAN, M. S., LEE, S. H., AND GARRETT, Chem. Soc., 69, 113 (1947). (11) o., "The . . SHAH,B. v., BUEHLER,R. J., AND KEMPTHORNE, ~ e t h a dof Pzrallel ~ a n p n t (PARTAN) e for Finding an Optimum," Technical Report No. 2, Office of Naval Research Contract Nonr-530 (05), Iowa State University Statistical Laboratory, Ames (Apr. 1961, rev. Aug. 1962). (12) WILDE,D. J., "Optimum Seeking Methods," Prentiee Hall, Inc., Englewood Cliffs, New Jersey, 1964.
NOTEADDEDIN PROOF: The author recently learned that Professor L. G. Sillen and co-workers in Stockholm have previously developed z graphical Least squares method similar to and in fact ronsiderahly more comprehensive than the technique described here. References to some of the recent work of the Swedish group include: N., Acta C h m . Scand., 16, 173 SILLEV,L. S., d c t a Chem. Scand., 18, 1085 (1964); INGRI, (1962); SILLEN,L. G., Acla C h m . Sand., 16, 159 (1962). Volume 42, Number 1 1, November 1965
/
607
