distributions.h

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

   * \see rng::rbeta(double a, double b)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  betafn(double a, double b)
  {
    SCYTHE_CHECK_10(a <= 0 || b <= 0, scythe_invalid_arg, "a or b < 0");

    if (a + b < 171.61447887182298) /* ~= 171.61 for IEEE */
      return gammafn(a) * gammafn(b) / gammafn(a+b);

    double val = lnbetafn(a, b);
    SCYTHE_CHECK_10(val < -708.39641853226412, scythe_range_error,
        "Underflow");
    
    return std::exp(val);
  }

  /* The natural log of the beta function */
  /*! \brief The natural log of the beta function.
   *
   * Computes the natural log of the beta function, 
   * evaluated at (\a a, \a b).
   *
   * \param a The first parameter.
   * \param b The second parameter.
   *
   * \see betafn(double a, double b)
   * \see pbeta(double x, double a, double b)
   * \see dbeta(double x, double a, double b)
   * \see rng::rbeta(double a, double b)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  lnbetafn (double a, double b)
  {
    double p, q;

    p = q = a;
    if(b < p) p = b;/* := min(a,b) */
    if(b > q) q = b;/* := max(a,b) */

    SCYTHE_CHECK_10(p <= 0 || q <= 0,scythe_invalid_arg, "a or b <= 0");

    if (p >= 10) {
      /* p and q are big. */
      double corr = lngammacor(p) + lngammacor(q) - lngammacor(p + q);
      return std::log(q) * -0.5 + M_LN_SQRT_2PI + corr
        + (p - 0.5) * std::log(p / (p + q)) + q * std::log(1 + (-p / (p + q)));
    } else if (q >= 10) {
      /* p is small, but q is big. */
      double corr = lngammacor(q) - lngammacor(p + q);
      return lngammafn(p) + corr + p - p * std::log(p + q)
        + (q - 0.5) * std::log(1 + (-p / (p + q)));
    }
    
    /* p and q are small: p <= q > 10. */
    return std::log(gammafn(p) * (gammafn(q) / gammafn(p + q)));
  }

  /* Compute the factorial of a non-negative integer */
  /*! \brief The factorial function.
   *
   * Computes the factorial of \a n.
   *
   * \param n The non-negative integer value to compute the factorial of.
   *
   * \see lnfactorial(unsigned int n)
   *
   */
  inline int
  factorial (unsigned int n)
  {
    if (n == 0)
      return 1;

    return n * factorial(n - 1);
  }

  /* Compute the natural log of the factorial of a non-negative
   * integer
   */
  /*! \brief The log of the factorial function.
   *
   * Computes the natural log of the factorial of \a n.
   *
   * \param n The non-negative integer value to compute the natural log of the factorial of.
   *
   * \see factorial(unsigned int n)
   *
   */
  inline double
  lnfactorial (unsigned int n)
  {
    double x = n+1;
    double cof[6] = {
      76.18009172947146, -86.50532032941677,
      24.01409824083091, -1.231739572450155,
      0.1208650973866179e-2, -0.5395239384953e-5
    };
    double y = x;
    double tmp = x + 5.5 - (x + 0.5) * std::log(x + 5.5);
    double ser = 1.000000000190015;
    for (int j = 0; j <= 5; j++) {
      ser += (cof[j] / ++y);
    }
    return(std::log(2.5066282746310005 * ser / x) - tmp);
  }

  /*********************************
   * Fully Specified Distributions * 
   *********************************/

  /* These macros provide a nice shorthand for the matrix versions of
   * the pdf and cdf functions.
   */
 
#define SCYTHE_ARGSET(...) __VA_ARGS__

