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