Viewing contents of file '../idllib/astron/contrib/freudenreich/boot_xyfit.pro'
FUNCTION BOOT_XYFIT, X,Y, N_SAMPLE,YINT,SLOP,SYINT,SSLOP,NUMOUT,ROBUST=whatever
;
;+
; NAME:
; BOOT_XYFIT
;
; PURPOSE:
; Bootstrap linear fit to data in which there are errors in both Y and X,
; but the measurement errors are not available. Calculates the bisector
; of the Y vs X and X vs Y regression.
;
; CALLING SEQUENCE:
; COEFF = BOOT_XYFIT( X,Y, NSAMPLE, YINT, SLOP, SYINT, SSLOP, NUMOUT,
; [ /ROBUST ] )
;
; INPUT ARGUMENT:
; X = x vector
; Y = y vector
; N_SAMPLE = the number of bootstrap samples
;
; RETURNS:
; COEFF = array of coefficients. Dimensions=(NDEG+1,N_SAMPLE)
; First dimension in the order A0, A1
;
; OPTIONAL OUTPUT ARGUMENT:
; YINT = Y intercept (average of COEFF(0,*))
; SLOP = slope (average of (COEFF(1,*))
; SYINT = standard error of Y intercept (std. dev. of COEFF(0,*))
; SSLOP = standard error of slope (std. dev. of COEFF(1,*))
; NUMOUT = actual number of samples returned. If there is an error in
; the fit to any sample, NUMOUT is decreased by one.
;
; OPTIONAL INPUT KEYWORD:
; ROBUST if present, an outlier-resistant fit is performed
;
; NOTES:
; This program randomly selects (x,y) points and fits a line to them. The
; line is the bisector of the Y vs X and X vs Y regressions. It does this
; N_SAMPLE times.
;
; WARNING:
; At least NDEG+1 points must be input. It is best to have at least twice
; that many for a robust fit.
;
; REVISION HISTORY:
; Written H.T. Freudenreich, HSTX, 3/16/94. Adapted from obsolete
; BOOT_LINE.
;-
IF KEYWORD_SET(WHATEVER) THEN ROBUST=1 ELSE ROBUST=0
N=N_ELEMENTS(X)
IF N LT 2 THEN BEGIN
PRINT,'BOOT_XYFIT: Too few points! Returning 0'
RETURN,0.
ENDIF
COEFF=FLTARR(2,N_SAMPLE)
; Get the random number seeds:
R0=SYSTIME(1)*2.+1.
SEEDS=RANDOMU(R0,N_SAMPLE)*1.0E6+1.
NGOOD = -1
FOR L=0,N_SAMPLE-1 DO BEGIN
R=SEEDS(L)
; Uniform random numbers between 0 and N-1:
PICK=RANDOMU(R,N)
R=FIX(PICK*N)
U=X(R) & V=Y(R)
IF ROBUST EQ 1 THEN BEGIN
CC=ROBUST_LINEFIT(U,V,/BISECT)
ENDIF ELSE BEGIN
CC1=POLY_FIT(U,V,1)
KK =POLY_FIT(V,U,1)
CC2=[-KK(0)/KK(1),1./KK(1)]
YYINT = CC1(0) & YSLOP = CC1(1)
XYINT = CC2(0) & XSLOP = CC2(1)
IF YSLOP GT XSLOP THEN BEGIN
A1=YYINT & B1=YSLOP & R1=SQRT(1.+YSLOP^2)
A2=XYINT & B2=XSLOP & R2=SQRT(1.+XSLOP^2)
ENDIF ELSE BEGIN
A2=YYINT & B2=YSLOP & R2=SQRT(1.+YSLOP^2)
A1=XYINT & B1=XSLOP & R1=SQRT(1.+XSLOP^2)
ENDELSE
YINT=(R1*A2+R2*A1)/(R1+R2)
SLOP=(R1*B2+R2*B1)/(R1+R2)
CC=[YINT,SLOP]
ENDELSE
IF (N_ELEMENTS(CC) EQ 2) THEN BEGIN
NGOOD=NGOOD+1
COEFF(*,NGOOD)=CC
ENDIF
ENDFOR
NGOOD = NGOOD+1
IF NGOOD LT N_SAMPLE THEN BEGIN
PRINT,'BOOT_XYFIT: Some samples rejected.'
COEFF=COEFF(*,0:NGOOD-1)
ENDIF
YINT=TOTAL(COEFF(0,*))/NGOOD
SLOP=TOTAL(COEFF(1,*))/NGOOD
SYINT=STDEV(COEFF(0,*))
SSLOP=STDEV(COEFF(1,*))
RETURN,COEFF
END