subroutine aptalsl (xl, xr, pl, pr, nbins, nb, np, xran, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTALSL
c
c call aptalsl (xl, xr, pl, pr, nbins, nb, np, xran, nerr)
c
c Version: aptalsl Updated 1990 November 27 14:00.
c aptalsl Originated 1990 July 19 13:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np values of xran from nbins bins.
c Each bin has a relative probability distribution function
c varying linearly from pl(n) at xl(n) on the left, to pr(n)
c at xr(n) on the right. Each value xran(n) is to be sampled
c from bin nb(n), where nb must be in the range from 1 to nbins.
c Flag nerr indicates any input error.
c
c Note: Use subroutines aptalsb, aptalst, aptcums, aptcumt,
c apteqsb or apteqxs to randomly sample the bin indices nb
c from any tabulated aliased, cumulative, or equal-probability
c probability distribution function.
c Use subroutine aptalsh instead of aptalsl when the distribution
c in each bin is uniform instead of linear.
c
c Input: xl, xr, pl, pr, nbins, nb, np.
c
c Output: xran, nerr.
c
c Glossary:
c
c nb Input Bin number for which a value of xran is to be randomly
c sampled. Must be in the range from 1 to nbins.
c Size np.
c
c nbins Input Size of arrays xl, xr, pl, pr.
c
c nerr Output Indicates an input error, if not 0.
c 1 if nbins is not positive.
c 2 if np is not positive.
c
c np Input Number of random samples xran. Size of nb, xran.
c
c pl, pr Input Relative probabilities of random variable xran at
c xl and xr, respectively. Size nbins.
c Need not be normalized, but may not be negative.
c Relative probability varies linearly across bin.
c To find pl and pr, when is given:
c must be in the central third of the zone,
c or between xma and xmb, where
c xma = (2.0 * xl + xr) / 3.0,
c xmb = (xl + 2.0 * xr) / 3.0.
c Then, for any pl, or for any pr,
c ( - xma) * pl = (xmb - ) * pr.
c
c xl, xr Input Values of random variable xran at the left and right
c boundaries of the bin, respectively. Size nbins.
c
c xran Output Randomly sampled value from bin nb(n). Size np.
c in the range from xl(nb(n)) to xr(nb(n)).
c In a bin, the expected value of xran is
c = fa * xl + fb * xr,
c where fa = (2.0 * pl + pr) / (3.0 * (pl + pr))
c and fb = (pl + 2.0 * pr) / (3.0 * (pl + pr)).
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Index of bin to be sampled.
dimension nb (1)
c---- Probability at xl.
dimension pl (1)
c---- Probability at xr.
dimension pr (1)
c---- Value of xran with probability pl.
dimension xl (1)
c---- Value of xran with probability pr.
dimension xr (1)
c---- Randomly sampled value of xran.
dimension xran (1)
c.... Local variables.
c---- Index in local array.
common /laptalsl/ n
c---- First index of subset of data.
common /laptalsl/ n1
c---- Last index of subset of data.
common /laptalsl/ n2
c---- Index of a bin.
common /laptalsl/ nbin
c---- Index in external array.
common /laptalsl/ nn
c---- Size of current subset of data.
common /laptalsl/ ns
c---- Probability of x1.
common /laptalsl/ p1
c---- Sum of pl and pr.
common /laptalsl/ psum
c---- First random number.
common /laptalsl/ ranf1 (64)
c---- Second random number.
common /laptalsl/ ranf2 (64)
c---- First random number ranf1.
common /laptalsl/ u1
c---- Value 1.0 - u1.
common /laptalsl/ um
c---- First random value of xran.
common /laptalsl/ x1
c---- Second random value of xran.
common /laptalsl/ x2
cbugc***DEBUG begins.
cbugc---- Mean value of xran.
cbug common /laptalsl/ xmean
cbugc---- Standard deviation from mean xran.
cbug common /laptalsl/ xdev
cbug 9901 format (/ 'aptalsl sampling from a linear distribution.' /
cbug & (i3,' xl,xr=',1p2e22.14 /
cbug & ' pl,pr=',1p2e22.14))
cbug 9902 format (/ 'np=',i5,' nb=' / (20i4))
cbug write ( 3, 9901) (n, xl(n), xr(n), pl(n), pr(n), n = 1, nbins)
cbug write ( 3, 9902) np, (nb(n), n = 1, np)
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for input errors.
if (nbins .le. 0) then
nerr = 1
go to 210
endif
if (np .le. 0) then
nerr = 2
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Generate ns pairs of random numbers.
c---- Loop over local arrays.
do 120 n = 1, ns
ranf1(n) = ranf( )
c---- End of loop over local arrays.
120 continue
c---- Loop over local arrays.
do 130 n = 1, ns
ranf2(n) = ranf( )
c---- End of loop over local arrays.
130 continue
c.... Randomly sample values of xran from the bins.
c---- Loop over local arrays.
do 140 n = 1, ns
nn = n + n1 - 1
nbin = nb(nn)
u1 = ranf1(n)
um = 1.0 - u1
x1 = um * xl(nbin) + u1 * xr(nbin)
x2 = u1 * xl(nbin) + um * xr(nbin)
p1 = um * pl(nbin) + u1 * pr(nbin)
psum = pl(nbin) + pr(nbin)
if (p1 .ge. ranf2(n) * psum) then
xran(nn) = x1
else
xran(nn) = x2
endif
c---- End of loop over local arrays.
140 continue
c.... See if all data subsets are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptalsl results: xran')
cbug write ( 3, 9903)
cbug call aptmean (xran, np, 1.e-11, xmean, xdev, nerr)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptalsl. (+1 line.)
end
UCRL-WEB-209832