distributions.h

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

   *
   */
  inline double
  pchisq(double x, double df)
  {
    return pgamma(x, df/2.0, 2.0);
  }

  SCYTHE_DISTFUN_MATRIX(pchisq, double, df, double df)

  /* PDFs */
  /*! \brief The \f$\chi^2\f$ density function.
   *
   * Computes the value of the \f$\chi^2\f$ probability density
   * 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 pchisq(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)
   *
   */
  inline double
  dchisq(double x, double df)
  {
    return dgamma(x, df / 2.0, 2.0);
  }

  SCYTHE_DISTFUN_MATRIX(dchisq, double, df, double df)

  /**** The Exponential Distribution ****/

  /* CDFs */
  /*! \brief The exponential distribution function.
   *
   * Computes the value of the exponential cumulative distribution
   * function with given \a scale, 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 scale The positive scale of the function.
   *
   * \see dexp(double x, double scale)
   * \see rng::rexp(double scale)
   *
   * \throw scythe_invalid_arg (Level 1)
   */
  inline double
  pexp(double x, double scale)
  {
    SCYTHE_CHECK_10(scale <= 0, scythe_invalid_arg, "scale <= 0");
    
    if (x <= 0)
      return 0;
    
    return (1 - std::exp(-x*scale));
  }

  SCYTHE_DISTFUN_MATRIX(pexp, double, scale, double scale)

  /* PDFs */
  /*! \brief The exponential density function.
   *
   * Computes the value of the exponential probability density
   * function with given \a scale, 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 scale The positive scale of the function.
   *
   * \see pexp(double x, double scale)
   * \see rng::rexp(double scale)
   *
   * \throw scythe_invalid_arg (Level 1)
   */
  inline double
  dexp(double x, double scale)
  {
    SCYTHE_CHECK_10(scale <= 0, scythe_invalid_arg, "scale <= 0");
    
    if (x < 0)
      return 0;
      
    return std::exp(-x * scale) * scale;
  }

  SCYTHE_DISTFUN_MATRIX(dexp, double, scale, double scale)

  /**** The f Distribution ****/

  /* CDFs */
  /*! \brief The F distribution function.
   *
   * Computes the value of the F cumulative distribution function with
   * \a df1 and \a df2 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 df1 The non-negative degrees of freedom for the
   * \f$\chi^2\f$ variate in the nominator of the F statistic.
   * \param df2 The non-negative degrees of freedom for the
   * \f$\chi^2\f$ variate in the denominator of the F statistic.
   *
   *
   * \see df(double x, double df1, double df2)
   * \see rng::rf(double df1, double df2)
   * 
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   * \throw scythe_convergence_error (Level 1)
   */
  inline double
  pf(double x, double df1, double df2)
  {
    SCYTHE_CHECK_10(df1 <= 0 || df2 <= 0, scythe_invalid_arg, 
        "df1 or df2 <= 0");
  
    if (x <= 0)
      return 0;
    
    if (df2 > 4e5)
      return pchisq(x*df1,df1);
    if (df1 > 4e5)
      return 1-pchisq(df2/x,df2);
    
    return (1-pbeta(df2 / (df2 + df1 * x), df2 / 2.0, df1 / 2.0));
  }

  SCYTHE_DISTFUN_MATRIX(pf, double, SCYTHE_ARGSET(df1, df2), double df1, double df2)

  /* PDFs */

  /*! \brief The F density function.
   *
   * Computes the value of the F probability density function with
   * \a df1 and \a df2 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 df1 The non-negative degrees of freedom for the
   * \f$\chi^2\f$ variate in the nominator of the F statistic.
   * \param df2 The non-negative degrees of freedom for the
   * \f$\chi^2\f$ variate in the denominator of the F statistic.
   *
   * \see df(double x, double df1, double df2)
   * \see rng::rf(double df1, double df2)
   * 
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   * \throw scythe_convergence_error (Level 1)
   */
  inline double
  df(double x, double df1, double df2)
  {
    double dens;
    
    SCYTHE_CHECK_10(df1 <= 0 || df2 <= 0, scythe_invalid_arg, 
        "df1 or df2 <= 0");
    
    if (x <= 0)
      return 0;
      
    double f = 1 / (df2 + x * df1);
    double q = df2 * f;
    double p = x * df1 * f;
    
    if (df1 >= 2) {
      f = df1 * q / 2;
      dens = dbinom_raw((df1 - 2) / 2,(df1 + df2 - 2) / 2, p, q);
    } else {
      f = (df1 * df1 * q) /(2 * p * (df1 + df2));
      dens = dbinom_raw(df1 / 2,(df1 + df2)/ 2, p, q);
    }
    
    return f*dens;
  }

  SCYTHE_DISTFUN_MATRIX(df, double, SCYTHE_ARGSET(df1, df2), double df1, double df2)

  /**** The Gamma Distribution ****/

  /* CDFs */
  /*! \brief The gamma distribution function.
   *
   * Computes the value of the gamma cumulative distribution
   * function with given \a shape and \a scale, 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 shape The non-negative shape of the distribution.
   * \param scale The non-negative scale of the distribution.
   *
   * \see dgamma(double x, double shape, double scale)
   * \see rng::rgamma(double shape, double scale)
   * \see gammafn(double x)
   * \see lngammafn(double x)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   * \throw scythe_convergence_error (Level 1)
   */
  inline double
  pgamma (double x, double shape, double scale)
  {
    const double xbig = 1.0e+8, xlarge = 1.0e+37, 
      alphlimit = 1000.;/* normal approx. for shape > alphlimit */
      
    int lower_tail = 1;

    double pn1, pn2, pn3, pn4, pn5, pn6, arg, a, b, c, an, osum, sum;
    long n;
    int pearson;

    /* check that we have valid values for x and shape */

    SCYTHE_CHECK_10(shape <= 0. || scale <= 0., scythe_invalid_arg,
        "shape or scale <= 0");

    x /= scale;
    
    if (x <= 0.)
      return 0.0;

    /* use a normal approximation if shape > alphlimit */

    if (shape > alphlimit) {
      pn1 = std::sqrt(shape) * 3. * (std::pow(x/shape, 1./3.) + 1.
            / (9. * shape) - 1.);
      return pnorm(pn1, 0., 1.);
    }

    /* if x is extremely large __compared to shape__ then return 1 */

    if (x > xbig * shape)
      return 1.0;

    if (x <= 1. || x < shape) {
      pearson = 1;/* use pearson's series expansion. */
      arg = shape * std::log(x) - x - lngammafn(shape + 1.);
      c = 1.;
      sum = 1.;
      a = shape;
      do {
        a += 1.;
        c *= x / a;
        sum += c;
      } while (c > DBL_EPSILON);
      arg += std::log(sum);
    }
    else { /* x >= max( 1, shape) */
      pearson = 0;/* use a continued fraction expansion */
      arg = shape * std::log(x) - x - lngammafn(shape);
      a = 1. - shape;
      b = a + x + 1.;
      pn1 = 1.;
      pn2 = x;
      pn3 = x + 1.;
      pn4 = x * b;
      sum = pn3 / pn4;
      for (n = 1; ; n++) {
        a += 1.;/* =   n+1 -shape */
        b += 2.;/* = 2(n+1)-shape+x */
        an = a * n;
        pn5 = b * pn3 - an * pn1;
        pn6 = b * pn4 - an * pn2;
        if (std::fabs(pn6) > 0.) {
          osum = sum;
          sum = pn5 / pn6;
          if (std::fabs(osum - sum) <= DBL_EPSILON * std::min(1., sum))
            break;
        }
        pn1 = pn3;
        pn2 = pn4;
        pn3 = pn5;
        pn4 = pn6;
        if (std::fabs(pn5) >= xlarge) {
          /* re-scale terms in continued fraction if they are large */
          pn1 /= xlarge;
          pn2 /= xlarge;
          pn3 /= xlarge;
          pn4 /= xlarge;
        }
      }
      arg += std::log(sum);
    }

    lower_tail = (lower_tail == pearson);

    sum = std::exp(arg);

    return (lower_tail) ? sum : 1 - sum;
  }

  SCYTHE_DISTFUN_MATRIX(pgamma, double, SCYTHE_ARGSET(shape, scale), double shape, double scale)

  /* PDFs */
  /*! \brief The gamma density function.
   *
   * Computes the value of the gamma probability density
   * function with given \a shape and \a scale, 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 shape The non-negative shape of the distribution.
   * \param scale The non-negative scale of the distribution.
   *
   * \see pgamma(double x, double shape, double scale)
   * \see rng::rgamma(double shape, double scale)
   * \see gammafn(double x)
   * \see lngammafn(double x)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   * \throw scythe_convergence_error (Level 1)
   */
  inline double
  dgamma(double x, double shape, double scale)
  {
    SCYTHE_CHECK_10(shape <= 0 || scale <= 0,scythe_invalid_arg,
        "shape or scale <= 0");

    if (x < 0)
      return 0.0;
    
    if (x == 0) {
      SCYTHE_CHECK_10(shape < 1,scythe_invalid_arg,