;$Id: ladfit.pro,v 1.9 1997/01/15 03:11:50 ali Exp $
;
; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       LADFIT
;
; PURPOSE:
;       This function fits the paired data {X(i), Y(i)} to the linear model,
;       y = A + Bx, using a "robust" least absolute deviation method. The 
;       result is a two-element vector containing the model parameters, A 
;       and B.
;
; CATEGORY:
;       Statistics.
;
; CALLING SEQUENCE:
;       Result = LADFIT(X, Y)
;
; INPUTS:
;       X:    An n-element vector of type integer, float or double.
;
;       Y:    An n-element vector of type integer, float or double.
;
; KEYWORD PARAMETERS:
;  ABSDEV:    Use this keyword to specify a named variable which returns the
;             mean absolute deviation for each data-point in the y-direction.
;
;  DOUBLE:    If set to a non-zero value, computations are done in double 
;             precision arithmetic.
;
; EXAMPLE:
;       Define two n-element vectors of paired data.
;         x = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89, -0.12, 1.41, $
;               2.95, 2.18,  3.72, 5.26]
;         y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98, -2.87, -1.66, $
;              -0.78, -2.61,  0.31,  1.74]
;       Compute the model parameters, A and B.
;         result = ladfit(x, y, absdev = absdev)
;       The result should be the two-element vector:
;         [-3.15301, 0.930440]
;       The keyword parameter should be returned as:
;         absdev = 0.636851
;
; REFERENCE:
;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
;       Cambridge University Press
;       ISBN 0-521-43108-5
;
; MODIFICATION HISTORY:
;       Written by:  GGS, RSI, September 1994
;       Modified:    GGS, RSI, July 1995
;                    Corrected an infinite loop condition that occured when
;                    the X input parameter contained mostly negative data.
;       Modified:    GGS, RSI, October 1996
;                    If least-absolute-deviation convergence condition is not
;                    satisfied, the algorithm switches to a chi-squared model.
;                    Modified keyword checking and use of double precision.
;       Modified:    GGS, RSI, November 1996
;                    Fixed an error in the computation of the median with
;                    even-length input data. See EVEN keyword to MEDIAN.
;-

FUNCTION Sign, z, d
  RETURN, ABS(z) * d / ABS(d)
END

FUNCTION MDfunc, b, x, y, arr, aa, absdev, nX, Double = Double
  ON_ERROR, 2
  eps = 1.0e-7
  arr = y - b*x
  if nX MOD 2 eq 0 then $ ;X is of even length.
    ;j = nX / 2
    ;Average Kth and K-1st medians.
    aa = MEDIAN(arr, /EVEN) $
  else aa = MEDIAN(arr)
  sum = 0
  absdev = 0

  ;for k = 0L, nX-1 do begin
  ;  d = y(k) - (b * x(k) + aa)
  ;  absdev = absdev + abs(d)
  ;  if y(k) ne 0 then d = d / abs(y(k))
  ;  if abs(d) gt eps then sum = sum + Sign(x(k), d)
  ;endfor

  d = y - (b * x + aa)
  absdev = TOTAL(abs(d), Double = Double)
  zeros = where(y ne 0)
  d = d[zeros] / abs(y[zeros])  
  zeros = where(abs(d) gt eps, nzeros)
  if nzeros ne 0 then $
    sum = TOTAL(x[zeros]*Sign(REPLICATE(1.0, nzeros), d[zeros]), $
                                                 Double = Double)
  if Double eq 0 then RETURN, FLOAT(sum) else RETURN, sum
END

FUNCTION MDfunc2, x, y, nX, Double = Double
  ON_ERROR, 2
  ss = nX + 0.0
  sx = TOTAL(x, Double = Double)
  sy = TOTAL(y, Double = Double)
  t = x - sx/ss
  st2 = TOTAL(t^2, Double = Double)
  b = TOTAL(t * y, Double = Double)
  b = b / st2
  a = (sy - sx * b) / ss
  sdeva = sqrt((1.0 + sx * sx / (ss * st2)) / ss)
  sdevb = sqrt(1.0 / st2)
  chisqr = TOTAL( (y - a - b * x)^2, Double = Double )
  prob = 1.0
  sdevdat = sqrt(chisqr / (nX-2))
  sdeva = sdeva * sdevdat
  sdevb = sdevb * sdevdat
  sig_ab = [sdeva, sdevb]
  if Double eq 0 then RETURN, FLOAT([a, b]) else RETURN, [a, b]
END

FUNCTION LadFit, x, y, absdev = absdev, Double = Double
  ON_ERROR, 2

  TypeX = SIZE(X)
  TypeY = SIZE(Y)
  nX = TypeX[TypeX[0]+2]
  nY = TypeY[TypeY[0]+2]

  if nX ne nY then $
    MESSAGE, "X and Y must be vectors of equal length."

  ;If the DOUBLE keyword is not set then the internal precision and
  ;result are identical to the type of input.
  if N_ELEMENTS(Double) eq 0 then $
    Double = (TypeX[TypeX[0]+1] eq 5 or TypeY[TypeY[0]+1] eq 5)

  sx = TOTAL(x, Double = Double)
  sy = TOTAL(y, Double = Double)
  sxy = TOTAL(x*y, Double = Double)
  sxx = TOTAL(x*x, Double = Double)
  del = nX * sxx - sx^2
  aa = (sxx * sy - sx * sxy) / del
  bb = (nX * sxy - sx * sy) / del
  chisqr = TOTAL((y - (aa + bb*x))^2, Double = Double)
  sigb = sqrt(chisqr / del)
  b1 = bb
  f1 = MDfunc(b1, x, y, arr, aa, absdev, nX, Double = Double)
  if f1 eq 0 then RETURN, MDfunc2(x, y, nX, Double = Double)
  b2 = bb + Sign(3.0 * sigb, f1) 
  f2 = MDfunc(b2, x, y, arr, aa, absdev, nX, Double = Double)
  while f1*f2 gt 0 do begin
    bb = 2.0 * b2 - b1
    b1 = b2
    f1 = f2
    b2 = bb
    f2 = MDfunc(b2, x, y, arr, aa, absdev, nX, Double = Double)
  endwhile
  sigb = 0.01 * sigb
  while abs(b2-b1) gt sigb do begin
    bb = 0.5 * (b1 + b2)
    if bb eq b1 or bb eq b2 then begin
      absdev = absdev / nX
      if Double eq 0 then RETURN, FLOAT([aa, bb]) else RETURN, [aa, bb]
    endif else begin
      f = MDfunc(bb, x, y, arr, aa, absdev, nX, Double = Double)
      if f*f1 ge 0 then begin
        f1 = f
        b1 = bb
      endif else begin
        f2 = f
        b2 = bb
      endelse
    endelse
  endwhile
  absdev = absdev / nX
  if Double eq 0 then RETURN, FLOAT([aa, bb]) else RETURN, [aa, bb]
END