distributions.h

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

    dbinom_raw (double x, double n, double p, double q)
    { 
      double f, lc;

      if (p == 0)
        return((x == 0) ? 1.0 : 0.0);
      if (q == 0)
        return((x == n) ? 1.0 : 0.0);

      if (x == 0) { 
        if(n == 0)
          return 1.0;
        
        lc = (p < 0.1) ? -bd0(n, n * q) - n * p : n * std::log(q);
        return(std::exp(lc));
      }
      if (x == n) { 
        lc = (q < 0.1) ? -bd0(n,n * p) - n * q : n * std::log(p);
        return(std::exp(lc));
      }

      if (x < 0 || x > n)
        return 0.0;

      lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) -
        bd0(n - x, n * q);
      
      f = (M_2PI * x * (n-x)) / n;

      return (std::exp(lc) / std::sqrt(f));
    }

    /* The normal probability density function implementation. */

#define SIXTEN 16
#define do_del(X)              \
    xsq = trunc(X * SIXTEN) / SIXTEN;        \
    del = (X - xsq) * (X + xsq);          \
    if(log_p) {              \
        *cum = (-xsq * xsq * 0.5) + (-del * 0.5) + std::log(temp);  \
        if((lower && x > 0.) || (upper && x <= 0.))      \
        *ccum = ::log1p(-std::exp(-xsq * xsq * 0.5) *     \
          std::exp(-del * 0.5) * temp);    \
    }                \
    else {                \
        *cum = std::exp(-xsq * xsq * 0.5) * std::exp(-del * 0.5) * temp;  \
        *ccum = 1.0 - *cum;            \
    }

#define swap_tail            \
    if (x > 0.) {/* swap  ccum <--> cum */      \
        temp = *cum; if(lower) *cum = *ccum; *ccum = temp;  \
    }

    void
    pnorm_both(double x, double *cum, double *ccum, int i_tail,
                bool log_p)
    {
      const double a[5] = {
        2.2352520354606839287,
        161.02823106855587881,
        1067.6894854603709582,
        18154.981253343561249,
        0.065682337918207449113
      };
      const double b[4] = {
        47.20258190468824187,
        976.09855173777669322,
        10260.932208618978205,
        45507.789335026729956
      };
      const double c[9] = {
        0.39894151208813466764,
        8.8831497943883759412,
        93.506656132177855979,
        597.27027639480026226,
        2494.5375852903726711,
        6848.1904505362823326,
        11602.651437647350124,
        9842.7148383839780218,
        1.0765576773720192317e-8
      };
      const double d[8] = {
        22.266688044328115691,
        235.38790178262499861,
        1519.377599407554805,
        6485.558298266760755,
        18615.571640885098091,
        34900.952721145977266,
        38912.003286093271411,
        19685.429676859990727
      };
      const double p[6] = {
        0.21589853405795699,
        0.1274011611602473639,
        0.022235277870649807,
        0.001421619193227893466,
        2.9112874951168792e-5,
        0.02307344176494017303
      };
      const double q[5] = {
        1.28426009614491121,
        0.468238212480865118,
        0.0659881378689285515,
        0.00378239633202758244,
        7.29751555083966205e-5
      };
      
      double xden, xnum, temp, del, eps, xsq, y;
      int i, lower, upper;

      /* Consider changing these : */
      eps = DBL_EPSILON * 0.5;

      /* i_tail in {0,1,2} =^= {lower, upper, both} */
      lower = i_tail != 1;
      upper = i_tail != 0;

      y = std::fabs(x);
      if (y <= 0.67448975) {
        /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */
        if (y > eps) {
          xsq = x * x;
          xnum = a[4] * xsq;
          xden = xsq;
          for (i = 0; i < 3; ++i) {
            xnum = (xnum + a[i]) * xsq;
            xden = (xden + b[i]) * xsq;
          }
        } else xnum = xden = 0.0;
        
        temp = x * (xnum + a[3]) / (xden + b[3]);
        if(lower)  *cum = 0.5 + temp;
        if(upper) *ccum = 0.5 - temp;
        if(log_p) {
          if(lower)  *cum = std::log(*cum);
          if(upper) *ccum = std::log(*ccum);
        }
      } else if (y <= M_SQRT_32) {
        /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32) 
         * ~= 5.657 */

        xnum = c[8] * y;
        xden = y;
        for (i = 0; i < 7; ++i) {
          xnum = (xnum + c[i]) * y;
          xden = (xden + d[i]) * y;
        }
        temp = (xnum + c[7]) / (xden + d[7]);
        do_del(y);
        swap_tail;
      } else if (log_p
                || (lower && -37.5193 < x && x < 8.2924)
                || (upper && -8.2929 < x && x < 37.5193)
          ) {
        /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
        xsq = 1.0 / (x * x);
        xnum = p[5] * xsq;
        xden = xsq;
        for (i = 0; i < 4; ++i) {
          xnum = (xnum + p[i]) * xsq;
          xden = (xden + q[i]) * xsq;
        }
        temp = xsq * (xnum + p[4]) / (xden + q[4]);
        temp = (M_1_SQRT_2PI - temp) / y;
        do_del(x);
        swap_tail;
      } else {
        if (x > 0) {
          *cum = 1.;
          *ccum = 0.;
        } else {
          *cum = 0.;
          *ccum = 1.;
        }
        //XXX commented out for debug-on testing of ordfactanal
        //(and perhaps others) since they tend to throw on the first
        //iteration
        //SCYTHE_THROW_10(scythe_convergence_error, "Did not converge");
      }

      return;
    }
