subroutine apteqin (noptx, xl, xr, pl, pr, nbins, nep, nbep,
& pmin, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTEQIN
c
c call apteqin (noptx, xl, xr, pl, pr, nbins, nep, nbep,
c & pmin, nerr)
c
c Version: apteqin Updated 1990 August 13 13:50.
c apteqin Originated 1990 August 7 10:00.
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 corresponding to a
c piecewise linear tabulated probability distribution function
c (PDF). The PDF is defined by either of two methods:
c
c For noptx = 0, the PDF is a histogram with nbins bins.
c In each bin, pl(n) = pr(n), and xl and xr are not used.
c The random variable to be sampled is the bin index nbep.
c Random sampling of the initial bin index n is done as follows:
c
c nn = 1 + ranf( ) * nep
c n = nbep(nn)
c
c For noptx = 1, the PDF is a piecewise linear tabulated function
c of p(x) vs x, with nbins x intervals, where x is the random
c variable to be sampled. In the n'th interval, the relative
c probability varies linearly from pl(n) at xl(n) on the left,
c to pr(n) at xr(n) on the right. For a histogram PDF,
c pl(n) = pr(n). For a continuous PDF, pl(n+1) = pr(n) and
c xl(n+1) = xr(n). There may be gaps in x, but no overlaps.
c For each of the nep EP bins, nbep will be calculated.
c
c Also, to find the minimum normalized probability pmin of any
c PDF bin.
c Flag nerr indicates any input error.
c
c Random sampling of the random variable x is done as follows:
c
c nn = 1 + ranf( ) * nep
c n = nbep(nn)
c Randomly sample a value of x in the range from xl(n)
c to xr(n), based on the linear probability distribution
c from pl(n) to pr(n) (uniform if pl(n) = pr(n)).
c
c Note: Use apteqxn to find the x values bounding a set of EP bins.
c Use apteqsb to sample bin numbers from an EP distribution
c of bin numbers. Use aptalsh (uniform) or aptalsl (linear) to
c sample x values from the sampled bin numbers.
c Use apteqxs to sample x values from an EP distribution of
c x intervals.
c
c Input: noptx, xl, xr, pl, pr, nbins, nep.
c
c Output: nbep, pmin, nerr.
c
c Glossary:
c
c nbep Output The initial PDF bin index corresponding to an
c EP bin. Size nep.
c
c nbins Input Number of initial PDF bins or intervals.
c Size of arrays pl, pr, and if noptx = 1,
c xl and xr.
c
c nep Output Number of final EP bins. Size of array nbeq.
c
c nerr Output Indicates an input error, if not 0.
c 1 if nbins is not positive.
c 2 if noptx is not 0 or 1.
c 3 if any pl(n) or pr(n) is negative.
c 4 if any xr(n) is less than xl(n).
c 5 if any xr(n) is greater than xl(n+1).
c 6 if nep is not positive.
c
c noptx Input Indicates the type of PDF:
c 0 if the bin index n is the random variable,
c xl and xr are not used, and pl(n) = pr(n)
c (a histogram).
c 1 if xl(n) and xr(n) are the values of the random
c variable with relative probabilities pl(n) and
c pr(n), respectively.
c
c pl Input For initial bin n, the unnormalized relative
c probability of the bin (noptx = 0), or of xl(n),
c per unit x interval (noptx = 1). Size nbins.
c 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 initial bin n, the unnormalized relative
c probability of the bin (noptx = 0), or of xr(n),
c per unit x interval (noptx = 1). Size nbins.
c 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 xl Input For initial bin n, the value of the random variable
c x at the left boundary of the bin (noptx = 1).
c For any n, xl(n) must not be greater than xr(n),
c and must be equal or greater than xr(n-1).
c If the PDF includes a continuous range of x values,
c in monotonically increasing order, with no gaps, then
c xl(n) = xr(n-1), for n = 2, nbins.
c Size nbins, if noptx = 1. Otherwise, not used.
c
c xr Input For initial bin n, the value of the random variable
c x at the right boundary of the bin (noptx = 1).
c For any n, xr(n) must not be less than xl(n),
c and must be equal or less than xl(n+1).
c If the PDF includes a continuous range of x values,
c in monotonically increasing order, with no gaps, then
c xr(n) = xl(n+1), for n = 1, nbins - 1.
c Size nbins, if noptx = 1. Otherwise, not used.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Initial bin index for an EP bin.
dimension nbep (1)
c---- Left-hand probability in initial bin.
dimension pl (1)
c---- Right-hand probability in initial bin.
dimension pr (1)
c---- Left-hand value of x in initial bin.
dimension xl (1)
c---- Right-hand value of x in initial bin.
dimension xr (1)
c.... Local variables.
c---- Bin width, 1.0 or xr(n) - xl(n).
common /lapteqin/ dx
c---- Index of PDF or EP bin.
common /lapteqin/ n
c---- Index of initial PDF bin.
common /lapteqin/ nbin
c---- Index of max contributing PDF bin.
common /lapteqin/ nbinmx
c---- Relative probability of a PDF bin.
common /lapteqin/ pbin
c---- Max contributing probability.
common /lapteqin/ pbinmx
c---- Cumulative relative probability.
common /lapteqin/ pcum
c---- Cum rel prob at left of PDF bin.
common /lapteqin/ pcuml
c---- Cum rel prob at left of EP bin.
common /lapteqin/ pcumle
c---- Cum rel prob at right of PDF bin.
common /lapteqin/ pcumr
c---- Cum rel prob at right of EP bin.
common /lapteqin/ pcumre
c---- PDF bin contributing probability.
common /lapteqin/ ptry
cbugc***DEBUG begins.
cbugc---- Minimum nep needed.
cbug common /lapteqxn/ nepneed
cbugc---- First index with same nbep.
cbug common /lapteqin/ nfirst
cbugc---- Last index with same nbep.
cbug common /lapteqin/ nlast
cbug 9901 format (/ 'apteqin finding equal-probability bins.',
cbug & ' noptx=',i3,' nep=',i5)
cbug 9902 format (i5,' pl,pr=',1p2e22.14)
cbug 9903 format (i5,' xl,xr=',1p2e22.14 / ' pl,pr=',1p2e22.14)
cbug write ( 3, 9901) noptx, nep
cbug if (noptx .eq. 0) then
cbug write (3, 9902) (n, pl(n), pr(n), n = 1, nbins)
cbug elseif (noptx .eq. 1) then
cbug write (3, 9903) (n, xl(n), xr(n), pl(n), pr(n), n = 1, nbins)
cbug endif
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for errors in scalar input.
if (nbins .le. 0) then
nerr = 1
go to 210
endif
if ((noptx .ne. 0) .and. (noptx .ne. 1)) then
nerr = 2
go to 210
endif
if (nep .le. 0) then
nerr = 6
go to 210
endif
c.... Test for errors in input arrays, and find total unnormalized probability.
pcum = 0.0
c---- Loop over bins.
do 120 n = 1, nbins
if (noptx .eq. 0) then
dx = 1.0
else
dx = (xr(n) - xl(n))
if (dx .lt. 0.0) then
nerr = 4
go to 210
endif
if ((n .lt. nbins) .and. (xr(n) .gt. xl(n+1))) then
nerr = 5
go to 210
endif
endif
if ((pl(n) .lt. 0.0) .or. (pr(n) .lt. 0.0)) then
nerr = 3
go to 210
endif
pbin = 0.5 * (pl(n) + pr(n)) * dx
pcum = pcum + pbin
if (n .eq. 1) then
pmin = pbin
else
pmin = amin1 (pmin, pbin)
endif
c---- End of loop over 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
c---- Loop over output EP bins.
do 140 n = 1, nep
pcumle = pcumre
pcumre = n * pcum / nep
nbinmx = nbin
pbinmx = pcumr - pcumle
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
nbep(n) = nbinmx
go to 140
endif
nbin = nbin + 1
pcuml = pcumr
if (noptx .eq. 0) then
dx = 1.0
else
dx = xr(nbin) - xl(nbin)
endif
pbin = 0.5 * (pl(nbin) + pr(nbin)) * dx
pcumr = pcuml + pbin
ptry = amin1 (pbin, pcumre - pcuml)
c---- A bigger probability increment.
if (ptry .gt. pbinmx) then
pbinmx = ptry
nbinmx = nbin
c---- Tested ptry.
endif
go to 130
c---- Tested current input PDF bin.
endif
nbep(n) = nbinmx
c---- End of loop over output EP bins.
140 continue
cbugc***DEBUG begins.
cbug 9905 format (/ 'apteqin results. pmin=',1pe22.14,' nepneed=',i5)
cbug 9906 format (i5,' to',i5,' nbep=',i5,' pcumre=',1pe22.14)
cbug nepneed = 1.0 / (pmin + 1.e-99)
cbug write ( 3, 9905) pmin, nepneed
cbug nfirst = 1
cbug do 170 n = 1, nep
cbug if (nbep(n) .ne. nbep(nfirst)) then
cbug nlast = n - 1
cbug pcumre = nlast * pcum / nep
cbug write ( 3, 9906) nfirst, nlast, nbep(nfirst), pcumre
cbug nfirst = n
cbug endif
cbug 170 continue
cbug nlast = nep
cbug pcumre = pcum
cbug write ( 3, 9906) nfirst, nlast, nbep(nfirst), pcumre
cbugc***DEBUG ends.
210 return
c.... End of subroutine apteqin. (+1 line.)
end
UCRL-WEB-209832