distributions.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/distributions.h
 *
 */

/*!  \file distributions.h
 *
 * \brief Definitions for probability density functions
 * (PDFs), cumulative distribution functions (CDFs), and some common
 * functions (gamma, beta, etc).
 *
 * This file provides functions that evaluate the PDFs and CDFs of a
 * number of probability distributions.  In addition, it includes
 * definitions for another of related functions, such as the gamma
 * and beta functions.
 *
 * The various distribution functions in this file operate on both
 * scalar quantiles and matrices of quantiles and the
 * definitions of both forms of these functions appear below.  We
 * provide explicit documentation only for the scalar versions of the
 * these functions and describe the Matrix versions in the scalar
 * calls' documents.  Much like the operators in matrix.h, we
 * implement these overloaded versions of the distribution functions
 * in terms of both generalized and default templates to allow for
 * explicit control over the template type of the returned Matrix.
 *
 * \note Doxygen does not correctly expand the macro definitions we use
 * to generate the Matrix versions of the various distribution
 * functions.  Therefore, it incorrectly substitutes the macro
 * variable
 * __VA_ARGS__ for the actual parameter values in the parameter lists
 * of each of these functions.  For example, the definitions of the
 * Matrix versions of pbeta are listed as
 * \code
 * template<matrix_order RO, matrix_style RS, matrix_order PO, matrix_style PS>
 * Matrix<double, RO, RS> scythe::pbeta (const Matrix<double, PO, PS> &X, __VA_ARGS__)
 *
 * template<matrix_order O, matrix_style S>
 * Matrix<double, O, Concrete> scythe::pbeta (const Matrix<double, O, S> &X, __VA_ARGS__)
 * \endcode
 * when they should be
 * \code
 * template<matrix_order RO, matrix_style RS, matrix_order PO, matrix_style PS>
 * Matrix<double, RO, RS> scythe::pbeta (const Matrix<double, PO, PS> &X, double a, double b)
 *
 * template<matrix_order O, matrix_style S>
 * Matrix<double, O, Concrete> scythe::pbeta (const Matrix<double, O, S> &X, double a, double b)
 * \endcode
 *
 * \par
 * Furthermore, Doxygen erroneously lists a number of variables at the
 * end of this document that are not, in fact, declared in
 * distributions.h.  Again, this error is a result of Doxygen's macro
 * parsing capabilities.
 * 
 */

/* TODO: Figure out how to get doxygen to stop listing erroneous
 * variables at the end of the doc for this file.  They stem from it
 * misreading the nested macro calls used to generate matrix procs.
 */

/* TODO: Full R-style versions of these function including arbitrary
 * recycling of matrix arguments.  This is going to have to wait for
 * variadic templates to be doable without a complete mess.  There is
 * currently a variadic patch available for g++, perhaps we can add a
 * SCYTHE_VARIADIC flag and include these as option until they become
 * part of the standard in 2009.  Something to come back to after 1.0.
 */

#ifndef SCYTHE_DISTRIBUTIONS_H
#define SCYTHE_DISTRIBUTIONS_H

#include <iostream>
#include <cmath>
#include <cfloat>
#include <climits>
#include <algorithm>
#include <limits>

#ifdef HAVE_IEEEFP_H
#include <ieeefp.h>
#endif

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "ide.h"
#include "error.h"
#else
#include "scythestat/matrix.h"
#include "scythestat/ide.h"
#include "scythestat/error.h"
#endif

/* Fill in some defs from R that aren't in math.h */
#ifndef M_PI
#define M_PI 3.141592653589793238462643383280
#endif
#define M_LN_SQRT_2PI 0.918938533204672741780329736406
#define M_LN_SQRT_PId2  0.225791352644727432363097614947
#define M_1_SQRT_2PI  0.39894228040143267793994605993
#define M_2PI   6.28318530717958647692528676655
#define M_SQRT_32 5.656854249492380195206754896838

#ifndef HAVE_TRUNC
/*! @cond */
inline double trunc(double x) throw ()
{
    if (x >= 0) 
      return std::floor(x);
    else
      return std::ceil(x);
}
/*! @endcond */
#endif

/* Many random number generators, pdfs, cdfs, and functions (gamma,
 * etc) in this file are based on code from 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
 */

