rng.h

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

      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);

      /*! \brief Generate a Weibull distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * Weibull distribution with given \a shape and \a scale.
       *
       * \param shape The strictly positive shape of the distribution.
			 * \param scale The strictly positive scale of the distribution.
			 * 
			 * \see pweibull(double x, double shape, double scale)
			 * \see dweibull(double x, double shape, double scale)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

			/*! \brief Generate an inverse \f$\chi^2\f$ distributed random
			 * variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * inverse \f$\chi^2\f$ distribution with \a nu degress of
			 * freedom.
       *
       * \param nu The degrees of freedom.
			 *
			 * \see rchisq(double df)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

      /*! \brief Generate an inverse gamma distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * inverse gamma distribution with given \a shape and \a scale.
       *
       * \param shape The strictly positive shape of the distribution.
			 * \param scale The strictly positive scale of the distribution.
			 * 
			 * \see rgamma(double alpha, double beta)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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 */

			/*! \brief Generate a truncated normally distributed random
			 * variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * normal distribution with given \a mean and \a variance,
			 * truncated both above and below.  It uses the inverse CDF
			 * method.
       *
       * \param mean The mean of the distribution.
			 * \param variance The variance of the distribution.
			 * \param below The lower truncation point of the distribution.
			 * \param above The upper truncation point of the distribution.
			 * 
			 * \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 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 
      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);

			/*! \brief Generate a truncated normally distributed random
			 * variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * normal distribution with given \a mean and \a variance,
			 * truncated both above and below.  It uses a combination of
			 * rejection sampling (when \a below <= mean <= \a above)
			 * sampling method of Robert and Casella (1999), pp. 288-289
			 * (when \a meam < \a below or \a mean > \a above).
       *
       * \param mean The mean of the distribution.
			 * \param variance The variance of the distribution.
			 * \param below The lower truncation point of the distribution.
			 * \param above The upper truncation point of the distribution.
			 * 
			 * \see rtnorm(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 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
      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);

			/*! \brief Generate a normally distributed random variate,
			 * truncated below.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * normal distribution with given \a mean and \a variance,
			 * truncated below.  It uses the slice sampling method of
			 * Robert and Casella (1999), pp. 288-289.
       *
			 * \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 use.
			 * 
			 * \see rtnorm(double mean, double variance, double below, double above)
			 * \see rtnorm_combo(double mean, double variance, double below, double above)
			 * \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 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_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);

			/*! \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 the slice sampling method of Robert
			 * and Casella (1999), pp. 288-289.
       *
       * \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 use.
			 * 
			 * \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 rtbnorm_combo(double mean, double variance, double below, 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
      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);

			/*! \brief Generate a normally distributed random
			 * variate, truncated below.
       *