*
* @ranNormalTrunc(options) draw
*
* Parameters:
* draw (output) is the generated (single-real value) draw
*
* Options:
* mu = mean of Normal distribution [0.0]
* sigma = standard error of Normal distribution [1.0]
* lower = lower bound of interval [-infinity]
* upper = upper bound of interval [+infinity]
*
* This uses the rejection method for draws from a truncated Normal. The comparison
* distribution is a truncated logistic. This could be made more efficient by
* grabbing a whole vector of draws at one time, since everything prior to the
* "loop" instruction is the same each time this is executed with a single set of
* parameters.
*
* Revision schedule:
* 09/1998 Written by Thomas Doan, Estima
*
proc ranNormalTrunc draw
type real *draw
*
option real mu 0.0
option real sigma 1.0
option real lower
option real upper
*
local real xlower xupper ylower yupper
local real xdraw udraw logp
local real sfac xmax xmaxs logc
*
if .not.%defined(draw) {
display "Syntax: @RanNormalTrunc(options) >>draw<<"
return
}
*
* The ratio of the (standard) Normal density to the logistic has two global maxima
* at +/- xmax, where xmax depends upon the scale parameter in the logistic. We
* choose the scale parameter which minimizes the maximum value.
*
compute sfac=.648
compute xmax=1.000
*
* Map endpoints. Normalize lower and upper to xlower and xupper. Map xlower and
* xupper to their fractiles in the reference logistic distribution.
*
if %defined(lower)
compute xlower=(lower-mu)/sigma , ylower=1.0/(1+exp(-xlower/sfac))
else
compute xlower=-1.e+10 , ylower=0.0
*
if %defined(upper)
compute xupper=(upper-mu)/sigma , yupper=1.0/(1+exp(-xupper/sfac))
else
compute xupper=1.e+10 , yupper=1.0
*
* We know what the global max of the ratio of f/g is. We here figure out what the
* maximum is over the interval requested.
*
* 1. If the support interval is contained in (-xmax,+xmax), the
* maximum will be at the larger endpoint in absolute value.
* 2. If the support interval is contained in (xmax,inf), the
* maximum will be at xlower.
* 3. If the support interval is contained in (-inf,-xmax), the
* maximum will be at xupper.
* 4. Otherwise, we have one of the two global maxes. Since we
* need the value only, we don't care which it is.
*
if xlower>-xmax.and.xupperxmax
compute xmaxs=xlower
else
if xupper<-xmax
compute xmaxs=-xupper
else
compute xmaxs=1.0
*
* Compute the log of the reciprocal of the ratio of f/g at the maximum. We're
* ignoring integrating factors in the two density functions
*
compute logc=+.5*xmaxs**2-xmaxs/sfac-2*log(1+exp(-xmaxs/sfac))
*
* To draw a truncated logistic, we first draw a uniform from (ylower,yupper). Back
* transform to a (standardized) logistic. Use the rejection method. When we break
* the loop, we'll have a truncated standardized Normal.
*
loop
compute udraw=%uniform(ylower,yupper)
compute xdraw=-sfac*log((1-udraw)/udraw)
compute logp =logc + $
-.5*xdraw**2+xdraw/sfac+2*log(1+exp(-xdraw/sfac))
if log(%uniform(0.0,1.0))<= logp
break
end loop
*
* Translate to the desired mean and variance
*
compute draw=xdraw*sigma+mu
end ranNormalTrunc