rng.h

Go to the documentation of this file.

/* 
 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
 * and Daniel Pemstein.  All Rights Reserved.
 *
 * This program is free software; you can redistribute it and/or
 * modify under the terms of the GNU General Public License as
 * published by Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.  See the text files
 * COPYING and LICENSE, distributed with this source code, for further
 * information.
 * --------------------------------------------------------------------
 *  scythestat/rng.h
 *
 * The code for many of the RNGs defined in this file and implemented
 * in rng.cc is based on that in the R project, version 1.6.0-1.7.1.
 * This code is available under the terms of the GNU GPL.  Original
 * copyright:
 * 
 * Copyright (C) 1998      Ross Ihaka
 * Copyright (C) 2000-2002 The R Development Core Team
 * Copyright (C) 2003      The R Foundation
 */

/* Doxygen doesn't deal well with the macros that we use to make
 * matrix versions of rngs easy to define.
 */

#ifndef SCYTHE_RNG_H
#define SCYTHE_RNG_H

#include <iostream>
#include <cmath>

#ifdef HAVE_IEEEFP_H
#include <ieeefp.h>
#endif

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "error.h"
#include "algorithm.h"
#include "distributions.h"
#include "ide.h"
#include "la.h"
#else
#include "scythestat/matrix.h"
#include "scythestat/error.h"
#include "scythestat/algorithm.h"
#include "scythestat/distributions.h"
#include "scythestat/ide.h"
#include "scythestat/la.h"
#endif

namespace scythe {

/* Shorthand for the matrix versions of the various distributions'
 * random number generators.
 */
  
#define SCYTHE_RNGMETH_MATRIX(NAME, RTYPE, ARGNAMES, ...)             \
  template <matrix_order O, matrix_style S>                           \
  Matrix<RTYPE, O, S>                                                 \
  NAME (unsigned int rows, unsigned int cols, __VA_ARGS__)            \
  {                                                                   \
    Matrix<RTYPE, O, Concrete> ret(rows, cols, false);                \
    typename Matrix<RTYPE,O,Concrete>::forward_iterator it;           \
    typename Matrix<RTYPE,O,Concrete>::forward_iterator last          \
      = ret.end_f();                                                  \
    for (it = ret.begin_f(); it != last; ++it)                        \
      *it = NAME (ARGNAMES);                                          \
    SCYTHE_VIEW_RETURN(RTYPE, O, S, ret)                              \
  }                                                                   \
                                                                      \
  Matrix<RTYPE, Col, Concrete>                                        \
  NAME (unsigned int rows, unsigned int cols, __VA_ARGS__)            \
  {                                                                   \
    return NAME <Col,Concrete> (rows, cols, ARGNAMES);                \
  }

  template <class RNGTYPE>
  class rng
  {
    public:

      /* This declaration allows users to treat rng objects like
       * functors that generate random uniform numbers.  This can be
       * quite convenient.
       */
      double operator() ()
      {
        return runif();
      }

      /* Returns random uniform numbers on (0, 1).  This function must
       * be implemented by extending classes */
      double runif ()
      {
        return as_derived().runif();
      }


      /* No point in declaring these virtual because we have to
       * override them anyway because C++ isn't too bright.  Also, it
       * is illegal to make template methods virtual
       */
      template <matrix_order O, matrix_style S>
      Matrix<double,O,S> runif(unsigned int rows, 
                               unsigned int cols) 
      {
        Matrix<double, O, S> ret(rows, cols, false);
        typename Matrix<double,O,S>::forward_iterator it;
        typename Matrix<double,O,S>::forward_iterator last=ret.end_f();
        for (it = ret.begin_f(); it != last; ++it)
          *it = runif();

        return ret;
      }

      Matrix<double,Col,Concrete> runif(unsigned int rows, 
                                        unsigned int cols)
      {
        return runif<Col,Concrete>(rows, cols);
      }

      double
      rbeta (double alpha, double beta)
      {
        double report;
        double xalpha, xbeta;
        
        // Check for allowable parameters
        SCYTHE_CHECK_10(alpha <= 0, scythe_invalid_arg, "alpha <= 0");
        SCYTHE_CHECK_10(beta <= 0, scythe_invalid_arg, "beta <= 0");
        
        xalpha = rchisq (2 * alpha);
        xbeta = rchisq (2 * beta);
        report = xalpha / (xalpha + xbeta);
        
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rbeta, double, SCYTHE_ARGSET(alpha, beta),
          double alpha, double beta);

