subroutine aptcump (noptx, xl, xr, pl, pr, nbins, pcum, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCUMP
c
c call aptcump (noptx, xl, xr, pl, pr, nbins, pcum, nerr)
c
c Version: aptcump Updated 1990 July 19 10:20.
c aptcump Originated 1990 June 18 11:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the cumulative probability distribution function (CPDF)
c corresponding to a piecewise linear tabulated probability
c distribution function (PDF). The PDF is defined by either of
c 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 n.
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
c The CPDF will be pcum(n), n = 1, nbins.
c
c Values of n and (if noptx = 1) x may be randomly sampled from
c the CPDF as follows:
c
c (1) Select a random number u1 uniformly between 0.0 and 1.0,
c and find the product pran = u1 * pcum(nbins).
c (2) Find the table interval n, between 1 and nbins, for which
c pcum(n-1) < pran < pcum(n).
c (3) Randomly sample a value of x in the range from xl(n) to
c xr(n), based on the linear probability distribution from
c pl(n) to pr(n) (uniform if pl(n) = pr(n)).
c
c Flag nerr indicates any input error, if not zero.
c
c Note: Use subroutine aptcums or aptcumt to sample bin indices from
c one or more cumulative PDFs.
c Sample x values from the resulting set of bin indices
c with aptalsh (uniform bins) or aptalsl (linear bins).
c
c Input: noptx, xl, xr, pl, pr, nbins.
c
c Output: pcum, nerr.
c
c Glossary:
c
c nbins Input Number of PDF bins or intervals. Size of arrays xl,
c xr, pl, pr, pcum.
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
c noptx Input Indicates the type of PDF:
c 0 if the bin index is the random variable, xl and xr
c are not used, and pl(n) = pr(n) (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, in the n'th interval.
c
c pcum In/Out The unnormalized relative cumulative probability of
c bins or x intervals 1 through n. Size nbins.
c For n = 1, nbins, pcum(n) is obtained by adding
c 0.5 * (pl(n) + pr(n)) (noptx = 0) or
c 0.5 * (pl(n) + pr(n)) * (xr(n) - xl(n)) (noptx = 1)
c to pcum(n-1) (substituting 0.0 for pcum(0)).
c
c pl Input For bin n, the unnormalized relative probability of the
c bin (noptx = 0) or of xl(n), per unit x interval
c (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 pr Input For bin n, the unnormalized relative probability of the
c bin (noptx = 0) or of xr(n), per unit x interval
c (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 bin n, the value of the random variable x at the
c 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 bin n, the value of the random variable x at the
c 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---- Cumulative probability, bins 1 to n.
dimension pcum (1)
c---- Left-hand probability in bin.
dimension pl (1)
c---- Right-hand probability in bin.
dimension pr (1)
c---- Left-hand value of x in bin.
dimension xl (1)
c---- Right-hand value of x in bin.
dimension xr (1)
c.... Local variables.
c---- Bin width, 1.0 or xr(n) - xl(n).
common /laptcump/ dx
c---- Index in pcum, pl, pr, xl, xr.
common /laptcump/ n
c---- Bin relative probability.
common /laptcump/ pbin
c---- Cumulative relative probability.
common /laptcump/ psum
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcump finding cumulative probability. noptx=',i3)
cbug 9902 format (i3,' pl,pr=',1p2e22.14)
cbug 9903 format (i3,' xl,xr=',1p2e22.14 / ' pl,pr=',1p2e22.14)
cbug write ( 3, 9901) noptx
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
c.... Test for errors in input arrays, and find the cumulative relative
c.... probability for each bin.
psum = 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
psum = psum + pbin
pcum(n) = psum
c---- End of loop over bins.
120 continue
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptcump results.')
cbug 9905 format (i3,' pbin=',1pe22.14,' pcum=',1pe22.14)
cbug write ( 3, 9904)
cbug n = 1
cbug write ( 3, 9905) n, pcum(1), pcum(1)
cbug if (nbins .ge. 2) then
cbug do 130 n = 2, nbins
cbug pbin = pcum(n) - pcum(n-1)
cbug write ( 3, 9905) n, pbin, pcum(n)
cbug 130 continue
cbug endif
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptcump. (+1 line.)
end
UCRL-WEB-209832