fitexy.pro

basic median filter simulation

;+
; 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
;
; 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
;-
function chisq_fitexy, B_angle
;
; NAME:
;       chisq_fitexy
; PURPOSE:
;       Function minimized by fitexy  (computes chi-square of linear fit).
;       It is called by minimization procedures during execution of fitexy.
; CALLING SEQUENCE:
;               chisq = chisq_fitexy( B_angle )
; INPUTS:
;       B_angle = arc-tangent of B_slope of linear fit.
; OUTPUTS:
;       Result of function = chi_square - offs  (offs is in COMMON).
; COMMON:
;       common fitexy, communicates the data from pro fitexy.
; PROCEDURE:
;       From "Numerical Recipes" column: Computer in Physics Vol.6 No.3
; MODIFICATION HISTORY:
;       Written, Frank Varosi NASA/GSFC 1992.

  common fitexy, xx, yy, sigx, sigy, ww, Ai, offs

        B_slope = tan( B_angle )
        ww = 1/( ( (B_slope * sigx)^2 + sigy^2 ) > 1.e-30 )
        if N_elements( ww ) EQ 1 then sumw = ww * N_elements( xx ) $
                                 else sumw = total( ww )
        y_Bx = yy - B_slope * xx
        Ai = total( ww * y_Bx )/sumw

return, total( ww * (y_Bx - Ai)^2 ) - offs
end
;-------------------------------------------------------------------------------
pro fitexy, x, y, A_intercept, B_slope, sigma_A_B, chi_sq, q, TOLERANCE=Tol, $
                                        X_SIGMA=x_sigma, Y_SIGMA=y_sigma
  compile_opt idl2					
  common fitexy, xx, yy, sigx, sigy, ww, Ai, offs

  if N_params() LT 4 then begin
     print,'Syntax -  fitexy, x, y, A, B, X_SIG=sigx, Y_SIG=sigy,' 
     print,'                  [sigma_A_B, chi_sq, q, TOLERANCE = ]'
     return
  endif
  
; Normalize data before running fitexy

  xm = (MOMENT(x, SDEV = xs, /DOUBLE))[0]
  ym = (MOMENT(y, SDEV = ys, /DOUBLE))[0]
  xx = (x - xm) / xs
  yy = (y - ym) / ys
  sigx = x_sigma / xs
  sigy = y_sigma / ys
 
   
;Compute first guess for B_slope using standard 1-D Linear Least-squares fit,
; where the non-linear term involving errors in x are ignored.
; (note that Tx is a transform to reduce roundoff errors)

        ww = sigx^2 + sigy^2
        if N_elements( ww ) EQ 1 then sumw = ww * N_elements( xx ) $
                                 else sumw = total( ww )
        Sx = total( xx * ww )
        Tx = xx - Sx/sumw
        B = total( ww * yy * Tx ) / total( ww * Tx^2 )

;Find the minimum chi-sq while including the non-linear term (B * sigx)^2
; involving variance in x data (computed by function chisq_fitexy):
; using minf_bracket (=MNBRAK) and minf_parabolic (=BRENT)
        offs = 0
        ang = [ 0, atan( B ), 1.571 ]
        chi = fltarr( 3 )
        for j=0,2 do chi[j] = chisq_fitexy( ang[j] )    ;this is for later...
        if N_elements( Tol ) NE 1 then Tol=1.e-3
        a0 = ang[0]
        a1 = ang[1]
        minf_bracket, a0,a1,a2, c0,c1,c2, FUNC="chisq_fitexy"
        minf_parabolic, a0,a1,a2, Bang, chi_sq, FUNC="chisq_fitexy", TOL=Tol

        if N_params() EQ 7 then q = 1 - chisqr_pdf( chi_sq, N_elements(x) - 2 )
        A_intercept = Ai        ;computed in function chisq_fitexy
        ang = [a0,a1,a2,ang]
        chi = [c0,c1,c2,chi]

;Now compute the variances of estimated parameters,
; by finding roots of ( (chi_sq + 1) - chisq_fitexy ).
;Note: ww, Ai are computed in function chisq_fitexy.

        offs = chi_sq + 1
        wc = where( chi GT offs, nc )

        if (nc GT 0) then begin

                angw = [ang[wc]]
                d1 = abs( angw - Bang ) MOD !PI
                d2 = !PI - d1
                wa = where( angw LT Bang, na )

                if (na GT 0) then begin
                        d = d1[wa]
                        d1[wa] = d2[wa]
                        d2[wa] = d
                   endif

                Bmax = zbrent( Bang,Bang+max(d1),F="chisq_fitexy",T=Tol ) -Bang
                Amax = Ai - A_intercept
                Bmin = zbrent( Bang,Bang-min(d2),F="chisq_fitexy",T=Tol ) -Bang
                Amin = Ai - A_intercept

                if N_elements( ww ) EQ 1 then r2 = 2/( ww * N_elements( x ) ) $
                                         else r2 = 2/total( ww )

                sigma_A_B = [ Amin^2 + Amax^2 + r2 , Bmin^2 + Bmax^2 ]
                sig_A_B = sqrt( sigma_A_B/2 ) / ([1,cos(Bang)^2])

          endif 

;Finally, transform parameters back to orignal units.


        B_slope = tan( Bang ) *ys /xs
        A_intercept = A_intercept*ys - tan(Bang) * ys / xs *xm + ym
        if Nc GT 0 then sigma_A_B = [SQRT( (sig_A_B[0] * ys)^2 +  $
                    (sig_A_B[1] * ys / xs * xm)^2 ), sig_A_B[1] * ys / xs] $
                else sigma_A_B = [1.e33,1.e33]    

return
end