procedure wfractiles work workwt ndraws desired fractiles
*
* wfractiles returns a collection of fractiles of a weighted set of
* data.
*
* Syntax:
*
* @wfractiles work workwt ndraws desired fractiles
*
* Parameters:
* work series of sample values. Should be defined from entries 1 to ndraws
* workwt corresponding weights for the entries of work. This is altered within
* the procedure, so be sure to send a copy if you need it for other purposes.
* ndraws number of draws, ending entry of work and workwt
* desired VECTOR of percentiles that you want
* fractiles output VECTOR of computed fractiles
*
* Revision Schedule
* 11/02 Written by Estima
*
type series work workwt
type integer ndraws
type vector desired *fractiles
*
local integer lb ub trial
local integer n i
local real fractile
*
* Order the work, workwt pair on the values of work. Compute a normalized
* accumulation of the weights.
*
order work 1 ndraws workwt
acc(scale) workwt 1 ndraws
compute n=%rows(desired)
dim fractiles(n)
do i=1,n
compute fractile=desired(i)
*
* Bracket the desired fractile
*
compute lb=0,ub=ndraws+1
compute trial=fix(ndraws*fractile)
loop
if workwt(trial)=ndraws
compute ub=ndraws
else
if trial<=0
compute lb=1
if lb>0.and.ub<=ndraws
break
end loop
if workwt(lb)>=fractile {
compute fractiles(i)=work(lb)
next
}
if workwt(ub)<=fractile {
compute fractiles(i)=work(ub)
next
}
*
* Now pinch it by bisection, until (lb,ub) includes
* at most two points.
*
loop
compute trial=(lb+ub)/2
if workwt(trial)=ub-1
break
end loop
compute fractiles(i)=work(ub)
end do i
end