;+ ; NAME: ; get_asymmetric_uncertainties ; ; PURPOSE: ; Given a distribution of parameter values, this routine returns the ; left/right uncertainties. ; ; CALLING SEQUENCE: ; errors = get_asymmetric_uncertainties(dst, num_sig=num_sig) ; ; INPUTS: ; dst -- distribution from which to extract uncertainties ; ; KEYWORD PARAMETERS: ; num_sig -- If the distribution is Gaussian about a mean value, ; this parameter represents the number of sigma away from that mean. ; sign -- If set, the routine will return a negative left error ; ; EXAMPLE: ; x = randomn(seed, 1001) ; errors = get_asymmetric_uncertainties(x, num_sig=1.) ; ;Should both be close to 1 ; print, errors ; ; MODIFICATION HISTORY: ; 2013 Jun 10 -- Written by Brian Jackson (decaelus@gmail.com) ; 2013 Jul 31 -- Dealt with the case that there are NOT elements to the left or ; right of the mean value by returning min/max values ; 2013 Jul 31 -- Added sign keyword ;- function get_asymmetric_uncertainties, dst, num_sig=num_sig, sign=sign if(~keyword_set(num_sig)) then num_sig = 1. ;The histogram value that corresponds to the number of sigma requested val = (gauss_pdf(num_sig) - gauss_pdf(-num_sig)) ;condition sign if(keyword_set(sign)) then begin sign = -1. end else sign = 1. men = mean(dst) ;left uncertainty ind = where(dst le men) ;cumulative histogram srt = sort(men - dst[ind]) srt_data = men - dst[ind[srt]] uniq_data = srt_data[uniq(srt_data)] left = sign*(men - min(dst)) ;Routine doesn't work if there's only one value in cum if(n_elements(uniq_data) gt 1) then begin cum = 0 for i = 0L, n_elements(uniq_data)-1 do \$ push, cum, double(n_elements(where(srt_data le uniq_data[i])))/\$ double(n_elements(srt_data)) left = sign*uniq_data[value_locate(cum, val)] end;if(n_elements(uniq_data) gt 1) ;right uncertainty ind = where(dst ge men) ;cumulative histogram srt = sort(dst[ind] - men) srt_data = dst[ind[srt]] - men uniq_data = srt_data[uniq(srt_data)] right = (max(dst) - men) ;Routine doesn't work if there's only one value in cum if(n_elements(uniq_data) gt 1) then begin cum = 0 for i = 0L, n_elements(uniq_data)-1 do \$ push, cum, double(n_elements(where(srt_data le uniq_data[i])))/\$ double(n_elements(srt_data)) right = uniq_data[value_locate(cum, val)] end;if(n_elements(uniq_data) gt 1) return, [left, right] end;function