idlastro / Math and Statistics: FITEXY

[Source code]

NAME
FITEXY
PURPOSE
Best straight-line fit to data with errors in both coordinates
EXPLANATION
Linear Least-squares approximation in one-dimension (y = a + b*x),
when both x and y data have errors   Users might be interested in 
Michael Williams MPFITEXY routines which include a number of 
enhancements to FITEXY. 
( http://user.astro.columbia.edu/~williams/mpfitexy/ )
CALLING EXAMPLE
FITEXY, x, y, A, B, X_SIG= , Y_SIG= , [sigma_A_B, chi_sq, q, TOL=]
INPUTS
x = array of values for independent variable.
y = array of data values assumed to be linearly dependent on x.
REQUIRED INPUT KEYWORDS
X_SIGMA = scalar or array specifying the standard deviation of x data.
Y_SIGMA = scalar or array specifying the standard deviation of y data.
OPTIONAL INPUT KEYWORD
TOLERANCE = desired accuracy of minimum & zero location, default=1.e-3.
OUTPUTS
A_intercept = constant parameter result of linear fit,
B_slope = slope parameter, so that:
                ( A_intercept + B_slope * x ) approximates the y data.
OPTIONAL OUTPUT
sigma_A_B = two element array giving standard deviation of 
         A_intercept and B_slope parameters, respectively.
         The standard deviations are not meaningful if (i) the
         fit is poor (see parameter q), or (ii) b is so large that
         the data are consistent with a vertical (infinite b) line.
         If the data are consistent with *all* values of b, then
         sigma_A_B = [1e33,e33]  
chi_sq = resulting minimum Chi-Square of Linear fit, scalar
q - chi-sq probability, scalar (0-1) giving the probability that
       a correct model would give a value equal or larger than the
       observed chi squared.   A small value of q indicates a poor
       fit, perhaps because the errors are underestimated.   As 
       discussed by Tremaine et al. (2002, ApJ, 574, 740) an 
       underestimate of the errors (e.g. due to an intrinsic dispersion)
       can lead to a bias in the derived slope, and it may be worth
       enlarging the error bars to get a reduced chi_sq ~ 1
COMMON
common fitexy, communicates the data for computation of chi-square.
PROCEDURE CALLS
CHISQ_FITEXY()            ;Included in this file
MINF_BRACKET, MINF_PARABOLIC, ZBRENT    ;In IDL Astronomy Library 
MOMENT(), CHISQR_PDF()     ;In standard IDL distribution
PROCEDURE
From "Numerical Recipes" column by Press and Teukolsky: 
in "Computer in Physics",  May, 1992 Vol.6 No.3
Also see the 2nd edition of the book "Numerical Recipes" by Press et al.
In order to avoid  problems with data sets where X and Y are of very 
different order of magnitude the data are normalized before the fitting
process is started.     The following normalization is used:
xx = (x - xm) / xs    and    sigx = x_sigma / xs    
                      where xm = MEAN(x) and xs = STDDEV(x)
yy = (y - ym) / ys    and    sigy = y_sigma / ys    
                      where ym = MEAN(y) and ys = STDDEV(y)
MODIFICATION HISTORY
Written, Frank Varosi NASA/GSFC  September 1992.
Now returns q rather than 1-q   W. Landsman  December 1992
Use CHISQR_PDF, MOMENT instead of STDEV,CHI_SQR1 W. Landsman April 1998
Fixed typo for initial guess of slope, this error was nearly
      always insignificant          W. Landsman   March 2000
Normalize X,Y before calculation (from F. Holland) W. Landsman Nov 2006