*
* %ranTruncate(mu,sigma,lower,upper) (returns one real value)
*
* mu = mean of Normal distribution
* sigma = standard error of Normal distribution
* lower = lower bound of interval (use %NA for no truncation below)
* upper = upper bound of interval (use %NA for no truncation above)
*
* 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:
* 05/2009 Written by Thomas Doan, Estima
*
function %ranTruncate mu sigma lower upper
type real %ranTruncate mu sigma lower upper
*
local real xlower xupper ylower yupper
local real xdraw udraw logp
local real sfac xmax xmaxs logc
*
* 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 percentiles in the reference logistic
* distribution.
*
if %valid(lower)
compute xlower=(lower-mu)/sigma , ylower=1.0/(1+exp(-xlower/sfac))
else
compute xlower=-1.e+10 , ylower=0.0
*
if %valid(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 %rantruncate=xdraw*sigma+mu
end %rantruncate