rng.h

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

/* 
 * 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
 */

/*!
 * \file rng.h
 *
 * \brief The definition of the random number generator base class.
 *
 */

/* 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);                \
  }

   /*! \brief Random number generator.
    *
    * This class provides objects capable of generating random numbers
    * from a variety of probability distributions.  This
    * abstract class forms the foundation of random number generation in
    * Scythe.  Specific random number generators should extend this class
    * and implement the virtual void function runif(); this function
    * should take no arguments and return uniformly distributed random
    * numbers on the interval (0, 1).  The rng class provides no
    * interface for seed-setting or initialization, allowing for maximal
    * flexibility in underlying implementation.  This class does provide
    * implementations of functions that return random numbers from a wide
    * variety of commonly (and not-so-commonly) used distributions, by
    * manipulating the uniform variates returned by runif().  See
    * rng/mersenne.h and rng/lecuyer.h for the rng implementations
    * offered by Scythe.
    *
    * Each univariate distribution is represented by three overloaded
    * versions of the same method.  The first is a simple method
    * returning a single value.  The remaining method versions return
    * Matrix values and are equivalent to calling the single-valued
    * method multiple times to fill a Matrix object.  They each take
    * two arguments describing the number of rows and columns in the
    * returned Matrix object and as many subsequent arguments as is
    * necessary to describe the distribution.  As is the case
    * throughout the library, the Matrix-returning versions of the
    * method include both a general and default template.  We
    * explicitly document only the single-valued versions of the
    * univariate methods.  For matrix-valued distributions we provide
    * only a single method per distribution.
    *
    * \note Doxygen incorrectly parses the macros we use to
    * automatically generate the Matrix returning versions of the
    * various univariate methods in this class.  Whenever you see the
    * macro variable __VA_ARGS__ in the public member function list
    * below, simply substitute in the arguments in the explicitly
    * documented single-valued version of the method.
    *
    */
  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.
       */
       /*! \brief Generate uniformly distributed random variates.
        *
        * This operator acts as an alias for runif() and generates
        * pseudo-random variates from the uniform distribution on the
        * interval (0, 1).  We include this operator to allow rng
        * objects to behave as function objects.
        */
      double operator() ()
      {
        return runif();
      }

      /* Returns random uniform numbers on (0, 1).  This function must
       * be implemented by extending classes */
       /*! \brief Generate uniformly distributed random variates.
        *
        * This method generates pseudo-random variates from the
        * uniform distribution on the interval (0, 1).
        *
        * This function is pure virtual and is implemented by
        * extending concrete classes, like scythe::mersenne and
        * scythe::lecuyer.
        */
      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);
      }

      /*! \brief Generate a beta distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * beta distribution described by the shape parameters \a a and
			 * \a b.
       *
       * \param alpha The first positive beta shape parameter.
       * \param beta the second positive beta shape parameter.
			 * 
			 * \see pbeta(double x, double a, double b)
			 * \see dbeta(double x, double a, double b)
			 * \see betafn(double a, double b)
			 * \see lnbetafn(double a, double b)
			 *
       * \throw scythe_invalid_arg (Level 1)
       */
      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);

      /*! \brief Generate a non-central hypergeometric disributed
			 * random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * non-centrial hypergeometric distribution described by the
       * number of positive outcomes \a m1, the two group size
       * parameters \a n1 and \a n2, and the odds ratio \a psi.
       *
       * \param m1 The number of positive outcomes in both groups.
       * \param n1 The size of group one.
       * \param n2 The size of group two.
       * \param psi The odds ratio
       * \param delta The precision.
       *
       * \throw scythe_convergence_error (Level 0)
       */
      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);

      /*! \brief Generate a Bernoulli distributed random variate.
       *
			 * This function returns a pseudo-random variate drawn from the
			 * Bernoulli distribution with probability of success \a p.
       *
       * \param p The probability of success on a trial.
			 * 
       * \throw scythe_invalid_arg (Level 1)
       */
      unsigned int
      rbern (double p)
      {
        unsigned int report;
        double unif;