subroutine rwill(n, x, y, u, v, bstart, mxiter, delmax, 1 output, xres, yres, work, ier) --------------------------------------------------------------------- Package: SLRPACK Version: October, 1985 ---------------------------------------------------------------------- PURPOSE ------- This user-callable subroutine estimates simple linear regression coefficients when both variables are subject to errors which are not necessarily homogeneous in variance. The algorithm for the estimation is due to Williamson (1968) (see reference below). DESCRIPTION ----------- 1. The input data are observations (x(i),y(i)), i = 1,...,n, and variances u(i) and v(i) associated with observations x(i) and y(i) respectively, for i = 1,...,n. 2. The function S = sum over i ((x(i)-X(i))**2/u(i) + (y(i)-Y(i))**2/v(i)) is minimized with respect to the fitted values X(i) and Y(i), i = 1,...,n. The set of points (X(i),Y(i)) i = 1,...,n, are constrained to be collinear. This subroutine straightforwardly implements the calculations described in Williamson (1968) (see references). 3. The slope estimate is obtained iteratively, with an initial estimate supplied by the user. One suggested initial estimate is the geometric mean of the OLS-y and OLS-x slopes obtained by executing SLRPACK subroutine RGM (see documentation of that subroutine). 4. At each iteration a linear equation is solved and the solution is selected as the slope estimate at that iteration. 5. The iterations terminate when 1) the angle between the lines est- imated in two successive iterations is less than user-set DELMAX, or 2) the maximum number of iterations (user-set MXITER) has been reached, whichever occurs first. 6. The equations for the standard deviations and the weighted averages are given in Williamson (1968). 7. Another algorithm for minimizing S, described by York (1966), is implemented in the SLRPACK subroutine RYORK. INPUT PARAMETERS ---------------- N Integer scalar (unchanged on output) Number of observations. N must be at least 3. X Real vector dimensioned at least N (unchanged on output) The X-observations. Y Real vector dimensioned at least N (unchanged on output) The Y-observations. U Real vector dimensioned at least N (unchanged on output) U(i) is the variance associated with X(i). U(i) > 0.0 for i = 1,...,N. V Real vector dimensioned at least N (unchanged on output) V(i) is the variance associated with Y(i). V(i) > 0.0 for i = 1,...,N. BSTART Real scalar (unchanged on output) Initial slope estimate. MXITER Integer scalar (unchanged on output) Maximum number of iterations allowed. DELMAX Real scalar (unchanged on output) Angular distance in radians such that execution terminates if the angular distance between lines estimated in two successive iterations is less than DELMAX. WORK Real vector dimensioned at least 5 * N Work vector. OUTPUT PARAMETERS ----------------- OUTPUT Real vector dimensioned at least 11 OUTPUT(1) = Slope estimate at final iteration. OUTPUT(2) = Intercept estimate at final iteration. OUTPUT(3) = Standard deviation of slope estimate. OUTPUT(4) = Standard deviation of intercept estimate. OUTPUT(5) = Standard deviation of the slope estimate multiplied by the chi-square test statistic weighted sum of squares divided by degrees of freedom. OUTPUT(6) = Standard deviation of the intercept estimate multiplied by the chi-square test statistic weighted sum of squares divided by degrees of freedom. OUTPUT(7) = Weighted average of x observations at final iteration. OUTPUT(8) = Weighted average of y observations at final iteration. OUTPUT(9) = Root mean square error in the x-direction. OUTPUT(10)= Root mean square error in the y-direction. OUTPUT(11)= Number of iterations. XRES Real vector dimensioned at least N The x residuals. YRES Real vector dimensioned at least N The y residuals. IER Integer scalar IER = 0 if there are no execution time errors or warnings IER = 1 Williamson technique regression line slope estimates cannot be calculated because the data set is too small Cannot compute OUTPUT(I) for I=1,...,11 and residuals IER = 2 Williamson technique regression line slope estimates cannot be calculated because all X-observations are equal Cannot compute OUTPUT(I) for I=1,...,11 and residuals IER = 3 Williamson technique regression line slope estimates cannot be calculated because all variances are not greater than 0.0 Cannot compute OUTPUT(I) for I=1,...,11 and residuals IER = 4 Williamson technique regression line slope estimates cannot be determined because the sum of the weights is zero (equation for W(i) is on page 1846 of Williamson). Cannot compute OUTPUT(I) for I=1,...,11 and residuals IER = 5 Uncertainties in slope and intercept cannot be calculated since the final slope is 0. Cannot compute OUTPUT(I) for I=3,4,5,6,9,10 and residuals IER = 6 Uncertainties in slope and intercept cannot be determined because certain calculations involve division by zero Cannot compute OUTPUT(I) for I=3,4,5,6,9,10 and residuals IER = 7 Williamson technique slope estimates have not converged in allowed number of iterations Have computed OUTPUT(I) for I=1,...,11 only for the last iteration. EXAMPLE ------- Input: N = 10 IRES = 1 MXITER = 20 BSTART = -.54600 DELMAX = .00010 I X(I) U(I) Y(I) V(I) 1 0.000 0.00100 5.900 1.00000 2 0.900 0.00100 5.400 0.55556 3 1.800 0.00200 4.400 0.25000 4 2.600 0.00125 4.600 0.12500 5 3.300 0.00500 3.500 0.05000 6 4.400 0.01250 3.700 0.05000 7 5.200 0.01667 2.800 0.01429 8 6.100 0.05000 2.800 0.01429 9 6.500 0.55556 2.400 0.01000 10 7.400 1.00000 1.500 0.00200 Call sequence: subroutine rwill(n, x, y, u, v, bstart, mxiter, delmax, 1 output, xres, yres, work, ier) Output: IER = 0 OUTPUT(1) = -.481 OUTPUT(2) = 5.480 OUTPUT(3) = .071 OUTPUT(4) = .362 OUTPUT(5) = .0576 OUTPUT(6) = .2919 OUTPUT(7) = 4.911 OUTPUT(8) = 3.120 OUTPUT(9) = .2888 OUTPUT(10)= .2776 OUTPUT(11)= 4. Weighted Residuals I XRES(I) YRES(I) 1 .000 -.420 2 .000 -.352 3 .001 .214 4 -.002 -.369 5 .019 .385 6 -.038 -.316 7 .080 .143 8 -.234 -.139 9 -.085 -.003 10 .874 .004 PRECISION --------- All calculations are done in single precision. OTHER SUBROUTINES ----------------- PORT subroutine R1MACH LANGUAGE -------- The routine is coded in standard Fortran 77. REFERENCE --------- Riggs, D. S., Guarnieri, J. A., and Addelman, S. (1978). Fitting straight lines when both variables are subject to error. Life Sciences, 22, 1305-1360. York, D. (1966). Least-square fitting of a straight line. Canadian Journal of Physics, 44, 1079-1086. Williamson, J. H. (1968). Least-squares fitting of a straight line. Canadian Journal of Physics, 46, 1845-1847. AUTHORS ------- R. Lindstrom Inorganic Analytical Reasearch Division Center for Analytical Chemistry Sally E. Howe Scientific Computing Division Center for Applied Mathematics National Bureau of Standards Gaithersburg, MD 20899 NBS CONTACT ----------- Sally E. Howe Scientific Computing Division ------------------------------------------------------------------------