rng.h

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

			 * This function returns a pseudo-random variate drawn from the
			 * normal distribution with given \a mean and \a variance,
			 * truncated below.  It uses a combination of
			 * rejection sampling (when \a mean >= \a below) and the slice
			 * sampling method of Robert and Casella (1999), pp. 288-289
			 * (when \a mean < \a below).
       *
       * \param mean The mean of the distribution.
			 * \param variance The variance of the distribution.
			 * \param below The lower truncation point of the distribution.
			 * \param iter The number of iterations to run the slice
			 * sampler.
			 * 
			 * \see rtnorm(double mean, double variance, double below, double above)
			 * \see rtnorm_combo(double mean, double variance, double below, double above)
			 * \see rtbnorm_slice(double mean, double variance, double below, unsigned int iter = 10)
			 * \see rtanorm_slice(double mean, double variance, double above, unsigned int iter = 10)
			 * \see rtanorm_combo(double mean, double variance, double above, unsigned int iter = 10)
			 * \see rnorm(double x, double mean, double sd)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

			/*! \brief Generate a normally distributed random variate,
			 * truncated above.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * normal distribution with given \a mean and \a variance,
			 * truncated above.  It uses a combination of rejection sampling
			 * (when \a mean <= \a above) and the slice sampling method of
			 * Robert and Casella (1999), pp. 288-289 (when \a mean > \a
			 * above).
       *
			 * \param mean The mean of the distribution.
			 * \param variance The variance of the distribution.
			 * \param above The upper truncation point of the distribution.
			 * \param iter The number of iterations to run the slice sampler.
			 * 
			 * \see rtnorm(double mean, double variance, double below, double above)
			 * \see rtnorm_combo(double mean, double variance, double below, double above)
			 * \see rtbnorm_slice(double mean, double variance, double below, unsigned int iter = 10)
			 * \see rtanorm_slice(double mean, double variance, double above, unsigned int iter = 10)
			 * \see rtbnorm_combo(double mean, double variance, double below, unsigned int iter = 10)
			 * \see rnorm(double x, double mean, double sd)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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 */
      
      /*! \brief Generate a Wishart distributed random variate Matrix.
       *
       * This function returns a pseudo-random matrix-valued variate
       * drawn from the Wishart disribution described by the scale
       * matrix \a Sigma, with \a v degrees of freedom.
       *
       * \param v The degrees of freedom of the distribution.
			 * \param Sigma The square scale matrix of the distribution.
			 *
       * \throw scythe_invalid_arg (Level 1)
       * \throw scythe_dimension_error (Level 1)
       */
      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 = A + (alpha * (t(alpha)));
        }

        return A;
      }

      /*! \brief Generate a Dirichlet distributed random variate Matrix.
       *
       * This function returns a pseudo-random matrix-valued variate
       * drawn from the Dirichlet disribution described by the vector
       * \a alpha.
       *
       * \param alpha A vector of non-negative reals.
			 *
       * \throw scythe_invalid_arg (Level 1)
       * \throw scythe_dimension_error (Level 1)
       */
      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;
          ++ait;
        }

        y /= ysum;

        return y;
      }

      /*! \brief Generate a multivariate normal distributed random
       * variate Matrix.
       *
       * This function returns a pseudo-random matrix-valued variate
       * drawn from the multivariate normal disribution with means \mu
       * and variance-covariance matrix \a sigma.
       *
       * \param mu A vector containing the distribution means.
       * \param sigma The distribution variance-covariance matrix.
			 *
       * \throw scythe_invalid_arg (Level 1)
       * \throw scythe_dimension_error (Level 1)
       */
      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));
      }

      /*! \brief Generate a multivariate Student t distributed random
       * variate Matrix.
       *
       * This function returns a pseudo-random matrix-valued variate
       * drawn from the multivariate Student t disribution with
       * and variance-covariance matrix \a sigma, and degrees of
       * freedom \a nu
       *
       * \param sigma The distribution variance-covariance matrix.
       * \param nu The strictly positive degrees of freedom.
			 *
       * \throw scythe_invalid_arg (Level 1)
       * \throw scythe_dimension_error (Level 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 */
      /*! \brief Default constructor
       *
       * Instantiate a random number generator
       */
      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 */