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