FUNCTION BIWEIGHT_MEAN,Y,SIGMA, WEIGHTs ; ;+ ; NAME: ; BIWEIGHT_MEAN ; ; PURPOSE: ; Calculate the center and dispersion (like mean and sigma) of a ; distribution using bisquare weighting. ; ; CALLING SEQUENCE: ; Mean = BIWEIGHT_MEAN( Vector, [ Sigma, Weights ] ) ; ; INPUTS: ; Vector = Distribution in vector form ; ; OUTPUT: ; Mean - The location of the center. ; ; OPTIONAL OUTPUT ARGUMENTS: ; ; Sigma = An outlier-resistant measure of the dispersion about the ; center, analogous to the standard deviation. ; ; Weights = The weights applied to the data in the last iteration, ; floating point vector ; ; NOTES: ; Since a sample mean scaled by sigma/sqrt(N), has a Student's T ; distribution, the half-width of the 95% confidence interval for ; the sample mean can be determined as follows: ; ABS( T_CVF( .975, .7*(N-1) )*SIGMA/SQRT(N) ) ; where N = number of points, and 0.975 = 1 - (1 - 0.95)/2. ; PROCEDURES USED: ; ROBUST_SIGMA() ; REVISION HISTORY ; Written, H. Freudenreich, STX, 12/89 ; Modified 2/94, H.T.F.: use a biweighted standard deviation rather than ; median absolute deviation. ; Modified 2/94, H.T.F.: use the fractional change in SIGMA as the ; convergence criterion rather than the change in center/SIGMA. ; Modified May 2002 Use MEDIAN(/EVEN) ; Modified October 2002, Faster computation of weights ; Corrected documentation on 95% confidence interval of mean ; P.Broos/W. Landsman July 2003 ;- ON_ERROR,2 maxit = 20 ; Allow 20 iterations, this should nearly always be sufficient eps = 1.0e-24 n = n_elements(y) close_enough =.03*sqrt(.5/(n-1)) ; compare to fractional change in width diff = 1.0e30 itnum = 0 ; As an initial estimate of the center, use the median: y0=median(y,/even) ; Calculate the weights: dev = y-y0 sigma = ROBUST_SIGMA( dev ) if sigma lt EPS then begin ; The median is IT. Do we need the weights? if arg_present(weights) then begin ; Flag any value away from the median: limit=3.*sigma weights = float(abs(dev) LE limit) endif diff = 0. ; (skip rest of routine) endif ; Repeat: while( (diff gt close_enough) and (itnum lt maxit) )do begin itnum = itnum + 1 uu = ( (y-y0)/(6.*sigma) )^2 uu = uu < 1. weights=(1.-uu)^2 & weights=weights/total(weights) y0 = total( weights*y ) dev = y-y0 prev_sigma = sigma & sigma = robust_sigma( dev,/zero ) if sigma gt eps then diff=abs(prev_sigma-sigma)/prev_sigma else diff=0. endwhile return,y0 end