#undef SIXTEN
#undef do_del
#undef swap_tail

    /* The standard normal distribution function */
    double
    pnorm1 (double x, bool lower_tail, bool log_p)
    {
      SCYTHE_CHECK_10(! finite(x), scythe_invalid_arg,
          "Quantile x is inifinte (+/-Inf) or NaN");

      double p, cp;
      pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);

      return (lower_tail ? p : cp);
    }
  } // anonymous namespace

  /*************
   * Functions *
   *************/
  
  /* The gamma function */

  /*! \brief The gamma function.
   *
   * Computes the gamma function, evaluated at \a x.
   *
   * \param x The value to compute gamma at.
   *
   * \see lngammafn(double x)
   * \see pgamma(double x, double shape, double scale)
   * \see dgamma(double x, double shape, double scale)
   * \see rng::rgamma(double shape, double scale)
   *
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double 
  gammafn (double x)
  {
    const double gamcs[22] = {
      +.8571195590989331421920062399942e-2,
      +.4415381324841006757191315771652e-2,
      +.5685043681599363378632664588789e-1,
      -.4219835396418560501012500186624e-2,
      +.1326808181212460220584006796352e-2,
      -.1893024529798880432523947023886e-3,
      +.3606925327441245256578082217225e-4,
      -.6056761904460864218485548290365e-5,
      +.1055829546302283344731823509093e-5,
      -.1811967365542384048291855891166e-6,
      +.3117724964715322277790254593169e-7,
      -.5354219639019687140874081024347e-8,
      +.9193275519859588946887786825940e-9,
      -.1577941280288339761767423273953e-9,
      +.2707980622934954543266540433089e-10,
      -.4646818653825730144081661058933e-11,
      +.7973350192007419656460767175359e-12,
      -.1368078209830916025799499172309e-12,
      +.2347319486563800657233471771688e-13,
      -.4027432614949066932766570534699e-14,
      +.6910051747372100912138336975257e-15,
      -.1185584500221992907052387126192e-15,
    };


    double y = std::fabs(x);

    if (y <= 10) {

      /* Compute gamma(x) for -10 <= x <= 10
       * Reduce the interval and find gamma(1 + y) for 0 <= y < 1
       * first of all. */

      int n = (int) x;
      if (x < 0)
        --n;
      
      y = x - n;/* n = floor(x)  ==>  y in [ 0, 1 ) */
      --n;
      double value = chebyshev_eval(y * 2 - 1, gamcs, 22) + .9375;
      
      if (n == 0)
        return value;/* x = 1.dddd = 1+y */

      if (n < 0) {
        /* compute gamma(x) for -10 <= x < 1 */

        /* If the argument is exactly zero or a negative integer */
        /* then return NaN. */
        SCYTHE_CHECK_10(x == 0 || (x < 0 && x == n + 2),
            scythe_range_error, "x is 0 or a negative integer");

        /* The answer is less than half precision */
        /* because x too near a negative integer. */
        SCYTHE_CHECK_10(x < -0.5 && 
            std::fabs(x - (int)(x - 0.5) / x) < 67108864.0,
            scythe_precision_error,
            "Answer < 1/2 precision because x is too near" <<
            " a negative integer");

        /* The argument is so close to 0 that the result
         * * would overflow. */
        SCYTHE_CHECK_10(y < 2.2474362225598545e-308, scythe_range_error,
            "x too close to 0");

        n = -n;

        for (int i = 0; i < n; i++)
          value /= (x + i);
        
        return value;
      } else {
        /* gamma(x) for 2 <= x <= 10 */

        for (int i = 1; i <= n; i++) {
          value *= (y + i);
        }
        return value;
      }
    } else {
      /* gamma(x) for   y = |x| > 10. */

      /* Overflow */
      SCYTHE_CHECK_10(x > 171.61447887182298, 
          scythe_range_error,"Overflow");

      /* Underflow */
      SCYTHE_CHECK_10(x < -170.5674972726612,
          scythe_range_error, "Underflow");

      double value = std::exp((y - 0.5) * std::log(y) - y 
          + M_LN_SQRT_2PI + lngammacor(y));

      if (x > 0)
        return value;

      SCYTHE_CHECK_10(std::fabs((x - (int)(x - 0.5))/x) < 67108864.0,
          scythe_precision_error, 
          "Answer < 1/2 precision because x is " <<
            "too near a negative integer");

      double sinpiy = std::sin(M_PI * y);

      /* Negative integer arg - overflow */
      SCYTHE_CHECK_10(sinpiy == 0, scythe_range_error, "Overflow");

      return -M_PI / (y * sinpiy * value);
    }
  }

  /* The natural log of the absolute value of the gamma function */
  /*! \brief The natural log of the absolute value of the gamma 
   * function.
   *
   * Computes the natural log of the absolute value of the gamma 
   * function, evaluated at \a x.
   *
   * \param x The value to compute log(abs(gamma())) at.
   *
   * \see gammafn(double x)
   * \see pgamma(double x, double shape, double scale)
   * \see dgamma(double x, double shape, double scale)
   * \see rng::rgamma(double shape, double scale)
   *
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  lngammafn(double x)
  {
    SCYTHE_CHECK_10(x <= 0 && x == (int) x, scythe_range_error,
        "x is 0 or a negative integer");

    double y = std::fabs(x);

    if (y <= 10)
      return std::log(std::fabs(gammafn(x)));

    SCYTHE_CHECK_10(y > 2.5327372760800758e+305, scythe_range_error,
        "Overflow");

    if (x > 0) /* i.e. y = x > 10 */
      return M_LN_SQRT_2PI + (x - 0.5) * std::log(x) - x
        + lngammacor(x);
    
    /* else: x < -10; y = -x */
    double sinpiy = std::fabs(std::sin(M_PI * y));

    if (sinpiy == 0) /* Negative integer argument */
      throw scythe_exception("UNEXPECTED ERROR",
           __FILE__, __func__, __LINE__,
           "ERROR:  Should never happen!");

    double ans = M_LN_SQRT_PId2 + (x - 0.5) * std::log(y) - x - std::log(sinpiy)
      - lngammacor(y);

    SCYTHE_CHECK_10(std::fabs((x - (int)(x - 0.5)) * ans / x) 
        < 1.490116119384765696e-8, scythe_precision_error, 
        "Answer < 1/2 precision because x is " 
        << "too near a negative integer");
    
    return ans;
  }

  /* The beta function */
  /*! \brief The beta function.
   *
   * Computes beta function, evaluated at (\a a, \a b).
   *
   * \param a The first parameter.
   * \param b The second parameter.
   *
   * \see lnbetafn(double a, double b)
   * \see pbeta(double x, double a, double b)
   * \see dbeta(double x, double a, double b)