Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/ChaMPlane/quantile/pre-release/quantile.pro Дата изменения: Wed Jul 26 06:40:36 2006 Дата индексирования: Tue Oct 2 08:21:26 2012 Кодировка: Поисковые слова: m 5

; JaeSub Hong, 2003-2005, version 1.7
; Please report any problem or suggestion at jaesub@head.cfa.harvard.edu
;
; Calculate quantiles from a distribution
; Refer to
; J. Hong, E.M. Schlegel & J.E. Grindlay,
; 2004 ApJ 614 p508 and references therein
;
; required routines : interpol.pro and ibeta.pro
; they are noramlly included in the standard idl distribution
;
; usage
; qt = quantile(frac, src, range=range, $
; bkg, ratio=ratio, err_Ex=err_Ex)
; examples
; range=[0.3,8.0]
; src=randomu(seed, 100) * (range[1]-range[0]) +range[0]
; qt = quantile(frac, src, range=range, err_Ex=err_Ex)
;
; bkg_on = randomu(seed, 20) * (range[1]-range[0]) +range[0]
; bkg = randomu(seed, 100) * (range[1]-range[0]) +range[0]
; src = [src,bkg_on] ; src and bkg photons in the src region
; qt = quantile(frac, src, range=range, bkg, ratio=0.2, err_Ex=err_Ex)
;
; required input:
; frac: list of quantile fractions
; default: 0.25,0.33,0.5,0.67,0.75
; src : value lists in the source region
; (e.g. energies of photons in the source region)
; range : the full range of values
; required input for bkgnd subtraction:
; bkg : value lists in the bkgnd region
; (e.g. energies of photons in the bkgnd region)
; should be sorted in a increasing order
; ratio : ration of the source to bkgnd region
; optional input
; nip : number of interplation points for bkgnd subtraction
; default 2000
; /nofixerror : unless requested, the following correction
; for error estimation is automatically applied
; 1. when the error estimation fails, it normally returns
; -1 for the error, but we can set the error to
; the whole range, which is reasonable.
; 2. when net<100, the error could be overestimated
; as much as 20-30% (See Fig.9 in Hong et al 2004)
; a relatively moderate correction
; x1./(1.+0.5/net)
; is applied to correct this.
; 3. correction for error due to bkgnd is applied
; sqrt(1+Nbkg/Nsrc) = sqrt(1.+(src-net)/net) (See Fig.9)
; 4. Error for QDy is likely overestimated
; due to correlation of Q25/Q75,
; so x0.8 applied to compensate this
; Refer to Hong et al 2004
; /nosort: src and bkg will be sorted in ascending order, but
; they are already sorted, you can skip the sorting procedure
;
; return value and optional output
; return value : quantiles Ex%
; err_Ex : error estimation of Ex%
; Qx : quantiles Qx = (Ex% - Elo) / (Eup - Elo)
; err_Qx : error restimation of Qx
; QDx : QCCD x value
; QDy : QCCD y value
; err_QDx : error of QDx
; err_QDy : error of QDy
;
; required routines
; interpol.pro
; ibeta.pro : use version 1.21 or higher
; to-do
; include Harrell-Davis tech
; routines to generate grid pattern

;------------------------------------------------------------------------
function qt_os, frac, values, range=range
; quantile estimation by order statistics

n_values=n_elements(values)
i = findgen(n_values)
ifrac = (i*2.+1.)/2./n_values

; force the boundary condition
values_ = [range[0],values,range[1]]
ifrac_ = [0.0,ifrac,1.0]
ans =interpol(values_,ifrac_, frac)
return, ans
end

function qt_err_mj, frac, values, range=range
; Error estimation by Maritz-Jarrett method

n=n_elements(values)
n_frac=n_elements(frac)

m=frac*n+0.5
a = m-1.
b = n-m

; ibeta is computational demanding in idl
i_src = findgen(n+1.)/n
ans=frac*0.0
for i=0,n_frac-1 do begin
if (a[i] gt 0.0) and (b[i] gt 0.0) then begin
beta = ibeta(a[i],b[i],i_src)
; print,i,beta
w=beta[1:*]-beta[0:n-1]
c = w*values
c1 = total(c)
c2 = total(c*values)
c0 = c2-c1^2
if c0 ge 0.0 then ans[i]=sqrt(c0) $
else ans[i]=-1.
endif else begin
ans[i]=-1.
endelse
endfor

return,ans
end

function qt_fix_err, Ex, error, n_net, n_src, range=range

maxerror = (range[1]-Ex)>(Ex-range[0])

w=where(error gt 0.0 and error lt maxerror,cw)
if cw gt 0 then $
error[w]= error[w] / (1.+0.5/(n_net>1.)) $
* sqrt(1.+(n_src-n_net)/n_net)

w=where(error le 0.0,cw)
if cw gt 0 then error[w] = maxerror[w]
w=where(error ge maxerror,cw)
if cw gt 0 then error[w] = maxerror[w]
return, error
end

;------------------------------------------------------------------------
function qt_simple, frac, values, range=range, $
err_Ex=err_Ex
ans=qt_os(frac, values, range=range)
err_Ex=qt_err_mj(frac,values,range=range)
return, ans
end

pro qt_qccd, range, $
frac=frac, $
Ex=Ex, err_Ex=err_Ex, $
Qx=Qx, err_Qx=err_Qx, $
QDx=QDx, QDy=QDy, $
err_QDx=err_QDx, err_QDy=err_QDy, $
nofixerror=nofixerror

rl = range[1]-range[0]
Qx = (Ex-range[0])/rl
err_Qx = err_Ex/rl

QDx = alog10(Qx[1]/(1.-Qx[1]))
m_l = Qx[1] - err_Qx[1]
m_u = Qx[1] + err_Qx[1]

bigN = 1.e3
; big number for divergence, is this big enough or too big

err_QDx_l = bigN
err_QDx_u = bigN
if (err_Qx[1] gt 0.0) then begin
if (m_l gt 0.0) then $
err_QDx_l = QDx-alog10(m_l/(1.-m_l))
if (m_u lt 1.0) then $
err_QDx_u = alog10(m_u/(1.-m_u))-QDx
endif
err_QDx = [-err_QDx_l, err_QDx_u]

if (Qx[2] le 0.0) then QDy = 3.0 $
else QDy = 3.*Qx[0]/Qx[2]
if (Qx[0] le 0.0 or err_Qx[0] lt 0.0 $
or Err_Qx[2] lt 0.0) then err_QDy = 3. $
else err_QDy = QDy * sqrt((err_Qx[0]/Qx[0])^2. +(err_Qx[2]/Qx[2])^2.);

; 0.8 is correction factor #3.
if not keyword_set(nofixerror) then ec3 = 0.8 $
else ec3 = 1.0

if (QDy - err_QDy le 0.0) then err_QDy_l = QDy $
else err_QDy_l = ec3 * err_QDy
if (QDy + err_QDy ge 3.0) then err_QDy_u = 3.-QDy $
else err_QDy_u = ec3 * err_QDy

err_QDy = [-err_QDy_l, err_QDy_u]

frac=frac[3:*]
Ex = Ex[3:*]
err_Ex = err_Ex[3:*]
Qx = Qx[3:*]
err_Qx = err_Qx[3:*]


end

function quantile, frac, src, range=range, $
bkg, ratio=ratio, $
err_Ex=err_Ex, $
Qx=Qx, err_Qx=err_Qx, $
QDx=QDx, QDy=QDy, $
err_QDx=err_QDx, err_QDy=err_QDy, $
Nip=Nip, $
nosort=nosort, $
nofixerror=nofixerror

; src_ph and bkg_ph should be sorted

if not keyword_set(nosort) then begin
order=sort(src)
src=src[order]
endif

frac=[0.25,0.50,0.75,frac]

if not keyword_set(bkg) then begin
ans = qt_simple(frac, src, range=range, err_Ex=err_Ex)
n_src = n_elements(src)
n_net = n_src
if not keyword_set(nofixerror) then $
err_Ex=qt_fix_err(ans,err_Ex,n_net,n_src,range=range)

qt_qccd, range, frac=frac, Ex=ans, err_Ex=err_Ex, $
Qx=Qx, err_Qx=err_Qx, $
QDx=QDx, QDy=QDy, $
err_QDx=err_QDx, err_QDy=err_QDy, $
nofixerror=nofixerror

return,ans
endif

if not keyword_set(nosort) then begin
order=sort(bkg)
bkg=bkg[order]
endif

n_src = n_elements(src)
n_bkg = n_elements(bkg)
if 0 eq n_elements(ratio) then ratio = 1.0

n_net = (n_src-ratio*n_bkg)>0

if n_net le 0.5 then begin
ans=interpol(range, [0.,1.0], frac)
err_Ex = (ans-range[0])>(range[1]-ans)
qt_qccd, range, frac=frac, Ex=ans, err_Ex=err_Ex, $
Qx=Qx, err_Qx=err_Qx, $
QDx=QDx, QDy=QDy, $
err_QDx=err_QDx, err_QDy=err_QDy, $
nofixerror=nofixerror

return, ans
endif

if not keyword_set(Nip) then Nip=2000 ; interpolation points
Narr = findgen(Nip)
ifrac=(Narr+0.5)/Nip
iE = Narr/Nip*(range[1]-range[0])+range[0]

iqt_src = qt_os(ifrac, src, range=range)
ic_src=interpol(n_src*ifrac, iqt_src, iE)

iqt_bkg = qt_os(ifrac, bkg, range=range)
ic_bkg=interpol(n_bkg*ifrac, iqt_bkg, iE)

; forward
ic_src = (ic_src>0.0) ic_bkg = (ic_bkg>0.0) ic_net = ic_src - ic_bkg*ratio
ic_net = (ic_net>0.0) for i=1, Nip-1 do $
if ic_net[i] lt ic_net[i-1] then ic_net[i]=ic_net[i-1]

;backward
ic_src_ = n_src-ic_src
ic_bkg_ = n_bkg-ic_bkg
ic_net_ = ic_src_ - ic_bkg_*ratio
ic_net_ = (ic_net_>0.0) for i=Nip-2,0,-1 do $
if ic_net_[i] lt ic_net_[i+1] then ic_net_[i]=ic_net_[i+1]

; average forward and backword
net_frac = (ic_net-ic_net_+n_net)/2./n_net
ans = interpol(iE, net_frac, frac)

w=where(finite(ans),cw)
n_frac=n_elements(frac)
if cw ne n_frac then print,'warning : some elements are not finite'

; now in order to estimate the error
; regenerate photons based on distribution
n_net_ = long(n_net+0.5)
;print,n_net,n_net_
tqt = (findgen(n_net_)*2+1.)/n_net_/2.
ip_src_ph = interpol(iE, net_frac, tqt)
err_Ex = qt_err_mj(frac, ip_src_ph, range=range)

if not keyword_set(nofixerror) then $
err_Ex=qt_fix_err(ans,err_Ex,n_net,n_src,range=range)
qt_qccd, range, frac=frac, Ex=ans, err_Ex=err_Ex, $
Qx=Qx, err_Qx=err_Qx, $
QDx=QDx, QDy=QDy, $
err_QDx=err_QDx, err_QDy=err_QDy, $
nofixerror=nofixerror

return, ans
end
;------------------------------------------------------------------------