      double 
      rnchypgeom(double m1, double n1, double n2, double psi,
                 double delta)
      {
        // Calculate mode of mass function
        double a = psi - 1;
        double b = -1 * ((n1+m1+2)*psi + n2 - m1);
        double c = psi * (n1+1) * (m1+1);
        double q = -0.5 * ( b + sgn(b) * 
            std::sqrt(std::pow(b,2) - 4*a*c));
        double root1 = c/q;
        double root2 = q/a;
        double el = std::max(0.0, m1-n2);
        double u = std::min(n1,m1);
        double mode = std::floor(root1);
        int exactcheck = 0;
        if (u<mode || mode<el) {
          mode = std::floor(root2);
          exactcheck = 1;
        }
     

        int size = static_cast<int>(u+1);

        double *fvec = new double[size];
        fvec[static_cast<int>(mode)] = 1.0;
        double s;
        // compute the mass function at y
        if (delta <= 0 || exactcheck==1){  //exact evaluation 
          // sum from mode to u
          double f = 1.0;
          s = 1.0;
          for (double i=(mode+1); i<=u; ++i){
            double r = ((n1-i+1)*(m1-i+1))/(i*(n2-m1+i)) * psi;
            f = f*r;
            s += f;
            fvec[static_cast<int>(i)] = f;
          }
         
          // sum from mode to el
          f = 1.0;
          for (double i=(mode-1); i>=el; --i){
            double r = ((n1-i)*(m1-i))/((i+1)*(n2-m1+i+1)) * psi;
            f = f/r;
            s += f;
            fvec[static_cast<int>(i)] = f;
          }
        } else { // approximation
          double epsilon = delta/10.0;
          // sum from mode to ustar
          double f = 1.0;
          s = 1.0;
          double i = mode+1;
          double r;
          do {
            if (i>u) break;
            r = ((n1-i+1)*(m1-i+1))/(i*(n2-m1+i)) * psi;
            f = f*r;
            s += f;
            fvec[static_cast<int>(i)] = f;
            ++i;
          } while(f>=epsilon || r>=5.0/6.0);
         
          // sum from mode to elstar
          f = 1.0;
          i = mode-1;
          do {
            if (i<el) break;
            r = ((n1-i)*(m1-i))/((i+1)*(n2-m1+i+1)) * psi;
            f = f/r;
            s += f;
            fvec[static_cast<int>(i)] = f;
            --i;
          } while(f>=epsilon || r <=6.0/5.0);         
        }

        double udraw = runif();
        double psum = fvec[static_cast<int>(mode)]/s;
        if (udraw<=psum)
          return mode;
        double lower = mode-1;
        double upper = mode+1;

        do{
          double fl;
          double fu;
          if (lower >= el)
            fl = fvec[static_cast<int>(lower)];
          else 
            fl = 0.0;

          if (upper <= u)
            fu = fvec[static_cast<int>(upper)];
          else
            fu = 0.0;

          if (fl > fu) {
            psum += fl/s;
            if (udraw<=psum)
              return lower;
            --lower;
          } else {
            psum += fu/s;
            if (udraw<=psum)
              return upper;
            ++upper;
          }
        } while(udraw>psum);
       
        delete [] fvec;
        SCYTHE_THROW(scythe_convergence_error,
          "Algorithm did not converge");
      }

      SCYTHE_RNGMETH_MATRIX(rnchypgeom, double,
          SCYTHE_ARGSET(m1, n1, n2, psi, delta), double m1, double n1,
          double n2, double psi, double delta);

      unsigned int
      rbern (double p)
      {
        unsigned int report;
        double unif;
          
        // Check for allowable paramters
        SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg,
            "p parameter not in[0,1]");
        
        unif = runif ();
        if (unif < p)
          report = 1;
        else
          report = 0;
        
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rbern, unsigned int, p, double p);
      
      unsigned int
      rbinom (unsigned int n, double p)
      {
        unsigned int report;
        unsigned int count = 0;
        double hold;
          
        // Check for allowable parameters
        SCYTHE_CHECK_10(n == 0, scythe_invalid_arg, "n == 0");
        SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg, 
            "p not in [0,1]");
          