namespace scythe {
  
  /*! @cond */
  
  /* Forward declarations */
  double gammafn (double);
  double lngammafn (double);
  double lnbetafn (double, double);
  double pgamma (double, double, double);
  double dgamma(double, double, double);
  double pnorm (double, double, double);

  /*! @endcond */

  /********************
   * Helper Functions *
   ********************/
  namespace {

    /* Evaluate a Chebysheve series at a given point */
    double
    chebyshev_eval (double x, const double *a, int n)
    {
      SCYTHE_CHECK_10(n < 1 || n > 1000, scythe_invalid_arg,
          "n not on [1, 1000]");
  
      SCYTHE_CHECK_10(x < -1.1 || x > 1.1, scythe_invalid_arg,
          "x not on [-1.1, 1.1]");
      
      double b0, b1, b2;
      b0 = b1 = b2 = 0;
  
      double twox = x * 2;
  
      for (int i = 1; i <= n; ++i) {
        b2 = b1;
        b1 = b0;
        b0 = twox * b1 - b2 + a[n - i];
      }
  
      return (b0 - b2) * 0.5;
    }

    /* Computes the log gamma correction factor for x >= 10 */
    double
    lngammacor(double x)
    {
      const double algmcs[15] = {
        +.1666389480451863247205729650822e+0,
        -.1384948176067563840732986059135e-4,
        +.9810825646924729426157171547487e-8,
        -.1809129475572494194263306266719e-10,
        +.6221098041892605227126015543416e-13,
      };
    
      SCYTHE_CHECK_10(x < 10, scythe_invalid_arg,
          "This function requires x >= 10");  
      SCYTHE_CHECK_10(x >= 3.745194030963158e306, scythe_range_error,
          "Underflow");
      
      if (x < 94906265.62425156) {
        double tmp = 10 / x;
        return chebyshev_eval(tmp * tmp * 2 - 1, algmcs, 5) / x;
      }
      
      return 1 / (x * 12);
    }

    /* Evaluates the "deviance part" */
    double
    bd0(double x, double np)
    {
      
      if(std::fabs(x - np) < 0.1 * (x + np)) {
        double v = (x - np) / (x + np);
        double s = (x - np) * v;
        double ej = 2 * x * v;
        v = v * v;
        for (int j = 1; ; j++) {
          ej *= v;
          double s1 = s + ej / ((j << 1) + 1);
          if (s1 == s)
            return s1;
          s = s1;
        }
      }
      
      return x * std::log(x / np) + np - x;
    }
  
    /* Computes the log of the error term in Stirling's formula */
    double
    stirlerr(double n)
    {
#define S0 0.083333333333333333333       /* 1/12 */
#define S1 0.00277777777777777777778     /* 1/360 */
#define S2 0.00079365079365079365079365  /* 1/1260 */
#define S3 0.000595238095238095238095238 /* 1/1680 */
#define S4 0.0008417508417508417508417508/* 1/1188 */
      
      /* error for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0 */
      const double sferr_halves[31] = {
        0.0, /* n=0 - wrong, place holder only */
        0.1534264097200273452913848,  /* 0.5 */
        0.0810614667953272582196702,  /* 1.0 */
        0.0548141210519176538961390,  /* 1.5 */
        0.0413406959554092940938221,  /* 2.0 */
        0.03316287351993628748511048, /* 2.5 */
        0.02767792568499833914878929, /* 3.0 */
        0.02374616365629749597132920, /* 3.5 */
        0.02079067210376509311152277, /* 4.0 */
        0.01848845053267318523077934, /* 4.5 */
        0.01664469118982119216319487, /* 5.0 */
        0.01513497322191737887351255, /* 5.5 */
        0.01387612882307074799874573, /* 6.0 */
        0.01281046524292022692424986, /* 6.5 */
        0.01189670994589177009505572, /* 7.0 */
        0.01110455975820691732662991, /* 7.5 */
        0.010411265261972096497478567, /* 8.0 */
        0.009799416126158803298389475, /* 8.5 */
        0.009255462182712732917728637, /* 9.0 */
        0.008768700134139385462952823, /* 9.5 */
        0.008330563433362871256469318, /* 10.0 */
        0.007934114564314020547248100, /* 10.5 */
        0.007573675487951840794972024, /* 11.0 */
        0.007244554301320383179543912, /* 11.5 */
        0.006942840107209529865664152, /* 12.0 */
        0.006665247032707682442354394, /* 12.5 */
        0.006408994188004207068439631, /* 13.0 */
        0.006171712263039457647532867, /* 13.5 */
        0.005951370112758847735624416, /* 14.0 */
        0.005746216513010115682023589, /* 14.5 */
        0.005554733551962801371038690  /* 15.0 */
      };
      double nn;
      
      if (n <= 15.0) {
        nn = n + n;
        if (nn == (int)nn)
          return(sferr_halves[(int)nn]);
        return (scythe::lngammafn(n + 1.) - (n + 0.5) * std::log(n) + n -
            std::log(std::sqrt(2 * M_PI)));
      }
      
      nn = n*n;
      if (n > 500)
        return((S0 - S1 / nn) / n);
      if (n > 80)
        return((S0 - (S1 - S2 / nn) / nn) / n);
      if (n > 35)
        return((S0 - (S1 - (S2 - S3 / nn) / nn) / nn) / n);
      /* 15 < n <= 35 : */
      return((S0 - (S1 - (S2 - (S3 - S4 / nn) / nn) / nn) / nn) / n);
#undef S1
#undef S2
#undef S3
#undef S4
    }


