The RATS Software Forum

Code: Select all

*
* %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.xupper<xmax
   compute xmaxs=%max(abs(xlower),abs(xupper))
else
if xlower>xmax
   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