#define SCYTHE_DISTFUN_MATRIX(NAME, XTYPE, ARGNAMES, ...)             \
  template <matrix_order RO, matrix_style RS,                         \
            matrix_order PO, matrix_style PS>                         \
  Matrix<double, RO, RS>                                              \
  NAME (const Matrix<XTYPE, PO, PS>& X, __VA_ARGS__)                  \
  {                                                                   \
    Matrix<double, RO, Concrete> ret(X.rows(), X.cols(), false);      \
    const_matrix_forward_iterator<XTYPE,RO,PO,PS> xit;                \
    const_matrix_forward_iterator<XTYPE,RO,PO,PS> xlast               \
      = X.template end_f<RO>();                                       \
    typename Matrix<double,RO,Concrete>::forward_iterator rit         \
      = ret.begin_f();                                                \
    for (xit = X.template begin_f<RO>(); xit != xlast; ++xit) {       \
      *rit = NAME (*xit, ARGNAMES);                                   \
      ++rit;                                                          \
    }                                                                 \
    SCYTHE_VIEW_RETURN(double, RO, RS, ret)                           \
  }                                                                   \
                                                                      \
  template <matrix_order O, matrix_style S>                           \
  Matrix<double, O, Concrete>                                         \
  NAME (const Matrix<XTYPE, O, S>& X, __VA_ARGS__)                    \
  {                                                                   \
    return NAME <O, Concrete, O, S> (X, ARGNAMES);                    \
  }

  /**** The Beta Distribution ****/

  /* CDFs */

  /*! \brief The beta distribution function.
   *
   * Computes the value of the beta cumulative distribution function
   * with shape parameters \a a and \a b at the desired quantile,
   * \a x.
   *
   * It is also possible to call this function with a Matrix of
   * doubles as its first argument.  In this case the function will
   * return a Matrix of doubles of the same dimension as \a x,
   * containing the result of evaluating this function at each value
   * in \a x, given the remaining fixed parameters.  By default, the
   * returned Matrix will be concrete and have the same matrix_order
   * as \a x, but you may invoke a generalized version of the function
   * with an explicit template call.
   *
   * \param x The desired quantile, between 0 and 1.
   * \param a The first non-negative beta shape parameter.
   * \param b The second non-negative beta shape parameter.
   *
   * \see dbeta(double x, double a, double b)
   * \see rng::rbeta(double a, double b)
   * \see betafn(double a, double b)
   * \see lnbetafn(double a, double b)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  pbeta(double x, double a, double b)
  {
    SCYTHE_CHECK_10(a <= 0 || b <= 0,scythe_invalid_arg, "a or b <= 0");
    
    if (x <= 0)
      return 0.;
    if (x >= 1)
      return 1.;
    
    return pbeta_raw(x,a,b);
  }

  SCYTHE_DISTFUN_MATRIX(pbeta, double, SCYTHE_ARGSET(a, b), double a, double b)

  /* PDFs */
  /*! \brief The beta density function.
   *
   * Computes the value of the beta probability density function
   * with shape parameters \a a and \a b at the desired quantile,
   * \a x.
   *
   * It is also possible to call this function with a Matrix of
   * doubles as its first argument.  In this case the function will
   * return a Matrix of doubles of the same dimension as \a x,
   * containing the result of evaluating this function at each value
   * in \a x, given the remaining fixed parameters.  By default, the
   * returned Matrix will be concrete and have the same matrix_order
   * as \a x, but you may invoke a generalized version of the function
   * with an explicit template call.
   *
   * \param x The desired quantile, between 0 and 1.
   * \param a The first non-negative beta shape parameter.
   * \param b The second non-negative beta shape parameter.
   *
   * \see pbeta(double x, double a, double b)
   * \see rng::rbeta(double a, double b)
   * \see betafn(double a, double b)
   * \see lnbetafn(double a, double b)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  dbeta(double x, double a, double b)
  {
    SCYTHE_CHECK_10((x < 0.0) || (x > 1.0), scythe_invalid_arg,
        "x not in [0,1]");
    SCYTHE_CHECK_10(a < 0.0, scythe_invalid_arg, "a < 0");
    SCYTHE_CHECK_10(b < 0.0, scythe_invalid_arg, "b < 0");

    return (std::pow(x, (a-1.0)) * std::pow((1.0-x), (b-1.0)) )
      / betafn(a,b);
  }

  SCYTHE_DISTFUN_MATRIX(dbeta, double, SCYTHE_ARGSET(a, b), double a, double b)

  /* Returns the natural log of the ordinate of the Beta density
   * evaluated at x with Shape1 a, and Shape2 b
   */

  /*! \brief The natural log of the ordinate of the beta density
   * function.
   *
   * Computes the value of the natural log of the ordinate of the beta
   * probability density function
   * with shape parameters \a a and \a b at the desired quantile,
   * \a x.
   *
   * It is also possible to call this function with a Matrix of
   * doubles as its first argument.  In this case the function will
   * return a Matrix of doubles of the same dimension as \a x,
   * containing the result of evaluating this function at each value
   * in \a x, given the remaining fixed parameters.  By default, the
   * returned Matrix will be concrete and have the same matrix_order
   * as \a x, but you may invoke a generalized version of the function
   * with an explicit template call.
   *
   * \param x The desired quantile, between 0 and 1.
   * \param a The first non-negative beta shape parameter.
   * \param b The second non-negative beta shape parameter.
   *
   * \see dbeta(double x, double a, double b)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  lndbeta1(double x, double a, double b)
  { 
    SCYTHE_CHECK_10((x < 0.0) || (x > 1.0), scythe_invalid_arg,
        "x not in [0,1]");
    SCYTHE_CHECK_10(a < 0.0, scythe_invalid_arg, "a < 0");
    SCYTHE_CHECK_10(b < 0.0, scythe_invalid_arg, "b < 0");

    return (a-1.0) * std::log(x) + (b-1) * std::log(1.0-x)
      - lnbetafn(a,b);
  }

  SCYTHE_DISTFUN_MATRIX(lndbeta1, double, SCYTHE_ARGSET(a, b), double a, double b)


  /**** The Binomial Distribution ****/

  /* CDFs */

  /*! \brief The binomial distribution function.
   *
   * Computes the value of the binomial cumulative distribution function
   * with \a n trials and \a p probability of success on each trial,
   * at the desired quantile \a x.
   *
   * It is also possible to call this function with a Matrix of
   * doubles as its first argument.  In this case the function will
   * return a Matrix of doubles of the same dimension as \a x,
   * containing the result of evaluating this function at each value
   * in \a x, given the remaining fixed parameters.  By default, the
   * returned Matrix will be concrete and have the same matrix_order
   * as \a x, but you may invoke a generalized version of the function
   * with an explicit template call.
   *
   * \param x The desired quantile.
   * \param n The number of trials.
   * \param p The probability of success on each trial.
   *
   * \see dbinom(double x, unsigned int n, double p)
   * \see rng::rbinom(unsigned int n, double p)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  pbinom(double x, unsigned int n, double p)
  {
      
    SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg, "p not in [0,1]");
    double X = std::floor(x);
      
    if (X < 0.0)
      return 0;
    
    if (n <= X)
      return 1;
      
    return pbeta(1 - p, n - X, X + 1);
  }

  SCYTHE_DISTFUN_MATRIX(pbinom, double, SCYTHE_ARGSET(n, p), unsigned int n, double p)

  /* PDFs */
  /*! \brief The binomial density function.
   *
   * Computes the value of the binomial probability density function
   * with \a n trials and \a p probability of success on each trial,
   * at the desired quantile \a x.
   *
   * It is also possible to call this function with a Matrix of
   * doubles as its first argument.  In this case the function will
   * return a Matrix of doubles of the same dimension as \a x,
   * containing the result of evaluating this function at each value
   * in \a x, given the remaining fixed parameters.  By default, the
   * returned Matrix will be concrete and have the same matrix_order
   * as \a x, but you may invoke a generalized version of the function
   * with an explicit template call.
   *
   * \param x The desired quantile.
   * \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 rng::rbinom(unsigned int n, double p)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  dbinom(double x, unsigned int n, double p)
  {
    SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg, "p not in [0, 1]");
    double X = std::floor(x + 0.5);
    return dbinom_raw(X, n, p, 1 - p);
  }

  SCYTHE_DISTFUN_MATRIX(dbinom, double, SCYTHE_ARGSET(n, p), unsigned int n, double p)

  /**** The Chi Squared Distribution ****/
  
  /* CDFs */
  /*! \brief The \f$\chi^2\f$ distribution function.
   *
   * Computes the value of the \f$\chi^2\f$ cumulative distribution
   * function with \a df degrees of freedom, at the desired quantile
   * \a x.
   *
   * It is also possible to call this function with a Matrix of
   * doubles as its first argument.  In this case the function will
   * return a Matrix of doubles of the same dimension as \a x,
   * containing the result of evaluating this function at each value
   * in \a x, given the remaining fixed parameters.  By default, the
   * returned Matrix will be concrete and have the same matrix_order
   * as \a x, but you may invoke a generalized version of the function
   * with an explicit template call.
   *
   * \param x The desired quantile.
   * \param df The degrees of freedom.

   * \see dchisq(double x, double df)
   * \see rng::rchisq(double df)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   * \throw scythe_convergence_error (Level 1)