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