pro estimator_simple,a,p,m,low,high,plot=plot ; ;+ ; From a given vector a of measurements and a second ; vector p with relative probabilities for each measure, we build the ; probability distribution for the measured quantity. Then we output just ; three quantities the mean and the 2-sigma low and high errors. ; ; IN: a - fltarr measurements ; p - fltarr relative prob. for each measurement ; (normalization NOT required) ; binsize - float binsize for grouping the measurements ; ; OUT: m - float mean ; low- float 2-sigma error on the left () ; high-float 2-sigma error on the right side ; ; ; limit is the one-side Gaussian integrated probability to be reached ; for 1-sigma limits p=0.68 and limit=0.34, for 3-sigma limits p=0.997 and ; limit=0.499, for 2-sigma limits p=0.955 and limit=p/2. ; in general n-sigma=gauss_cvf((1.-p)/2) ; ; C. Allende Prieto, UT, Dec 2002 ;- limit=0.478 binsize=(max(a)-min(a))/20. h=histograma(a,x=x,binsize=binsize) prob=fltarr(n_elements(x)) for i=0,n_elements(x)-1 do begin www=where(a ge x(i)-binsize/2. and a lt x(i)+binsize/2.) if (max(www) gt -1) then prob(i)=total(p(www)) endfor prob=prob/total(prob) m=total(x*prob)/total(prob) hotpixel=median(where(abs(x-m) eq min(abs(x-m)))) leftarea=0.0 & rightarea=0.0 low=x(0) high=x(n_elements(x)-1) i=hotpixel while (i lt n_elements(prob) and high eq x(n_elements(x)-1)) do begin rightarea=rightarea+prob(i) if (rightarea ge limit) then high=x(i) i=i+1 endwhile i=hotpixel while (i gt -1 and low eq x(0)) do begin leftarea=leftarea+prob(i) if (leftarea ge limit) then low=x(i) i=i-1 endwhile if keyword_set(plot) then begin print,m,low,high plot,x,prob,psy=10,xr=[min(x)*0.9,max(x)*1.1],yr=[0.,max(prob)*1.1],charsi=2.0 oplot,[m,m],[0,1],thick=2 oplot,[low,low],[0,1],linestyle=2 oplot,[high,high],[0,1],linestyle=2 endif end