idlastro / Robust Statistics procedures: ROBUST_LINEFIT

[Source code]

NAME
ROBUST_LINEFIT
PURPOSE
An outlier-resistant two-variable linear regression. 
EXPLANATION
Either Y on X or, for the case in which there is no true independent 
variable, the bisecting line of Y vs X and X vs Y is calculated. No 
knowledge of the errors of the input points is assumed.
CALLING SEQUENCE
COEFF = ROBUST_LINEFIT( X, Y, YFIT, SIG, COEF_SIG, [ /BISECT,
                BiSquare_Limit = , Close_factor = , NumIT = ] )
INPUTS
X = Independent variable vector, floating-point or double-precision
Y = Dependent variable vector
OUTPUTS
Function result = coefficient vector. 
If = 0.0 (scalar), no fit was possible.
If vector has more than 2 elements (the last=0) then the fit is dubious.
OPTIONAL OUTPUT PARAMETERS
YFIT = Vector of calculated y's
SIG  = The "standard deviation" of the fit's residuals. If BISECTOR 
        is set, this will be smaller by ~ sqrt(2).
COEF_SIG  = The estimated standard deviations of the coefficients. If 
        BISECTOR is set, however, this becomes the vector of fit 
        residuals measured orthogonal to the line.
OPTIONAL INPUT KEYWORDS
NUMIT = the number of iterations allowed. Default = 25
BISECT  if set, the bisector of the "Y vs X" and "X vs Y" fits is 
        determined.  The distance PERPENDICULAR to this line is used 
        in calculating weights. This is better when the uncertainties 
        in X and Y are comparable, so there is no true independent 
        variable.  Bisquare_Limit  Limit used for calculation of 
        bisquare weights. In units of outlier-resistant standard 
        deviations. Default: 6.
        Smaller limit ==>more resistant, less efficient
Close_Factor - Factor used to determine when the calculation has converged.
Convergence if the computed standard deviation changes by less 
than Close_Factor * ( uncertainty of the std dev of a normal
distribution ). Default: 0.03.
SUBROUTINE CALLS
ROB_CHECKFIT
ROBUST_SIGMA, to calculate a robust analog to the std. deviation
PROCEDURE
For the initial estimate, the data is sorted by X and broken into 2
groups. A line is fitted to the x and y medians of each group.
Bisquare ("Tukey's Biweight") weights are then calculated, using the 
a limit of 6 outlier-resistant standard deviations.
This is done iteratively until the standard deviation changes by less 
than CLOSE_ENOUGH = CLOSE_FACTOR * {uncertainty of the standard 
        deviation of a normal distribution}
REVISION HISTORY
Written, H. Freudenreich, STX, 4/91.
4/13/93 to return more realistic SS's HF
2/94 --more error-checking, changed convergence criterion HF
5/94 --added BISECT option. HF.
8/94 --added Close_Factor and Bisquare_Limit options  Jack Saba.
4/02 --V5.0 version, use MEDIAN(/EVEN)  W. Landsman
6/18 -- Fix error in computation of coeff_sig[0], W. Landsman/F. Navarate