distributions.h

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

          "x == 0 and shape < 1");
      
      if (shape > 1)
        return 0.0;
      
      return 1 / scale;
    }
    
    if (shape < 1) { 
      double pr = dpois_raw(shape, x/scale);
      return pr * shape / x;
    }
    
    /* else  shape >= 1 */
    double pr = dpois_raw(shape - 1, x / scale);
    return pr / scale;
  }

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

  /**** The Logistic Distribution ****/

  /* CDFs */
  /*! \brief The logistic distribution function.
   *
   * Computes the value of the logistic cumulative distribution
   * function with given \a location 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 location The location of the distribution.
   * \param scale The positive scale of the distribution.
   *
   * \see dlogis(double x, double location, double scale)
   * \see rng::rlogis(double location, double scale)
   * 
   * \throw scythe_invalid_arg (Level 1)
   */
  inline double
  plogis (double x, double location, double scale)
  {
    SCYTHE_CHECK_10(scale <= 0.0, scythe_invalid_arg, "scale <= 0");
    
    double X = (x-location) / scale;
      
    X = std::exp(-X);
      
    return 1 / (1+X);
  }

  SCYTHE_DISTFUN_MATRIX(plogis, double, SCYTHE_ARGSET(location, scale), double location, double scale)

  /* PDFs */
  /*! \brief The logistic density function.
   *
   * Computes the value of the logistic probability density
   * function with given \a location 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 location The location of the distribution.
   * \param scale The positive scale of the distribution.
   *
   * \see plogis(double x, double location, double scale)
   * \see rng::rlogis(double location, double scale)
   * 
   * \throw scythe_invalid_arg (Level 1)
   */
  inline double
  dlogis(double x, double location, double scale)
  {
    SCYTHE_CHECK_10(scale <= 0.0, scythe_invalid_arg, "scale <= 0");
    
    double X = (x - location) / scale;
    double e = std::exp(-X);
    double f = 1.0 + e;
      
    return e / (scale * f * f);
  }

  SCYTHE_DISTFUN_MATRIX(dlogis, double, SCYTHE_ARGSET(location, scale), double location, double scale)

  /**** The Log Normal Distribution ****/

  /* CDFs */
  /*! \brief The log-normal distribution function.
   *
   * Computes the value of the log-normal cumulative distribution
   * function with mean \a logmean and standard
   * deviation \a logsd, 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 logmean The mean of the distribution.
   * \param logsd The positive standard deviation of the distribution.
   *
   * \see dlnorm(double x, double logmean, double logsd)
   * \see rng::rlnorm(double logmean, double logsd)
   * \see pnorm(double x, double logmean, double logsd)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_convergence_error (Level 1)
   */
  inline double
  plnorm (double x, double logmean, double logsd)
  {
    SCYTHE_CHECK_10(logsd <= 0, scythe_invalid_arg, "logsd <= 0");
    
    if (x > 0)
      return pnorm(std::log(x), logmean, logsd);
    
    return 0;
  }

  SCYTHE_DISTFUN_MATRIX(plnorm, double, SCYTHE_ARGSET(logmean, logsd), double logmean, double logsd)

  /* PDFs */
  /*! \brief The log-normal density function.
   *
   * Computes the value of the log-normal probability density
   * function with mean \a logmean and standard
   * deviation \a logsd, 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 logmean The mean of the distribution.
   * \param logsd The positive standard deviation of the distribution.
   *
   * \see plnorm(double x, double logmean, double logsd)
   * \see rng::rlnorm(double logmean, double logsd)
   * \see dnorm(double x, double logmean, double logsd)
   *
   * \throw scythe_invalid_arg (Level 1)
   */
  inline double
  dlnorm(double x, double logmean, double logsd)
  {
    SCYTHE_CHECK_10(logsd <= 0, scythe_invalid_arg, "logsd <= 0");
    
    if (x == 0)
      return 0;
    
    double y = (std::log(x) - logmean) / logsd;
    
    return (1 / (std::sqrt(2 * M_PI))) * std::exp(-0.5 * y * y) / (x * logsd);
  }

  SCYTHE_DISTFUN_MATRIX(dlnorm, double, SCYTHE_ARGSET(logmean, logsd), double logmean, double logsd)

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

  /* CDFs */
  /*! \brief The negative binomial distribution function.
   *
   * Computes the value of the negative binomial cumulative distribution
   * function with \a n target number of successful 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 non-negative, integer, quantile.
   * \param n The positive target number of successful trials
   * (dispersion parameter).
   * \param p The probability of success on each trial.
   *
   * \see dnbinom(unsigned int x, double n, double p)
   * \see rng::rnbinom(double n, double p)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  pnbinom(unsigned int x, double n, double p)
  {
    SCYTHE_CHECK_10(n == 0 || p <= 0 || p >= 1, scythe_invalid_arg,
        "n == 0 or p not in (0,1)");
    
    return pbeta(p, n, x + 1);
  }

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

  /* PDFs */
  /*! \brief The negative binomial density function.
   *
   * Computes the value of the negative binomial probability density
   * function with \a n target number of successful 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 non-negative, integer, quantile.
   * \param n The positive target number of successful trials
   * (dispersion parameter).
   * \param p The probability of success on each trial.
   *
   * \see dnbinom(unsigned int x, double n, double p)
   * \see rng::rnbinom(double n, double p)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_range_error (Level 1)
   * \throw scythe_precision_error (Level 1)
   */
  inline double
  dnbinom(unsigned int x, double n, double p)
  {
    SCYTHE_CHECK_10(n == 0 || p <= 0 || p >= 1, scythe_invalid_arg,
        "n == 0 or p not in (0,1)");
    
    double prob = dbinom_raw(n, x + n, p, 1 - p);
    double P = (double) n / (n + x);
    
    return P * prob;
  }

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

  /**** The Normal Distribution ****/
  
  /* CDFs */
  /*! \brief The normal distribution function.
   *
   * Computes the value of the normal cumulative distribution
   * function with given \a mean and standard deviation \a sd, 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 mean The mean of the distribution.
   * \param sd The positive standard deviation of the distribution.
   *
   * \see dnorm(double x, double mean, double sd)
   * \see rng::rnorm(double mean, double sd)
   *
   * \throw scythe_invalid_arg (Level 1)
   * \throw scythe_convergence_error (Level 1)
   */
  inline double
  pnorm (double x, double mean, double sd)
  
  {
    SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg,
        "negative standard deviation");

    return pnorm1((x - mean) / sd, true, false);
  }

  SCYTHE_DISTFUN_MATRIX(pnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd)
  

  /* PDFs */
  /*! \brief The normal density function.
   *
   * Computes the value of the normal probability density
   * function with given \a mean and standard deviation \a sd, 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 mean The mean of the distribution.
   * \param sd The positive standard deviation of the distribution.
   *
   * \see pnorm(double x, double mean, double sd)
   * \see rng::rnorm(double mean, double sd)
   *
   * \throw scythe_invalid_arg (Level 1)
   */
  inline double
  dnorm(double x, double mean, double sd)
  {
    SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg,
        "negative standard deviation");
    
    double X = (x - mean) / sd;
    
    return (M_1_SQRT_2PI * std::exp(-0.5 * X * X) / sd);
  }

  SCYTHE_DISTFUN_MATRIX(dnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd)
 

  /* Return the natural log of the normal PDF */
  /*! \brief The natural log of normal density function.
   *
   * Computes the value of the natural log of the normal probability
   * density function with given \a mean and standard deviation \a sd,
   * 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 mean The mean of the distribution.
   * \param sd The positive standard deviation of the distribution.
   *
   * \see dnorm(double x, double mean, double sd)
   * \see pnorm(double x, double mean, double sd)
   * \see rng::rnorm(double mean, double sd)
   *
   * \throw scythe_invalid_arg (Level 1)
   */
  inline double
  lndnorm (double x, double mean, double sd)
  {
    SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg,
        "negative standard deviation");
    
    double X = (x - mean) / sd;
    
    return -(M_LN_SQRT_2PI  +  0.5 * X * X + std::log(sd));
  }

  SCYTHE_DISTFUN_MATRIX(lndnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd)

  /* Quantile functions */
  /*! \brief The standard normal quantile function.
   *
   * Computes the value of the standard normal quantile function
   * at the desired probability \a in_p.
   *
   * 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.