    /* Helper for dpois and dgamma */
    double
    dpois_raw (double x, double lambda)
    {
      if (lambda == 0)
        return ( (x == 0) ? 1.0 : 0.0);

      if (x == 0)
        return std::exp(-lambda);

      if (x < 0)
        return 0.0;

      return std::exp(-stirlerr(x) - bd0(x, lambda))
        / std::sqrt(2 * M_PI * x);
    }

  
    /* helper for pbeta */
    double
    pbeta_raw(double x, double pin, double qin)
    {
      double ans, c, finsum, p, ps, p1, q, term, xb, xi, y;
      int n, i, ib, swap_tail;
      
      const double eps = .5 * DBL_EPSILON;
      const double sml = DBL_MIN;
      const double lneps = std::log(eps);
      const double lnsml = std::log(eps);
      
      if (pin / (pin + qin) < x) {
        swap_tail = 1;
        y = 1 - x;
        p = qin;
        q = pin;
      } else {
        swap_tail=0;
        y = x;
        p = pin;
        q = qin;
      }
      
      if ((p + q) * y / (p + 1) < eps) {
        ans = 0;
        xb = p * std::log(std::max(y,sml)) - std::log(p) - lnbetafn(p,q);
        if (xb > lnsml && y != 0)
          ans = std::exp(xb);
        if (swap_tail)
          ans = 1-ans;
      } else {
        ps = q - std::floor(q);
        if (ps == 0)
          ps = 1;
        xb = p * std::log(y) - lnbetafn(ps, p) - std::log(p);
        ans = 0;
        if (xb >= lnsml) {
          ans = std::exp(xb);
          term = ans * p;
          if (ps != 1) {
            n = (int)std::max(lneps/std::log(y), 4.0);
            for(i = 1; i <= n; i++){
              xi = i;
              term *= (xi-ps)*y/xi;
              ans += term/(p+xi);
            }
          }
        }
        if (q > 1) {
          xb = p * std::log(y) + q * std::log(1 - y)
            - lnbetafn(p, q) - std::log(q);
          ib = (int) std::max(xb / lnsml, 0.0);
          term = std::exp(xb - ib * lnsml);
          c = 1 / (1 - y);
          p1 = q * c / (p + q - 1);
              
          finsum = 0;
          n = (int) q;
          if(q == n)
            n--;
          for (i = 1; i <= n; i++) {
            if(p1 <= 1 && term / eps <= finsum)
              break;
            xi = i;
            term = (q -xi + 1) * c * term / (p + q - xi);
            if (term > 1) {
              ib--;
              term *= sml;
            }
            if (ib == 0)
              finsum += term;
          }
          ans += finsum;
        }
        
        if(swap_tail)
          ans = 1-ans;
        ans = std::max(std::min(ans,1.),0.);
      }
      return ans;
    }
  
   /* Helper for dbinom */
    double