subroutine aptalin (noptx, xl, xr, pl, pr, nbins, pbin, fbin,
& nalb, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTALIN
c
c call aptalin (noptx, xl, xr, pl, pr, nbins, pbin, fbin,
c & nalb, nerr)
c
c Version: aptalin Updated 1990 July 19 10:20.
c aptalin Originated 1990 June 28 10:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To convert a piecewise linear tabulated probability distribution
c function (PFD) to the form required for random sampling by the
c aliasing method. 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 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 After the conversion, each bin will have a probability fbin(n)
c of its own index n, and a probability 1.0 - fbin(n) of the
c aliased index, nalb(n), and a total relative probability pbin(n)
c of sampling either bin n or bin nalb(n).
c
c After aliasing, values of n and x may be randomly sampled as
c follows:
c
c (1) Select a random number u1 uniformly between 0.0 and 1.0.
c (2) Randomly sample the initial bin index n:
c xn = 1 + u1 * nbins
c n = xn
c (3) Choose either bin index n or the alias bin index, nalb(n):
c u2 = xn - n (an independent random number)
c if (u2 .gt. fbin(n)) n = nalb(n)
c (4) 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 aptalsb to sample bin indices from the aliased
c bins. Sample x values from the resulting set of bin indices
c with aptalsh (uniform bins) or aptalsl (linear bins).
c
c History: 1990 July 17 10:40. Changed argument list, added option noptx.
c
c Input: noptx, xl, xr, pl, pr, nbins.
c
c Output: pbin, fbin, nalb, nerr.
c
c Glossary:
c
c fbin Output The probability of selecting bin n instead of the alias
c bin nalb(n). Size nbins.
c
c nalb Output The index of the aliased bin associated with bin n.
c Size nbins.
c
c nbins Input Number of PDF bins or intervals. Size of arrays xl,
c xr, pl, pr, pbin, fbin, nalb.
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 pbin In/Out The unnormalized relative probability of the n'th bin.
c Initially 0.5 * (pl(n) + pr(n)) (noptx = 0) or
c 0.5 * (pl(n) + pr(n)) * (xr(n) - xl(n)) (noptx = 1).
c After aliasing, the output values of pbin will each
c be equal to the mean of all of the initial values.
c Size nbins.
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---- Probability of n in aliased bin n.
dimension fbin (1)
c---- Aliased bin index for bin n.
dimension nalb (1)
c---- Relative prob of n and nalb in bin n.
dimension pbin (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 /laptalin/ dx
c---- Index in fbin, nalb or pbin.
common /laptalin/ n
c---- Index of bin with max pbin value.
common /laptalin/ nmax
c---- Index of bin with min pbin value.
common /laptalin/ nmin
c---- Number of bins aliased.
common /laptalin/ nn
c---- Maximum value of pbin.
common /laptalin/ pmax
c---- Mean value of pbin.
common /laptalin/ pmean
c---- Minimum value of pbin.
common /laptalin/ pmin
c---- Sum of pbin values.
common /laptalin/ psum
cbugc***DEBUG begins.
cbugc---- Number of occurrences of a bin index.
cbug common /laptalin/ nbin
cbug 9901 format (/ 'aptalin converting a PDF for aliasing. 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 relative
c.... probability of each bin, and the average for all the bins.
psum = 0.0
c---- Loop over pbin.
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(n) = 0.5 * (pl(n) + pr(n)) * dx
psum = psum + pbin(n)
c---- End of loop over pbin.
120 continue
pmean = psum / nbins
c.... Initialize.
c---- Loop over bins.
do 130 n = 1, nbins
fbin(n) = 1.0
nalb(n) = n
c---- End of loop over bins.
130 continue
if (nbins .gt. 1) then
c.... Alias the least probable bin to the most probable, until done.
c---- Loop over min and max pbin.
do 140 nn = 1, nbins - 1
c.... Find the min and max pbin values, and their indices.
call aminmx (pbin, 1, nbins, 1, pmin, pmax, nmin, nmax)
c++++ Alias min bin to max bin.
if (nalb(nmin) .eq. nmin) then
fbin(nmin) = pmin / pmean
nalb(nmin) = nmax
pbin(nmin) = pmean
pbin(nmax) = pmin + pmax - pmean
endif
c---- End of loop over min and max pbin.
140 continue
c---- Tested nbins.
endif
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptalin results. pmean=',1pe22.14 /
cbug & (i3,' nalb=',i3,' fbin=',1pe22.14, ' pbin=',1pe22.14))
cbug write ( 3, 9904) pmean, (n, nalb(n), fbin(n), pbin(n),
cbug & n = 1, nbins)
cbug
cbug 9905 format (/ ' check total probabilities.')
cbug 9906 format (' nn=',i3,' nbin=',i3,' psum=',1pe22.14)
cbug write ( 3, 9905)
cbug do 160 nn = 1, nbins
cbug nbin = 1
cbug psum = fbin(nn) * pbin(nn)
cbug do 150 n = 1, nbins
cbug if (nalb(n) .eq. nn) then
cbug psum = psum + (1.0 - fbin(n)) * pbin(n)
cbug if (n .ne. nn) then
cbug nbin = nbin + 1
cbug endif
cbug endif
cbug 150 continue
cbug write ( 3, 9906) nn, nbin, psum
cbug 160 continue
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptalin. (+1 line.)
end
UCRL-WEB-209832