        // Loop and count successes
        for (unsigned int i = 0; i < n; i++) {
          hold = runif ();
          if (hold < p)
            ++count;
        }
        report = count;
        
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rbinom, unsigned int, SCYTHE_ARGSET(n, p),
          unsigned int n, double p);

      double
      rchisq (double df)
      {
        double report;
          
        // Check for allowable paramter
        SCYTHE_CHECK_10(df <= 0, scythe_invalid_arg,
            "Degrees of freedom <= 0");
      
        // Return Gamma(nu/2, 1/2) variate
        report = rgamma (df / 2, .5);
        
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rchisq, double, df, double df);

      double
      rexp (double invscale)
      {
        double report;
        
        // Check for allowable parameter
        SCYTHE_CHECK_10(invscale <= 0, scythe_invalid_arg,
            "Inverse scale parameter <= 0");
        
        report = -std::log (runif ()) / invscale;
        
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rexp, double, invscale, double invscale);
    
      double
      rf (double df1, double df2)
      {
        SCYTHE_CHECK_10(df1 <= 0 || df2 <= 0, scythe_invalid_arg,
            "n1 or n2 <= 0");

        return ((rchisq(df1) / df1) / (rchisq(df2) / df2));
      }

      SCYTHE_RNGMETH_MATRIX(rf, double, SCYTHE_ARGSET(df1, df2),
          double df1, double df2);

      double
      rgamma (double shape, double rate)
      {
        double report;

        // Check for allowable parameters
        SCYTHE_CHECK_10(shape <= 0, scythe_invalid_arg, "shape <= 0");
        SCYTHE_CHECK_10(rate <= 0, scythe_invalid_arg, "rate <= 0");

        if (shape > 1)
          report = rgamma1 (shape) / rate;
        else if (shape == 1)
          report = -std::log (runif ()) / rate;
        else
          report = rgamma1 (shape + 1) 
            * std::pow (runif (), 1 / shape) / rate;

        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rgamma, double, SCYTHE_ARGSET(shape, rate),
          double shape, double rate);

      double
      rlogis (double location, double scale)
      {
        double report;
        double unif;
          
        // Check for allowable paramters
        SCYTHE_CHECK_10(scale <= 0, scythe_invalid_arg, "scale <= 0");
        
        unif = runif ();
        report = location + scale * std::log (unif / (1 - unif));
        
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rlogis, double, 
          SCYTHE_ARGSET(location, scale),
          double location, double scale);

      double
      rlnorm (double logmean, double logsd)
      {
        SCYTHE_CHECK_10(logsd < 0.0, scythe_invalid_arg,
            "standard deviation < 0");

        return std::exp(rnorm(logmean, logsd));
      }

      SCYTHE_RNGMETH_MATRIX(rlnorm, double, 
          SCYTHE_ARGSET(logmean, logsd),
          double logmean, double logsd);

      unsigned int
      rnbinom (double n, double p)
      {
        SCYTHE_CHECK_10(n == 0 || p <= 0 || p > 1, scythe_invalid_arg,
            "n == 0, p <= 0, or p > 1");

        return rpois(rgamma(n, (1 - p) / p));
      }

      SCYTHE_RNGMETH_MATRIX(rnbinom, unsigned int,
          SCYTHE_ARGSET(n, p), double n, double p);

      double
      rnorm (double mean = 0, double sd = 1)
      {
        SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg, 
            "Negative standard deviation");
        
        return (mean + rnorm1 () * sd);
      }

      SCYTHE_RNGMETH_MATRIX(rnorm, double, SCYTHE_ARGSET(mean, sd),
          double mean, double sd);

      unsigned int
      rpois(double lambda)
      {
        SCYTHE_CHECK_10(lambda <= 0, scythe_invalid_arg, "lambda <= 0");
        unsigned int n;
        
        if (lambda < 33) {
          double cutoff = std::exp(-lambda);
          n = -1;
          double t = 1.0;
          do {
            ++n;
            t *= runif();
          } while (t > cutoff);    
        } else {
          bool accept = false;
          double c = 0.767 - 3.36/lambda;
          double beta = M_PI/std::sqrt(3*lambda);
          double alpha = lambda*beta;
          double k = std::log(c) - lambda - std::log(beta);
            
          while (! accept){
            double u1 = runif();
            double x = (alpha - std::log((1-u1)/u1))/beta;
            while (x <= -0.5){
              u1 = runif();
              x = (alpha - std::log((1-u1)/u1))/beta;
            } 
            n = static_cast<int>(x + 0.5);
            double u2 = runif();
            double lhs = alpha - beta*x +
              std::log(u2/std::pow(1+std::exp(alpha-beta*x),2));
            double rhs = k + n*std::log(lambda) - lnfactorial(n);
            if (lhs <= rhs)
              accept = true;
          }
        }
        
        return n;
      }

      SCYTHE_RNGMETH_MATRIX(rpois, unsigned int, lambda, double lambda);

      /* There is a naming issue here, with respect to the p- and d-
       * functions in distributions.  This is really analagous to rt1-
       * and dt1- XXX Clear up.  Also, we should probably have a
       * random number generator for both versions of the student t.
       */

      double
      rt (double mu, double sigma2, double nu)
      {
        double report;
        double x, z;
          
        // Check for allowable paramters
        SCYTHE_CHECK_10(sigma2 <= 0, scythe_invalid_arg,
            "Variance parameter sigma2 <= 0");
        SCYTHE_CHECK_10(nu <= 0, scythe_invalid_arg,
            "D.O.F parameter nu <= 0");
        
        z = rnorm1 ();
        x = rchisq (nu);
        report = mu + std::sqrt (sigma2) * z 
          * std::sqrt (nu) / std::sqrt (x);
        
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rt1, double, SCYTHE_ARGSET(mu, sigma2, nu),
          double mu, double sigma2, double nu);

      double
      rweibull (double shape, double scale)
      {
        SCYTHE_CHECK_10(shape <= 0 || scale <= 0, scythe_invalid_arg,
            "shape or scale <= 0");

        return scale * std::pow(-std::log(runif()), 1.0 / shape);
      }

      SCYTHE_RNGMETH_MATRIX(rweibull, double,
          SCYTHE_ARGSET(shape, scale), double shape, double scale);

      double
      richisq (double nu)
      {
        double report;
          
        // Check for allowable parameter
        SCYTHE_CHECK_10(nu <= 0, scythe_invalid_arg,
            "Degrees of freedom <= 0");
          
        // Return Inverse-Gamma(nu/2, 1/2) variate
        report = rigamma (nu / 2, .5);
        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(richisq, double, nu, double nu);

      double
      rigamma (double alpha, double beta)
      {
        double report;
        
        // Check for allowable parameters
        SCYTHE_CHECK_10(alpha <= 0, scythe_invalid_arg, "alpha <= 0");
        SCYTHE_CHECK_10(beta <= 0, scythe_invalid_arg, "beta <= 0");

        // Return reciprocal of gamma variate
        report = std::pow (rgamma (alpha, beta), -1);

        return (report);
      }

      SCYTHE_RNGMETH_MATRIX(rigamma, double, SCYTHE_ARGSET(alpha, beta),
          double alpha, double beta);

      /* Truncated Distributions */

      double 
      rtnorm(double mean, double variance, double below, double above)
      {  
        SCYTHE_CHECK_10(below >= above, scythe_invalid_arg,
            "Truncation bound not logically consistent");
        SCYTHE_CHECK_10(variance <= 0, scythe_invalid_arg,
            "Variance <= 0");
        
        double sd = std::sqrt(variance);
        double FA = 0.0;
        double FB = 0.0;
        if ((std::fabs((above-mean)/sd) < 8.2) 
            && (std::fabs((below-mean)/sd) < 8.2)){
          FA = pnorm1((above-mean)/sd, true, false);
          FB = pnorm1((below-mean)/sd, true, false);
        }
        if ((((above-mean)/sd) < 8.2)  && (((below-mean)/sd) <= -8.2) ){ 
          FA = pnorm1((above-mean)/sd, true, false);
          FB = 0.0;
        }
        if ( (((above-mean)/sd) >= 8.2)  && (((below-mean)/sd) > -8.2) ){ 
          FA = 1.0;
          FB = pnorm1((below-mean)/sd, true, false);
        } 
        if ( (((above-mean)/sd) >= 8.2) && (((below-mean)/sd) <= -8.2)){
          FA = 1.0;
          FB = 0.0;
        }
        double term = runif()*(FA-FB)+FB;
        if (term < 5.6e-17)
          term = 5.6e-17;
        if (term > (1 - 5.6e-17))
          term = 1 - 5.6e-17;
        double draw = mean + sd * qnorm1(term);
        if (draw > above)
          draw = above;
        if (draw < below)
          draw = below;
         
        return draw;
      }

      SCYTHE_RNGMETH_MATRIX(rtnorm, double, 
          SCYTHE_ARGSET(mean, variance, above, below), double mean,
          double variance, double above, double below);

      double
      rtnorm_combo(double mean, double variance, double below, 
                   double above)
      {
        SCYTHE_CHECK_10(below >= above, scythe_invalid_arg,
            "Truncation bound not logically consistent");
        SCYTHE_CHECK_10(variance <= 0, scythe_invalid_arg,
            "Variance <= 0");
        
        double sd = std::sqrt(variance);
        if ((((above-mean)/sd > 0.5) && ((mean-below)/sd > 0.5))
            ||
            (((above-mean)/sd > 2.0) && ((below-mean)/sd < 0.25))
            ||
            (((mean-below)/sd > 2.0) && ((above-mean)/sd > -0.25))) { 
          double x = rnorm(mean, sd);
          while ((x > above) || (x < below))
            x = rnorm(mean,sd);
          return x;
        } else {
          // use the inverse cdf method
          double FA = 0.0;
          double FB = 0.0;
          if ((std::fabs((above-mean)/sd) < 8.2) 
              && (std::fabs((below-mean)/sd) < 8.2)){
            FA = pnorm1((above-mean)/sd, true, false);
            FB = pnorm1((below-mean)/sd, true, false);
          }
          if ((((above-mean)/sd) < 8.2)  && (((below-mean)/sd) <= -8.2) ){ 
            FA = pnorm1((above-mean)/sd, true, false);
            FB = 0.0;
          }
          if ( (((above-mean)/sd) >= 8.2)  && (((below-mean)/sd) > -8.2) ){ 
            FA = 1.0;
            FB = pnorm1((below-mean)/sd, true, false);
          } 
          if ( (((above-mean)/sd) >= 8.2) && (((below-mean)/sd) <= -8.2)){
            FA = 1.0;
            FB = 0.0;
          }
          double term = runif()*(FA-FB)+FB;
          if (term < 5.6e-17)
            term = 5.6e-17;
          if (term > (1 - 5.6e-17))
            term = 1 - 5.6e-17;
          double x = mean + sd * qnorm1(term);
          if (x > above)
            x = above;
          if (x < below)
            x = below;
          return x;
        }    
      }

      SCYTHE_RNGMETH_MATRIX(rtnorm_combo, double, 
          SCYTHE_ARGSET(mean, variance, above, below), double mean,
          double variance, double above, double below);

      double
      rtbnorm_slice (double mean, double variance, double below,
                     unsigned int iter = 10)
      {
        SCYTHE_CHECK_10(below < mean, scythe_invalid_arg,
            "Truncation point < mean");
        SCYTHE_CHECK_10(variance <= 0, scythe_invalid_arg,
            "Variance <= 0");
         
        double z = 0;
        double x = below + .00001;
         
        for (unsigned int i=0; i<iter; ++i){
          z = runif()*std::exp(-1*std::pow((x-mean),2)/(2*variance));
          x = runif()*
            ((mean + std::sqrt(-2*variance*std::log(z))) - below) + below;
        }

        if (! finite(x)) {
          SCYTHE_WARN("Mean extremely far from truncation point. "
              << "Returning truncation point");
          return below; 
        }

        return x;
      }

      SCYTHE_RNGMETH_MATRIX(rtbnorm_slice, double, 
          SCYTHE_ARGSET(mean, variance, below, iter), double mean, 
          double variance, double below, unsigned int iter = 10);

      double
      rtanorm_slice (double mean, double variance, double above, 
          unsigned int iter = 10)
      {
        SCYTHE_CHECK_10(above > mean, scythe_invalid_arg,
            "Truncation point > mean");
        SCYTHE_CHECK_10(variance <= 0, scythe_invalid_arg,
            "Variance <= 0");
      
        double below = -1*above;
        double newmu = -1*mean;
        double z = 0;
        double x = below + .00001;
         
        for (unsigned int i=0; i<iter; ++i){
          z = runif()*std::exp(-1*std::pow((x-newmu),2)
              /(2*variance));
          x = runif()
            *( (newmu + std::sqrt(-2*variance*std::log(z))) - below) 
            + below;
        }
        if (! finite(x)) {
          SCYTHE_WARN("Mean extremely far from truncation point. "
              << "Returning truncation point");
          return above; 
        }
        
        return -1*x;
      }

      SCYTHE_RNGMETH_MATRIX(rtanorm_slice, double, 
          SCYTHE_ARGSET(mean, variance, above, iter), double mean, 
          double variance, double above, unsigned int iter = 10);

      double
      rtbnorm_combo (double mean, double variance, double below, 
          unsigned int iter = 10)
      {
        SCYTHE_CHECK_10(variance <= 0, scythe_invalid_arg,
            "Variance <= 0");
        
        double s = std::sqrt(variance);
        // do rejection sampling and return value
        //if (m >= below){
        if ((mean/s - below/s ) > -0.5){
          double x = rnorm(mean, s);
          while (x < below)
            x = rnorm(mean,s);
          return x; 
        } else if ((mean/s - below/s ) > -5.0 ){
          // use the inverse cdf method
          double above =  std::numeric_limits<double>::infinity();
          double x = rtnorm(mean, variance, below, above);
          return x;
        } else {
          // do slice sampling and return value
          double z = 0;
          double x = below + .00001;
          for (unsigned int i=0; i<iter; ++i){
            z = runif() * std::exp(-1 * std::pow((x - mean), 2)
                / (2 * variance));
            x = runif() 
              * ((mean + std::sqrt(-2 * variance * std::log(z))) 
                - below) + below;
          }
          if (! finite(x)) {
            SCYTHE_WARN("Mean extremely far from truncation point. "
                << "Returning truncation point");
            return below; 
          }
          return x;
        }
      }

      SCYTHE_RNGMETH_MATRIX(rtbnorm_combo, double, 
          SCYTHE_ARGSET(mean, variance, below, iter), double mean, 
          double variance, double below, unsigned int iter = 10);

      double
      rtanorm_combo (double mean, double variance, double above, 
          const unsigned int iter = 10)
      {
        SCYTHE_CHECK_10(variance <= 0, scythe_invalid_arg,
            "Variance <= 0");

        double s = std::sqrt(variance);
        // do rejection sampling and return value
        if ((mean/s - above/s ) < 0.5){ 
          double x = rnorm(mean, s);
          while (x > above)
            x = rnorm(mean,s);
          return x;
        } else if ((mean/s - above/s ) < 5.0 ){
          // use the inverse cdf method
          double below =  -std::numeric_limits<double>::infinity();
          double x = rtnorm(mean, variance, below, above);
          return x;
        } else {
          // do slice sampling and return value
          double below = -1*above;
          double newmu = -1*mean;
          double z = 0;
          double x = below + .00001;
             
          for (unsigned int i=0; i<iter; ++i){
            z = runif() * std::exp(-1 * std::pow((x-newmu), 2)
                /(2 * variance));
            x = runif() 
              * ((newmu + std::sqrt(-2 * variance * std::log(z)))
                  - below) + below;
          }
          if (! finite(x)) {
            SCYTHE_WARN("Mean extremely far from truncation point. "
                << "Returning truncation point");
            return above; 
          }
          return -1*x;
        }
      }

      SCYTHE_RNGMETH_MATRIX(rtanorm_combo, double, 
          SCYTHE_ARGSET(mean, variance, above, iter), double mean, 
          double variance, double above, unsigned int iter = 10);

      /* Multivariate Distributions */
      
      template <matrix_order O, matrix_style S>
      Matrix<double, O, Concrete>
      rwish(unsigned int v, const Matrix<double, O, S> &Sigma)
      {
        SCYTHE_CHECK_10(! Sigma.isSquare(), scythe_dimension_error,
            "Sigma not square");
        SCYTHE_CHECK_10(v < Sigma.rows(), scythe_invalid_arg, 
            "v < Sigma.rows()");
          
        Matrix<double,O,Concrete> 
          A(Sigma.rows(), Sigma.rows());
        Matrix<double,O,Concrete> C = cholesky<O,Concrete>(Sigma);
        Matrix<double,O,Concrete> alpha;
          
        for (unsigned int i = 0; i < v; ++i) {
          alpha = C * rnorm(Sigma.rows(), 1, 0, 1);
          A += (alpha * (t(alpha)));
        }

        return A;
      }

      template <matrix_order O, matrix_style S>
      Matrix<double, O, Concrete>
      rdirich(const Matrix<double, O, S>& alpha) 
      { 
        // Check for allowable parameters
        SCYTHE_CHECK_10(std::min(alpha) <= 0, scythe_invalid_arg,
            "alpha has elements < 0");
        SCYTHE_CHECK_10(! alpha.isColVector(), scythe_dimension_error,
            "alpha not column vector");
     
        Matrix<double, O, Concrete> y(alpha.rows(), 1);
        double ysum = 0;

        // We would use std::transform here but rgamma is a function
        // and wouldn't get inlined.
        const_matrix_forward_iterator<double,O,O,S> ait;
        const_matrix_forward_iterator<double,O,O,S> alast
          = alpha.template end_f();
        typename Matrix<double,O,Concrete>::forward_iterator yit 
          = y.begin_f();
        for (ait = alpha.begin_f(); ait != alast; ++ait) {
          *yit = rgamma(*ait, 1);
          ysum += *yit;
          ++yit;
        }

        y /= ysum;

        return y;
      }

      template <matrix_order PO1, matrix_style PS1,
                matrix_order PO2, matrix_style PS2>
      Matrix<double, PO1, Concrete>
      rmvnorm(const Matrix<double, PO1, PS1>& mu, 
              const Matrix<double, PO2, PS2>& sigma)
      {  
        unsigned int dim = mu.rows();
        SCYTHE_CHECK_10(! mu.isColVector(), scythe_dimension_error,
            "mu not column vector");
        SCYTHE_CHECK_10(! sigma.isSquare(), scythe_dimension_error,
            "sigma not square");
        SCYTHE_CHECK_10(sigma.rows() != dim, scythe_conformation_error,
            "mu and sigma not conformable");
        
        return(mu + cholesky(sigma) * rnorm(dim, 1, 0, 1));
      }

      template <matrix_order O, matrix_style S>
      Matrix<double, O, Concrete>
      rmvt (const Matrix<double, O, S>& sigma, double nu)
      {
        Matrix<double, O, Concrete> result;
        SCYTHE_CHECK_10(nu <= 0, scythe_invalid_arg,
            "D.O.F parameter nu <= 0");

        result = 
          rmvnorm(Matrix<double, O>(sigma.rows(), 1, true, 0), sigma);
        result /= std::sqrt(rchisq(nu) / nu);
        return result;
      }

    protected:
      /* Default (and only) constructor */
      rng() 
        : rnorm_count_ (1) // Initialize the normal counter
      {}

      /* For Barton and Nackman trick. */
      RNGTYPE& as_derived()
      {
        return static_cast<RNGTYPE&>(*this);
      }


      /* Generate Standard Normal variates */

      /* These instance variables were static in the old
       * implementation.  Making them instance variables provides
       * thread safety, as long as two threads don't access the same
       * rng at the same time w/out precautions.  Fixes possible
       * previous issues with lecuyer.  See the similar approach in
       * rgamma1 below.
       */
      int rnorm_count_;
      double x2_;

      double
      rnorm1 ()
      {
        double nu1, nu2, rsquared, sqrt_term;
        if (rnorm_count_ == 1){ // odd numbered passses
          do {
            nu1 = -1 +2*runif();
            nu2 = -1 +2*runif();
            rsquared = ::pow(nu1,2) + ::pow(nu2,2);
          } while (rsquared >= 1 || rsquared == 0.0);
          sqrt_term = std::sqrt(-2*std::log(rsquared)/rsquared);
          x2_ = nu2*sqrt_term;
          rnorm_count_ = 2;
          return nu1*sqrt_term;
        } else { // even numbered passes
          rnorm_count_ = 1;
          return x2_;
        } 
      }

      /* Generate standard gamma variates */
      double accept_;

      double
      rgamma1 (double alpha)
      {
        int test;
        double u, v, w, x, y, z, b, c;

        // Check for allowable parameters
        SCYTHE_CHECK_10(alpha <= 1, scythe_invalid_arg, "alpha <= 1");

        // Implement Best's (1978) simulator
        b = alpha - 1;
        c = 3 * alpha - 0.75;
        test = 0;
        while (test == 0) {
          u = runif ();
          v = runif ();

          w = u * (1 - u);
          y = std::sqrt (c / w) * (u - .5);
          x = b + y;

          if (x > 0) {
            z = 64 * std::pow (v, 2) * std::pow (w, 3);
            if (z <= (1 - (2 * std::pow (y, 2) / x))) {
              test = 1;
              accept_ = x;
            } else if ((2 * (b * std::log (x / b) - y)) >= ::log (z)) {
              test = 1;
              accept_ = x;
            } else {
              test = 0;
            }
          }
        }
        
        return (accept_);
      }

  };

  
} // end namespace scythe    
#endif /* RNG_H */