subroutine apteqxn (x, pl, pr, nbins, nep, xep, pmin, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTEQXN
c
c call apteqxn (x, pl, pr, nbins, nep, xep, pmin, nerr)
c
c Version: apteqxn Updated 1990 August 13 13:50.
c apteqxn Originated 1990 August 13 13:50.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the nep equal probability (EP) bins, each bounded by
c xep(n) on the left and xep(n+1) on the right, corresponding to
c the nbins probability distribution function (PDF) bins, each
c with relative probability pl(n) at x(n) on the left, and
c relative probability pr(n) at x(n+1) on the right, varying
c linearly. For a histogram PDF, pl(n) = pr(n). Also, to find
c the minimum normalized probability pmin of any PDF bin.
c Flag nerr indicates any input error.
c
c Random sampling of the random variable x is done as follows:
c
c n = 1 + ranf( ) * nep
c x = xep(n) + ranf( ) * (xep(n+1) - xep(n))
c
c Input: x, pl, pr, nbins, nep.
c
c Output: xep, pmin, nerr.
c
c Glossary:
c
c nbins Input Number of input PDF bins or intervals.
c
c nep Output Number of output EP bins.
c
c nerr Output Indicates an input error, if not 0.
c 1 if nbins is not positive.
c 2 if nep is not positive.
c 3 if any pl(n) or pr(n) is negative.
c 4 if any x(n+1) is less than x(n).
c
c pl Input For the PDF bin n, the unnormalized relative
c probability of x(n), per unit x interval.
c Size nbins. May not be negative for any n.
c If the PDF is a histogram, pl(n) = pr(n).
c If the PDF is continuous, pl(n) = pr(n-1), for
c n = 2, nbins.
c
c pmin Output The minimum normalized probability of any of the
c nbins PDF bins. The reciprocal of pmin is the
c minimum value of nep needed to guarantee that no
c EP bin contains more than one PDF bin boundary.
c
c pr Input For the PDF bin n, the unnormalized relative
c probability of x(n+1), per unit x interval.
c Size nbins. May not be negative for any n.
c If the PDF is a histogram, pr(n) = pl(n).
c If the PDF is continuous, pr(n) = pl(n+1), for
c n = 1, nbins - 1.
c
c x Input The input PDF bin with index n is bounded at the left
c by x(n) and at the right by x(n+1). Size nbins + 1.
c For any n, x(n) must not be greater than x(n+1).
c
c xep Output The output EP bin with index n is bounded at the left
c by xep(n) and at the right by xep(n+1).
c Size nep + 1.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Left-hand probability in PDF bin.
dimension pl (1)
c---- Right-hand probability in PDF bin.
dimension pr (1)
c---- Bounding value of x in PDF bin.
dimension x (1)
c---- Bounding value of x in EP bin.
dimension xep (1)
cbugc***DEBUG begins.
cbugc---- Final PDF bin used by EP bin.
cbug dimension nbep (1000)
cbugc***DEBUG ends.
c.... Local variables.
c---- Probability increment in EP bin.
common /lapteqxn/ dpcum
c---- A very small number.
common /lapteqxn/ fuz
c---- Index of PDF or EP bin.
common /lapteqxn/ n
c---- Index of a PDF bin.
common /lapteqxn/ nbin
c---- Relative probability of a PDF bin.
common /lapteqxn/ pbin
c---- Cumulative relative probability.
common /lapteqxn/ pcum
c---- Cum rel prob at left of PDF bin.
common /lapteqxn/ pcuml
c---- Cum rel prob at left of EP bin.
common /lapteqxn/ pcumle
c---- Maximum of pcuml, pcumle.
common /lapteqxn/ pcumll
c---- Cum rel prob at right of PDF bin.
common /lapteqxn/ pcumr
c---- Cum rel prob at right of EP bin.
common /lapteqxn/ pcumre
c---- Relative probability at xll.
common /lapteqxn/ pll
c---- Slope (pr - pl) / (xR - xL).
common /lapteqxn/ slope
c---- Maximum of left-side x, xep.
common /lapteqxn/ xll
cbugc***DEBUG begins.
cbugc---- Minimum nep needed.
cbug common /lapteqxn/ nepneed
cbug 9901 format (/ 'apteqxn finding equal-probability bins.',
cbug & ' nep=',i5)
cbug 9902 format (i5,' xL,xR=',1p2e22.14 / ' pl,pr=',1p2e22.14)
cbug write ( 3, 9901) nep
cbug write ( 3, 9902) (n, x(n), x(n+1), pl(n), pr(n), n = 1, nbins)
cbugc***DEBUG ends.
c.... initialize.
c---- A very small number.
fuz = 1.e-99
nerr = 0
c.... Test for errors in scalar input.
if (nbins .le. 0) then
nerr = 1
go to 210
endif
if (nep .le. 0) then
nerr = 2
go to 210
endif
c.... Test for errors in input arrays, and find total unnormalized probability.
pcum = 0.0
c---- Loop over input PDF bins.
do 120 n = 1, nbins
if ((pl(n) .lt. 0.0) .or. (pr(n) .lt. 0.0)) then
nerr = 3
go to 210
endif
if (x(n+1) .lt. x(n)) then
nerr = 4
go to 210
endif
pbin = 0.5 * (pl(n) + pr(n)) * (x(n+1) - x(n))
pcum = pcum + pbin
if (n .eq. 1) then
pmin = pbin
else
pmin = amin1 (pmin, pbin)
endif
c---- End of loop over input PDF bins.
120 continue
pmin = pmin / pcum
c.... Integrate over input PDF bins to find output EP bins.
nbin = 0
pcumr = 0.0
pcumre = 0.0
xep(1) = x(1)
c---- Loop over output EP bins.
do 140 n = 1, nep
pcumle = pcumre
pcumre = n * pcum / nep
c---- Add one or more input PDF bins.
130 if (pcumre .gt. pcumr) then
c---- Reached end of input PDF table.
if (nbin .eq. nbins) then
xep(n+1) = x(nbins+1)
cbugc***DEBUG begins.
cbug nbep(n) = nbins
cbugc***DEBUG ends.
go to 140
endif
nbin = nbin + 1
pcuml = pcumr
pbin = 0.5 * (pl(nbin) + pr(nbin)) * (x(nbin+1) - x(nbin))
pcumr = pcuml + pbin
go to 130
c---- Tested current input PDF bin.
endif
c.... Right side of output EP bin is in this input PDF bin.
slope = (pr(nbin) - pl(nbin)) / (x(nbin+1) - x(nbin) + fuz)
pcumll = amax1 (pcuml, pcumle)
dpcum = pcumre - pcumll
xll = amax1 (x(nbin), xep(n))
pll = pl(nbin) + slope * (xll - x(nbin))
xep(n+1) = xll + 2.0 * dpcum /
& (pll + sqrt (pll**2 + 2.0 * slope * dpcum))
cbugc***DEBUG begins.
cbug nbep(n) = nbin
cbugc###fx01 format (' n= ',i5,' xepL= ',1pe22.14,' xepR= ',1pe22.14 /
cbugc### & 13x,'pcumle=',1pe22.14,' pcumre=',1pe22.14 /
cbugc### & 13x,'dpcum= ',1pe22.14,' xll= ',1pe22.14)
cbugc###fx02 format (' nbin=',i5,' xL= ',1pe22.14,' xR= ',1pe22.14 /
cbugc### & 13x,'pcuml= ',1pe22.14,' pcumr= ',1pe22.14)
cbugc### write ( 3, fx01) n, xep(n), xep(n+1), pcumle, pcumre, dpcum, xll
cbugc### write ( 3, fx02) nbin, x(nbin), x(nbin+1), pcuml, pcumr
cbugc***DEBUG ends.
c---- End of loop over output EP bins.
140 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'apteqxn results. pmin=',1pe22.14,' nepneed=',i5)
cbug 9904 format (i5,' xepL=',1pe22.14)
cbug 9905 format (i5,' xepR=',1pe22.14,' pcumre=',1pe22.14)
cbug nepneed = 1.0 / (pmin + fuz)
cbug write ( 3, 9903) pmin, nepneed
cbug n = 1
cbug write ( 3, 9904) n, xep(n)
cbug do 170 n = 1, nep
cbug pcumre = n * pcum / nep
cbug write ( 3, 9905) n, xep(n+1), pcumre
cbug 170 continue
cbugc***DEBUG ends.
210 return
c.... End of subroutine apteqxn. (+1 line.)
end
UCRL-WEB-209832