 
      subroutine ryork (n, x, y, wx, wy, bstart, mxiter, delmax,
     1                  ires, maxd, work, output, res, 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 York (1966) (see reference below).
 
 
      DESCRIPTION
      -----------
 
      1. The input data are observations (x(i),y(i)), i = 1,...,n, and
         weights wx(i) and wy(i) associated with observations x(i) and
         y(i) respectively, for i = 1,...,n.
 
      2. The function
 
         S = sum over i (wx(i)*(x(i)-X(i))**2 + wy(i)*(y(i)-Y(i))**2)
 
         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 York (1966) (see references below).
 
      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 cubic equation is solved and the "best" of
         the up to three possible real-valued solutions is selected as
         the slope estimate at that iteration.
 
         There are between one and three solutions of the cubic equation
         at a given iteration. If there is one real solution (and two com-
         plex solutions), the iterations proceed. If there is more than
         one real solution, the solution closest to the previous solution
         is found and the iterations proceed.
 
      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 weighted averages of the x- and y-
         observations are given by equations (19) in the York reference,
         and the equations for the two standard deviations are given
         on pages 1084 and 1085 of the same reference.
 
      7. Another algorithm for minimizing S, described by Williamson
         (1968), is implemented in the SLRPACK subroutine RWILL.
 
 
      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.
 
      WX      Real vector dimensioned at least N (unchanged on output)
              Weights associated with x-observations (commonly the
              multiplicative inverses of the variances associated with
              the x-observations).
 
              WX(i) > 0   for i = 1,...,N
 
      WY      Real vector dimensioned at least N (unchanged on output)
              Weights associated with y-observations (commonly the
              multiplicative inverses of the variances associated with
              the y-observations).
 
              WY(i) > 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.
 
      IRES    Integer scalar (unchanged on output)
 
              IRES = 1 if residuals are to be calculated
              IRES = 0 otherwise
 
      MAXD    Integer scalar (unchanged on output)
              Row dimension of RES in the calling program.
 
              MAXD must be at least N if IRES = 1, otherwise MAXD
              may be 1.
 
      WORK    Real vector dimensioned at least 10 + (9*N)
              Work vector.
 
      OUTPUT PARAMETERS
      -----------------
 
      OUTPUT  Real vector dimensioned at least 7
 
              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) = weighted average of x observations at final
                          iteration
 
              OUTPUT(6) = weighted average of y observations at final
                          iteration
 
              OUTPUT(7) = number of iterations
 
      RES     Real matrix dimensioned at least MAXD by 2
              The i-th x and y residuals are in RES(i,1) and  RES(i,2)
              respectively.
 
      IER     Integer scalar
 
              IER =  0  if there are no execution time errors or
                           warnings
 
              IER =  1  York technique regression line slope estimates
                           cannot be calculated because the data set is
                           too small
 
                        Cannot compute OUTPUT(I) for I=1,...,7
                           and cannot compute the residuals
 
              IER =  2  York technique regression line slope estimates
                           cannot be calculated because all x-observations
                           are equal
 
                        Cannot compute OUTPUT(I) for I=1,2,3,4,7
                           and cannot compute the residuals
 
              IER =  3  York technique regression line slope estimates
                           cannot be calculated because all of the
                           weights are not positive
 
                        Cannot compute OUTPUT(I) for I=1,...,7
                           and cannot compute the residuals
 
              IER =  4  York technique slope estimates have not converged
                           in allowed number of iterations
 
                        Have computed OUTPUT(I) for I=1,...,7 and residuals
                           for the final iteration
 
      EXAMPLE
      -------
 
    INPUT:
 
               N = 10
            IRES =  1
          MXITER = 20
          BSTART =  -.54600
          DELMAX =   .00010
 
              I        X(I)      WX(I)       Y(I)      WY(I)
              1        .000   1000.000      5.900      1.000
              2        .900   1000.000      5.400      1.800
              3       1.800    500.000      4.400      4.000
              4       2.600    800.000      4.600      8.000
              5       3.300    200.000      3.500     20.000
              6       4.400     80.000      3.700     20.000
              7       5.200     60.000      2.800     70.000
              8       6.100     20.000      2.800     70.000
              9       6.500      1.800      2.400    100.000
             10       7.400      1.000      1.500    500.000
 
    CALL SEQUENCE:
 
      call ryork (n, x, y, wx, wy, bstart, mxiter, delmax,
     1            ires, maxd, work, output, res, ier)
 
    OUTPUT:
 
                    IER = 0
 
                    OUTPUT(1) =  -.481
                    OUTPUT(2) =  5.480
                    OUTPUT(3) =   .071
                    OUTPUT(4) =   .362
                    OUTPUT(5) =  4.911
                    OUTPUT(6) =  3.120
                    OUTPUT(7) =  5.
 
               Weighted Residuals
            I         RES(I,1)    RES(I,2)
            1           .000       -.420
            2           .000       -.352
            3           .001        .215
            4          -.002       -.369
            5           .019        .385
            6          -.038       -.316
            7           .080        .143
            8          -.234       -.139
            9          -.084       -.003
           10           .875        .004
 
      Note:  Given initial slope bstart = -0.546, after one "iteration"
             of the York technique, the estimated slope is -0.477, in
             agreement with York (1966).
 
      PRECISION
      ---------
 
      All calculations are done in single precision.
 
      LANGUAGE
      --------
 
      The routine is written in standard Fortran 77.
 
      OTHER SUBROUTINES
      -----------------
 
      SLRPACK  subroutines  WYORK, RYBAR, RYARNG
 
      PORT     subroutines  R1MACH
 
      REFERENCES
      ----------
 
      Madansky, A. (1959). The fitting of straight lines when both
         variables are subject to error. JASA, 54, 173-205.
 
      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.
 
      CONTRIBUTORS
      ------------
 
      Sally E. Howe
      Katherine Rensenbrink
      Amy DelGiorno
      Gregory Rhoads
      Scientific Computing Division
      Center for Applied Mathematics
      National Bureau of Standards
      Gaithersburg, MD  20899
 
      NBS CONTACT
      -----------
 
      Sally E. Howe
      Center for Applied Mathematics
 
-------------------------------------------------------------------------
 
