;+
; NAME:
;  madscale
; PURPOSE: (one line)
;  minimize the average deviation (MAD) for y = a * x.
; DESCRIPTION:
;  This routine estimates the scaling factor, a, such that
;  the average deviation, total(abs(y-a*x) ), is minimized.
;
;  The statistics are robust in that the result is insensitive
;  to outliers.  
;
;  Note that for y[i] = a * x[i], x can be thought of as
;  either a coordinate, or a template (e.g., a line profile).
;
; CATEGORY:
;  Statistics
; CALLING SEQUENCE:
;  a = madscale(x, y)
; INPUTS:
;  x   - dependent variable.
;  y   - data, proportional to x.
; OPTIONAL INPUT PARAMETERS:
; INPUT KEYWORD PARAMETERS:
; OUTPUT KEYWORD PARAMETERS:
; OUTPUTS:
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  None.
; RESTRICTIONS:
;  None.
; PROCEDURE:
;
; The goal is to minimize f = total(abs(y-a*x)) = total( abs(x) * abs(y/x - a) )
; This is equivalent to finding the weighted median of y/x, with
; abs(x) as the weights.
;
; MODIFICATION HISTORY:
;  Written by Leslie A. Young, Soutwest Research Institute, 2002 Oct 14.
;-
function madscale, x, y
  
  gd = where(x ne 0, ngd)
  if ngd eq 0 then return, 0
  return, medianw(y[gd]/x[gd], abs(x[gd]))
  
end

pro madscaleTEST

  x = [-2.,-1.,0.,1.,2.]
  a = -1
  y = a * x
  yerr = [0.432511,    -0.690602,     -1.03700,      2.23537,     0.505393]*0.2
  ydat=y+yerr
  ydat[1] = 3.5
  mada = madscale(x,ydat)
  lsa = total(x*ydat)/total(x^2)
  print, a, lsa, mada
  plot, x, y, xr=[-3,3], yr=[-4,4]
  oplot, x, ydat, ps=4
  oplot, x, lsa*x, line=2
  oplot, x, mada*x, line=2, thick=2
end