rng.h

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

          
        // 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);
      
      /*! \brief Generate a binomial distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * binomial distribution with \a n trials and \p probability of
			 * success on each trial.
       *
       * \param n The number of trials.
       * \param p The probability of success on each trial.
			 * 
			 * \see pbinom(double x, unsigned int n, double p)
			 * \see dbinom(double x, unsigned int n, double p)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

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

      /*! \brief Generate an exponentially distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * exponential distribution described by the inverse scale
			 * parameter \a invscale.
       *
       * \param invscale The inverse scale parameter.
			 * 
			 * \see pexp(double x, double scale)
			 * \see dexp(double x, double scale)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);
    
      /*! \brief Generate an F distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * F distribution with degress of freedom \a df1 and \a df2.
       *
       * \param df1 The positive degrees of freedom for the
			 * \f$chi^2\f$ variate in the nominator of the F statistic.
       * \param df2 The positive degrees of freedom for the
			 * \f$chi^2\f$ variate in the denominator of the F statistic.
			 *
			 * \see pf(double x, double df1, double df2)
			 * \see df(double x, double df1, double df2)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

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

      /*! \brief Generate a logistically distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * logistic distribution described by the given \a location and
			 * \a scale variables.
       *
       * \param location The location of the distribution.
			 * \param scale The scale of the distribution.
			 * 
			 * \see plogis(double x, double location, double scale)
			 * \see dlogis(double x, double location, double scale)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

      /*! \brief Generate a log-normal distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * log-normal distribution with given logged mean and standard
			 * deviation.
       *
       * \param logmean The logged mean of the distribtion.
			 * \param logsd The strictly positive logged standard deviation
			 * of the distribution.
			 * 
			 * \see plnorm(double x, double logmean, double logsd)
			 * \see dlnorm(double x, double logmean, double logsd)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

			/*! \brief Generate a negative binomial distributed random
			 * variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * negative binomial distribution with given dispersion
			 * parameter and probability of success on each trial.
       *
       * \param n The strictly positive target number of successful
			 * trials (dispersion parameters).
			 * \param p The probability of success on each trial.
			 * 
			 * \see pnbinom(unsigned int x, double n, double p)
			 * \see dnbinom(unsigned int x, double n, double p)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

      /*! \brief Generate a normally distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * normal distribution with given \a mean and \a standard
			 * distribution.
       *
       * \param mean The mean of the distribution.
			 * \param sd The standard deviation of the distribution.
			 * 
			 * \see pnorm(double x, double mean, double sd)
			 * \see dnorm(double x, double mean, double sd)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

      /*! \brief Generate a Poisson distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * Poisson distribution with expected number of occurrences \a
			 * lambda.
       *
       * \param lambda The strictly positive expected number of
			 * occurrences.
			 * 
			 * \see ppois(double x, double lambda)
			 * \see dpois(double x, double lambda)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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.
       */

      /*! \brief Generate a Student t distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * Student's t distribution with given mean \a mu, variance \a
			 * sigma2, and degrees of freedom \a nu
       *
       * \param mu The mean of the distribution.
			 * \param sigma2 The variance of the distribution.
			 * \param nu The degrees of freedom of the distribution.
			 * 
			 * \see dt1(double x, double mu, double sigma2, double nu)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */