robust_mean.pro

function robust_mean,Y,CUT,Sigma,Num_Rej,GoodInd=GoodInd
;+
; NAME:
;    Robust_Mean 
;
; PURPOSE:
;    Outlier-resistant determination of the mean and standard deviation.
;
; EXPLANATION:
;    Robust_Mean trims away outliers using the median and the median
;    absolute deviation.    An approximation formula is used to correct for
;    the trunction caused by trimming away outliers
;
; CALLING SEQUENCE:
;    mean = Robust_Mean( VECTOR, Sigma_CUT, Sigma_Mean, Num_RejECTED)
;
; INPUT ARGUMENT:
;       VECTOR    = Vector to average
;       Sigma_CUT = Data more than this number of standard deviations from the
;               median is ignored. Suggested values: 2.0 and up.
;
; OUTPUT ARGUMENT:
;       Mean  = the mean of the input vector, numeric scalar
;
; KEYWORDS:
;
;       GoodInd = The indices of the values not rejected
;
; OPTIONAL OUTPUTS:
;Sigma_Mean = the approximate standard deviation of the mean, numeric
;            scalar.  This is the Sigma of the distribution divided by sqrt(N-1)
;            where N is the number of unrejected points. The larger
;            SIGMA_CUT, the more accurate. It will tend to underestimate the
;            true uncertainty of the mean, and this may become significant for
;            cuts of 2.0 or less.
;       Num_RejECTED = the number of points trimmed, integer scalar
;
; EXAMPLE:
;       IDL> a = randomn(seed, 10000)    ;Normal distribution with 10000 pts
;       IDL> Robust_Mean,a, 3, mean, meansig, num    ;3 Sigma clipping   
;       IDL> print, mean, meansig,num
;
;       The mean should be near 0, and meansig should be near 0.01 ( =
;        1/sqrt(10000) ).    
; PROCEDURES USED:
;       AVG() - compute simple mean
; REVISION HISTORY:
;       Written, H. Freudenreich, STX, 1989; Second iteration added 5/91.
;       Use MEDIAN(/EVEN)    W. Landsman   April 2002
;       Correct conditional test, higher order truncation correction formula
;                R. Arendt/W. Landsman   June 2002
;       New truncation formula for sigma H. Freudenriech  July 2002
;-  
On_Error,2
 if N_params() LT 2 then begin
     print,'Syntax - Robust_Mean(Vector, Sigma_cut, [ Sigma_mean, '
     print,'                                  Num_Rejected ])'
     return,1
 endif  
 Npts    = N_Elements(Y)
 YMed    = MEDIAN(Y,/EVEN)
 AbsDev  = ABS(Y-YMed)
 MedAbsDev = MEDIAN(AbsDev,/EVEN)/0.6745
 IF MedAbsDev LT 1.0E-24 THEN MedAbsDev = AVG(AbsDev)/.8  
 Cutoff    = Cut*MedAbsDev  
 GoodInd = WHERE( AbsDev LE Cutoff, Num_Good )
 GoodPts = Y[ GoodInd ]
 Mean    = AVG( GoodPts )
 Sigma   = SQRT( TOTAL((GoodPts-Mean)^2)/Num_Good )
 Num_Rej = Npts - Num_Good ; Compenate Sigma for truncation (formula by HF):
 SC = Cut > 1.0
 IF SC LE 4.50 THEN $
    SIGMA=SIGMA/(-0.15405+0.90723*SC-0.23584*SC^2+0.020142*SC^3)  
 Cutoff = Cut*Sigma  
 GoodInd = WHERE( AbsDev LE Cutoff, Num_Good )
 GoodPts = Y[ GoodInd ]
 mean    = AVG( GoodPts )
 Sigma   = SQRT( TOTAL((GoodPts-mean)^2)/Num_Good )
 Num_Rej = Npts - Num_Good ; Fixed bug (should check for SC not Sigma) & add higher order correction
 SC = Cut > 1.0
 IF SC LE 4.50 THEN $
    SIGMA=SIGMA/(-0.15405+0.90723*SC-0.23584*SC^2+0.020142*SC^3) 
; Now the standard deviation of the mean:
 Sigma = Sigma/SQRT(Npts-1.)  
return,mean
 END