one draw from truncated normal distribution (mu,sigma^2)

// ---------------------------------------------
// Returns one draw from truncated normal distribution (mu,sigma^2) with range
// (bound,+inf) if lb=TRUE and
// (-inf,bound) if lb=FALSE
// from http://athens.src.uchicago.edu/jenni/econ319_2003/lecture.html
// ---------------------------------------------

double truncn(double bound, bool lb, double mu, double sigma){

  double c, z, w;

  // 1. standardised cut-off c for truncation from below or above
  if(lb == TRUE){
    c = (bound-mu)/sigma;
  } else{
    c = -(bound-mu)/sigma;
  }

  // 2. standardised draw using Geweke (1991)
  if(c < 0.45){ // normal rejection sampling
    z = ::Rf_rnorm(0.0,1.0);
    while(z < c){
      z = ::Rf_rnorm(0.0,1.0);
    }
  } else{ // exponential rejection sampling
    z = -log(1-::Rf_runif(0.0,1.0))/c;
    w = ::Rf_runif(0.0,1.0);
    while(w > exp(-0.5*pow(z,2))){
      z = -log(1-::Rf_runif(0.0,1.0))/c;
      w = ::Rf_runif(0.0,1.0);
    }
    z = z+c;
  }

  // 3. reverse standardisation
  if(lb == TRUE){
    return mu + sigma*z;
  } else{
    return mu - sigma*z;
  }
}