diff options
Diffstat (limited to 'src/third_party/boost-1.56.0/boost/math/distributions')
23 files changed, 9029 insertions, 0 deletions
diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/beta.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/beta.hpp new file mode 100644 index 00000000000..5ecf902d990 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/beta.hpp @@ -0,0 +1,541 @@ +// boost\math\distributions\beta.hpp + +// Copyright John Maddock 2006. +// Copyright Paul A. Bristow 2006. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// http://en.wikipedia.org/wiki/Beta_distribution +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm +// http://mathworld.wolfram.com/BetaDistribution.html + +// The Beta Distribution is a continuous probability distribution. +// The beta distribution is used to model events which are constrained to take place +// within an interval defined by maxima and minima, +// so is used extensively in PERT and other project management systems +// to describe the time to completion. +// The cdf of the beta distribution is used as a convenient way +// of obtaining the sum over a set of binomial outcomes. +// The beta distribution is also used in Bayesian statistics. + +#ifndef BOOST_MATH_DIST_BETA_HPP +#define BOOST_MATH_DIST_BETA_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/beta.hpp> // for beta. +#include <boost/math/distributions/complement.hpp> // complements. +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks +#include <boost/math/special_functions/fpclassify.hpp> // isnan. +#include <boost/math/tools/roots.hpp> // for root finding. + +#if defined (BOOST_MSVC) +# pragma warning(push) +# pragma warning(disable: 4702) // unreachable code +// in domain_error_imp in error_handling +#endif + +#include <utility> + +namespace boost +{ + namespace math + { + namespace beta_detail + { + // Common error checking routines for beta distribution functions: + template <class RealType, class Policy> + inline bool check_alpha(const char* function, const RealType& alpha, RealType* result, const Policy& pol) + { + if(!(boost::math::isfinite)(alpha) || (alpha <= 0)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Alpha argument is %1%, but must be > 0 !", alpha, pol); + return false; + } + return true; + } // bool check_alpha + + template <class RealType, class Policy> + inline bool check_beta(const char* function, const RealType& beta, RealType* result, const Policy& pol) + { + if(!(boost::math::isfinite)(beta) || (beta <= 0)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Beta argument is %1%, but must be > 0 !", beta, pol); + return false; + } + return true; + } // bool check_beta + + template <class RealType, class Policy> + inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) + { + if((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); + return false; + } + return true; + } // bool check_prob + + template <class RealType, class Policy> + inline bool check_x(const char* function, const RealType& x, RealType* result, const Policy& pol) + { + if(!(boost::math::isfinite)(x) || (x < 0) || (x > 1)) + { + *result = policies::raise_domain_error<RealType>( + function, + "x argument is %1%, but must be >= 0 and <= 1 !", x, pol); + return false; + } + return true; + } // bool check_x + + template <class RealType, class Policy> + inline bool check_dist(const char* function, const RealType& alpha, const RealType& beta, RealType* result, const Policy& pol) + { // Check both alpha and beta. + return check_alpha(function, alpha, result, pol) + && check_beta(function, beta, result, pol); + } // bool check_dist + + template <class RealType, class Policy> + inline bool check_dist_and_x(const char* function, const RealType& alpha, const RealType& beta, RealType x, RealType* result, const Policy& pol) + { + return check_dist(function, alpha, beta, result, pol) + && beta_detail::check_x(function, x, result, pol); + } // bool check_dist_and_x + + template <class RealType, class Policy> + inline bool check_dist_and_prob(const char* function, const RealType& alpha, const RealType& beta, RealType p, RealType* result, const Policy& pol) + { + return check_dist(function, alpha, beta, result, pol) + && check_prob(function, p, result, pol); + } // bool check_dist_and_prob + + template <class RealType, class Policy> + inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol) + { + if(!(boost::math::isfinite)(mean) || (mean <= 0)) + { + *result = policies::raise_domain_error<RealType>( + function, + "mean argument is %1%, but must be > 0 !", mean, pol); + return false; + } + return true; + } // bool check_mean + template <class RealType, class Policy> + inline bool check_variance(const char* function, const RealType& variance, RealType* result, const Policy& pol) + { + if(!(boost::math::isfinite)(variance) || (variance <= 0)) + { + *result = policies::raise_domain_error<RealType>( + function, + "variance argument is %1%, but must be > 0 !", variance, pol); + return false; + } + return true; + } // bool check_variance + } // namespace beta_detail + + // typedef beta_distribution<double> beta; + // is deliberately NOT included to avoid a name clash with the beta function. + // Use beta_distribution<> mybeta(...) to construct type double. + + template <class RealType = double, class Policy = policies::policy<> > + class beta_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + beta_distribution(RealType l_alpha = 1, RealType l_beta = 1) : m_alpha(l_alpha), m_beta(l_beta) + { + RealType result; + beta_detail::check_dist( + "boost::math::beta_distribution<%1%>::beta_distribution", + m_alpha, + m_beta, + &result, Policy()); + } // beta_distribution constructor. + // Accessor functions: + RealType alpha() const + { + return m_alpha; + } + RealType beta() const + { // . + return m_beta; + } + + // Estimation of the alpha & beta parameters. + // http://en.wikipedia.org/wiki/Beta_distribution + // gives formulae in section on parameter estimation. + // Also NIST EDA page 3 & 4 give the same. + // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm + // http://www.epi.ucdavis.edu/diagnostictests/betabuster.html + + static RealType find_alpha( + RealType mean, // Expected value of mean. + RealType variance) // Expected value of variance. + { + static const char* function = "boost::math::beta_distribution<%1%>::find_alpha"; + RealType result = 0; // of error checks. + if(false == + ( + beta_detail::check_mean(function, mean, &result, Policy()) + && beta_detail::check_variance(function, variance, &result, Policy()) + ) + ) + { + return result; + } + return mean * (( (mean * (1 - mean)) / variance)- 1); + } // RealType find_alpha + + static RealType find_beta( + RealType mean, // Expected value of mean. + RealType variance) // Expected value of variance. + { + static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; + RealType result = 0; // of error checks. + if(false == + ( + beta_detail::check_mean(function, mean, &result, Policy()) + && + beta_detail::check_variance(function, variance, &result, Policy()) + ) + ) + { + return result; + } + return (1 - mean) * (((mean * (1 - mean)) /variance)-1); + } // RealType find_beta + + // Estimate alpha & beta from either alpha or beta, and x and probability. + // Uses for these parameter estimators are unclear. + + static RealType find_alpha( + RealType beta, // from beta. + RealType x, // x. + RealType probability) // cdf + { + static const char* function = "boost::math::beta_distribution<%1%>::find_alpha"; + RealType result = 0; // of error checks. + if(false == + ( + beta_detail::check_prob(function, probability, &result, Policy()) + && + beta_detail::check_beta(function, beta, &result, Policy()) + && + beta_detail::check_x(function, x, &result, Policy()) + ) + ) + { + return result; + } + return ibeta_inva(beta, x, probability, Policy()); + } // RealType find_alpha(beta, a, probability) + + static RealType find_beta( + // ibeta_invb(T b, T x, T p); (alpha, x, cdf,) + RealType alpha, // alpha. + RealType x, // probability x. + RealType probability) // probability cdf. + { + static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; + RealType result = 0; // of error checks. + if(false == + ( + beta_detail::check_prob(function, probability, &result, Policy()) + && + beta_detail::check_alpha(function, alpha, &result, Policy()) + && + beta_detail::check_x(function, x, &result, Policy()) + ) + ) + { + return result; + } + return ibeta_invb(alpha, x, probability, Policy()); + } // RealType find_beta(alpha, x, probability) + + private: + RealType m_alpha; // Two parameters of the beta distribution. + RealType m_beta; + }; // template <class RealType, class Policy> class beta_distribution + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const beta_distribution<RealType, Policy>& /* dist */) + { // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const beta_distribution<RealType, Policy>& /* dist */) + { // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); + } + + template <class RealType, class Policy> + inline RealType mean(const beta_distribution<RealType, Policy>& dist) + { // Mean of beta distribution = np. + return dist.alpha() / (dist.alpha() + dist.beta()); + } // mean + + template <class RealType, class Policy> + inline RealType variance(const beta_distribution<RealType, Policy>& dist) + { // Variance of beta distribution = np(1-p). + RealType a = dist.alpha(); + RealType b = dist.beta(); + return (a * b) / ((a + b ) * (a + b) * (a + b + 1)); + } // variance + + template <class RealType, class Policy> + inline RealType mode(const beta_distribution<RealType, Policy>& dist) + { + static const char* function = "boost::math::mode(beta_distribution<%1%> const&)"; + + RealType result; + if ((dist.alpha() <= 1)) + { + result = policies::raise_domain_error<RealType>( + function, + "mode undefined for alpha = %1%, must be > 1!", dist.alpha(), Policy()); + return result; + } + + if ((dist.beta() <= 1)) + { + result = policies::raise_domain_error<RealType>( + function, + "mode undefined for beta = %1%, must be > 1!", dist.beta(), Policy()); + return result; + } + RealType a = dist.alpha(); + RealType b = dist.beta(); + return (a-1) / (a + b - 2); + } // mode + + //template <class RealType, class Policy> + //inline RealType median(const beta_distribution<RealType, Policy>& dist) + //{ // Median of beta distribution is not defined. + // return tools::domain_error<RealType>(function, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); + //} // median + + //But WILL be provided by the derived accessor as quantile(0.5). + + template <class RealType, class Policy> + inline RealType skewness(const beta_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // ADL of std functions. + RealType a = dist.alpha(); + RealType b = dist.beta(); + return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b)); + } // skewness + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const beta_distribution<RealType, Policy>& dist) + { + RealType a = dist.alpha(); + RealType b = dist.beta(); + RealType a_2 = a * a; + RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2)); + RealType d = a * b * (a + b + 2) * (a + b + 3); + return n / d; + } // kurtosis_excess + + template <class RealType, class Policy> + inline RealType kurtosis(const beta_distribution<RealType, Policy>& dist) + { + return 3 + kurtosis_excess(dist); + } // kurtosis + + template <class RealType, class Policy> + inline RealType pdf(const beta_distribution<RealType, Policy>& dist, const RealType& x) + { // Probability Density/Mass Function. + BOOST_FPU_EXCEPTION_GUARD + + static const char* function = "boost::math::pdf(beta_distribution<%1%> const&, %1%)"; + + BOOST_MATH_STD_USING // for ADL of std functions + + RealType a = dist.alpha(); + RealType b = dist.beta(); + + // Argument checks: + RealType result = 0; + if(false == beta_detail::check_dist_and_x( + function, + a, b, x, + &result, Policy())) + { + return result; + } + using boost::math::beta; + return ibeta_derivative(a, b, x, Policy()); + } // pdf + + template <class RealType, class Policy> + inline RealType cdf(const beta_distribution<RealType, Policy>& dist, const RealType& x) + { // Cumulative Distribution Function beta. + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; + + RealType a = dist.alpha(); + RealType b = dist.beta(); + + // Argument checks: + RealType result = 0; + if(false == beta_detail::check_dist_and_x( + function, + a, b, x, + &result, Policy())) + { + return result; + } + // Special cases: + if (x == 0) + { + return 0; + } + else if (x == 1) + { + return 1; + } + return ibeta(a, b, x, Policy()); + } // beta cdf + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c) + { // Complemented Cumulative Distribution Function beta. + + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; + + RealType const& x = c.param; + beta_distribution<RealType, Policy> const& dist = c.dist; + RealType a = dist.alpha(); + RealType b = dist.beta(); + + // Argument checks: + RealType result = 0; + if(false == beta_detail::check_dist_and_x( + function, + a, b, x, + &result, Policy())) + { + return result; + } + if (x == 0) + { + return 1; + } + else if (x == 1) + { + return 0; + } + // Calculate cdf beta using the incomplete beta function. + // Use of ibeta here prevents cancellation errors in calculating + // 1 - x if x is very small, perhaps smaller than machine epsilon. + return ibetac(a, b, x, Policy()); + } // beta cdf + + template <class RealType, class Policy> + inline RealType quantile(const beta_distribution<RealType, Policy>& dist, const RealType& p) + { // Quantile or Percent Point beta function or + // Inverse Cumulative probability distribution function CDF. + // Return x (0 <= x <= 1), + // for a given probability p (0 <= p <= 1). + // These functions take a probability as an argument + // and return a value such that the probability that a random variable x + // will be less than or equal to that value + // is whatever probability you supplied as an argument. + + static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; + + RealType result = 0; // of argument checks: + RealType a = dist.alpha(); + RealType b = dist.beta(); + if(false == beta_detail::check_dist_and_prob( + function, + a, b, p, + &result, Policy())) + { + return result; + } + // Special cases: + if (p == 0) + { + return 0; + } + if (p == 1) + { + return 1; + } + return ibeta_inv(a, b, p, static_cast<RealType*>(0), Policy()); + } // quantile + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c) + { // Complement Quantile or Percent Point beta function . + // Return the number of expected x for a given + // complement of the probability q. + + static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; + + // + // Error checks: + RealType q = c.param; + const beta_distribution<RealType, Policy>& dist = c.dist; + RealType result = 0; + RealType a = dist.alpha(); + RealType b = dist.beta(); + if(false == beta_detail::check_dist_and_prob( + function, + a, + b, + q, + &result, Policy())) + { + return result; + } + // Special cases: + if(q == 1) + { + return 0; + } + if(q == 0) + { + return 1; + } + + return ibetac_inv(a, b, q, static_cast<RealType*>(0), Policy()); + } // Quantile Complement + + } // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#if defined (BOOST_MSVC) +# pragma warning(pop) +#endif + +#endif // BOOST_MATH_DIST_BETA_HPP + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/binomial.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/binomial.hpp new file mode 100644 index 00000000000..a48c89c5b93 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/binomial.hpp @@ -0,0 +1,728 @@ +// boost\math\distributions\binomial.hpp + +// Copyright John Maddock 2006. +// Copyright Paul A. Bristow 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// http://en.wikipedia.org/wiki/binomial_distribution + +// Binomial distribution is the discrete probability distribution of +// the number (k) of successes, in a sequence of +// n independent (yes or no, success or failure) Bernoulli trials. + +// It expresses the probability of a number of events occurring in a fixed time +// if these events occur with a known average rate (probability of success), +// and are independent of the time since the last event. + +// The number of cars that pass through a certain point on a road during a given period of time. +// The number of spelling mistakes a secretary makes while typing a single page. +// The number of phone calls at a call center per minute. +// The number of times a web server is accessed per minute. +// The number of light bulbs that burn out in a certain amount of time. +// The number of roadkill found per unit length of road + +// http://en.wikipedia.org/wiki/binomial_distribution + +// Given a sample of N measured values k[i], +// we wish to estimate the value of the parameter x (mean) +// of the binomial population from which the sample was drawn. +// To calculate the maximum likelihood value = 1/N sum i = 1 to N of k[i] + +// Also may want a function for EXACTLY k. + +// And probability that there are EXACTLY k occurrences is +// exp(-x) * pow(x, k) / factorial(k) +// where x is expected occurrences (mean) during the given interval. +// For example, if events occur, on average, every 4 min, +// and we are interested in number of events occurring in 10 min, +// then x = 10/4 = 2.5 + +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366i.htm + +// The binomial distribution is used when there are +// exactly two mutually exclusive outcomes of a trial. +// These outcomes are appropriately labeled "success" and "failure". +// The binomial distribution is used to obtain +// the probability of observing x successes in N trials, +// with the probability of success on a single trial denoted by p. +// The binomial distribution assumes that p is fixed for all trials. + +// P(x, p, n) = n!/(x! * (n-x)!) * p^x * (1-p)^(n-x) + +// http://mathworld.wolfram.com/BinomialCoefficient.html + +// The binomial coefficient (n; k) is the number of ways of picking +// k unordered outcomes from n possibilities, +// also known as a combination or combinatorial number. +// The symbols _nC_k and (n; k) are used to denote a binomial coefficient, +// and are sometimes read as "n choose k." +// (n; k) therefore gives the number of k-subsets possible out of a set of n distinct items. + +// For example: +// The 2-subsets of {1,2,3,4} are the six pairs {1,2}, {1,3}, {1,4}, {2,3}, {2,4}, and {3,4}, so (4; 2)==6. + +// http://functions.wolfram.com/GammaBetaErf/Binomial/ for evaluation. + +// But note that the binomial distribution +// (like others including the poisson, negative binomial & Bernoulli) +// is strictly defined as a discrete function: only integral values of k are envisaged. +// However because of the method of calculation using a continuous gamma function, +// it is convenient to treat it as if a continous function, +// and permit non-integral values of k. +// To enforce the strict mathematical model, users should use floor or ceil functions +// on k outside this function to ensure that k is integral. + +#ifndef BOOST_MATH_SPECIAL_BINOMIAL_HPP +#define BOOST_MATH_SPECIAL_BINOMIAL_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/beta.hpp> // for incomplete beta. +#include <boost/math/distributions/complement.hpp> // complements +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks +#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> // error checks +#include <boost/math/special_functions/fpclassify.hpp> // isnan. +#include <boost/math/tools/roots.hpp> // for root finding. + +#include <utility> + +namespace boost +{ + namespace math + { + + template <class RealType, class Policy> + class binomial_distribution; + + namespace binomial_detail{ + // common error checking routines for binomial distribution functions: + template <class RealType, class Policy> + inline bool check_N(const char* function, const RealType& N, RealType* result, const Policy& pol) + { + if((N < 0) || !(boost::math::isfinite)(N)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Number of Trials argument is %1%, but must be >= 0 !", N, pol); + return false; + } + return true; + } + template <class RealType, class Policy> + inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol) + { + if((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol); + return false; + } + return true; + } + template <class RealType, class Policy> + inline bool check_dist(const char* function, const RealType& N, const RealType& p, RealType* result, const Policy& pol) + { + return check_success_fraction( + function, p, result, pol) + && check_N( + function, N, result, pol); + } + template <class RealType, class Policy> + inline bool check_dist_and_k(const char* function, const RealType& N, const RealType& p, RealType k, RealType* result, const Policy& pol) + { + if(check_dist(function, N, p, result, pol) == false) + return false; + if((k < 0) || !(boost::math::isfinite)(k)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Number of Successes argument is %1%, but must be >= 0 !", k, pol); + return false; + } + if(k > N) + { + *result = policies::raise_domain_error<RealType>( + function, + "Number of Successes argument is %1%, but must be <= Number of Trials !", k, pol); + return false; + } + return true; + } + template <class RealType, class Policy> + inline bool check_dist_and_prob(const char* function, const RealType& N, RealType p, RealType prob, RealType* result, const Policy& pol) + { + if(check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol) == false) + return false; + return true; + } + + template <class T, class Policy> + T inverse_binomial_cornish_fisher(T n, T sf, T p, T q, const Policy& pol) + { + BOOST_MATH_STD_USING + // mean: + T m = n * sf; + // standard deviation: + T sigma = sqrt(n * sf * (1 - sf)); + // skewness + T sk = (1 - 2 * sf) / sigma; + // kurtosis: + // T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf)); + // Get the inverse of a std normal distribution: + T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>(); + // Set the sign: + if(p < 0.5) + x = -x; + T x2 = x * x; + // w is correction term due to skewness + T w = x + sk * (x2 - 1) / 6; + /* + // Add on correction due to kurtosis. + // Disabled for now, seems to make things worse? + // + if(n >= 10) + w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36; + */ + w = m + sigma * w; + if(w < tools::min_value<T>()) + return sqrt(tools::min_value<T>()); + if(w > n) + return n; + return w; + } + + template <class RealType, class Policy> + RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q, bool comp) + { // Quantile or Percent Point Binomial function. + // Return the number of expected successes k, + // for a given probability p. + // + // Error checks: + BOOST_MATH_STD_USING // ADL of std names + RealType result = 0; + RealType trials = dist.trials(); + RealType success_fraction = dist.success_fraction(); + if(false == binomial_detail::check_dist_and_prob( + "boost::math::quantile(binomial_distribution<%1%> const&, %1%)", + trials, + success_fraction, + p, + &result, Policy())) + { + return result; + } + + // Special cases: + // + if(p == 0) + { // There may actually be no answer to this question, + // since the probability of zero successes may be non-zero, + // but zero is the best we can do: + return 0; + } + if(p == 1) + { // Probability of n or fewer successes is always one, + // so n is the most sensible answer here: + return trials; + } + if (p <= pow(1 - success_fraction, trials)) + { // p <= pdf(dist, 0) == cdf(dist, 0) + return 0; // So the only reasonable result is zero. + } // And root finder would fail otherwise. + if(success_fraction == 1) + { // our formulae break down in this case: + return p > 0.5f ? trials : 0; + } + + // Solve for quantile numerically: + // + RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, success_fraction, p, q, Policy()); + RealType factor = 8; + if(trials > 100) + factor = 1.01f; // guess is pretty accurate + else if((trials > 10) && (trials - 1 > guess) && (guess > 3)) + factor = 1.15f; // less accurate but OK. + else if(trials < 10) + { + // pretty inaccurate guess in this area: + if(guess > trials / 64) + { + guess = trials / 4; + factor = 2; + } + else + guess = trials / 1024; + } + else + factor = 2; // trials largish, but in far tails. + + typedef typename Policy::discrete_quantile_type discrete_quantile_type; + boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); + return detail::inverse_discrete_quantile( + dist, + comp ? q : p, + comp, + guess, + factor, + RealType(1), + discrete_quantile_type(), + max_iter); + } // quantile + + } + + template <class RealType = double, class Policy = policies::policy<> > + class binomial_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + binomial_distribution(RealType n = 1, RealType p = 0.5) : m_n(n), m_p(p) + { // Default n = 1 is the Bernoulli distribution + // with equal probability of 'heads' or 'tails. + RealType r; + binomial_detail::check_dist( + "boost::math::binomial_distribution<%1%>::binomial_distribution", + m_n, + m_p, + &r, Policy()); + } // binomial_distribution constructor. + + RealType success_fraction() const + { // Probability. + return m_p; + } + RealType trials() const + { // Total number of trials. + return m_n; + } + + enum interval_type{ + clopper_pearson_exact_interval, + jeffreys_prior_interval + }; + + // + // Estimation of the success fraction parameter. + // The best estimate is actually simply successes/trials, + // these functions are used + // to obtain confidence intervals for the success fraction. + // + static RealType find_lower_bound_on_p( + RealType trials, + RealType successes, + RealType probability, + interval_type t = clopper_pearson_exact_interval) + { + static const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p"; + // Error checks: + RealType result = 0; + if(false == binomial_detail::check_dist_and_k( + function, trials, RealType(0), successes, &result, Policy()) + && + binomial_detail::check_dist_and_prob( + function, trials, RealType(0), probability, &result, Policy())) + { return result; } + + if(successes == 0) + return 0; + + // NOTE!!! The Clopper Pearson formula uses "successes" not + // "successes+1" as usual to get the lower bound, + // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm + return (t == clopper_pearson_exact_interval) ? ibeta_inv(successes, trials - successes + 1, probability, static_cast<RealType*>(0), Policy()) + : ibeta_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(0), Policy()); + } + static RealType find_upper_bound_on_p( + RealType trials, + RealType successes, + RealType probability, + interval_type t = clopper_pearson_exact_interval) + { + static const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p"; + // Error checks: + RealType result = 0; + if(false == binomial_detail::check_dist_and_k( + function, trials, RealType(0), successes, &result, Policy()) + && + binomial_detail::check_dist_and_prob( + function, trials, RealType(0), probability, &result, Policy())) + { return result; } + + if(trials == successes) + return 1; + + return (t == clopper_pearson_exact_interval) ? ibetac_inv(successes + 1, trials - successes, probability, static_cast<RealType*>(0), Policy()) + : ibetac_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(0), Policy()); + } + // Estimate number of trials parameter: + // + // "How many trials do I need to be P% sure of seeing k events?" + // or + // "How many trials can I have to be P% sure of seeing fewer than k events?" + // + static RealType find_minimum_number_of_trials( + RealType k, // number of events + RealType p, // success fraction + RealType alpha) // risk level + { + static const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials"; + // Error checks: + RealType result = 0; + if(false == binomial_detail::check_dist_and_k( + function, k, p, k, &result, Policy()) + && + binomial_detail::check_dist_and_prob( + function, k, p, alpha, &result, Policy())) + { return result; } + + result = ibetac_invb(k + 1, p, alpha, Policy()); // returns n - k + return result + k; + } + + static RealType find_maximum_number_of_trials( + RealType k, // number of events + RealType p, // success fraction + RealType alpha) // risk level + { + static const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials"; + // Error checks: + RealType result = 0; + if(false == binomial_detail::check_dist_and_k( + function, k, p, k, &result, Policy()) + && + binomial_detail::check_dist_and_prob( + function, k, p, alpha, &result, Policy())) + { return result; } + + result = ibeta_invb(k + 1, p, alpha, Policy()); // returns n - k + return result + k; + } + + private: + RealType m_n; // Not sure if this shouldn't be an int? + RealType m_p; // success_fraction + }; // template <class RealType, class Policy> class binomial_distribution + + typedef binomial_distribution<> binomial; + // typedef binomial_distribution<double> binomial; + // IS now included since no longer a name clash with function binomial. + //typedef binomial_distribution<double> binomial; // Reserved name of type double. + + template <class RealType, class Policy> + const std::pair<RealType, RealType> range(const binomial_distribution<RealType, Policy>& dist) + { // Range of permissible values for random variable k. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials()); + } + + template <class RealType, class Policy> + const std::pair<RealType, RealType> support(const binomial_distribution<RealType, Policy>& dist) + { // Range of supported values for random variable k. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + return std::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials()); + } + + template <class RealType, class Policy> + inline RealType mean(const binomial_distribution<RealType, Policy>& dist) + { // Mean of Binomial distribution = np. + return dist.trials() * dist.success_fraction(); + } // mean + + template <class RealType, class Policy> + inline RealType variance(const binomial_distribution<RealType, Policy>& dist) + { // Variance of Binomial distribution = np(1-p). + return dist.trials() * dist.success_fraction() * (1 - dist.success_fraction()); + } // variance + + template <class RealType, class Policy> + RealType pdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k) + { // Probability Density/Mass Function. + BOOST_FPU_EXCEPTION_GUARD + + BOOST_MATH_STD_USING // for ADL of std functions + + RealType n = dist.trials(); + + // Error check: + RealType result = 0; // initialization silences some compiler warnings + if(false == binomial_detail::check_dist_and_k( + "boost::math::pdf(binomial_distribution<%1%> const&, %1%)", + n, + dist.success_fraction(), + k, + &result, Policy())) + { + return result; + } + + // Special cases of success_fraction, regardless of k successes and regardless of n trials. + if (dist.success_fraction() == 0) + { // probability of zero successes is 1: + return static_cast<RealType>(k == 0 ? 1 : 0); + } + if (dist.success_fraction() == 1) + { // probability of n successes is 1: + return static_cast<RealType>(k == n ? 1 : 0); + } + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + if (n == 0) + { + return 1; // Probability = 1 = certainty. + } + if (k == 0) + { // binomial coeffic (n 0) = 1, + // n ^ 0 = 1 + return pow(1 - dist.success_fraction(), n); + } + if (k == n) + { // binomial coeffic (n n) = 1, + // n ^ 0 = 1 + return pow(dist.success_fraction(), k); // * pow((1 - dist.success_fraction()), (n - k)) = 1 + } + + // Probability of getting exactly k successes + // if C(n, k) is the binomial coefficient then: + // + // f(k; n,p) = C(n, k) * p^k * (1-p)^(n-k) + // = (n!/(k!(n-k)!)) * p^k * (1-p)^(n-k) + // = (tgamma(n+1) / (tgamma(k+1)*tgamma(n-k+1))) * p^k * (1-p)^(n-k) + // = p^k (1-p)^(n-k) / (beta(k+1, n-k+1) * (n+1)) + // = ibeta_derivative(k+1, n-k+1, p) / (n+1) + // + using boost::math::ibeta_derivative; // a, b, x + return ibeta_derivative(k+1, n-k+1, dist.success_fraction(), Policy()) / (n+1); + + } // pdf + + template <class RealType, class Policy> + inline RealType cdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k) + { // Cumulative Distribution Function Binomial. + // The random variate k is the number of successes in n trials. + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + + // Returns the sum of the terms 0 through k of the Binomial Probability Density/Mass: + // + // i=k + // -- ( n ) i n-i + // > | | p (1-p) + // -- ( i ) + // i=0 + + // The terms are not summed directly instead + // the incomplete beta integral is employed, + // according to the formula: + // P = I[1-p]( n-k, k+1). + // = 1 - I[p](k + 1, n - k) + + BOOST_MATH_STD_USING // for ADL of std functions + + RealType n = dist.trials(); + RealType p = dist.success_fraction(); + + // Error check: + RealType result = 0; + if(false == binomial_detail::check_dist_and_k( + "boost::math::cdf(binomial_distribution<%1%> const&, %1%)", + n, + p, + k, + &result, Policy())) + { + return result; + } + if (k == n) + { + return 1; + } + + // Special cases, regardless of k. + if (p == 0) + { // This need explanation: + // the pdf is zero for all cases except when k == 0. + // For zero p the probability of zero successes is one. + // Therefore the cdf is always 1: + // the probability of k or *fewer* successes is always 1 + // if there are never any successes! + return 1; + } + if (p == 1) + { // This is correct but needs explanation: + // when k = 1 + // all the cdf and pdf values are zero *except* when k == n, + // and that case has been handled above already. + return 0; + } + // + // P = I[1-p](n - k, k + 1) + // = 1 - I[p](k + 1, n - k) + // Use of ibetac here prevents cancellation errors in calculating + // 1-p if p is very small, perhaps smaller than machine epsilon. + // + // Note that we do not use a finite sum here, since the incomplete + // beta uses a finite sum internally for integer arguments, so + // we'll just let it take care of the necessary logic. + // + return ibetac(k + 1, n - k, p, Policy()); + } // binomial cdf + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c) + { // Complemented Cumulative Distribution Function Binomial. + // The random variate k is the number of successes in n trials. + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + + // Returns the sum of the terms k+1 through n of the Binomial Probability Density/Mass: + // + // i=n + // -- ( n ) i n-i + // > | | p (1-p) + // -- ( i ) + // i=k+1 + + // The terms are not summed directly instead + // the incomplete beta integral is employed, + // according to the formula: + // Q = 1 -I[1-p]( n-k, k+1). + // = I[p](k + 1, n - k) + + BOOST_MATH_STD_USING // for ADL of std functions + + RealType const& k = c.param; + binomial_distribution<RealType, Policy> const& dist = c.dist; + RealType n = dist.trials(); + RealType p = dist.success_fraction(); + + // Error checks: + RealType result = 0; + if(false == binomial_detail::check_dist_and_k( + "boost::math::cdf(binomial_distribution<%1%> const&, %1%)", + n, + p, + k, + &result, Policy())) + { + return result; + } + + if (k == n) + { // Probability of greater than n successes is necessarily zero: + return 0; + } + + // Special cases, regardless of k. + if (p == 0) + { + // This need explanation: the pdf is zero for all + // cases except when k == 0. For zero p the probability + // of zero successes is one. Therefore the cdf is always + // 1: the probability of *more than* k successes is always 0 + // if there are never any successes! + return 0; + } + if (p == 1) + { + // This needs explanation, when p = 1 + // we always have n successes, so the probability + // of more than k successes is 1 as long as k < n. + // The k == n case has already been handled above. + return 1; + } + // + // Calculate cdf binomial using the incomplete beta function. + // Q = 1 -I[1-p](n - k, k + 1) + // = I[p](k + 1, n - k) + // Use of ibeta here prevents cancellation errors in calculating + // 1-p if p is very small, perhaps smaller than machine epsilon. + // + // Note that we do not use a finite sum here, since the incomplete + // beta uses a finite sum internally for integer arguments, so + // we'll just let it take care of the necessary logic. + // + return ibeta(k + 1, n - k, p, Policy()); + } // binomial cdf + + template <class RealType, class Policy> + inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p) + { + return binomial_detail::quantile_imp(dist, p, RealType(1-p), false); + } // quantile + + template <class RealType, class Policy> + RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c) + { + return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true); + } // quantile + + template <class RealType, class Policy> + inline RealType mode(const binomial_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // ADL of std functions. + RealType p = dist.success_fraction(); + RealType n = dist.trials(); + return floor(p * (n + 1)); + } + + template <class RealType, class Policy> + inline RealType median(const binomial_distribution<RealType, Policy>& dist) + { // Bounds for the median of the negative binomial distribution + // VAN DE VEN R. ; WEBER N. C. ; + // Univ. Sydney, school mathematics statistics, Sydney N.S.W. 2006, AUSTRALIE + // Metrika (Metrika) ISSN 0026-1335 CODEN MTRKA8 + // 1993, vol. 40, no3-4, pp. 185-189 (4 ref.) + + // Bounds for median and 50 percetage point of binomial and negative binomial distribution + // Metrika, ISSN 0026-1335 (Print) 1435-926X (Online) + // Volume 41, Number 1 / December, 1994, DOI 10.1007/BF01895303 + BOOST_MATH_STD_USING // ADL of std functions. + RealType p = dist.success_fraction(); + RealType n = dist.trials(); + // Wikipedia says one of floor(np) -1, floor (np), floor(np) +1 + return floor(p * n); // Chose the middle value. + } + + template <class RealType, class Policy> + inline RealType skewness(const binomial_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // ADL of std functions. + RealType p = dist.success_fraction(); + RealType n = dist.trials(); + return (1 - 2 * p) / sqrt(n * p * (1 - p)); + } + + template <class RealType, class Policy> + inline RealType kurtosis(const binomial_distribution<RealType, Policy>& dist) + { + RealType p = dist.success_fraction(); + RealType n = dist.trials(); + return 3 - 6 / n + 1 / (n * p * (1 - p)); + } + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const binomial_distribution<RealType, Policy>& dist) + { + RealType p = dist.success_fraction(); + RealType q = 1 - p; + RealType n = dist.trials(); + return (1 - 6 * p * q) / (n * p * q); + } + + } // namespace math + } // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_MATH_SPECIAL_BINOMIAL_HPP + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/cauchy.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/cauchy.hpp new file mode 100644 index 00000000000..5a3a64f0f2c --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/cauchy.hpp @@ -0,0 +1,362 @@ +// Copyright John Maddock 2006, 2007. +// Copyright Paul A. Bristow 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_CAUCHY_HPP +#define BOOST_STATS_CAUCHY_HPP + +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable : 4127) // conditional expression is constant +#endif + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/distributions/complement.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/config/no_tr1/cmath.hpp> + +#include <utility> + +namespace boost{ namespace math +{ + +template <class RealType, class Policy> +class cauchy_distribution; + +namespace detail +{ + +template <class RealType, class Policy> +RealType cdf_imp(const cauchy_distribution<RealType, Policy>& dist, const RealType& x, bool complement) +{ + // + // This calculates the cdf of the Cauchy distribution and/or its complement. + // + // The usual formula for the Cauchy cdf is: + // + // cdf = 0.5 + atan(x)/pi + // + // But that suffers from cancellation error as x -> -INF. + // + // Recall that for x < 0: + // + // atan(x) = -pi/2 - atan(1/x) + // + // Substituting into the above we get: + // + // CDF = -atan(1/x) ; x < 0 + // + // So the proceedure is to calculate the cdf for -fabs(x) + // using the above formula, and then subtract from 1 when required + // to get the result. + // + BOOST_MATH_STD_USING // for ADL of std functions + static const char* function = "boost::math::cdf(cauchy<%1%>&, %1%)"; + RealType result = 0; + RealType location = dist.location(); + RealType scale = dist.scale(); + if(false == detail::check_location(function, location, &result, Policy())) + { + return result; + } + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) + { // cdf +infinity is unity. + return static_cast<RealType>((complement) ? 0 : 1); + } + if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) + { // cdf -infinity is zero. + return static_cast<RealType>((complement) ? 1 : 0); + } + if(false == detail::check_x(function, x, &result, Policy())) + { // Catches x == NaN + return result; + } + RealType mx = -fabs((x - location) / scale); // scale is > 0 + if(mx > -tools::epsilon<RealType>() / 8) + { // special case first: x extremely close to location. + return 0.5; + } + result = -atan(1 / mx) / constants::pi<RealType>(); + return (((x > location) != complement) ? 1 - result : result); +} // cdf + +template <class RealType, class Policy> +RealType quantile_imp( + const cauchy_distribution<RealType, Policy>& dist, + const RealType& p, + bool complement) +{ + // This routine implements the quantile for the Cauchy distribution, + // the value p may be the probability, or its complement if complement=true. + // + // The procedure first performs argument reduction on p to avoid error + // when calculating the tangent, then calulates the distance from the + // mid-point of the distribution. This is either added or subtracted + // from the location parameter depending on whether `complement` is true. + // + static const char* function = "boost::math::quantile(cauchy<%1%>&, %1%)"; + BOOST_MATH_STD_USING // for ADL of std functions + + RealType result = 0; + RealType location = dist.location(); + RealType scale = dist.scale(); + if(false == detail::check_location(function, location, &result, Policy())) + { + return result; + } + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(false == detail::check_probability(function, p, &result, Policy())) + { + return result; + } + // Special cases: + if(p == 1) + { + return (complement ? -1 : 1) * policies::raise_overflow_error<RealType>(function, 0, Policy()); + } + if(p == 0) + { + return (complement ? 1 : -1) * policies::raise_overflow_error<RealType>(function, 0, Policy()); + } + + RealType P = p - floor(p); // argument reduction of p: + if(P > 0.5) + { + P = P - 1; + } + if(P == 0.5) // special case: + { + return location; + } + result = -scale / tan(constants::pi<RealType>() * P); + return complement ? RealType(location - result) : RealType(location + result); +} // quantile + +} // namespace detail + +template <class RealType = double, class Policy = policies::policy<> > +class cauchy_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + cauchy_distribution(RealType l_location = 0, RealType l_scale = 1) + : m_a(l_location), m_hg(l_scale) + { + static const char* function = "boost::math::cauchy_distribution<%1%>::cauchy_distribution"; + RealType result; + detail::check_location(function, l_location, &result, Policy()); + detail::check_scale(function, l_scale, &result, Policy()); + } // cauchy_distribution + + RealType location()const + { + return m_a; + } + RealType scale()const + { + return m_hg; + } + +private: + RealType m_a; // The location, this is the median of the distribution. + RealType m_hg; // The scale )or shape), this is the half width at half height. +}; + +typedef cauchy_distribution<double> cauchy; + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const cauchy_distribution<RealType, Policy>&) +{ // Range of permissible values for random variable x. + if (std::numeric_limits<RealType>::has_infinity) + { + return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity. + } + else + { // Can only use max_value. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max. + } +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const cauchy_distribution<RealType, Policy>& ) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + if (std::numeric_limits<RealType>::has_infinity) + { + return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity. + } + else + { // Can only use max_value. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-tools::max_value<RealType>(), max_value<RealType>()); // - to + max. + } +} + +template <class RealType, class Policy> +inline RealType pdf(const cauchy_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::pdf(cauchy<%1%>&, %1%)"; + RealType result = 0; + RealType location = dist.location(); + RealType scale = dist.scale(); + if(false == detail::check_scale("boost::math::pdf(cauchy<%1%>&, %1%)", scale, &result, Policy())) + { + return result; + } + if(false == detail::check_location("boost::math::pdf(cauchy<%1%>&, %1%)", location, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + return 0; // pdf + and - infinity is zero. + } + // These produce MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity()) + //{ // pdf + and - infinity is zero. + // return 0; + //} + + if(false == detail::check_x(function, x, &result, Policy())) + { // Catches x = NaN + return result; + } + + RealType xs = (x - location) / scale; + result = 1 / (constants::pi<RealType>() * scale * (1 + xs * xs)); + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const cauchy_distribution<RealType, Policy>& dist, const RealType& x) +{ + return detail::cdf_imp(dist, x, false); +} // cdf + +template <class RealType, class Policy> +inline RealType quantile(const cauchy_distribution<RealType, Policy>& dist, const RealType& p) +{ + return detail::quantile_imp(dist, p, false); +} // quantile + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<cauchy_distribution<RealType, Policy>, RealType>& c) +{ + return detail::cdf_imp(c.dist, c.param, true); +} // cdf complement + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<cauchy_distribution<RealType, Policy>, RealType>& c) +{ + return detail::quantile_imp(c.dist, c.param, true); +} // quantile complement + +template <class RealType, class Policy> +inline RealType mean(const cauchy_distribution<RealType, Policy>&) +{ // There is no mean: + typedef typename Policy::assert_undefined_type assert_type; + BOOST_STATIC_ASSERT(assert_type::value == 0); + + return policies::raise_domain_error<RealType>( + "boost::math::mean(cauchy<%1%>&)", + "The Cauchy distribution does not have a mean: " + "the only possible return value is %1%.", + std::numeric_limits<RealType>::quiet_NaN(), Policy()); +} + +template <class RealType, class Policy> +inline RealType variance(const cauchy_distribution<RealType, Policy>& /*dist*/) +{ + // There is no variance: + typedef typename Policy::assert_undefined_type assert_type; + BOOST_STATIC_ASSERT(assert_type::value == 0); + + return policies::raise_domain_error<RealType>( + "boost::math::variance(cauchy<%1%>&)", + "The Cauchy distribution does not have a variance: " + "the only possible return value is %1%.", + std::numeric_limits<RealType>::quiet_NaN(), Policy()); +} + +template <class RealType, class Policy> +inline RealType mode(const cauchy_distribution<RealType, Policy>& dist) +{ + return dist.location(); +} + +template <class RealType, class Policy> +inline RealType median(const cauchy_distribution<RealType, Policy>& dist) +{ + return dist.location(); +} +template <class RealType, class Policy> +inline RealType skewness(const cauchy_distribution<RealType, Policy>& /*dist*/) +{ + // There is no skewness: + typedef typename Policy::assert_undefined_type assert_type; + BOOST_STATIC_ASSERT(assert_type::value == 0); + + return policies::raise_domain_error<RealType>( + "boost::math::skewness(cauchy<%1%>&)", + "The Cauchy distribution does not have a skewness: " + "the only possible return value is %1%.", + std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? +} + +template <class RealType, class Policy> +inline RealType kurtosis(const cauchy_distribution<RealType, Policy>& /*dist*/) +{ + // There is no kurtosis: + typedef typename Policy::assert_undefined_type assert_type; + BOOST_STATIC_ASSERT(assert_type::value == 0); + + return policies::raise_domain_error<RealType>( + "boost::math::kurtosis(cauchy<%1%>&)", + "The Cauchy distribution does not have a kurtosis: " + "the only possible return value is %1%.", + std::numeric_limits<RealType>::quiet_NaN(), Policy()); +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const cauchy_distribution<RealType, Policy>& /*dist*/) +{ + // There is no kurtosis excess: + typedef typename Policy::assert_undefined_type assert_type; + BOOST_STATIC_ASSERT(assert_type::value == 0); + + return policies::raise_domain_error<RealType>( + "boost::math::kurtosis_excess(cauchy<%1%>&)", + "The Cauchy distribution does not have a kurtosis: " + "the only possible return value is %1%.", + std::numeric_limits<RealType>::quiet_NaN(), Policy()); +} + +} // namespace math +} // namespace boost + +#ifdef _MSC_VER +#pragma warning(pop) +#endif + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_CAUCHY_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/chi_squared.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/chi_squared.hpp new file mode 100644 index 00000000000..071c7756f49 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/chi_squared.hpp @@ -0,0 +1,364 @@ +// Copyright John Maddock 2006, 2007. +// Copyright Paul A. Bristow 2008, 2010. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP +#define BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/gamma.hpp> // for incomplete beta. +#include <boost/math/distributions/complement.hpp> // complements +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks +#include <boost/math/special_functions/fpclassify.hpp> + +#include <utility> + +namespace boost{ namespace math{ + +template <class RealType = double, class Policy = policies::policy<> > +class chi_squared_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + chi_squared_distribution(RealType i) : m_df(i) + { + RealType result; + detail::check_df( + "boost::math::chi_squared_distribution<%1%>::chi_squared_distribution", m_df, &result, Policy()); + } // chi_squared_distribution + + RealType degrees_of_freedom()const + { + return m_df; + } + + // Parameter estimation: + static RealType find_degrees_of_freedom( + RealType difference_from_variance, + RealType alpha, + RealType beta, + RealType variance, + RealType hint = 100); + +private: + // + // Data member: + // + RealType m_df; // degrees of freedom is a positive real number. +}; // class chi_squared_distribution + +typedef chi_squared_distribution<double> chi_squared; + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4127) +#endif + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const chi_squared_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x. + if (std::numeric_limits<RealType>::has_infinity) + { + return std::pair<RealType, RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()); // 0 to + infinity. + } + else + { + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + max. + } +} + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const chi_squared_distribution<RealType, Policy>& /*dist*/) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity. +} + +template <class RealType, class Policy> +RealType pdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square) +{ + BOOST_MATH_STD_USING // for ADL of std functions + RealType degrees_of_freedom = dist.degrees_of_freedom(); + // Error check: + RealType error_result; + + static const char* function = "boost::math::pdf(const chi_squared_distribution<%1%>&, %1%)"; + + if(false == detail::check_df( + function, degrees_of_freedom, &error_result, Policy())) + return error_result; + + if((chi_square < 0) || !(boost::math::isfinite)(chi_square)) + { + return policies::raise_domain_error<RealType>( + function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy()); + } + + if(chi_square == 0) + { + // Handle special cases: + if(degrees_of_freedom < 2) + { + return policies::raise_overflow_error<RealType>( + function, 0, Policy()); + } + else if(degrees_of_freedom == 2) + { + return 0.5f; + } + else + { + return 0; + } + } + + return gamma_p_derivative(degrees_of_freedom / 2, chi_square / 2, Policy()) / 2; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square) +{ + RealType degrees_of_freedom = dist.degrees_of_freedom(); + // Error check: + RealType error_result; + static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; + + if(false == detail::check_df( + function, degrees_of_freedom, &error_result, Policy())) + return error_result; + + if((chi_square < 0) || !(boost::math::isfinite)(chi_square)) + { + return policies::raise_domain_error<RealType>( + function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy()); + } + + return boost::math::gamma_p(degrees_of_freedom / 2, chi_square / 2, Policy()); +} // cdf + +template <class RealType, class Policy> +inline RealType quantile(const chi_squared_distribution<RealType, Policy>& dist, const RealType& p) +{ + RealType degrees_of_freedom = dist.degrees_of_freedom(); + static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; + // Error check: + RealType error_result; + if(false == + ( + detail::check_df(function, degrees_of_freedom, &error_result, Policy()) + && detail::check_probability(function, p, &error_result, Policy())) + ) + return error_result; + + return 2 * boost::math::gamma_p_inv(degrees_of_freedom / 2, p, Policy()); +} // quantile + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c) +{ + RealType const& degrees_of_freedom = c.dist.degrees_of_freedom(); + RealType const& chi_square = c.param; + static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; + // Error check: + RealType error_result; + if(false == detail::check_df( + function, degrees_of_freedom, &error_result, Policy())) + return error_result; + + if((chi_square < 0) || !(boost::math::isfinite)(chi_square)) + { + return policies::raise_domain_error<RealType>( + function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy()); + } + + return boost::math::gamma_q(degrees_of_freedom / 2, chi_square / 2, Policy()); +} + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c) +{ + RealType const& degrees_of_freedom = c.dist.degrees_of_freedom(); + RealType const& q = c.param; + static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; + // Error check: + RealType error_result; + if(false == ( + detail::check_df(function, degrees_of_freedom, &error_result, Policy()) + && detail::check_probability(function, q, &error_result, Policy())) + ) + return error_result; + + return 2 * boost::math::gamma_q_inv(degrees_of_freedom / 2, q, Policy()); +} + +template <class RealType, class Policy> +inline RealType mean(const chi_squared_distribution<RealType, Policy>& dist) +{ // Mean of Chi-Squared distribution = v. + return dist.degrees_of_freedom(); +} // mean + +template <class RealType, class Policy> +inline RealType variance(const chi_squared_distribution<RealType, Policy>& dist) +{ // Variance of Chi-Squared distribution = 2v. + return 2 * dist.degrees_of_freedom(); +} // variance + +template <class RealType, class Policy> +inline RealType mode(const chi_squared_distribution<RealType, Policy>& dist) +{ + RealType df = dist.degrees_of_freedom(); + static const char* function = "boost::math::mode(const chi_squared_distribution<%1%>&)"; + // Most sources only define mode for df >= 2, + // but for 0 <= df <= 2, the pdf maximum actually occurs at random variate = 0; + // So one could extend the definition of mode thus: + //if(df < 0) + //{ + // return policies::raise_domain_error<RealType>( + // function, + // "Chi-Squared distribution only has a mode for degrees of freedom >= 0, but got degrees of freedom = %1%.", + // df, Policy()); + //} + //return (df <= 2) ? 0 : df - 2; + + if(df < 2) + return policies::raise_domain_error<RealType>( + function, + "Chi-Squared distribution only has a mode for degrees of freedom >= 2, but got degrees of freedom = %1%.", + df, Policy()); + return df - 2; +} + +//template <class RealType, class Policy> +//inline RealType median(const chi_squared_distribution<RealType, Policy>& dist) +//{ // Median is given by Quantile[dist, 1/2] +// RealType df = dist.degrees_of_freedom(); +// if(df <= 1) +// return tools::domain_error<RealType>( +// BOOST_CURRENT_FUNCTION, +// "The Chi-Squared distribution only has a mode for degrees of freedom >= 2, but got degrees of freedom = %1%.", +// df); +// return df - RealType(2)/3; +//} +// Now implemented via quantile(half) in derived accessors. + +template <class RealType, class Policy> +inline RealType skewness(const chi_squared_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // For ADL + RealType df = dist.degrees_of_freedom(); + return sqrt (8 / df); // == 2 * sqrt(2 / df); +} + +template <class RealType, class Policy> +inline RealType kurtosis(const chi_squared_distribution<RealType, Policy>& dist) +{ + RealType df = dist.degrees_of_freedom(); + return 3 + 12 / df; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const chi_squared_distribution<RealType, Policy>& dist) +{ + RealType df = dist.degrees_of_freedom(); + return 12 / df; +} + +// +// Parameter estimation comes last: +// +namespace detail +{ + +template <class RealType, class Policy> +struct df_estimator +{ + df_estimator(RealType a, RealType b, RealType variance, RealType delta) + : alpha(a), beta(b), ratio(delta/variance) + { // Constructor + } + + RealType operator()(const RealType& df) + { + if(df <= tools::min_value<RealType>()) + return 1; + chi_squared_distribution<RealType, Policy> cs(df); + + RealType result; + if(ratio > 0) + { + RealType r = 1 + ratio; + result = cdf(cs, quantile(complement(cs, alpha)) / r) - beta; + } + else + { // ratio <= 0 + RealType r = 1 + ratio; + result = cdf(complement(cs, quantile(cs, alpha) / r)) - beta; + } + return result; + } +private: + RealType alpha; + RealType beta; + RealType ratio; // Difference from variance / variance, so fractional. +}; + +} // namespace detail + +template <class RealType, class Policy> +RealType chi_squared_distribution<RealType, Policy>::find_degrees_of_freedom( + RealType difference_from_variance, + RealType alpha, + RealType beta, + RealType variance, + RealType hint) +{ + static const char* function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)"; + // Check for domain errors: + RealType error_result; + if(false == + detail::check_probability(function, alpha, &error_result, Policy()) + && detail::check_probability(function, beta, &error_result, Policy())) + { // Either probability is outside 0 to 1. + return error_result; + } + + if(hint <= 0) + { // No hint given, so guess df = 1. + hint = 1; + } + + detail::df_estimator<RealType, Policy> f(alpha, beta, variance, difference_from_variance); + tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); + boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); + std::pair<RealType, RealType> r = + tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy()); + RealType result = r.first + (r.second - r.first) / 2; + if(max_iter >= policies::get_max_root_iterations<Policy>()) + { + policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" + " either there is no answer to how many degrees of freedom are required" + " or the answer is infinite. Current best guess is %1%", result, Policy()); + } + return result; +} + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/complement.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/complement.hpp new file mode 100644 index 00000000000..26d0d49e6de --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/complement.hpp @@ -0,0 +1,195 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Paul A. Bristow 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_COMPLEMENT_HPP +#define BOOST_STATS_COMPLEMENT_HPP + +// +// This code really defines our own tuple type. +// It would be nice to reuse boost::math::tuple +// while retaining our own type safety, but it's +// not clear if that's possible. In any case this +// code is *very* lightweight. +// +namespace boost{ namespace math{ + +template <class Dist, class RealType> +struct complemented2_type +{ + complemented2_type( + const Dist& d, + const RealType& p1) + : dist(d), + param(p1) {} + + const Dist& dist; + const RealType& param; + +private: + complemented2_type& operator=(const complemented2_type&); +}; + +template <class Dist, class RealType1, class RealType2> +struct complemented3_type +{ + complemented3_type( + const Dist& d, + const RealType1& p1, + const RealType2& p2) + : dist(d), + param1(p1), + param2(p2) {} + + const Dist& dist; + const RealType1& param1; + const RealType2& param2; +private: + complemented3_type& operator=(const complemented3_type&); +}; + +template <class Dist, class RealType1, class RealType2, class RealType3> +struct complemented4_type +{ + complemented4_type( + const Dist& d, + const RealType1& p1, + const RealType2& p2, + const RealType3& p3) + : dist(d), + param1(p1), + param2(p2), + param3(p3) {} + + const Dist& dist; + const RealType1& param1; + const RealType2& param2; + const RealType3& param3; +private: + complemented4_type& operator=(const complemented4_type&); +}; + +template <class Dist, class RealType1, class RealType2, class RealType3, class RealType4> +struct complemented5_type +{ + complemented5_type( + const Dist& d, + const RealType1& p1, + const RealType2& p2, + const RealType3& p3, + const RealType4& p4) + : dist(d), + param1(p1), + param2(p2), + param3(p3), + param4(p4) {} + + const Dist& dist; + const RealType1& param1; + const RealType2& param2; + const RealType3& param3; + const RealType4& param4; +private: + complemented5_type& operator=(const complemented5_type&); +}; + +template <class Dist, class RealType1, class RealType2, class RealType3, class RealType4, class RealType5> +struct complemented6_type +{ + complemented6_type( + const Dist& d, + const RealType1& p1, + const RealType2& p2, + const RealType3& p3, + const RealType4& p4, + const RealType5& p5) + : dist(d), + param1(p1), + param2(p2), + param3(p3), + param4(p4), + param5(p5) {} + + const Dist& dist; + const RealType1& param1; + const RealType2& param2; + const RealType3& param3; + const RealType4& param4; + const RealType5& param5; +private: + complemented6_type& operator=(const complemented6_type&); +}; + +template <class Dist, class RealType1, class RealType2, class RealType3, class RealType4, class RealType5, class RealType6> +struct complemented7_type +{ + complemented7_type( + const Dist& d, + const RealType1& p1, + const RealType2& p2, + const RealType3& p3, + const RealType4& p4, + const RealType5& p5, + const RealType6& p6) + : dist(d), + param1(p1), + param2(p2), + param3(p3), + param4(p4), + param5(p5), + param6(p6) {} + + const Dist& dist; + const RealType1& param1; + const RealType2& param2; + const RealType3& param3; + const RealType4& param4; + const RealType5& param5; + const RealType6& param6; +private: + complemented7_type& operator=(const complemented7_type&); +}; + +template <class Dist, class RealType> +inline complemented2_type<Dist, RealType> complement(const Dist& d, const RealType& r) +{ + return complemented2_type<Dist, RealType>(d, r); +} + +template <class Dist, class RealType1, class RealType2> +inline complemented3_type<Dist, RealType1, RealType2> complement(const Dist& d, const RealType1& r1, const RealType2& r2) +{ + return complemented3_type<Dist, RealType1, RealType2>(d, r1, r2); +} + +template <class Dist, class RealType1, class RealType2, class RealType3> +inline complemented4_type<Dist, RealType1, RealType2, RealType3> complement(const Dist& d, const RealType1& r1, const RealType2& r2, const RealType3& r3) +{ + return complemented4_type<Dist, RealType1, RealType2, RealType3>(d, r1, r2, r3); +} + +template <class Dist, class RealType1, class RealType2, class RealType3, class RealType4> +inline complemented5_type<Dist, RealType1, RealType2, RealType3, RealType4> complement(const Dist& d, const RealType1& r1, const RealType2& r2, const RealType3& r3, const RealType4& r4) +{ + return complemented5_type<Dist, RealType1, RealType2, RealType3, RealType4>(d, r1, r2, r3, r4); +} + +template <class Dist, class RealType1, class RealType2, class RealType3, class RealType4, class RealType5> +inline complemented6_type<Dist, RealType1, RealType2, RealType3, RealType4, RealType5> complement(const Dist& d, const RealType1& r1, const RealType2& r2, const RealType3& r3, const RealType4& r4, const RealType5& r5) +{ + return complemented6_type<Dist, RealType1, RealType2, RealType3, RealType4, RealType5>(d, r1, r2, r3, r4, r5); +} + +template <class Dist, class RealType1, class RealType2, class RealType3, class RealType4, class RealType5, class RealType6> +inline complemented7_type<Dist, RealType1, RealType2, RealType3, RealType4, RealType5, RealType6> complement(const Dist& d, const RealType1& r1, const RealType2& r2, const RealType3& r3, const RealType4& r4, const RealType5& r5, const RealType6& r6) +{ + return complemented7_type<Dist, RealType1, RealType2, RealType3, RealType4, RealType5, RealType6>(d, r1, r2, r3, r4, r5, r6); +} + +} // namespace math +} // namespace boost + +#endif // BOOST_STATS_COMPLEMENT_HPP + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/detail/common_error_handling.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/detail/common_error_handling.hpp new file mode 100644 index 00000000000..71a771e5a14 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/detail/common_error_handling.hpp @@ -0,0 +1,194 @@ +// Copyright John Maddock 2006, 2007. +// Copyright Paul A. Bristow 2006, 2007, 2012. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DISTRIBUTIONS_COMMON_ERROR_HANDLING_HPP +#define BOOST_MATH_DISTRIBUTIONS_COMMON_ERROR_HANDLING_HPP + +#include <boost/math/policies/error_handling.hpp> +#include <boost/math/special_functions/fpclassify.hpp> +// using boost::math::isfinite; +// using boost::math::isnan; + +namespace boost{ namespace math{ namespace detail +{ + +template <class RealType, class Policy> +inline bool check_probability(const char* function, RealType const& prob, RealType* result, const Policy& pol) +{ + if((prob < 0) || (prob > 1) || !(boost::math::isfinite)(prob)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Probability argument is %1%, but must be >= 0 and <= 1 !", prob, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_df(const char* function, RealType const& df, RealType* result, const Policy& pol) +{ // df > 0 but NOT +infinity allowed. + if((df <= 0) || !(boost::math::isfinite)(df)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Degrees of freedom argument is %1%, but must be > 0 !", df, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_df_gt0_to_inf(const char* function, RealType const& df, RealType* result, const Policy& pol) +{ // df > 0 or +infinity are allowed. + if( (df <= 0) || (boost::math::isnan)(df) ) + { // is bad df <= 0 or NaN or -infinity. + *result = policies::raise_domain_error<RealType>( + function, + "Degrees of freedom argument is %1%, but must be > 0 !", df, pol); + return false; + } + return true; +} // check_df_gt0_to_inf + + +template <class RealType, class Policy> +inline bool check_scale( + const char* function, + RealType scale, + RealType* result, + const Policy& pol) +{ + if((scale <= 0) || !(boost::math::isfinite)(scale)) + { // Assume scale == 0 is NOT valid for any distribution. + *result = policies::raise_domain_error<RealType>( + function, + "Scale parameter is %1%, but must be > 0 !", scale, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_location( + const char* function, + RealType location, + RealType* result, + const Policy& pol) +{ + if(!(boost::math::isfinite)(location)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Location parameter is %1%, but must be finite!", location, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_x( + const char* function, + RealType x, + RealType* result, + const Policy& pol) +{ + // Note that this test catches both infinity and NaN. + // Some distributions permit x to be infinite, so these must be tested 1st and return, + // leaving this test to catch any NaNs. + // See Normal, Logistic, Laplace and Cauchy for example. + if(!(boost::math::isfinite)(x)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Random variate x is %1%, but must be finite!", x, pol); + return false; + } + return true; +} // bool check_x + +template <class RealType, class Policy> +inline bool check_x_gt0( + const char* function, + RealType x, + RealType* result, + const Policy& pol) +{ + if(x <= 0) + { + *result = policies::raise_domain_error<RealType>( + function, + "Random variate x is %1%, but must be > 0!", x, pol); + return false; + } + + return true; + // Note that this test catches both infinity and NaN. + // Some special cases permit x to be infinite, so these must be tested 1st, + // leaving this test to catch any NaNs. See Normal and cauchy for example. +} // bool check_x_gt0 + +template <class RealType, class Policy> +inline bool check_positive_x( + const char* function, + RealType x, + RealType* result, + const Policy& pol) +{ + if(!(boost::math::isfinite)(x) || (x < 0)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Random variate x is %1%, but must be finite and >= 0!", x, pol); + return false; + } + return true; + // Note that this test catches both infinity and NaN. + // Some special cases permit x to be infinite, so these must be tested 1st, + // leaving this test to catch any NaNs. see Normal and cauchy for example. +} + +template <class RealType, class Policy> +inline bool check_non_centrality( + const char* function, + RealType ncp, + RealType* result, + const Policy& pol) +{ + if((ncp < 0) || !(boost::math::isfinite)(ncp)) + { // Assume scale == 0 is NOT valid for any distribution. + *result = policies::raise_domain_error<RealType>( + function, + "Non centrality parameter is %1%, but must be > 0 !", ncp, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_finite( + const char* function, + RealType x, + RealType* result, + const Policy& pol) +{ + if(!(boost::math::isfinite)(x)) + { // Assume scale == 0 is NOT valid for any distribution. + *result = policies::raise_domain_error<RealType>( + function, + "Parameter is %1%, but must be finite !", x, pol); + return false; + } + return true; +} + +} // namespace detail +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_DISTRIBUTIONS_COMMON_ERROR_HANDLING_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/detail/derived_accessors.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/detail/derived_accessors.hpp new file mode 100644 index 00000000000..00f5a93258c --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/detail/derived_accessors.hpp @@ -0,0 +1,163 @@ +// Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_DERIVED_HPP +#define BOOST_STATS_DERIVED_HPP + +// This file implements various common properties of distributions +// that can be implemented in terms of other properties: +// variance OR standard deviation (see note below), +// hazard, cumulative hazard (chf), coefficient_of_variation. +// +// Note that while both variance and standard_deviation are provided +// here, each distribution MUST SPECIALIZE AT LEAST ONE OF THESE +// otherwise these two versions will just call each other over and over +// until stack space runs out ... + +// Of course there may be more efficient means of implementing these +// that are specific to a particular distribution, but these generic +// versions give these properties "for free" with most distributions. +// +// In order to make use of this header, it must be included AT THE END +// of the distribution header, AFTER the distribution and its core +// property accessors have been defined: this is so that compilers +// that implement 2-phase lookup and early-type-checking of templates +// can find the definitions refered to herein. +// + +#include <boost/type_traits/is_same.hpp> +#include <boost/static_assert.hpp> + +#ifdef BOOST_MSVC +# pragma warning(push) +# pragma warning(disable: 4723) // potential divide by 0 +// Suppressing spurious warning in coefficient_of_variation +#endif + +namespace boost{ namespace math{ + +template <class Distribution> +typename Distribution::value_type variance(const Distribution& dist); + +template <class Distribution> +inline typename Distribution::value_type standard_deviation(const Distribution& dist) +{ + BOOST_MATH_STD_USING // ADL of sqrt. + return sqrt(variance(dist)); +} + +template <class Distribution> +inline typename Distribution::value_type variance(const Distribution& dist) +{ + typename Distribution::value_type result = standard_deviation(dist); + return result * result; +} + +template <class Distribution, class RealType> +inline typename Distribution::value_type hazard(const Distribution& dist, const RealType& x) +{ // hazard function + // http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ + typedef typename Distribution::value_type value_type; + typedef typename Distribution::policy_type policy_type; + value_type p = cdf(complement(dist, x)); + value_type d = pdf(dist, x); + if(d > p * tools::max_value<value_type>()) + return policies::raise_overflow_error<value_type>( + "boost::math::hazard(const Distribution&, %1%)", 0, policy_type()); + if(d == 0) + { + // This protects against 0/0, but is it the right thing to do? + return 0; + } + return d / p; +} + +template <class Distribution, class RealType> +inline typename Distribution::value_type chf(const Distribution& dist, const RealType& x) +{ // cumulative hazard function. + // http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ + BOOST_MATH_STD_USING + return -log(cdf(complement(dist, x))); +} + +template <class Distribution> +inline typename Distribution::value_type coefficient_of_variation(const Distribution& dist) +{ + typedef typename Distribution::value_type value_type; + typedef typename Distribution::policy_type policy_type; + + using std::abs; + + value_type m = mean(dist); + value_type d = standard_deviation(dist); + if((abs(m) < 1) && (d > abs(m) * tools::max_value<value_type>())) + { // Checks too that m is not zero, + return policies::raise_overflow_error<value_type>("boost::math::coefficient_of_variation(const Distribution&, %1%)", 0, policy_type()); + } + return d / m; // so MSVC warning on zerodivide is spurious, and suppressed. +} +// +// Next follow overloads of some of the standard accessors with mixed +// argument types. We just use a typecast to forward on to the "real" +// implementation with all arguments of the same type: +// +template <class Distribution, class RealType> +inline typename Distribution::value_type pdf(const Distribution& dist, const RealType& x) +{ + typedef typename Distribution::value_type value_type; + return pdf(dist, static_cast<value_type>(x)); +} +template <class Distribution, class RealType> +inline typename Distribution::value_type cdf(const Distribution& dist, const RealType& x) +{ + typedef typename Distribution::value_type value_type; + return cdf(dist, static_cast<value_type>(x)); +} +template <class Distribution, class RealType> +inline typename Distribution::value_type quantile(const Distribution& dist, const RealType& x) +{ + typedef typename Distribution::value_type value_type; + return quantile(dist, static_cast<value_type>(x)); +} +/* +template <class Distribution, class RealType> +inline typename Distribution::value_type chf(const Distribution& dist, const RealType& x) +{ + typedef typename Distribution::value_type value_type; + return chf(dist, static_cast<value_type>(x)); +} +*/ +template <class Distribution, class RealType> +inline typename Distribution::value_type cdf(const complemented2_type<Distribution, RealType>& c) +{ + typedef typename Distribution::value_type value_type; + return cdf(complement(c.dist, static_cast<value_type>(c.param))); +} + +template <class Distribution, class RealType> +inline typename Distribution::value_type quantile(const complemented2_type<Distribution, RealType>& c) +{ + typedef typename Distribution::value_type value_type; + return quantile(complement(c.dist, static_cast<value_type>(c.param))); +} + +template <class Dist> +inline typename Dist::value_type median(const Dist& d) +{ // median - default definition for those distributions for which a + // simple closed form is not known, + // and for which a domain_error and/or NaN generating function is NOT defined. + typedef typename Dist::value_type value_type; + return quantile(d, static_cast<value_type>(0.5f)); +} + +} // namespace math +} // namespace boost + + +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif + +#endif // BOOST_STATS_DERIVED_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/detail/inv_discrete_quantile.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/detail/inv_discrete_quantile.hpp new file mode 100644 index 00000000000..44b0a29d15d --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/detail/inv_discrete_quantile.hpp @@ -0,0 +1,571 @@ +// Copyright John Maddock 2007. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE +#define BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE + +#include <algorithm> + +namespace boost{ namespace math{ namespace detail{ + +// +// Functor for root finding algorithm: +// +template <class Dist> +struct distribution_quantile_finder +{ + typedef typename Dist::value_type value_type; + typedef typename Dist::policy_type policy_type; + + distribution_quantile_finder(const Dist d, value_type p, bool c) + : dist(d), target(p), comp(c) {} + + value_type operator()(value_type const& x) + { + return comp ? value_type(target - cdf(complement(dist, x))) : value_type(cdf(dist, x) - target); + } + +private: + Dist dist; + value_type target; + bool comp; +}; +// +// The purpose of adjust_bounds, is to toggle the last bit of the +// range so that both ends round to the same integer, if possible. +// If they do both round the same then we terminate the search +// for the root *very* quickly when finding an integer result. +// At the point that this function is called we know that "a" is +// below the root and "b" above it, so this change can not result +// in the root no longer being bracketed. +// +template <class Real, class Tol> +void adjust_bounds(Real& /* a */, Real& /* b */, Tol const& /* tol */){} + +template <class Real> +void adjust_bounds(Real& /* a */, Real& b, tools::equal_floor const& /* tol */) +{ + BOOST_MATH_STD_USING + b -= tools::epsilon<Real>() * b; +} + +template <class Real> +void adjust_bounds(Real& a, Real& /* b */, tools::equal_ceil const& /* tol */) +{ + BOOST_MATH_STD_USING + a += tools::epsilon<Real>() * a; +} + +template <class Real> +void adjust_bounds(Real& a, Real& b, tools::equal_nearest_integer const& /* tol */) +{ + BOOST_MATH_STD_USING + a += tools::epsilon<Real>() * a; + b -= tools::epsilon<Real>() * b; +} +// +// This is where all the work is done: +// +template <class Dist, class Tolerance> +typename Dist::value_type + do_inverse_discrete_quantile( + const Dist& dist, + const typename Dist::value_type& p, + bool comp, + typename Dist::value_type guess, + const typename Dist::value_type& multiplier, + typename Dist::value_type adder, + const Tolerance& tol, + boost::uintmax_t& max_iter) +{ + typedef typename Dist::value_type value_type; + typedef typename Dist::policy_type policy_type; + + static const char* function = "boost::math::do_inverse_discrete_quantile<%1%>"; + + BOOST_MATH_STD_USING + + distribution_quantile_finder<Dist> f(dist, p, comp); + // + // Max bounds of the distribution: + // + value_type min_bound, max_bound; + boost::math::tie(min_bound, max_bound) = support(dist); + + if(guess > max_bound) + guess = max_bound; + if(guess < min_bound) + guess = min_bound; + + value_type fa = f(guess); + boost::uintmax_t count = max_iter - 1; + value_type fb(fa), a(guess), b =0; // Compiler warning C4701: potentially uninitialized local variable 'b' used + + if(fa == 0) + return guess; + + // + // For small expected results, just use a linear search: + // + if(guess < 10) + { + b = a; + while((a < 10) && (fa * fb >= 0)) + { + if(fb <= 0) + { + a = b; + b = a + 1; + if(b > max_bound) + b = max_bound; + fb = f(b); + --count; + if(fb == 0) + return b; + if(a == b) + return b; // can't go any higher! + } + else + { + b = a; + a = (std::max)(value_type(b - 1), value_type(0)); + if(a < min_bound) + a = min_bound; + fa = f(a); + --count; + if(fa == 0) + return a; + if(a == b) + return a; // We can't go any lower than this! + } + } + } + // + // Try and bracket using a couple of additions first, + // we're assuming that "guess" is likely to be accurate + // to the nearest int or so: + // + else if(adder != 0) + { + // + // If we're looking for a large result, then bump "adder" up + // by a bit to increase our chances of bracketing the root: + // + //adder = (std::max)(adder, 0.001f * guess); + if(fa < 0) + { + b = a + adder; + if(b > max_bound) + b = max_bound; + } + else + { + b = (std::max)(value_type(a - adder), value_type(0)); + if(b < min_bound) + b = min_bound; + } + fb = f(b); + --count; + if(fb == 0) + return b; + if(count && (fa * fb >= 0)) + { + // + // We didn't bracket the root, try + // once more: + // + a = b; + fa = fb; + if(fa < 0) + { + b = a + adder; + if(b > max_bound) + b = max_bound; + } + else + { + b = (std::max)(value_type(a - adder), value_type(0)); + if(b < min_bound) + b = min_bound; + } + fb = f(b); + --count; + } + if(a > b) + { + using std::swap; + swap(a, b); + swap(fa, fb); + } + } + // + // If the root hasn't been bracketed yet, try again + // using the multiplier this time: + // + if((boost::math::sign)(fb) == (boost::math::sign)(fa)) + { + if(fa < 0) + { + // + // Zero is to the right of x2, so walk upwards + // until we find it: + // + while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b)) + { + if(count == 0) + return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); + a = b; + fa = fb; + b *= multiplier; + if(b > max_bound) + b = max_bound; + fb = f(b); + --count; + BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); + } + } + else + { + // + // Zero is to the left of a, so walk downwards + // until we find it: + // + while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b)) + { + if(fabs(a) < tools::min_value<value_type>()) + { + // Escape route just in case the answer is zero! + max_iter -= count; + max_iter += 1; + return 0; + } + if(count == 0) + return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type()); + b = a; + fb = fa; + a /= multiplier; + if(a < min_bound) + a = min_bound; + fa = f(a); + --count; + BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); + } + } + } + max_iter -= count; + if(fa == 0) + return a; + if(fb == 0) + return b; + if(a == b) + return b; // Ran out of bounds trying to bracket - there is no answer! + // + // Adjust bounds so that if we're looking for an integer + // result, then both ends round the same way: + // + adjust_bounds(a, b, tol); + // + // We don't want zero or denorm lower bounds: + // + if(a < tools::min_value<value_type>()) + a = tools::min_value<value_type>(); + // + // Go ahead and find the root: + // + std::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type()); + max_iter += count; + BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); + return (r.first + r.second) / 2; +} +// +// Some special routine for rounding up and down: +// We want to check and see if we are very close to an integer, and if so test to see if +// that integer is an exact root of the cdf. We do this because our root finder only +// guarantees to find *a root*, and there can sometimes be many consecutive floating +// point values which are all roots. This is especially true if the target probability +// is very close 1. +// +template <class Dist> +inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c) +{ + BOOST_MATH_STD_USING + typename Dist::value_type cc = ceil(result); + typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1; + if(pp == p) + result = cc; + else + result = floor(result); + // + // Now find the smallest integer <= result for which we get an exact root: + // + while(result != 0) + { + cc = result - 1; + if(cc < support(d).first) + break; + typename Dist::value_type pp = c ? cdf(complement(d, cc)) : cdf(d, cc); + if(pp == p) + result = cc; + else if(c ? pp > p : pp < p) + break; + result -= 1; + } + + return result; +} + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4127) +#endif + +template <class Dist> +inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c) +{ + BOOST_MATH_STD_USING + typename Dist::value_type cc = floor(result); + typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0; + if(pp == p) + result = cc; + else + result = ceil(result); + // + // Now find the largest integer >= result for which we get an exact root: + // + while(true) + { + cc = result + 1; + if(cc > support(d).second) + break; + typename Dist::value_type pp = c ? cdf(complement(d, cc)) : cdf(d, cc); + if(pp == p) + result = cc; + else if(c ? pp < p : pp > p) + break; + result += 1; + } + + return result; +} + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif +// +// Now finally are the public API functions. +// There is one overload for each policy, +// each one is responsible for selecting the correct +// termination condition, and rounding the result +// to an int where required. +// +template <class Dist> +inline typename Dist::value_type + inverse_discrete_quantile( + const Dist& dist, + typename Dist::value_type p, + bool c, + const typename Dist::value_type& guess, + const typename Dist::value_type& multiplier, + const typename Dist::value_type& adder, + const policies::discrete_quantile<policies::real>&, + boost::uintmax_t& max_iter) +{ + if(p > 0.5) + { + p = 1 - p; + c = !c; + } + typename Dist::value_type pp = c ? 1 - p : p; + if(pp <= pdf(dist, 0)) + return 0; + return do_inverse_discrete_quantile( + dist, + p, + c, + guess, + multiplier, + adder, + tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()), + max_iter); +} + +template <class Dist> +inline typename Dist::value_type + inverse_discrete_quantile( + const Dist& dist, + const typename Dist::value_type& p, + bool c, + const typename Dist::value_type& guess, + const typename Dist::value_type& multiplier, + const typename Dist::value_type& adder, + const policies::discrete_quantile<policies::integer_round_outwards>&, + boost::uintmax_t& max_iter) +{ + typedef typename Dist::value_type value_type; + BOOST_MATH_STD_USING + typename Dist::value_type pp = c ? 1 - p : p; + if(pp <= pdf(dist, 0)) + return 0; + // + // What happens next depends on whether we're looking for an + // upper or lower quantile: + // + if(pp < 0.5f) + return round_to_floor(dist, do_inverse_discrete_quantile( + dist, + p, + c, + (guess < 1 ? value_type(1) : (value_type)floor(guess)), + multiplier, + adder, + tools::equal_floor(), + max_iter), p, c); + // else: + return round_to_ceil(dist, do_inverse_discrete_quantile( + dist, + p, + c, + (value_type)ceil(guess), + multiplier, + adder, + tools::equal_ceil(), + max_iter), p, c); +} + +template <class Dist> +inline typename Dist::value_type + inverse_discrete_quantile( + const Dist& dist, + const typename Dist::value_type& p, + bool c, + const typename Dist::value_type& guess, + const typename Dist::value_type& multiplier, + const typename Dist::value_type& adder, + const policies::discrete_quantile<policies::integer_round_inwards>&, + boost::uintmax_t& max_iter) +{ + typedef typename Dist::value_type value_type; + BOOST_MATH_STD_USING + typename Dist::value_type pp = c ? 1 - p : p; + if(pp <= pdf(dist, 0)) + return 0; + // + // What happens next depends on whether we're looking for an + // upper or lower quantile: + // + if(pp < 0.5f) + return round_to_ceil(dist, do_inverse_discrete_quantile( + dist, + p, + c, + ceil(guess), + multiplier, + adder, + tools::equal_ceil(), + max_iter), p, c); + // else: + return round_to_floor(dist, do_inverse_discrete_quantile( + dist, + p, + c, + (guess < 1 ? value_type(1) : floor(guess)), + multiplier, + adder, + tools::equal_floor(), + max_iter), p, c); +} + +template <class Dist> +inline typename Dist::value_type + inverse_discrete_quantile( + const Dist& dist, + const typename Dist::value_type& p, + bool c, + const typename Dist::value_type& guess, + const typename Dist::value_type& multiplier, + const typename Dist::value_type& adder, + const policies::discrete_quantile<policies::integer_round_down>&, + boost::uintmax_t& max_iter) +{ + typedef typename Dist::value_type value_type; + BOOST_MATH_STD_USING + typename Dist::value_type pp = c ? 1 - p : p; + if(pp <= pdf(dist, 0)) + return 0; + return round_to_floor(dist, do_inverse_discrete_quantile( + dist, + p, + c, + (guess < 1 ? value_type(1) : floor(guess)), + multiplier, + adder, + tools::equal_floor(), + max_iter), p, c); +} + +template <class Dist> +inline typename Dist::value_type + inverse_discrete_quantile( + const Dist& dist, + const typename Dist::value_type& p, + bool c, + const typename Dist::value_type& guess, + const typename Dist::value_type& multiplier, + const typename Dist::value_type& adder, + const policies::discrete_quantile<policies::integer_round_up>&, + boost::uintmax_t& max_iter) +{ + BOOST_MATH_STD_USING + typename Dist::value_type pp = c ? 1 - p : p; + if(pp <= pdf(dist, 0)) + return 0; + return round_to_ceil(dist, do_inverse_discrete_quantile( + dist, + p, + c, + ceil(guess), + multiplier, + adder, + tools::equal_ceil(), + max_iter), p, c); +} + +template <class Dist> +inline typename Dist::value_type + inverse_discrete_quantile( + const Dist& dist, + const typename Dist::value_type& p, + bool c, + const typename Dist::value_type& guess, + const typename Dist::value_type& multiplier, + const typename Dist::value_type& adder, + const policies::discrete_quantile<policies::integer_round_nearest>&, + boost::uintmax_t& max_iter) +{ + typedef typename Dist::value_type value_type; + BOOST_MATH_STD_USING + typename Dist::value_type pp = c ? 1 - p : p; + if(pp <= pdf(dist, 0)) + return 0; + // + // Note that we adjust the guess to the nearest half-integer: + // this increase the chances that we will bracket the root + // with two results that both round to the same integer quickly. + // + return round_to_floor(dist, do_inverse_discrete_quantile( + dist, + p, + c, + (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f), + multiplier, + adder, + tools::equal_nearest_integer(), + max_iter) + 0.5f, p, c); +} + +}}} // namespaces + +#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/exponential.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/exponential.hpp new file mode 100644 index 00000000000..bfe7e6b4ac8 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/exponential.hpp @@ -0,0 +1,275 @@ +// Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_EXPONENTIAL_HPP +#define BOOST_STATS_EXPONENTIAL_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/special_functions/log1p.hpp> +#include <boost/math/special_functions/expm1.hpp> +#include <boost/math/distributions/complement.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/config/no_tr1/cmath.hpp> + +#ifdef BOOST_MSVC +# pragma warning(push) +# pragma warning(disable: 4127) // conditional expression is constant +# pragma warning(disable: 4702) // unreachable code (return after domain_error throw). +#endif + +#include <utility> + +namespace boost{ namespace math{ + +namespace detail{ +// +// Error check: +// +template <class RealType, class Policy> +inline bool verify_lambda(const char* function, RealType l, RealType* presult, const Policy& pol) +{ + if((l <= 0) || !(boost::math::isfinite)(l)) + { + *presult = policies::raise_domain_error<RealType>( + function, + "The scale parameter \"lambda\" must be > 0, but was: %1%.", l, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool verify_exp_x(const char* function, RealType x, RealType* presult, const Policy& pol) +{ + if((x < 0) || (boost::math::isnan)(x)) + { + *presult = policies::raise_domain_error<RealType>( + function, + "The random variable must be >= 0, but was: %1%.", x, pol); + return false; + } + return true; +} + +} // namespace detail + +template <class RealType = double, class Policy = policies::policy<> > +class exponential_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + exponential_distribution(RealType l_lambda = 1) + : m_lambda(l_lambda) + { + RealType err; + detail::verify_lambda("boost::math::exponential_distribution<%1%>::exponential_distribution", l_lambda, &err, Policy()); + } // exponential_distribution + + RealType lambda()const { return m_lambda; } + +private: + RealType m_lambda; +}; + +typedef exponential_distribution<double> exponential; + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const exponential_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x. + if (std::numeric_limits<RealType>::has_infinity) + { + return std::pair<RealType, RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()); // 0 to + infinity. + } + else + { + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + max + } +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const exponential_distribution<RealType, Policy>& /*dist*/) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + using boost::math::tools::min_value; + return std::pair<RealType, RealType>(min_value<RealType>(), max_value<RealType>()); + // min_value<RealType>() to avoid a discontinuity at x = 0. +} + +template <class RealType, class Policy> +inline RealType pdf(const exponential_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::pdf(const exponential_distribution<%1%>&, %1%)"; + + RealType lambda = dist.lambda(); + RealType result = 0; + if(0 == detail::verify_lambda(function, lambda, &result, Policy())) + return result; + if(0 == detail::verify_exp_x(function, x, &result, Policy())) + return result; + // Workaround for VC11/12 bug: + if ((boost::math::isinf)(x)) + return 0; + result = lambda * exp(-lambda * x); + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const exponential_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const exponential_distribution<%1%>&, %1%)"; + + RealType result = 0; + RealType lambda = dist.lambda(); + if(0 == detail::verify_lambda(function, lambda, &result, Policy())) + return result; + if(0 == detail::verify_exp_x(function, x, &result, Policy())) + return result; + result = -boost::math::expm1(-x * lambda, Policy()); + + return result; +} // cdf + +template <class RealType, class Policy> +inline RealType quantile(const exponential_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const exponential_distribution<%1%>&, %1%)"; + + RealType result = 0; + RealType lambda = dist.lambda(); + if(0 == detail::verify_lambda(function, lambda, &result, Policy())) + return result; + if(0 == detail::check_probability(function, p, &result, Policy())) + return result; + + if(p == 0) + return 0; + if(p == 1) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = -boost::math::log1p(-p, Policy()) / lambda; + return result; +} // quantile + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<exponential_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const exponential_distribution<%1%>&, %1%)"; + + RealType result = 0; + RealType lambda = c.dist.lambda(); + if(0 == detail::verify_lambda(function, lambda, &result, Policy())) + return result; + if(0 == detail::verify_exp_x(function, c.param, &result, Policy())) + return result; + // Workaround for VC11/12 bug: + if (c.param >= tools::max_value<RealType>()) + return 0; + result = exp(-c.param * lambda); + + return result; +} + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<exponential_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const exponential_distribution<%1%>&, %1%)"; + + RealType result = 0; + RealType lambda = c.dist.lambda(); + if(0 == detail::verify_lambda(function, lambda, &result, Policy())) + return result; + + RealType q = c.param; + if(0 == detail::check_probability(function, q, &result, Policy())) + return result; + + if(q == 1) + return 0; + if(q == 0) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = -log(q) / lambda; + return result; +} + +template <class RealType, class Policy> +inline RealType mean(const exponential_distribution<RealType, Policy>& dist) +{ + RealType result = 0; + RealType lambda = dist.lambda(); + if(0 == detail::verify_lambda("boost::math::mean(const exponential_distribution<%1%>&)", lambda, &result, Policy())) + return result; + return 1 / lambda; +} + +template <class RealType, class Policy> +inline RealType standard_deviation(const exponential_distribution<RealType, Policy>& dist) +{ + RealType result = 0; + RealType lambda = dist.lambda(); + if(0 == detail::verify_lambda("boost::math::standard_deviation(const exponential_distribution<%1%>&)", lambda, &result, Policy())) + return result; + return 1 / lambda; +} + +template <class RealType, class Policy> +inline RealType mode(const exponential_distribution<RealType, Policy>& /*dist*/) +{ + return 0; +} + +template <class RealType, class Policy> +inline RealType median(const exponential_distribution<RealType, Policy>& dist) +{ + using boost::math::constants::ln_two; + return ln_two<RealType>() / dist.lambda(); // ln(2) / lambda +} + +template <class RealType, class Policy> +inline RealType skewness(const exponential_distribution<RealType, Policy>& /*dist*/) +{ + return 2; +} + +template <class RealType, class Policy> +inline RealType kurtosis(const exponential_distribution<RealType, Policy>& /*dist*/) +{ + return 9; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const exponential_distribution<RealType, Policy>& /*dist*/) +{ + return 6; +} + +} // namespace math +} // namespace boost + +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_EXPONENTIAL_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/extreme_value.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/extreme_value.hpp new file mode 100644 index 00000000000..ef9fbe817da --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/extreme_value.hpp @@ -0,0 +1,297 @@ +// Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_EXTREME_VALUE_HPP +#define BOOST_STATS_EXTREME_VALUE_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/special_functions/log1p.hpp> +#include <boost/math/special_functions/expm1.hpp> +#include <boost/math/distributions/complement.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/config/no_tr1/cmath.hpp> + +// +// This is the maximum extreme value distribution, see +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366g.htm +// and http://mathworld.wolfram.com/ExtremeValueDistribution.html +// Also known as a Fisher-Tippett distribution, a log-Weibull +// distribution or a Gumbel distribution. + +#include <utility> + +#ifdef BOOST_MSVC +# pragma warning(push) +# pragma warning(disable: 4702) // unreachable code (return after domain_error throw). +#endif + +namespace boost{ namespace math{ + +namespace detail{ +// +// Error check: +// +template <class RealType, class Policy> +inline bool verify_scale_b(const char* function, RealType b, RealType* presult, const Policy& pol) +{ + if((b <= 0) || !(boost::math::isfinite)(b)) + { + *presult = policies::raise_domain_error<RealType>( + function, + "The scale parameter \"b\" must be finite and > 0, but was: %1%.", b, pol); + return false; + } + return true; +} + +} // namespace detail + +template <class RealType = double, class Policy = policies::policy<> > +class extreme_value_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + extreme_value_distribution(RealType a = 0, RealType b = 1) + : m_a(a), m_b(b) + { + RealType err; + detail::verify_scale_b("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", b, &err, Policy()); + detail::check_finite("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", a, &err, Policy()); + } // extreme_value_distribution + + RealType location()const { return m_a; } + RealType scale()const { return m_b; } + +private: + RealType m_a, m_b; +}; + +typedef extreme_value_distribution<double> extreme_value; + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const extreme_value_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>( + std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(), + std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>()); +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const extreme_value_distribution<RealType, Policy>& /*dist*/) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline RealType pdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::pdf(const extreme_value_distribution<%1%>&, %1%)"; + + RealType a = dist.location(); + RealType b = dist.scale(); + RealType result = 0; + if((boost::math::isinf)(x)) + return 0.0f; + if(0 == detail::verify_scale_b(function, b, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if(0 == detail::check_x(function, x, &result, Policy())) + return result; + result = exp((a-x)/b) * exp(-exp((a-x)/b)) / b; + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)"; + + if((boost::math::isinf)(x)) + return x < 0 ? 0.0f : 1.0f; + RealType a = dist.location(); + RealType b = dist.scale(); + RealType result = 0; + if(0 == detail::verify_scale_b(function, b, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if(0 == detail::check_x("boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy())) + return result; + + result = exp(-exp((a-x)/b)); + + return result; +} // cdf + +template <class RealType, class Policy> +RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)"; + + RealType a = dist.location(); + RealType b = dist.scale(); + RealType result = 0; + if(0 == detail::verify_scale_b(function, b, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if(0 == detail::check_probability(function, p, &result, Policy())) + return result; + + if(p == 0) + return -policies::raise_overflow_error<RealType>(function, 0, Policy()); + if(p == 1) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = a - log(-log(p)) * b; + + return result; +} // quantile + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)"; + + if((boost::math::isinf)(c.param)) + return c.param < 0 ? 1.0f : 0.0f; + RealType a = c.dist.location(); + RealType b = c.dist.scale(); + RealType result = 0; + if(0 == detail::verify_scale_b(function, b, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if(0 == detail::check_x(function, c.param, &result, Policy())) + return result; + + result = -boost::math::expm1(-exp((a-c.param)/b), Policy()); + + return result; +} + +template <class RealType, class Policy> +RealType quantile(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)"; + + RealType a = c.dist.location(); + RealType b = c.dist.scale(); + RealType q = c.param; + RealType result = 0; + if(0 == detail::verify_scale_b(function, b, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if(0 == detail::check_probability(function, q, &result, Policy())) + return result; + + if(q == 0) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + if(q == 1) + return -policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = a - log(-boost::math::log1p(-q, Policy())) * b; + + return result; +} + +template <class RealType, class Policy> +inline RealType mean(const extreme_value_distribution<RealType, Policy>& dist) +{ + RealType a = dist.location(); + RealType b = dist.scale(); + RealType result = 0; + if(0 == detail::verify_scale_b("boost::math::mean(const extreme_value_distribution<%1%>&)", b, &result, Policy())) + return result; + if(0 == detail::check_scale("boost::math::mean(const extreme_value_distribution<%1%>&)", a, &result, Policy())) + return result; + return a + constants::euler<RealType>() * b; +} + +template <class RealType, class Policy> +inline RealType standard_deviation(const extreme_value_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions. + + RealType b = dist.scale(); + RealType result = 0; + if(0 == detail::verify_scale_b("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", b, &result, Policy())) + return result; + if(0 == detail::check_scale("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", dist.location(), &result, Policy())) + return result; + return constants::pi<RealType>() * b / sqrt(static_cast<RealType>(6)); +} + +template <class RealType, class Policy> +inline RealType mode(const extreme_value_distribution<RealType, Policy>& dist) +{ + return dist.location(); +} + +template <class RealType, class Policy> +inline RealType median(const extreme_value_distribution<RealType, Policy>& dist) +{ + using constants::ln_ln_two; + return dist.location() - dist.scale() * ln_ln_two<RealType>(); +} + +template <class RealType, class Policy> +inline RealType skewness(const extreme_value_distribution<RealType, Policy>& /*dist*/) +{ + // + // This is 12 * sqrt(6) * zeta(3) / pi^3: + // See http://mathworld.wolfram.com/ExtremeValueDistribution.html + // + return static_cast<RealType>(1.1395470994046486574927930193898461120875997958366L); +} + +template <class RealType, class Policy> +inline RealType kurtosis(const extreme_value_distribution<RealType, Policy>& /*dist*/) +{ + // See http://mathworld.wolfram.com/ExtremeValueDistribution.html + return RealType(27) / 5; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const extreme_value_distribution<RealType, Policy>& /*dist*/) +{ + // See http://mathworld.wolfram.com/ExtremeValueDistribution.html + return RealType(12) / 5; +} + + +} // namespace math +} // namespace boost + +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_EXTREME_VALUE_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/fisher_f.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/fisher_f.hpp new file mode 100644 index 00000000000..9e259bcc96c --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/fisher_f.hpp @@ -0,0 +1,387 @@ +// Copyright John Maddock 2006. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP +#define BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/beta.hpp> // for incomplete beta. +#include <boost/math/distributions/complement.hpp> // complements +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks +#include <boost/math/special_functions/fpclassify.hpp> + +#include <utility> + +namespace boost{ namespace math{ + +template <class RealType = double, class Policy = policies::policy<> > +class fisher_f_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + fisher_f_distribution(const RealType& i, const RealType& j) : m_df1(i), m_df2(j) + { + static const char* function = "fisher_f_distribution<%1%>::fisher_f_distribution"; + RealType result; + detail::check_df( + function, m_df1, &result, Policy()); + detail::check_df( + function, m_df2, &result, Policy()); + } // fisher_f_distribution + + RealType degrees_of_freedom1()const + { + return m_df1; + } + RealType degrees_of_freedom2()const + { + return m_df2; + } + +private: + // + // Data members: + // + RealType m_df1; // degrees of freedom are a real number. + RealType m_df2; // degrees of freedom are a real number. +}; + +typedef fisher_f_distribution<double> fisher_f; + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const fisher_f_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const fisher_f_distribution<RealType, Policy>& /*dist*/) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); +} + +template <class RealType, class Policy> +RealType pdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + static const char* function = "boost::math::pdf(fisher_f_distribution<%1%> const&, %1%)"; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + + if((x < 0) || !(boost::math::isfinite)(x)) + { + return policies::raise_domain_error<RealType>( + function, "Random variable parameter was %1%, but must be > 0 !", x, Policy()); + } + + if(x == 0) + { + // special cases: + if(df1 < 2) + return policies::raise_overflow_error<RealType>( + function, 0, Policy()); + else if(df1 == 2) + return 1; + else + return 0; + } + + // + // You reach this formula by direct differentiation of the + // cdf expressed in terms of the incomplete beta. + // + // There are two versions so we don't pass a value of z + // that is very close to 1 to ibeta_derivative: for some values + // of df1 and df2, all the change takes place in this area. + // + RealType v1x = df1 * x; + RealType result; + if(v1x > df2) + { + result = (df2 * df1) / ((df2 + v1x) * (df2 + v1x)); + result *= ibeta_derivative(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy()); + } + else + { + result = df2 + df1 * x; + result = (result * df1 - x * df1 * df1) / (result * result); + result *= ibeta_derivative(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy()); + } + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x) +{ + static const char* function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)"; + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + + if((x < 0) || !(boost::math::isfinite)(x)) + { + return policies::raise_domain_error<RealType>( + function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy()); + } + + RealType v1x = df1 * x; + // + // There are two equivalent formulas used here, the aim is + // to prevent the final argument to the incomplete beta + // from being too close to 1: for some values of df1 and df2 + // the rate of change can be arbitrarily large in this area, + // whilst the value we're passing will have lost information + // content as a result of being 0.999999something. Better + // to switch things around so we're passing 1-z instead. + // + return v1x > df2 + ? boost::math::ibetac(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy()) + : boost::math::ibeta(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy()); +} // cdf + +template <class RealType, class Policy> +inline RealType quantile(const fisher_f_distribution<RealType, Policy>& dist, const RealType& p) +{ + static const char* function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)"; + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + if(false == (detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy()) + && detail::check_probability( + function, p, &error_result, Policy()))) + return error_result; + + // With optimizations turned on, gcc wrongly warns about y being used + // uninitializated unless we initialize it to something: + RealType x, y(0); + + x = boost::math::ibeta_inv(df1 / 2, df2 / 2, p, &y, Policy()); + + return df2 * x / (df1 * y); +} // quantile + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c) +{ + static const char* function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)"; + RealType df1 = c.dist.degrees_of_freedom1(); + RealType df2 = c.dist.degrees_of_freedom2(); + RealType x = c.param; + // Error check: + RealType error_result = 0; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + + if((x < 0) || !(boost::math::isfinite)(x)) + { + return policies::raise_domain_error<RealType>( + function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy()); + } + + RealType v1x = df1 * x; + // + // There are two equivalent formulas used here, the aim is + // to prevent the final argument to the incomplete beta + // from being too close to 1: for some values of df1 and df2 + // the rate of change can be arbitrarily large in this area, + // whilst the value we're passing will have lost information + // content as a result of being 0.999999something. Better + // to switch things around so we're passing 1-z instead. + // + return v1x > df2 + ? boost::math::ibeta(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy()) + : boost::math::ibetac(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy()); +} + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c) +{ + static const char* function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)"; + RealType df1 = c.dist.degrees_of_freedom1(); + RealType df2 = c.dist.degrees_of_freedom2(); + RealType p = c.param; + // Error check: + RealType error_result = 0; + if(false == (detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy()) + && detail::check_probability( + function, p, &error_result, Policy()))) + return error_result; + + RealType x, y; + + x = boost::math::ibetac_inv(df1 / 2, df2 / 2, p, &y, Policy()); + + return df2 * x / (df1 * y); +} + +template <class RealType, class Policy> +inline RealType mean(const fisher_f_distribution<RealType, Policy>& dist) +{ // Mean of F distribution = v. + static const char* function = "boost::math::mean(fisher_f_distribution<%1%> const&)"; + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + if(df2 <= 2) + { + return policies::raise_domain_error<RealType>( + function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mean.", df2, Policy()); + } + return df2 / (df2 - 2); +} // mean + +template <class RealType, class Policy> +inline RealType variance(const fisher_f_distribution<RealType, Policy>& dist) +{ // Variance of F distribution. + static const char* function = "boost::math::variance(fisher_f_distribution<%1%> const&)"; + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + if(df2 <= 4) + { + return policies::raise_domain_error<RealType>( + function, "Second degree of freedom was %1% but must be > 4 in order for the distribution to have a valid variance.", df2, Policy()); + } + return 2 * df2 * df2 * (df1 + df2 - 2) / (df1 * (df2 - 2) * (df2 - 2) * (df2 - 4)); +} // variance + +template <class RealType, class Policy> +inline RealType mode(const fisher_f_distribution<RealType, Policy>& dist) +{ + static const char* function = "boost::math::mode(fisher_f_distribution<%1%> const&)"; + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + if(df2 <= 2) + { + return policies::raise_domain_error<RealType>( + function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mode.", df2, Policy()); + } + return df2 * (df1 - 2) / (df1 * (df2 + 2)); +} + +//template <class RealType, class Policy> +//inline RealType median(const fisher_f_distribution<RealType, Policy>& dist) +//{ // Median of Fisher F distribution is not defined. +// return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); +// } // median + +// Now implemented via quantile(half) in derived accessors. + +template <class RealType, class Policy> +inline RealType skewness(const fisher_f_distribution<RealType, Policy>& dist) +{ + static const char* function = "boost::math::skewness(fisher_f_distribution<%1%> const&)"; + BOOST_MATH_STD_USING // ADL of std names + // See http://mathworld.wolfram.com/F-Distribution.html + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + if(df2 <= 6) + { + return policies::raise_domain_error<RealType>( + function, "Second degree of freedom was %1% but must be > 6 in order for the distribution to have a skewness.", df2, Policy()); + } + return 2 * (df2 + 2 * df1 - 2) * sqrt((2 * df2 - 8) / (df1 * (df2 + df1 - 2))) / (df2 - 6); +} + +template <class RealType, class Policy> +RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist); + +template <class RealType, class Policy> +inline RealType kurtosis(const fisher_f_distribution<RealType, Policy>& dist) +{ + return 3 + kurtosis_excess(dist); +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist) +{ + static const char* function = "boost::math::kurtosis_excess(fisher_f_distribution<%1%> const&)"; + // See http://mathworld.wolfram.com/F-Distribution.html + RealType df1 = dist.degrees_of_freedom1(); + RealType df2 = dist.degrees_of_freedom2(); + // Error check: + RealType error_result = 0; + if(false == detail::check_df( + function, df1, &error_result, Policy()) + && detail::check_df( + function, df2, &error_result, Policy())) + return error_result; + if(df2 <= 8) + { + return policies::raise_domain_error<RealType>( + function, "Second degree of freedom was %1% but must be > 8 in order for the distribution to have a kutosis.", df2, Policy()); + } + RealType df2_2 = df2 * df2; + RealType df1_2 = df1 * df1; + RealType n = -16 + 20 * df2 - 8 * df2_2 + df2_2 * df2 + 44 * df1 - 32 * df2 * df1 + 5 * df2_2 * df1 - 22 * df1_2 + 5 * df2 * df1_2; + n *= 12; + RealType d = df1 * (df2 - 6) * (df2 - 8) * (df1 + df2 - 2); + return n / d; +} + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/fwd.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/fwd.hpp new file mode 100644 index 00000000000..609dc3b563e --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/fwd.hpp @@ -0,0 +1,146 @@ +// fwd.hpp Forward declarations of Boost.Math distributions. + +// Copyright Paul A. Bristow 2007, 2010, 2012. +// Copyright John Maddock 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DISTRIBUTIONS_FWD_HPP +#define BOOST_MATH_DISTRIBUTIONS_FWD_HPP + +// 31 distributions at Boost 1.52 + +namespace boost{ namespace math{ + +template <class RealType, class Policy> +class bernoulli_distribution; + +template <class RealType, class Policy> +class beta_distribution; + +template <class RealType, class Policy> +class binomial_distribution; + +template <class RealType, class Policy> +class cauchy_distribution; + +template <class RealType, class Policy> +class chi_squared_distribution; + +template <class RealType, class Policy> +class exponential_distribution; + +template <class RealType, class Policy> +class extreme_value_distribution; + +template <class RealType, class Policy> +class fisher_f_distribution; + +template <class RealType, class Policy> +class gamma_distribution; + +template <class RealType, class Policy> +class geometric_distribution; + +template <class RealType, class Policy> +class hypergeometric_distribution; + +template <class RealType, class Policy> +class inverse_chi_squared_distribution; + +template <class RealType, class Policy> +class inverse_gamma_distribution; + +template <class RealType, class Policy> +class inverse_gaussian_distribution; + +template <class RealType, class Policy> +class laplace_distribution; + +template <class RealType, class Policy> +class logistic_distribution; + +template <class RealType, class Policy> +class lognormal_distribution; + +template <class RealType, class Policy> +class negative_binomial_distribution; + +template <class RealType, class Policy> +class non_central_beta_distribution; + +template <class RealType, class Policy> +class non_central_chi_squared_distribution; + +template <class RealType, class Policy> +class non_central_f_distribution; + +template <class RealType, class Policy> +class non_central_t_distribution; + +template <class RealType, class Policy> +class normal_distribution; + +template <class RealType, class Policy> +class pareto_distribution; + +template <class RealType, class Policy> +class poisson_distribution; + +template <class RealType, class Policy> +class rayleigh_distribution; + +template <class RealType, class Policy> +class skew_normal_distribution; + +template <class RealType, class Policy> +class students_t_distribution; + +template <class RealType, class Policy> +class triangular_distribution; + +template <class RealType, class Policy> +class uniform_distribution; + +template <class RealType, class Policy> +class weibull_distribution; + +}} // namespaces + +#define BOOST_MATH_DECLARE_DISTRIBUTIONS(Type, Policy)\ + typedef boost::math::bernoulli_distribution<Type, Policy> bernoulli;\ + typedef boost::math::beta_distribution<Type, Policy> beta;\ + typedef boost::math::binomial_distribution<Type, Policy> binomial;\ + typedef boost::math::cauchy_distribution<Type, Policy> cauchy;\ + typedef boost::math::chi_squared_distribution<Type, Policy> chi_squared;\ + typedef boost::math::exponential_distribution<Type, Policy> exponential;\ + typedef boost::math::extreme_value_distribution<Type, Policy> extreme_value;\ + typedef boost::math::fisher_f_distribution<Type, Policy> fisher_f;\ + typedef boost::math::gamma_distribution<Type, Policy> gamma;\ + typedef boost::math::geometric_distribution<Type, Policy> geometric;\ + typedef boost::math::hypergeometric_distribution<Type, Policy> hypergeometric;\ + typedef boost::math::inverse_chi_squared_distribution<Type, Policy> inverse_chi_squared;\ + typedef boost::math::inverse_gaussian_distribution<Type, Policy> inverse_gaussian;\ + typedef boost::math::inverse_gamma_distribution<Type, Policy> inverse_gamma;\ + typedef boost::math::laplace_distribution<Type, Policy> laplace;\ + typedef boost::math::logistic_distribution<Type, Policy> logistic;\ + typedef boost::math::lognormal_distribution<Type, Policy> lognormal;\ + typedef boost::math::negative_binomial_distribution<Type, Policy> negative_binomial;\ + typedef boost::math::non_central_beta_distribution<Type, Policy> non_central_beta;\ + typedef boost::math::non_central_chi_squared_distribution<Type, Policy> non_central_chi_squared;\ + typedef boost::math::non_central_f_distribution<Type, Policy> non_central_f;\ + typedef boost::math::non_central_t_distribution<Type, Policy> non_central_t;\ + typedef boost::math::normal_distribution<Type, Policy> normal;\ + typedef boost::math::pareto_distribution<Type, Policy> pareto;\ + typedef boost::math::poisson_distribution<Type, Policy> poisson;\ + typedef boost::math::rayleigh_distribution<Type, Policy> rayleigh;\ + typedef boost::math::skew_normal_distribution<Type, Policy> skew_normal;\ + typedef boost::math::students_t_distribution<Type, Policy> students_t;\ + typedef boost::math::triangular_distribution<Type, Policy> triangular;\ + typedef boost::math::uniform_distribution<Type, Policy> uniform;\ + typedef boost::math::weibull_distribution<Type, Policy> weibull; + +#endif // BOOST_MATH_DISTRIBUTIONS_FWD_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/gamma.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/gamma.hpp new file mode 100644 index 00000000000..9a9e2a4f524 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/gamma.hpp @@ -0,0 +1,349 @@ +// Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_GAMMA_HPP +#define BOOST_STATS_GAMMA_HPP + +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm +// http://mathworld.wolfram.com/GammaDistribution.html +// http://en.wikipedia.org/wiki/Gamma_distribution + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/gamma.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/distributions/complement.hpp> + +#include <utility> + +namespace boost{ namespace math +{ +namespace detail +{ + +template <class RealType, class Policy> +inline bool check_gamma_shape( + const char* function, + RealType shape, + RealType* result, const Policy& pol) +{ + if((shape <= 0) || !(boost::math::isfinite)(shape)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Shape parameter is %1%, but must be > 0 !", shape, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_gamma_x( + const char* function, + RealType const& x, + RealType* result, const Policy& pol) +{ + if((x < 0) || !(boost::math::isfinite)(x)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Random variate is %1% but must be >= 0 !", x, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_gamma( + const char* function, + RealType scale, + RealType shape, + RealType* result, const Policy& pol) +{ + return check_scale(function, scale, result, pol) && check_gamma_shape(function, shape, result, pol); +} + +} // namespace detail + +template <class RealType = double, class Policy = policies::policy<> > +class gamma_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + gamma_distribution(RealType l_shape, RealType l_scale = 1) + : m_shape(l_shape), m_scale(l_scale) + { + RealType result; + detail::check_gamma("boost::math::gamma_distribution<%1%>::gamma_distribution", l_scale, l_shape, &result, Policy()); + } + + RealType shape()const + { + return m_shape; + } + + RealType scale()const + { + return m_scale; + } +private: + // + // Data members: + // + RealType m_shape; // distribution shape + RealType m_scale; // distribution scale +}; + +// NO typedef because of clash with name of gamma function. + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const gamma_distribution<RealType, Policy>& /* dist */) +{ // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const gamma_distribution<RealType, Policy>& /* dist */) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + using boost::math::tools::min_value; + return std::pair<RealType, RealType>(min_value<RealType>(), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline RealType pdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::pdf(const gamma_distribution<%1%>&, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_gamma_x(function, x, &result, Policy())) + return result; + + if(x == 0) + { + return 0; + } + result = gamma_p_derivative(shape, x / scale, Policy()) / scale; + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const gamma_distribution<%1%>&, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_gamma_x(function, x, &result, Policy())) + return result; + + result = boost::math::gamma_p(shape, x / scale, Policy()); + return result; +} // cdf + +template <class RealType, class Policy> +inline RealType quantile(const gamma_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_probability(function, p, &result, Policy())) + return result; + + if(p == 1) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = gamma_p_inv(shape, p, Policy()) * scale; + + return result; +} + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)"; + + RealType shape = c.dist.shape(); + RealType scale = c.dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_gamma_x(function, c.param, &result, Policy())) + return result; + + result = gamma_q(shape, c.param / scale, Policy()); + + return result; +} + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)"; + + RealType shape = c.dist.shape(); + RealType scale = c.dist.scale(); + RealType q = c.param; + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_probability(function, q, &result, Policy())) + return result; + + if(q == 0) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = gamma_q_inv(shape, q, Policy()) * scale; + + return result; +} + +template <class RealType, class Policy> +inline RealType mean(const gamma_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::mean(const gamma_distribution<%1%>&)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + + result = shape * scale; + return result; +} + +template <class RealType, class Policy> +inline RealType variance(const gamma_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::variance(const gamma_distribution<%1%>&)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + + result = shape * scale * scale; + return result; +} + +template <class RealType, class Policy> +inline RealType mode(const gamma_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::mode(const gamma_distribution<%1%>&)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + + if(shape < 1) + return policies::raise_domain_error<RealType>( + function, + "The mode of the gamma distribution is only defined for values of the shape parameter >= 1, but got %1%.", + shape, Policy()); + + result = (shape - 1) * scale; + return result; +} + +//template <class RealType, class Policy> +//inline RealType median(const gamma_distribution<RealType, Policy>& dist) +//{ // Rely on default definition in derived accessors. +//} + +template <class RealType, class Policy> +inline RealType skewness(const gamma_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::skewness(const gamma_distribution<%1%>&)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + + result = 2 / sqrt(shape); + return result; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const gamma_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::kurtosis_excess(const gamma_distribution<%1%>&)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_gamma(function, scale, shape, &result, Policy())) + return result; + + result = 6 / shape; + return result; +} + +template <class RealType, class Policy> +inline RealType kurtosis(const gamma_distribution<RealType, Policy>& dist) +{ + return kurtosis_excess(dist) + 3; +} + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_GAMMA_HPP + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/geometric.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/geometric.hpp new file mode 100644 index 00000000000..88947d6c57a --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/geometric.hpp @@ -0,0 +1,516 @@ +// boost\math\distributions\geometric.hpp + +// Copyright John Maddock 2010. +// Copyright Paul A. Bristow 2010. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// geometric distribution is a discrete probability distribution. +// It expresses the probability distribution of the number (k) of +// events, occurrences, failures or arrivals before the first success. +// supported on the set {0, 1, 2, 3...} + +// Note that the set includes zero (unlike some definitions that start at one). + +// The random variate k is the number of events, occurrences or arrivals. +// k argument may be integral, signed, or unsigned, or floating point. +// If necessary, it has already been promoted from an integral type. + +// Note that the geometric distribution +// (like others including the binomial, geometric & Bernoulli) +// is strictly defined as a discrete function: +// only integral values of k are envisaged. +// However because the method of calculation uses a continuous gamma function, +// it is convenient to treat it as if a continous function, +// and permit non-integral values of k. +// To enforce the strict mathematical model, users should use floor or ceil functions +// on k outside this function to ensure that k is integral. + +// See http://en.wikipedia.org/wiki/geometric_distribution +// http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/Statistics/DiscreteDistributions.html +// http://mathworld.wolfram.com/GeometricDistribution.html + +#ifndef BOOST_MATH_SPECIAL_GEOMETRIC_HPP +#define BOOST_MATH_SPECIAL_GEOMETRIC_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x) == Ix(a, b). +#include <boost/math/distributions/complement.hpp> // complement. +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks domain_error & logic_error. +#include <boost/math/special_functions/fpclassify.hpp> // isnan. +#include <boost/math/tools/roots.hpp> // for root finding. +#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> + +#include <boost/type_traits/is_floating_point.hpp> +#include <boost/type_traits/is_integral.hpp> +#include <boost/type_traits/is_same.hpp> +#include <boost/mpl/if.hpp> + +#include <limits> // using std::numeric_limits; +#include <utility> + +#if defined (BOOST_MSVC) +# pragma warning(push) +// This believed not now necessary, so commented out. +//# pragma warning(disable: 4702) // unreachable code. +// in domain_error_imp in error_handling. +#endif + +namespace boost +{ + namespace math + { + namespace geometric_detail + { + // Common error checking routines for geometric distribution function: + template <class RealType, class Policy> + inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol) + { + if( !(boost::math::isfinite)(p) || (p < 0) || (p > 1) ) + { + *result = policies::raise_domain_error<RealType>( + function, + "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol); + return false; + } + return true; + } + + template <class RealType, class Policy> + inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& pol) + { + return check_success_fraction(function, p, result, pol); + } + + template <class RealType, class Policy> + inline bool check_dist_and_k(const char* function, const RealType& p, RealType k, RealType* result, const Policy& pol) + { + if(check_dist(function, p, result, pol) == false) + { + return false; + } + if( !(boost::math::isfinite)(k) || (k < 0) ) + { // Check k failures. + *result = policies::raise_domain_error<RealType>( + function, + "Number of failures argument is %1%, but must be >= 0 !", k, pol); + return false; + } + return true; + } // Check_dist_and_k + + template <class RealType, class Policy> + inline bool check_dist_and_prob(const char* function, RealType p, RealType prob, RealType* result, const Policy& pol) + { + if(check_dist(function, p, result, pol) && detail::check_probability(function, prob, result, pol) == false) + { + return false; + } + return true; + } // check_dist_and_prob + } // namespace geometric_detail + + template <class RealType = double, class Policy = policies::policy<> > + class geometric_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + geometric_distribution(RealType p) : m_p(p) + { // Constructor stores success_fraction p. + RealType result; + geometric_detail::check_dist( + "geometric_distribution<%1%>::geometric_distribution", + m_p, // Check success_fraction 0 <= p <= 1. + &result, Policy()); + } // geometric_distribution constructor. + + // Private data getter class member functions. + RealType success_fraction() const + { // Probability of success as fraction in range 0 to 1. + return m_p; + } + RealType successes() const + { // Total number of successes r = 1 (for compatibility with negative binomial?). + return 1; + } + + // Parameter estimation. + // (These are copies of negative_binomial distribution with successes = 1). + static RealType find_lower_bound_on_p( + RealType trials, + RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test. + { + static const char* function = "boost::math::geometric<%1%>::find_lower_bound_on_p"; + RealType result = 0; // of error checks. + RealType successes = 1; + RealType failures = trials - successes; + if(false == detail::check_probability(function, alpha, &result, Policy()) + && geometric_detail::check_dist_and_k( + function, RealType(0), failures, &result, Policy())) + { + return result; + } + // Use complement ibeta_inv function for lower bound. + // This is adapted from the corresponding binomial formula + // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm + // This is a Clopper-Pearson interval, and may be overly conservative, + // see also "A Simple Improved Inferential Method for Some + // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY + // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf + // + return ibeta_inv(successes, failures + 1, alpha, static_cast<RealType*>(0), Policy()); + } // find_lower_bound_on_p + + static RealType find_upper_bound_on_p( + RealType trials, + RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test. + { + static const char* function = "boost::math::geometric<%1%>::find_upper_bound_on_p"; + RealType result = 0; // of error checks. + RealType successes = 1; + RealType failures = trials - successes; + if(false == geometric_detail::check_dist_and_k( + function, RealType(0), failures, &result, Policy()) + && detail::check_probability(function, alpha, &result, Policy())) + { + return result; + } + if(failures == 0) + { + return 1; + }// Use complement ibetac_inv function for upper bound. + // Note adjusted failures value: *not* failures+1 as usual. + // This is adapted from the corresponding binomial formula + // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm + // This is a Clopper-Pearson interval, and may be overly conservative, + // see also "A Simple Improved Inferential Method for Some + // Discrete Distributions" Yong CAI and K. Krishnamoorthy + // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf + // + return ibetac_inv(successes, failures, alpha, static_cast<RealType*>(0), Policy()); + } // find_upper_bound_on_p + + // Estimate number of trials : + // "How many trials do I need to be P% sure of seeing k or fewer failures?" + + static RealType find_minimum_number_of_trials( + RealType k, // number of failures (k >= 0). + RealType p, // success fraction 0 <= p <= 1. + RealType alpha) // risk level threshold 0 <= alpha <= 1. + { + static const char* function = "boost::math::geometric<%1%>::find_minimum_number_of_trials"; + // Error checks: + RealType result = 0; + if(false == geometric_detail::check_dist_and_k( + function, p, k, &result, Policy()) + && detail::check_probability(function, alpha, &result, Policy())) + { + return result; + } + result = ibeta_inva(k + 1, p, alpha, Policy()); // returns n - k + return result + k; + } // RealType find_number_of_failures + + static RealType find_maximum_number_of_trials( + RealType k, // number of failures (k >= 0). + RealType p, // success fraction 0 <= p <= 1. + RealType alpha) // risk level threshold 0 <= alpha <= 1. + { + static const char* function = "boost::math::geometric<%1%>::find_maximum_number_of_trials"; + // Error checks: + RealType result = 0; + if(false == geometric_detail::check_dist_and_k( + function, p, k, &result, Policy()) + && detail::check_probability(function, alpha, &result, Policy())) + { + return result; + } + result = ibetac_inva(k + 1, p, alpha, Policy()); // returns n - k + return result + k; + } // RealType find_number_of_trials complemented + + private: + //RealType m_r; // successes fixed at unity. + RealType m_p; // success_fraction + }; // template <class RealType, class Policy> class geometric_distribution + + typedef geometric_distribution<double> geometric; // Reserved name of type double. + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const geometric_distribution<RealType, Policy>& /* dist */) + { // Range of permissible values for random variable k. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer? + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const geometric_distribution<RealType, Policy>& /* dist */) + { // Range of supported values for random variable k. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer? + } + + template <class RealType, class Policy> + inline RealType mean(const geometric_distribution<RealType, Policy>& dist) + { // Mean of geometric distribution = (1-p)/p. + return (1 - dist.success_fraction() ) / dist.success_fraction(); + } // mean + + // median implemented via quantile(half) in derived accessors. + + template <class RealType, class Policy> + inline RealType mode(const geometric_distribution<RealType, Policy>&) + { // Mode of geometric distribution = zero. + BOOST_MATH_STD_USING // ADL of std functions. + return 0; + } // mode + + template <class RealType, class Policy> + inline RealType variance(const geometric_distribution<RealType, Policy>& dist) + { // Variance of Binomial distribution = (1-p) / p^2. + return (1 - dist.success_fraction()) + / (dist.success_fraction() * dist.success_fraction()); + } // variance + + template <class RealType, class Policy> + inline RealType skewness(const geometric_distribution<RealType, Policy>& dist) + { // skewness of geometric distribution = 2-p / (sqrt(r(1-p)) + BOOST_MATH_STD_USING // ADL of std functions. + RealType p = dist.success_fraction(); + return (2 - p) / sqrt(1 - p); + } // skewness + + template <class RealType, class Policy> + inline RealType kurtosis(const geometric_distribution<RealType, Policy>& dist) + { // kurtosis of geometric distribution + // http://en.wikipedia.org/wiki/geometric is kurtosis_excess so add 3 + RealType p = dist.success_fraction(); + return 3 + (p*p - 6*p + 6) / (1 - p); + } // kurtosis + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const geometric_distribution<RealType, Policy>& dist) + { // kurtosis excess of geometric distribution + // http://mathworld.wolfram.com/Kurtosis.html table of kurtosis_excess + RealType p = dist.success_fraction(); + return (p*p - 6*p + 6) / (1 - p); + } // kurtosis_excess + + // RealType standard_deviation(const geometric_distribution<RealType, Policy>& dist) + // standard_deviation provided by derived accessors. + // RealType hazard(const geometric_distribution<RealType, Policy>& dist) + // hazard of geometric distribution provided by derived accessors. + // RealType chf(const geometric_distribution<RealType, Policy>& dist) + // chf of geometric distribution provided by derived accessors. + + template <class RealType, class Policy> + inline RealType pdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k) + { // Probability Density/Mass Function. + BOOST_FPU_EXCEPTION_GUARD + BOOST_MATH_STD_USING // For ADL of math functions. + static const char* function = "boost::math::pdf(const geometric_distribution<%1%>&, %1%)"; + + RealType p = dist.success_fraction(); + RealType result = 0; + if(false == geometric_detail::check_dist_and_k( + function, + p, + k, + &result, Policy())) + { + return result; + } + if (k == 0) + { + return p; // success_fraction + } + RealType q = 1 - p; // Inaccurate for small p? + // So try to avoid inaccuracy for large or small p. + // but has little effect > last significant bit. + //cout << "p * pow(q, k) " << result << endl; // seems best whatever p + //cout << "exp(p * k * log1p(-p)) " << p * exp(k * log1p(-p)) << endl; + //if (p < 0.5) + //{ + // result = p * pow(q, k); + //} + //else + //{ + // result = p * exp(k * log1p(-p)); + //} + result = p * pow(q, k); + return result; + } // geometric_pdf + + template <class RealType, class Policy> + inline RealType cdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k) + { // Cumulative Distribution Function of geometric. + static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)"; + + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + RealType p = dist.success_fraction(); + // Error check: + RealType result = 0; + if(false == geometric_detail::check_dist_and_k( + function, + p, + k, + &result, Policy())) + { + return result; + } + if(k == 0) + { + return p; // success_fraction + } + //RealType q = 1 - p; // Bad for small p + //RealType probability = 1 - std::pow(q, k+1); + + RealType z = boost::math::log1p(-p, Policy()) * (k + 1); + RealType probability = -boost::math::expm1(z, Policy()); + + return probability; + } // cdf Cumulative Distribution Function geometric. + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c) + { // Complemented Cumulative Distribution Function geometric. + BOOST_MATH_STD_USING + static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)"; + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + RealType const& k = c.param; + geometric_distribution<RealType, Policy> const& dist = c.dist; + RealType p = dist.success_fraction(); + // Error check: + RealType result = 0; + if(false == geometric_detail::check_dist_and_k( + function, + p, + k, + &result, Policy())) + { + return result; + } + RealType z = boost::math::log1p(-p, Policy()) * (k+1); + RealType probability = exp(z); + return probability; + } // cdf Complemented Cumulative Distribution Function geometric. + + template <class RealType, class Policy> + inline RealType quantile(const geometric_distribution<RealType, Policy>& dist, const RealType& x) + { // Quantile, percentile/100 or Percent Point geometric function. + // Return the number of expected failures k for a given probability p. + + // Inverse cumulative Distribution Function or Quantile (percentile / 100) of geometric Probability. + // k argument may be integral, signed, or unsigned, or floating point. + + static const char* function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)"; + BOOST_MATH_STD_USING // ADL of std functions. + + RealType success_fraction = dist.success_fraction(); + // Check dist and x. + RealType result = 0; + if(false == geometric_detail::check_dist_and_prob + (function, success_fraction, x, &result, Policy())) + { + return result; + } + + // Special cases. + if (x == 1) + { // Would need +infinity failures for total confidence. + result = policies::raise_overflow_error<RealType>( + function, + "Probability argument is 1, which implies infinite failures !", Policy()); + return result; + // usually means return +std::numeric_limits<RealType>::infinity(); + // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR + } + if (x == 0) + { // No failures are expected if P = 0. + return 0; // Total trials will be just dist.successes. + } + // if (P <= pow(dist.success_fraction(), 1)) + if (x <= success_fraction) + { // p <= pdf(dist, 0) == cdf(dist, 0) + return 0; + } + if (x == 1) + { + return 0; + } + + // log(1-x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small + result = boost::math::log1p(-x, Policy()) / boost::math::log1p(-success_fraction, Policy()) - 1; + // Subtract a few epsilons here too? + // to make sure it doesn't slip over, so ceil would be one too many. + return result; + } // RealType quantile(const geometric_distribution dist, p) + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c) + { // Quantile or Percent Point Binomial function. + // Return the number of expected failures k for a given + // complement of the probability Q = 1 - P. + static const char* function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)"; + BOOST_MATH_STD_USING + // Error checks: + RealType x = c.param; + const geometric_distribution<RealType, Policy>& dist = c.dist; + RealType success_fraction = dist.success_fraction(); + RealType result = 0; + if(false == geometric_detail::check_dist_and_prob( + function, + success_fraction, + x, + &result, Policy())) + { + return result; + } + + // Special cases: + if(x == 1) + { // There may actually be no answer to this question, + // since the probability of zero failures may be non-zero, + return 0; // but zero is the best we can do: + } + if (-x <= boost::math::powm1(dist.success_fraction(), dist.successes(), Policy())) + { // q <= cdf(complement(dist, 0)) == pdf(dist, 0) + return 0; // + } + if(x == 0) + { // Probability 1 - Q == 1 so infinite failures to achieve certainty. + // Would need +infinity failures for total confidence. + result = policies::raise_overflow_error<RealType>( + function, + "Probability argument complement is 0, which implies infinite failures !", Policy()); + return result; + // usually means return +std::numeric_limits<RealType>::infinity(); + // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR + } + // log(x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small + result = log(x) / boost::math::log1p(-success_fraction, Policy()) - 1; + return result; + + } // quantile complement + + } // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#if defined (BOOST_MSVC) +# pragma warning(pop) +#endif + +#endif // BOOST_MATH_SPECIAL_GEOMETRIC_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/laplace.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/laplace.hpp new file mode 100644 index 00000000000..09b24c868b5 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/laplace.hpp @@ -0,0 +1,350 @@ +// Copyright Thijs van den Berg, 2008. +// Copyright John Maddock 2008. +// Copyright Paul A. Bristow 2008, 2014. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +// This module implements the Laplace distribution. +// Weisstein, Eric W. "Laplace Distribution." From MathWorld--A Wolfram Web Resource. +// http://mathworld.wolfram.com/LaplaceDistribution.html +// http://en.wikipedia.org/wiki/Laplace_distribution +// +// Abramowitz and Stegun 1972, p 930 +// http://www.math.sfu.ca/~cbm/aands/page_930.htm + +#ifndef BOOST_STATS_LAPLACE_HPP +#define BOOST_STATS_LAPLACE_HPP + +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/distributions/complement.hpp> +#include <boost/math/constants/constants.hpp> +#include <limits> + +namespace boost{ namespace math{ + +#ifdef BOOST_MSVC +# pragma warning(push) +# pragma warning(disable:4127) // conditional expression is constant +#endif + +template <class RealType = double, class Policy = policies::policy<> > +class laplace_distribution +{ +public: + // ---------------------------------- + // public Types + // ---------------------------------- + typedef RealType value_type; + typedef Policy policy_type; + + // ---------------------------------- + // Constructor(s) + // ---------------------------------- + laplace_distribution(RealType l_location = 0, RealType l_scale = 1) + : m_location(l_location), m_scale(l_scale) + { + RealType result; + check_parameters("boost::math::laplace_distribution<%1%>::laplace_distribution()", &result); + } + + + // ---------------------------------- + // Public functions + // ---------------------------------- + + RealType location() const + { + return m_location; + } + + RealType scale() const + { + return m_scale; + } + + bool check_parameters(const char* function, RealType* result) const + { + if(false == detail::check_scale(function, m_scale, result, Policy())) return false; + if(false == detail::check_location(function, m_location, result, Policy())) return false; + return true; + } + +private: + RealType m_location; + RealType m_scale; +}; // class laplace_distribution + +// +// Convenient type synonym for double. +typedef laplace_distribution<double> laplace; + +// +// Non-member functions. +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const laplace_distribution<RealType, Policy>&) +{ + if (std::numeric_limits<RealType>::has_infinity) + { // Can use infinity. + return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity. + } + else + { // Can only use max_value. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. + } + +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const laplace_distribution<RealType, Policy>&) +{ + if (std::numeric_limits<RealType>::has_infinity) + { // Can Use infinity. + return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity. + } + else + { // Can only use max_value. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. + } +} + +template <class RealType, class Policy> +inline RealType pdf(const laplace_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + // Checking function argument + RealType result = 0; + const char* function = "boost::math::pdf(const laplace_distribution<%1%>&, %1%))"; + + // Check scale and location. + if (false == dist.check_parameters(function, &result)) return result; + // Special pdf values. + if((boost::math::isinf)(x)) + { + return 0; // pdf + and - infinity is zero. + } + if (false == detail::check_x(function, x, &result, Policy())) return result; + + // General case + RealType scale( dist.scale() ); + RealType location( dist.location() ); + + RealType exponent = x - location; + if (exponent>0) exponent = -exponent; + exponent /= scale; + + result = exp(exponent); + result /= 2 * scale; + + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const laplace_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // For ADL of std functions. + + RealType result = 0; + // Checking function argument. + const char* function = "boost::math::cdf(const laplace_distribution<%1%>&, %1%)"; + // Check scale and location. + if (false == dist.check_parameters(function, &result)) return result; + + // Special cdf values: + if((boost::math::isinf)(x)) + { + if(x < 0) return 0; // -infinity. + return 1; // + infinity. + } + if (false == detail::check_x(function, x, &result, Policy())) return result; + + // General cdf values + RealType scale( dist.scale() ); + RealType location( dist.location() ); + + if (x < location) + { + result = exp( (x-location)/scale )/2; + } + else + { + result = 1 - exp( (location-x)/scale )/2; + } + return result; +} // cdf + + +template <class RealType, class Policy> +inline RealType quantile(const laplace_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions. + + // Checking function argument + RealType result = 0; + const char* function = "boost::math::quantile(const laplace_distribution<%1%>&, %1%)"; + if (false == dist.check_parameters(function, &result)) return result; + if(false == detail::check_probability(function, p, &result, Policy())) return result; + + // Extreme values of p: + if(p == 0) + { + result = policies::raise_overflow_error<RealType>(function, + "probability parameter is 0, but must be > 0!", Policy()); + return -result; // -std::numeric_limits<RealType>::infinity(); + } + + if(p == 1) + { + result = policies::raise_overflow_error<RealType>(function, + "probability parameter is 1, but must be < 1!", Policy()); + return result; // std::numeric_limits<RealType>::infinity(); + } + // Calculate Quantile + RealType scale( dist.scale() ); + RealType location( dist.location() ); + + if (p - 0.5 < 0.0) + result = location + scale*log( static_cast<RealType>(p*2) ); + else + result = location - scale*log( static_cast<RealType>(-p*2 + 2) ); + + return result; +} // quantile + + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<laplace_distribution<RealType, Policy>, RealType>& c) +{ + // Calculate complement of cdf. + BOOST_MATH_STD_USING // for ADL of std functions + + RealType scale = c.dist.scale(); + RealType location = c.dist.location(); + RealType x = c.param; + RealType result = 0; + + // Checking function argument. + const char* function = "boost::math::cdf(const complemented2_type<laplace_distribution<%1%>, %1%>&)"; + + // Check scale and location. + //if(false == detail::check_scale(function, scale, result, Policy())) return false; + //if(false == detail::check_location(function, location, result, Policy())) return false; + if (false == c.dist.check_parameters(function, &result)) return result; + + // Special cdf values. + if((boost::math::isinf)(x)) + { + if(x < 0) return 1; // cdf complement -infinity is unity. + return 0; // cdf complement +infinity is zero. + } + if(false == detail::check_x(function, x, &result, Policy()))return result; + + // Cdf interval value. + if (-x < -location) + { + result = exp( (-x+location)/scale )/2; + } + else + { + result = 1 - exp( (-location+x)/scale )/2; + } + return result; +} // cdf complement + + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<laplace_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions. + + // Calculate quantile. + RealType scale = c.dist.scale(); + RealType location = c.dist.location(); + RealType q = c.param; + RealType result = 0; + + // Checking function argument. + const char* function = "quantile(const complemented2_type<laplace_distribution<%1%>, %1%>&)"; + if (false == c.dist.check_parameters(function, &result)) return result; + + // Extreme values. + if(q == 0) + { + return std::numeric_limits<RealType>::infinity(); + } + if(q == 1) + { + return -std::numeric_limits<RealType>::infinity(); + } + if(false == detail::check_probability(function, q, &result, Policy())) return result; + + if (0.5 - q < 0.0) + result = location + scale*log( static_cast<RealType>(-q*2 + 2) ); + else + result = location - scale*log( static_cast<RealType>(q*2) ); + + + return result; +} // quantile + +template <class RealType, class Policy> +inline RealType mean(const laplace_distribution<RealType, Policy>& dist) +{ + return dist.location(); +} + +template <class RealType, class Policy> +inline RealType standard_deviation(const laplace_distribution<RealType, Policy>& dist) +{ + return constants::root_two<RealType>() * dist.scale(); +} + +template <class RealType, class Policy> +inline RealType mode(const laplace_distribution<RealType, Policy>& dist) +{ + return dist.location(); +} + +template <class RealType, class Policy> +inline RealType median(const laplace_distribution<RealType, Policy>& dist) +{ + return dist.location(); +} + +template <class RealType, class Policy> +inline RealType skewness(const laplace_distribution<RealType, Policy>& /*dist*/) +{ + return 0; +} + +template <class RealType, class Policy> +inline RealType kurtosis(const laplace_distribution<RealType, Policy>& /*dist*/) +{ + return 6; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const laplace_distribution<RealType, Policy>& /*dist*/) +{ + return 3; +} + +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_LAPLACE_HPP + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/lognormal.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/lognormal.hpp new file mode 100644 index 00000000000..4e6c0610d4b --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/lognormal.hpp @@ -0,0 +1,341 @@ +// Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_LOGNORMAL_HPP +#define BOOST_STATS_LOGNORMAL_HPP + +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm +// http://mathworld.wolfram.com/LogNormalDistribution.html +// http://en.wikipedia.org/wiki/Lognormal_distribution + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/distributions/normal.hpp> +#include <boost/math/special_functions/expm1.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> + +#include <utility> + +namespace boost{ namespace math +{ +namespace detail +{ + + template <class RealType, class Policy> + inline bool check_lognormal_x( + const char* function, + RealType const& x, + RealType* result, const Policy& pol) + { + if((x < 0) || !(boost::math::isfinite)(x)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Random variate is %1% but must be >= 0 !", x, pol); + return false; + } + return true; + } + +} // namespace detail + + +template <class RealType = double, class Policy = policies::policy<> > +class lognormal_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + lognormal_distribution(RealType l_location = 0, RealType l_scale = 1) + : m_location(l_location), m_scale(l_scale) + { + RealType result; + detail::check_scale("boost::math::lognormal_distribution<%1%>::lognormal_distribution", l_scale, &result, Policy()); + detail::check_location("boost::math::lognormal_distribution<%1%>::lognormal_distribution", l_location, &result, Policy()); + } + + RealType location()const + { + return m_location; + } + + RealType scale()const + { + return m_scale; + } +private: + // + // Data members: + // + RealType m_location; // distribution location. + RealType m_scale; // distribution scale. +}; + +typedef lognormal_distribution<double> lognormal; + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const lognormal_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x is >0 to +infinity. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const lognormal_distribution<RealType, Policy>& /*dist*/) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); +} + +template <class RealType, class Policy> +RealType pdf(const lognormal_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType mu = dist.location(); + RealType sigma = dist.scale(); + + static const char* function = "boost::math::pdf(const lognormal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(0 == detail::check_scale(function, sigma, &result, Policy())) + return result; + if(0 == detail::check_location(function, mu, &result, Policy())) + return result; + if(0 == detail::check_lognormal_x(function, x, &result, Policy())) + return result; + + if(x == 0) + return 0; + + RealType exponent = log(x) - mu; + exponent *= -exponent; + exponent /= 2 * sigma * sigma; + + result = exp(exponent); + result /= sigma * sqrt(2 * constants::pi<RealType>()) * x; + + return result; +} + +template <class RealType, class Policy> +inline RealType cdf(const lognormal_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const lognormal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(0 == detail::check_scale(function, dist.scale(), &result, Policy())) + return result; + if(0 == detail::check_location(function, dist.location(), &result, Policy())) + return result; + if(0 == detail::check_lognormal_x(function, x, &result, Policy())) + return result; + + if(x == 0) + return 0; + + normal_distribution<RealType, Policy> norm(dist.location(), dist.scale()); + return cdf(norm, log(x)); +} + +template <class RealType, class Policy> +inline RealType quantile(const lognormal_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const lognormal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(0 == detail::check_scale(function, dist.scale(), &result, Policy())) + return result; + if(0 == detail::check_location(function, dist.location(), &result, Policy())) + return result; + if(0 == detail::check_probability(function, p, &result, Policy())) + return result; + + if(p == 0) + return 0; + if(p == 1) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + normal_distribution<RealType, Policy> norm(dist.location(), dist.scale()); + return exp(quantile(norm, p)); +} + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<lognormal_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const lognormal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(0 == detail::check_scale(function, c.dist.scale(), &result, Policy())) + return result; + if(0 == detail::check_location(function, c.dist.location(), &result, Policy())) + return result; + if(0 == detail::check_lognormal_x(function, c.param, &result, Policy())) + return result; + + if(c.param == 0) + return 1; + + normal_distribution<RealType, Policy> norm(c.dist.location(), c.dist.scale()); + return cdf(complement(norm, log(c.param))); +} + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<lognormal_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const lognormal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(0 == detail::check_scale(function, c.dist.scale(), &result, Policy())) + return result; + if(0 == detail::check_location(function, c.dist.location(), &result, Policy())) + return result; + if(0 == detail::check_probability(function, c.param, &result, Policy())) + return result; + + if(c.param == 1) + return 0; + if(c.param == 0) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + normal_distribution<RealType, Policy> norm(c.dist.location(), c.dist.scale()); + return exp(quantile(complement(norm, c.param))); +} + +template <class RealType, class Policy> +inline RealType mean(const lognormal_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType mu = dist.location(); + RealType sigma = dist.scale(); + + RealType result = 0; + if(0 == detail::check_scale("boost::math::mean(const lognormal_distribution<%1%>&)", sigma, &result, Policy())) + return result; + if(0 == detail::check_location("boost::math::mean(const lognormal_distribution<%1%>&)", mu, &result, Policy())) + return result; + + return exp(mu + sigma * sigma / 2); +} + +template <class RealType, class Policy> +inline RealType variance(const lognormal_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType mu = dist.location(); + RealType sigma = dist.scale(); + + RealType result = 0; + if(0 == detail::check_scale("boost::math::variance(const lognormal_distribution<%1%>&)", sigma, &result, Policy())) + return result; + if(0 == detail::check_location("boost::math::variance(const lognormal_distribution<%1%>&)", mu, &result, Policy())) + return result; + + return boost::math::expm1(sigma * sigma, Policy()) * exp(2 * mu + sigma * sigma); +} + +template <class RealType, class Policy> +inline RealType mode(const lognormal_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType mu = dist.location(); + RealType sigma = dist.scale(); + + RealType result = 0; + if(0 == detail::check_scale("boost::math::mode(const lognormal_distribution<%1%>&)", sigma, &result, Policy())) + return result; + if(0 == detail::check_location("boost::math::mode(const lognormal_distribution<%1%>&)", mu, &result, Policy())) + return result; + + return exp(mu - sigma * sigma); +} + +template <class RealType, class Policy> +inline RealType median(const lognormal_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + RealType mu = dist.location(); + return exp(mu); // e^mu +} + +template <class RealType, class Policy> +inline RealType skewness(const lognormal_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + //RealType mu = dist.location(); + RealType sigma = dist.scale(); + + RealType ss = sigma * sigma; + RealType ess = exp(ss); + + RealType result = 0; + if(0 == detail::check_scale("boost::math::skewness(const lognormal_distribution<%1%>&)", sigma, &result, Policy())) + return result; + if(0 == detail::check_location("boost::math::skewness(const lognormal_distribution<%1%>&)", dist.location(), &result, Policy())) + return result; + + return (ess + 2) * sqrt(boost::math::expm1(ss, Policy())); +} + +template <class RealType, class Policy> +inline RealType kurtosis(const lognormal_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + //RealType mu = dist.location(); + RealType sigma = dist.scale(); + RealType ss = sigma * sigma; + + RealType result = 0; + if(0 == detail::check_scale("boost::math::kurtosis(const lognormal_distribution<%1%>&)", sigma, &result, Policy())) + return result; + if(0 == detail::check_location("boost::math::kurtosis(const lognormal_distribution<%1%>&)", dist.location(), &result, Policy())) + return result; + + return exp(4 * ss) + 2 * exp(3 * ss) + 3 * exp(2 * ss) - 3; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const lognormal_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + // RealType mu = dist.location(); + RealType sigma = dist.scale(); + RealType ss = sigma * sigma; + + RealType result = 0; + if(0 == detail::check_scale("boost::math::kurtosis_excess(const lognormal_distribution<%1%>&)", sigma, &result, Policy())) + return result; + if(0 == detail::check_location("boost::math::kurtosis_excess(const lognormal_distribution<%1%>&)", dist.location(), &result, Policy())) + return result; + + return exp(4 * ss) + 2 * exp(3 * ss) + 3 * exp(2 * ss) - 6; +} + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_STUDENTS_T_HPP + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/negative_binomial.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/negative_binomial.hpp new file mode 100644 index 00000000000..ca5723fa7d4 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/negative_binomial.hpp @@ -0,0 +1,607 @@ +// boost\math\special_functions\negative_binomial.hpp + +// Copyright Paul A. Bristow 2007. +// Copyright John Maddock 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// http://en.wikipedia.org/wiki/negative_binomial_distribution +// http://mathworld.wolfram.com/NegativeBinomialDistribution.html +// http://documents.wolfram.com/teachersedition/Teacher/Statistics/DiscreteDistributions.html + +// The negative binomial distribution NegativeBinomialDistribution[n, p] +// is the distribution of the number (k) of failures that occur in a sequence of trials before +// r successes have occurred, where the probability of success in each trial is p. + +// In a sequence of Bernoulli trials or events +// (independent, yes or no, succeed or fail) with success_fraction probability p, +// negative_binomial is the probability that k or fewer failures +// preceed the r th trial's success. +// random variable k is the number of failures (NOT the probability). + +// Negative_binomial distribution is a discrete probability distribution. +// But note that the negative binomial distribution +// (like others including the binomial, Poisson & Bernoulli) +// is strictly defined as a discrete function: only integral values of k are envisaged. +// However because of the method of calculation using a continuous gamma function, +// it is convenient to treat it as if a continous function, +// and permit non-integral values of k. + +// However, by default the policy is to use discrete_quantile_policy. + +// To enforce the strict mathematical model, users should use conversion +// on k outside this function to ensure that k is integral. + +// MATHCAD cumulative negative binomial pnbinom(k, n, p) + +// Implementation note: much greater speed, and perhaps greater accuracy, +// might be achieved for extreme values by using a normal approximation. +// This is NOT been tested or implemented. + +#ifndef BOOST_MATH_SPECIAL_NEGATIVE_BINOMIAL_HPP +#define BOOST_MATH_SPECIAL_NEGATIVE_BINOMIAL_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x) == Ix(a, b). +#include <boost/math/distributions/complement.hpp> // complement. +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks domain_error & logic_error. +#include <boost/math/special_functions/fpclassify.hpp> // isnan. +#include <boost/math/tools/roots.hpp> // for root finding. +#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> + +#include <boost/type_traits/is_floating_point.hpp> +#include <boost/type_traits/is_integral.hpp> +#include <boost/type_traits/is_same.hpp> +#include <boost/mpl/if.hpp> + +#include <limits> // using std::numeric_limits; +#include <utility> + +#if defined (BOOST_MSVC) +# pragma warning(push) +// This believed not now necessary, so commented out. +//# pragma warning(disable: 4702) // unreachable code. +// in domain_error_imp in error_handling. +#endif + +namespace boost +{ + namespace math + { + namespace negative_binomial_detail + { + // Common error checking routines for negative binomial distribution functions: + template <class RealType, class Policy> + inline bool check_successes(const char* function, const RealType& r, RealType* result, const Policy& pol) + { + if( !(boost::math::isfinite)(r) || (r <= 0) ) + { + *result = policies::raise_domain_error<RealType>( + function, + "Number of successes argument is %1%, but must be > 0 !", r, pol); + return false; + } + return true; + } + template <class RealType, class Policy> + inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol) + { + if( !(boost::math::isfinite)(p) || (p < 0) || (p > 1) ) + { + *result = policies::raise_domain_error<RealType>( + function, + "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol); + return false; + } + return true; + } + template <class RealType, class Policy> + inline bool check_dist(const char* function, const RealType& r, const RealType& p, RealType* result, const Policy& pol) + { + return check_success_fraction(function, p, result, pol) + && check_successes(function, r, result, pol); + } + template <class RealType, class Policy> + inline bool check_dist_and_k(const char* function, const RealType& r, const RealType& p, RealType k, RealType* result, const Policy& pol) + { + if(check_dist(function, r, p, result, pol) == false) + { + return false; + } + if( !(boost::math::isfinite)(k) || (k < 0) ) + { // Check k failures. + *result = policies::raise_domain_error<RealType>( + function, + "Number of failures argument is %1%, but must be >= 0 !", k, pol); + return false; + } + return true; + } // Check_dist_and_k + + template <class RealType, class Policy> + inline bool check_dist_and_prob(const char* function, const RealType& r, RealType p, RealType prob, RealType* result, const Policy& pol) + { + if(check_dist(function, r, p, result, pol) && detail::check_probability(function, prob, result, pol) == false) + { + return false; + } + return true; + } // check_dist_and_prob + } // namespace negative_binomial_detail + + template <class RealType = double, class Policy = policies::policy<> > + class negative_binomial_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + negative_binomial_distribution(RealType r, RealType p) : m_r(r), m_p(p) + { // Constructor. + RealType result; + negative_binomial_detail::check_dist( + "negative_binomial_distribution<%1%>::negative_binomial_distribution", + m_r, // Check successes r > 0. + m_p, // Check success_fraction 0 <= p <= 1. + &result, Policy()); + } // negative_binomial_distribution constructor. + + // Private data getter class member functions. + RealType success_fraction() const + { // Probability of success as fraction in range 0 to 1. + return m_p; + } + RealType successes() const + { // Total number of successes r. + return m_r; + } + + static RealType find_lower_bound_on_p( + RealType trials, + RealType successes, + RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test. + { + static const char* function = "boost::math::negative_binomial<%1%>::find_lower_bound_on_p"; + RealType result = 0; // of error checks. + RealType failures = trials - successes; + if(false == detail::check_probability(function, alpha, &result, Policy()) + && negative_binomial_detail::check_dist_and_k( + function, successes, RealType(0), failures, &result, Policy())) + { + return result; + } + // Use complement ibeta_inv function for lower bound. + // This is adapted from the corresponding binomial formula + // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm + // This is a Clopper-Pearson interval, and may be overly conservative, + // see also "A Simple Improved Inferential Method for Some + // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY + // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf + // + return ibeta_inv(successes, failures + 1, alpha, static_cast<RealType*>(0), Policy()); + } // find_lower_bound_on_p + + static RealType find_upper_bound_on_p( + RealType trials, + RealType successes, + RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test. + { + static const char* function = "boost::math::negative_binomial<%1%>::find_upper_bound_on_p"; + RealType result = 0; // of error checks. + RealType failures = trials - successes; + if(false == negative_binomial_detail::check_dist_and_k( + function, successes, RealType(0), failures, &result, Policy()) + && detail::check_probability(function, alpha, &result, Policy())) + { + return result; + } + if(failures == 0) + return 1; + // Use complement ibetac_inv function for upper bound. + // Note adjusted failures value: *not* failures+1 as usual. + // This is adapted from the corresponding binomial formula + // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm + // This is a Clopper-Pearson interval, and may be overly conservative, + // see also "A Simple Improved Inferential Method for Some + // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY + // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf + // + return ibetac_inv(successes, failures, alpha, static_cast<RealType*>(0), Policy()); + } // find_upper_bound_on_p + + // Estimate number of trials : + // "How many trials do I need to be P% sure of seeing k or fewer failures?" + + static RealType find_minimum_number_of_trials( + RealType k, // number of failures (k >= 0). + RealType p, // success fraction 0 <= p <= 1. + RealType alpha) // risk level threshold 0 <= alpha <= 1. + { + static const char* function = "boost::math::negative_binomial<%1%>::find_minimum_number_of_trials"; + // Error checks: + RealType result = 0; + if(false == negative_binomial_detail::check_dist_and_k( + function, RealType(1), p, k, &result, Policy()) + && detail::check_probability(function, alpha, &result, Policy())) + { return result; } + + result = ibeta_inva(k + 1, p, alpha, Policy()); // returns n - k + return result + k; + } // RealType find_number_of_failures + + static RealType find_maximum_number_of_trials( + RealType k, // number of failures (k >= 0). + RealType p, // success fraction 0 <= p <= 1. + RealType alpha) // risk level threshold 0 <= alpha <= 1. + { + static const char* function = "boost::math::negative_binomial<%1%>::find_maximum_number_of_trials"; + // Error checks: + RealType result = 0; + if(false == negative_binomial_detail::check_dist_and_k( + function, RealType(1), p, k, &result, Policy()) + && detail::check_probability(function, alpha, &result, Policy())) + { return result; } + + result = ibetac_inva(k + 1, p, alpha, Policy()); // returns n - k + return result + k; + } // RealType find_number_of_trials complemented + + private: + RealType m_r; // successes. + RealType m_p; // success_fraction + }; // template <class RealType, class Policy> class negative_binomial_distribution + + typedef negative_binomial_distribution<double> negative_binomial; // Reserved name of type double. + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const negative_binomial_distribution<RealType, Policy>& /* dist */) + { // Range of permissible values for random variable k. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer? + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const negative_binomial_distribution<RealType, Policy>& /* dist */) + { // Range of supported values for random variable k. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer? + } + + template <class RealType, class Policy> + inline RealType mean(const negative_binomial_distribution<RealType, Policy>& dist) + { // Mean of Negative Binomial distribution = r(1-p)/p. + return dist.successes() * (1 - dist.success_fraction() ) / dist.success_fraction(); + } // mean + + //template <class RealType, class Policy> + //inline RealType median(const negative_binomial_distribution<RealType, Policy>& dist) + //{ // Median of negative_binomial_distribution is not defined. + // return policies::raise_domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); + //} // median + // Now implemented via quantile(half) in derived accessors. + + template <class RealType, class Policy> + inline RealType mode(const negative_binomial_distribution<RealType, Policy>& dist) + { // Mode of Negative Binomial distribution = floor[(r-1) * (1 - p)/p] + BOOST_MATH_STD_USING // ADL of std functions. + return floor((dist.successes() -1) * (1 - dist.success_fraction()) / dist.success_fraction()); + } // mode + + template <class RealType, class Policy> + inline RealType skewness(const negative_binomial_distribution<RealType, Policy>& dist) + { // skewness of Negative Binomial distribution = 2-p / (sqrt(r(1-p)) + BOOST_MATH_STD_USING // ADL of std functions. + RealType p = dist.success_fraction(); + RealType r = dist.successes(); + + return (2 - p) / + sqrt(r * (1 - p)); + } // skewness + + template <class RealType, class Policy> + inline RealType kurtosis(const negative_binomial_distribution<RealType, Policy>& dist) + { // kurtosis of Negative Binomial distribution + // http://en.wikipedia.org/wiki/Negative_binomial is kurtosis_excess so add 3 + RealType p = dist.success_fraction(); + RealType r = dist.successes(); + return 3 + (6 / r) + ((p * p) / (r * (1 - p))); + } // kurtosis + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const negative_binomial_distribution<RealType, Policy>& dist) + { // kurtosis excess of Negative Binomial distribution + // http://mathworld.wolfram.com/Kurtosis.html table of kurtosis_excess + RealType p = dist.success_fraction(); + RealType r = dist.successes(); + return (6 - p * (6-p)) / (r * (1-p)); + } // kurtosis_excess + + template <class RealType, class Policy> + inline RealType variance(const negative_binomial_distribution<RealType, Policy>& dist) + { // Variance of Binomial distribution = r (1-p) / p^2. + return dist.successes() * (1 - dist.success_fraction()) + / (dist.success_fraction() * dist.success_fraction()); + } // variance + + // RealType standard_deviation(const negative_binomial_distribution<RealType, Policy>& dist) + // standard_deviation provided by derived accessors. + // RealType hazard(const negative_binomial_distribution<RealType, Policy>& dist) + // hazard of Negative Binomial distribution provided by derived accessors. + // RealType chf(const negative_binomial_distribution<RealType, Policy>& dist) + // chf of Negative Binomial distribution provided by derived accessors. + + template <class RealType, class Policy> + inline RealType pdf(const negative_binomial_distribution<RealType, Policy>& dist, const RealType& k) + { // Probability Density/Mass Function. + BOOST_FPU_EXCEPTION_GUARD + + static const char* function = "boost::math::pdf(const negative_binomial_distribution<%1%>&, %1%)"; + + RealType r = dist.successes(); + RealType p = dist.success_fraction(); + RealType result = 0; + if(false == negative_binomial_detail::check_dist_and_k( + function, + r, + dist.success_fraction(), + k, + &result, Policy())) + { + return result; + } + + result = (p/(r + k)) * ibeta_derivative(r, static_cast<RealType>(k+1), p, Policy()); + // Equivalent to: + // return exp(lgamma(r + k) - lgamma(r) - lgamma(k+1)) * pow(p, r) * pow((1-p), k); + return result; + } // negative_binomial_pdf + + template <class RealType, class Policy> + inline RealType cdf(const negative_binomial_distribution<RealType, Policy>& dist, const RealType& k) + { // Cumulative Distribution Function of Negative Binomial. + static const char* function = "boost::math::cdf(const negative_binomial_distribution<%1%>&, %1%)"; + using boost::math::ibeta; // Regularized incomplete beta function. + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + RealType p = dist.success_fraction(); + RealType r = dist.successes(); + // Error check: + RealType result = 0; + if(false == negative_binomial_detail::check_dist_and_k( + function, + r, + dist.success_fraction(), + k, + &result, Policy())) + { + return result; + } + + RealType probability = ibeta(r, static_cast<RealType>(k+1), p, Policy()); + // Ip(r, k+1) = ibeta(r, k+1, p) + return probability; + } // cdf Cumulative Distribution Function Negative Binomial. + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<negative_binomial_distribution<RealType, Policy>, RealType>& c) + { // Complemented Cumulative Distribution Function Negative Binomial. + + static const char* function = "boost::math::cdf(const negative_binomial_distribution<%1%>&, %1%)"; + using boost::math::ibetac; // Regularized incomplete beta function complement. + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + RealType const& k = c.param; + negative_binomial_distribution<RealType, Policy> const& dist = c.dist; + RealType p = dist.success_fraction(); + RealType r = dist.successes(); + // Error check: + RealType result = 0; + if(false == negative_binomial_detail::check_dist_and_k( + function, + r, + p, + k, + &result, Policy())) + { + return result; + } + // Calculate cdf negative binomial using the incomplete beta function. + // Use of ibeta here prevents cancellation errors in calculating + // 1-p if p is very small, perhaps smaller than machine epsilon. + // Ip(k+1, r) = ibetac(r, k+1, p) + // constrain_probability here? + RealType probability = ibetac(r, static_cast<RealType>(k+1), p, Policy()); + // Numerical errors might cause probability to be slightly outside the range < 0 or > 1. + // This might cause trouble downstream, so warn, possibly throw exception, but constrain to the limits. + return probability; + } // cdf Cumulative Distribution Function Negative Binomial. + + template <class RealType, class Policy> + inline RealType quantile(const negative_binomial_distribution<RealType, Policy>& dist, const RealType& P) + { // Quantile, percentile/100 or Percent Point Negative Binomial function. + // Return the number of expected failures k for a given probability p. + + // Inverse cumulative Distribution Function or Quantile (percentile / 100) of negative_binomial Probability. + // MAthCAD pnbinom return smallest k such that negative_binomial(k, n, p) >= probability. + // k argument may be integral, signed, or unsigned, or floating point. + // BUT Cephes/CodeCogs says: finds argument p (0 to 1) such that cdf(k, n, p) = y + static const char* function = "boost::math::quantile(const negative_binomial_distribution<%1%>&, %1%)"; + BOOST_MATH_STD_USING // ADL of std functions. + + RealType p = dist.success_fraction(); + RealType r = dist.successes(); + // Check dist and P. + RealType result = 0; + if(false == negative_binomial_detail::check_dist_and_prob + (function, r, p, P, &result, Policy())) + { + return result; + } + + // Special cases. + if (P == 1) + { // Would need +infinity failures for total confidence. + result = policies::raise_overflow_error<RealType>( + function, + "Probability argument is 1, which implies infinite failures !", Policy()); + return result; + // usually means return +std::numeric_limits<RealType>::infinity(); + // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR + } + if (P == 0) + { // No failures are expected if P = 0. + return 0; // Total trials will be just dist.successes. + } + if (P <= pow(dist.success_fraction(), dist.successes())) + { // p <= pdf(dist, 0) == cdf(dist, 0) + return 0; + } + if(p == 0) + { // Would need +infinity failures for total confidence. + result = policies::raise_overflow_error<RealType>( + function, + "Success fraction is 0, which implies infinite failures !", Policy()); + return result; + // usually means return +std::numeric_limits<RealType>::infinity(); + // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR + } + /* + // Calculate quantile of negative_binomial using the inverse incomplete beta function. + using boost::math::ibeta_invb; + return ibeta_invb(r, p, P, Policy()) - 1; // + */ + RealType guess = 0; + RealType factor = 5; + if(r * r * r * P * p > 0.005) + guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), P, RealType(1-P), Policy()); + + if(guess < 10) + { + // + // Cornish-Fisher Negative binomial approximation not accurate in this area: + // + guess = (std::min)(RealType(r * 2), RealType(10)); + } + else + factor = (1-P < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f); + BOOST_MATH_INSTRUMENT_CODE("guess = " << guess); + // + // Max iterations permitted: + // + boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); + typedef typename Policy::discrete_quantile_type discrete_type; + return detail::inverse_discrete_quantile( + dist, + P, + false, + guess, + factor, + RealType(1), + discrete_type(), + max_iter); + } // RealType quantile(const negative_binomial_distribution dist, p) + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<negative_binomial_distribution<RealType, Policy>, RealType>& c) + { // Quantile or Percent Point Binomial function. + // Return the number of expected failures k for a given + // complement of the probability Q = 1 - P. + static const char* function = "boost::math::quantile(const negative_binomial_distribution<%1%>&, %1%)"; + BOOST_MATH_STD_USING + + // Error checks: + RealType Q = c.param; + const negative_binomial_distribution<RealType, Policy>& dist = c.dist; + RealType p = dist.success_fraction(); + RealType r = dist.successes(); + RealType result = 0; + if(false == negative_binomial_detail::check_dist_and_prob( + function, + r, + p, + Q, + &result, Policy())) + { + return result; + } + + // Special cases: + // + if(Q == 1) + { // There may actually be no answer to this question, + // since the probability of zero failures may be non-zero, + return 0; // but zero is the best we can do: + } + if(Q == 0) + { // Probability 1 - Q == 1 so infinite failures to achieve certainty. + // Would need +infinity failures for total confidence. + result = policies::raise_overflow_error<RealType>( + function, + "Probability argument complement is 0, which implies infinite failures !", Policy()); + return result; + // usually means return +std::numeric_limits<RealType>::infinity(); + // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR + } + if (-Q <= boost::math::powm1(dist.success_fraction(), dist.successes(), Policy())) + { // q <= cdf(complement(dist, 0)) == pdf(dist, 0) + return 0; // + } + if(p == 0) + { // Success fraction is 0 so infinite failures to achieve certainty. + // Would need +infinity failures for total confidence. + result = policies::raise_overflow_error<RealType>( + function, + "Success fraction is 0, which implies infinite failures !", Policy()); + return result; + // usually means return +std::numeric_limits<RealType>::infinity(); + // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR + } + //return ibetac_invb(r, p, Q, Policy()) -1; + RealType guess = 0; + RealType factor = 5; + if(r * r * r * (1-Q) * p > 0.005) + guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), RealType(1-Q), Q, Policy()); + + if(guess < 10) + { + // + // Cornish-Fisher Negative binomial approximation not accurate in this area: + // + guess = (std::min)(RealType(r * 2), RealType(10)); + } + else + factor = (Q < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f); + BOOST_MATH_INSTRUMENT_CODE("guess = " << guess); + // + // Max iterations permitted: + // + boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); + typedef typename Policy::discrete_quantile_type discrete_type; + return detail::inverse_discrete_quantile( + dist, + Q, + true, + guess, + factor, + RealType(1), + discrete_type(), + max_iter); + } // quantile complement + + } // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#if defined (BOOST_MSVC) +# pragma warning(pop) +#endif + +#endif // BOOST_MATH_SPECIAL_NEGATIVE_BINOMIAL_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/normal.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/normal.hpp new file mode 100644 index 00000000000..32cf66e3ef0 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/normal.hpp @@ -0,0 +1,329 @@ +// Copyright John Maddock 2006, 2007. +// Copyright Paul A. Bristow 2006, 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_NORMAL_HPP +#define BOOST_STATS_NORMAL_HPP + +// http://en.wikipedia.org/wiki/Normal_distribution +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm +// Also: +// Weisstein, Eric W. "Normal Distribution." +// From MathWorld--A Wolfram Web Resource. +// http://mathworld.wolfram.com/NormalDistribution.html + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/erf.hpp> // for erf/erfc. +#include <boost/math/distributions/complement.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> + +#include <utility> + +namespace boost{ namespace math{ + +template <class RealType = double, class Policy = policies::policy<> > +class normal_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + normal_distribution(RealType l_mean = 0, RealType sd = 1) + : m_mean(l_mean), m_sd(sd) + { // Default is a 'standard' normal distribution N01. + static const char* function = "boost::math::normal_distribution<%1%>::normal_distribution"; + + RealType result; + detail::check_scale(function, sd, &result, Policy()); + detail::check_location(function, l_mean, &result, Policy()); + } + + RealType mean()const + { // alias for location. + return m_mean; + } + + RealType standard_deviation()const + { // alias for scale. + return m_sd; + } + + // Synonyms, provided to allow generic use of find_location and find_scale. + RealType location()const + { // location. + return m_mean; + } + RealType scale()const + { // scale. + return m_sd; + } + +private: + // + // Data members: + // + RealType m_mean; // distribution mean or location. + RealType m_sd; // distribution standard deviation or scale. +}; // class normal_distribution + +typedef normal_distribution<double> normal; + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4127) +#endif + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const normal_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x. + if (std::numeric_limits<RealType>::has_infinity) + { + return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity. + } + else + { // Can only use max_value. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. + } +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const normal_distribution<RealType, Policy>& /*dist*/) +{ // This is range values for random variable x where cdf rises from 0 to 1, and outside it, the pdf is zero. + if (std::numeric_limits<RealType>::has_infinity) + { + return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity. + } + else + { // Can only use max_value. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. + } +} + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif + +template <class RealType, class Policy> +inline RealType pdf(const normal_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType sd = dist.standard_deviation(); + RealType mean = dist.mean(); + + static const char* function = "boost::math::pdf(const normal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(false == detail::check_scale(function, sd, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, mean, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + return 0; // pdf + and - infinity is zero. + } + // Below produces MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity()) + //{ // pdf + and - infinity is zero. + // return 0; + //} + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + RealType exponent = x - mean; + exponent *= -exponent; + exponent /= 2 * sd * sd; + + result = exp(exponent); + result /= sd * sqrt(2 * constants::pi<RealType>()); + + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const normal_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType sd = dist.standard_deviation(); + RealType mean = dist.mean(); + static const char* function = "boost::math::cdf(const normal_distribution<%1%>&, %1%)"; + RealType result = 0; + if(false == detail::check_scale(function, sd, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, mean, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + if(x < 0) return 0; // -infinity + return 1; // + infinity + } + // These produce MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) + //{ // cdf +infinity is unity. + // return 1; + //} + //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) + //{ // cdf -infinity is zero. + // return 0; + //} + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + RealType diff = (x - mean) / (sd * constants::root_two<RealType>()); + result = boost::math::erfc(-diff, Policy()) / 2; + return result; +} // cdf + +template <class RealType, class Policy> +inline RealType quantile(const normal_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType sd = dist.standard_deviation(); + RealType mean = dist.mean(); + static const char* function = "boost::math::quantile(const normal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(false == detail::check_scale(function, sd, &result, Policy())) + return result; + if(false == detail::check_location(function, mean, &result, Policy())) + return result; + if(false == detail::check_probability(function, p, &result, Policy())) + return result; + + result= boost::math::erfc_inv(2 * p, Policy()); + result = -result; + result *= sd * constants::root_two<RealType>(); + result += mean; + return result; +} // quantile + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<normal_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType sd = c.dist.standard_deviation(); + RealType mean = c.dist.mean(); + RealType x = c.param; + static const char* function = "boost::math::cdf(const complement(normal_distribution<%1%>&), %1%)"; + + RealType result = 0; + if(false == detail::check_scale(function, sd, &result, Policy())) + return result; + if(false == detail::check_location(function, mean, &result, Policy())) + return result; + if((boost::math::isinf)(x)) + { + if(x < 0) return 1; // cdf complement -infinity is unity. + return 0; // cdf complement +infinity is zero + } + // These produce MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) + //{ // cdf complement +infinity is zero. + // return 0; + //} + //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) + //{ // cdf complement -infinity is unity. + // return 1; + //} + if(false == detail::check_x(function, x, &result, Policy())) + return result; + + RealType diff = (x - mean) / (sd * constants::root_two<RealType>()); + result = boost::math::erfc(diff, Policy()) / 2; + return result; +} // cdf complement + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<normal_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + RealType sd = c.dist.standard_deviation(); + RealType mean = c.dist.mean(); + static const char* function = "boost::math::quantile(const complement(normal_distribution<%1%>&), %1%)"; + RealType result = 0; + if(false == detail::check_scale(function, sd, &result, Policy())) + return result; + if(false == detail::check_location(function, mean, &result, Policy())) + return result; + RealType q = c.param; + if(false == detail::check_probability(function, q, &result, Policy())) + return result; + result = boost::math::erfc_inv(2 * q, Policy()); + result *= sd * constants::root_two<RealType>(); + result += mean; + return result; +} // quantile + +template <class RealType, class Policy> +inline RealType mean(const normal_distribution<RealType, Policy>& dist) +{ + return dist.mean(); +} + +template <class RealType, class Policy> +inline RealType standard_deviation(const normal_distribution<RealType, Policy>& dist) +{ + return dist.standard_deviation(); +} + +template <class RealType, class Policy> +inline RealType mode(const normal_distribution<RealType, Policy>& dist) +{ + return dist.mean(); +} + +template <class RealType, class Policy> +inline RealType median(const normal_distribution<RealType, Policy>& dist) +{ + return dist.mean(); +} + +template <class RealType, class Policy> +inline RealType skewness(const normal_distribution<RealType, Policy>& /*dist*/) +{ + return 0; +} + +template <class RealType, class Policy> +inline RealType kurtosis(const normal_distribution<RealType, Policy>& /*dist*/) +{ + return 3; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const normal_distribution<RealType, Policy>& /*dist*/) +{ + return 0; +} + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_NORMAL_HPP + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/poisson.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/poisson.hpp new file mode 100644 index 00000000000..e4665bff69b --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/poisson.hpp @@ -0,0 +1,527 @@ +// boost\math\distributions\poisson.hpp + +// Copyright John Maddock 2006. +// Copyright Paul A. Bristow 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// Poisson distribution is a discrete probability distribution. +// It expresses the probability of a number (k) of +// events, occurrences, failures or arrivals occurring in a fixed time, +// assuming these events occur with a known average or mean rate (lambda) +// and are independent of the time since the last event. +// The distribution was discovered by Simeon-Denis Poisson (1781-1840). + +// Parameter lambda is the mean number of events in the given time interval. +// The random variate k is the number of events, occurrences or arrivals. +// k argument may be integral, signed, or unsigned, or floating point. +// If necessary, it has already been promoted from an integral type. + +// Note that the Poisson distribution +// (like others including the binomial, negative binomial & Bernoulli) +// is strictly defined as a discrete function: +// only integral values of k are envisaged. +// However because the method of calculation uses a continuous gamma function, +// it is convenient to treat it as if a continous function, +// and permit non-integral values of k. +// To enforce the strict mathematical model, users should use floor or ceil functions +// on k outside this function to ensure that k is integral. + +// See http://en.wikipedia.org/wiki/Poisson_distribution +// http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/Statistics/DiscreteDistributions.html + +#ifndef BOOST_MATH_SPECIAL_POISSON_HPP +#define BOOST_MATH_SPECIAL_POISSON_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q +#include <boost/math/special_functions/trunc.hpp> // for incomplete gamma. gamma_q +#include <boost/math/distributions/complement.hpp> // complements +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks +#include <boost/math/special_functions/fpclassify.hpp> // isnan. +#include <boost/math/special_functions/factorials.hpp> // factorials. +#include <boost/math/tools/roots.hpp> // for root finding. +#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> + +#include <utility> + +namespace boost +{ + namespace math + { + namespace poisson_detail + { + // Common error checking routines for Poisson distribution functions. + // These are convoluted, & apparently redundant, to try to ensure that + // checks are always performed, even if exceptions are not enabled. + + template <class RealType, class Policy> + inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol) + { + if(!(boost::math::isfinite)(mean) || (mean < 0)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Mean argument is %1%, but must be >= 0 !", mean, pol); + return false; + } + return true; + } // bool check_mean + + template <class RealType, class Policy> + inline bool check_mean_NZ(const char* function, const RealType& mean, RealType* result, const Policy& pol) + { // mean == 0 is considered an error. + if( !(boost::math::isfinite)(mean) || (mean <= 0)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Mean argument is %1%, but must be > 0 !", mean, pol); + return false; + } + return true; + } // bool check_mean_NZ + + template <class RealType, class Policy> + inline bool check_dist(const char* function, const RealType& mean, RealType* result, const Policy& pol) + { // Only one check, so this is redundant really but should be optimized away. + return check_mean_NZ(function, mean, result, pol); + } // bool check_dist + + template <class RealType, class Policy> + inline bool check_k(const char* function, const RealType& k, RealType* result, const Policy& pol) + { + if((k < 0) || !(boost::math::isfinite)(k)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Number of events k argument is %1%, but must be >= 0 !", k, pol); + return false; + } + return true; + } // bool check_k + + template <class RealType, class Policy> + inline bool check_dist_and_k(const char* function, RealType mean, RealType k, RealType* result, const Policy& pol) + { + if((check_dist(function, mean, result, pol) == false) || + (check_k(function, k, result, pol) == false)) + { + return false; + } + return true; + } // bool check_dist_and_k + + template <class RealType, class Policy> + inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) + { // Check 0 <= p <= 1 + if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); + return false; + } + return true; + } // bool check_prob + + template <class RealType, class Policy> + inline bool check_dist_and_prob(const char* function, RealType mean, RealType p, RealType* result, const Policy& pol) + { + if((check_dist(function, mean, result, pol) == false) || + (check_prob(function, p, result, pol) == false)) + { + return false; + } + return true; + } // bool check_dist_and_prob + + } // namespace poisson_detail + + template <class RealType = double, class Policy = policies::policy<> > + class poisson_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + poisson_distribution(RealType l_mean = 1) : m_l(l_mean) // mean (lambda). + { // Expected mean number of events that occur during the given interval. + RealType r; + poisson_detail::check_dist( + "boost::math::poisson_distribution<%1%>::poisson_distribution", + m_l, + &r, Policy()); + } // poisson_distribution constructor. + + RealType mean() const + { // Private data getter function. + return m_l; + } + private: + // Data member, initialized by constructor. + RealType m_l; // mean number of occurrences. + }; // template <class RealType, class Policy> class poisson_distribution + + typedef poisson_distribution<double> poisson; // Reserved name of type double. + + // Non-member functions to give properties of the distribution. + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const poisson_distribution<RealType, Policy>& /* dist */) + { // Range of permissible values for random variable k. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer? + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const poisson_distribution<RealType, Policy>& /* dist */) + { // Range of supported values for random variable k. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); + } + + template <class RealType, class Policy> + inline RealType mean(const poisson_distribution<RealType, Policy>& dist) + { // Mean of poisson distribution = lambda. + return dist.mean(); + } // mean + + template <class RealType, class Policy> + inline RealType mode(const poisson_distribution<RealType, Policy>& dist) + { // mode. + BOOST_MATH_STD_USING // ADL of std functions. + return floor(dist.mean()); + } + + //template <class RealType, class Policy> + //inline RealType median(const poisson_distribution<RealType, Policy>& dist) + //{ // median = approximately lambda + 1/3 - 0.2/lambda + // RealType l = dist.mean(); + // return dist.mean() + static_cast<RealType>(0.3333333333333333333333333333333333333333333333) + // - static_cast<RealType>(0.2) / l; + //} // BUT this formula appears to be out-by-one compared to quantile(half) + // Query posted on Wikipedia. + // Now implemented via quantile(half) in derived accessors. + + template <class RealType, class Policy> + inline RealType variance(const poisson_distribution<RealType, Policy>& dist) + { // variance. + return dist.mean(); + } + + // RealType standard_deviation(const poisson_distribution<RealType, Policy>& dist) + // standard_deviation provided by derived accessors. + + template <class RealType, class Policy> + inline RealType skewness(const poisson_distribution<RealType, Policy>& dist) + { // skewness = sqrt(l). + BOOST_MATH_STD_USING // ADL of std functions. + return 1 / sqrt(dist.mean()); + } + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const poisson_distribution<RealType, Policy>& dist) + { // skewness = sqrt(l). + return 1 / dist.mean(); // kurtosis_excess 1/mean from Wiki & MathWorld eq 31. + // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess + // is more convenient because the kurtosis excess of a normal distribution is zero + // whereas the true kurtosis is 3. + } // RealType kurtosis_excess + + template <class RealType, class Policy> + inline RealType kurtosis(const poisson_distribution<RealType, Policy>& dist) + { // kurtosis is 4th moment about the mean = u4 / sd ^ 4 + // http://en.wikipedia.org/wiki/Curtosis + // kurtosis can range from -2 (flat top) to +infinity (sharp peak & heavy tails). + // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm + return 3 + 1 / dist.mean(); // NIST. + // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess + // is more convenient because the kurtosis excess of a normal distribution is zero + // whereas the true kurtosis is 3. + } // RealType kurtosis + + template <class RealType, class Policy> + RealType pdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) + { // Probability Density/Mass Function. + // Probability that there are EXACTLY k occurrences (or arrivals). + BOOST_FPU_EXCEPTION_GUARD + + BOOST_MATH_STD_USING // for ADL of std functions. + + RealType mean = dist.mean(); + // Error check: + RealType result = 0; + if(false == poisson_detail::check_dist_and_k( + "boost::math::pdf(const poisson_distribution<%1%>&, %1%)", + mean, + k, + &result, Policy())) + { + return result; + } + + // Special case of mean zero, regardless of the number of events k. + if (mean == 0) + { // Probability for any k is zero. + return 0; + } + if (k == 0) + { // mean ^ k = 1, and k! = 1, so can simplify. + return exp(-mean); + } + return boost::math::gamma_p_derivative(k+1, mean, Policy()); + } // pdf + + template <class RealType, class Policy> + RealType cdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) + { // Cumulative Distribution Function Poisson. + // The random variate k is the number of occurrences(or arrivals) + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + // Returns the sum of the terms 0 through k of the Poisson Probability Density or Mass (pdf). + + // But note that the Poisson distribution + // (like others including the binomial, negative binomial & Bernoulli) + // is strictly defined as a discrete function: only integral values of k are envisaged. + // However because of the method of calculation using a continuous gamma function, + // it is convenient to treat it as if it is a continous function + // and permit non-integral values of k. + // To enforce the strict mathematical model, users should use floor or ceil functions + // outside this function to ensure that k is integral. + + // The terms are not summed directly (at least for larger k) + // instead the incomplete gamma integral is employed, + + BOOST_MATH_STD_USING // for ADL of std function exp. + + RealType mean = dist.mean(); + // Error checks: + RealType result = 0; + if(false == poisson_detail::check_dist_and_k( + "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", + mean, + k, + &result, Policy())) + { + return result; + } + // Special cases: + if (mean == 0) + { // Probability for any k is zero. + return 0; + } + if (k == 0) + { // return pdf(dist, static_cast<RealType>(0)); + // but mean (and k) have already been checked, + // so this avoids unnecessary repeated checks. + return exp(-mean); + } + // For small integral k could use a finite sum - + // it's cheaper than the gamma function. + // BUT this is now done efficiently by gamma_q function. + // Calculate poisson cdf using the gamma_q function. + return gamma_q(k+1, mean, Policy()); + } // binomial cdf + + template <class RealType, class Policy> + RealType cdf(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) + { // Complemented Cumulative Distribution Function Poisson + // The random variate k is the number of events, occurrences or arrivals. + // k argument may be integral, signed, or unsigned, or floating point. + // If necessary, it has already been promoted from an integral type. + // But note that the Poisson distribution + // (like others including the binomial, negative binomial & Bernoulli) + // is strictly defined as a discrete function: only integral values of k are envisaged. + // However because of the method of calculation using a continuous gamma function, + // it is convenient to treat it as is it is a continous function + // and permit non-integral values of k. + // To enforce the strict mathematical model, users should use floor or ceil functions + // outside this function to ensure that k is integral. + + // Returns the sum of the terms k+1 through inf of the Poisson Probability Density/Mass (pdf). + // The terms are not summed directly (at least for larger k) + // instead the incomplete gamma integral is employed, + + RealType const& k = c.param; + poisson_distribution<RealType, Policy> const& dist = c.dist; + + RealType mean = dist.mean(); + + // Error checks: + RealType result = 0; + if(false == poisson_detail::check_dist_and_k( + "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", + mean, + k, + &result, Policy())) + { + return result; + } + // Special case of mean, regardless of the number of events k. + if (mean == 0) + { // Probability for any k is unity, complement of zero. + return 1; + } + if (k == 0) + { // Avoid repeated checks on k and mean in gamma_p. + return -boost::math::expm1(-mean, Policy()); + } + // Unlike un-complemented cdf (sum from 0 to k), + // can't use finite sum from k+1 to infinity for small integral k, + // anyway it is now done efficiently by gamma_p. + return gamma_p(k + 1, mean, Policy()); // Calculate Poisson cdf using the gamma_p function. + // CCDF = gamma_p(k+1, lambda) + } // poisson ccdf + + template <class RealType, class Policy> + inline RealType quantile(const poisson_distribution<RealType, Policy>& dist, const RealType& p) + { // Quantile (or Percent Point) Poisson function. + // Return the number of expected events k for a given probability p. + static const char* function = "boost::math::quantile(const poisson_distribution<%1%>&, %1%)"; + RealType result = 0; // of Argument checks: + if(false == poisson_detail::check_prob( + function, + p, + &result, Policy())) + { + return result; + } + // Special case: + if (dist.mean() == 0) + { // if mean = 0 then p = 0, so k can be anything? + if (false == poisson_detail::check_mean_NZ( + function, + dist.mean(), + &result, Policy())) + { + return result; + } + } + if(p == 0) + { + return 0; // Exact result regardless of discrete-quantile Policy + } + if(p == 1) + { + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + } + typedef typename Policy::discrete_quantile_type discrete_type; + boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); + RealType guess, factor = 8; + RealType z = dist.mean(); + if(z < 1) + guess = z; + else + guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, RealType(1-p), Policy()); + if(z > 5) + { + if(z > 1000) + factor = 1.01f; + else if(z > 50) + factor = 1.1f; + else if(guess > 10) + factor = 1.25f; + else + factor = 2; + if(guess < 1.1) + factor = 8; + } + + return detail::inverse_discrete_quantile( + dist, + p, + false, + guess, + factor, + RealType(1), + discrete_type(), + max_iter); + } // quantile + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) + { // Quantile (or Percent Point) of Poisson function. + // Return the number of expected events k for a given + // complement of the probability q. + // + // Error checks: + static const char* function = "boost::math::quantile(complement(const poisson_distribution<%1%>&, %1%))"; + RealType q = c.param; + const poisson_distribution<RealType, Policy>& dist = c.dist; + RealType result = 0; // of argument checks. + if(false == poisson_detail::check_prob( + function, + q, + &result, Policy())) + { + return result; + } + // Special case: + if (dist.mean() == 0) + { // if mean = 0 then p = 0, so k can be anything? + if (false == poisson_detail::check_mean_NZ( + function, + dist.mean(), + &result, Policy())) + { + return result; + } + } + if(q == 0) + { + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + } + if(q == 1) + { + return 0; // Exact result regardless of discrete-quantile Policy + } + typedef typename Policy::discrete_quantile_type discrete_type; + boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); + RealType guess, factor = 8; + RealType z = dist.mean(); + if(z < 1) + guess = z; + else + guess = boost::math::detail::inverse_poisson_cornish_fisher(z, RealType(1-q), q, Policy()); + if(z > 5) + { + if(z > 1000) + factor = 1.01f; + else if(z > 50) + factor = 1.1f; + else if(guess > 10) + factor = 1.25f; + else + factor = 2; + if(guess < 1.1) + factor = 8; + } + + return detail::inverse_discrete_quantile( + dist, + q, + true, + guess, + factor, + RealType(1), + discrete_type(), + max_iter); + } // quantile complement. + + } // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> +#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> + +#endif // BOOST_MATH_SPECIAL_POISSON_HPP + + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/students_t.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/students_t.hpp new file mode 100644 index 00000000000..0d6a6466913 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/students_t.hpp @@ -0,0 +1,490 @@ +// Copyright John Maddock 2006. +// Copyright Paul A. Bristow 2006, 2012. +// Copyright Thomas Mang 2012. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_STUDENTS_T_HPP +#define BOOST_STATS_STUDENTS_T_HPP + +// http://en.wikipedia.org/wiki/Student%27s_t_distribution +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3664.htm + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x). +#include <boost/math/distributions/complement.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/distributions/normal.hpp> + +#include <utility> + +#ifdef BOOST_MSVC +# pragma warning(push) +# pragma warning(disable: 4702) // unreachable code (return after domain_error throw). +#endif + +namespace boost{ namespace math{ + +template <class RealType = double, class Policy = policies::policy<> > +class students_t_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + students_t_distribution(RealType df) : df_(df) + { // Constructor. + RealType result; + detail::check_df_gt0_to_inf( // Checks that df > 0 or df == inf. + "boost::math::students_t_distribution<%1%>::students_t_distribution", df_, &result, Policy()); + } // students_t_distribution + + RealType degrees_of_freedom()const + { + return df_; + } + + // Parameter estimation: + static RealType find_degrees_of_freedom( + RealType difference_from_mean, + RealType alpha, + RealType beta, + RealType sd, + RealType hint = 100); + +private: + // Data member: + RealType df_; // degrees of freedom is a real number or +infinity. +}; + +typedef students_t_distribution<double> students_t; // Convenience typedef for double version. + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const students_t_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x. + // NOT including infinity. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const students_t_distribution<RealType, Policy>& /*dist*/) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_FPU_EXCEPTION_GUARD + BOOST_MATH_STD_USING // for ADL of std functions. + + RealType error_result; + if(false == detail::check_x( + "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy())) + return error_result; + RealType df = dist.degrees_of_freedom(); + if(false == detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity. + "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy())) + return error_result; + + RealType result; + if ((boost::math::isinf)(x)) + { // +infinity. + normal_distribution<RealType, Policy> n(0, 1); + result = pdf(n, x); + return result; + } + RealType limit = policies::get_epsilon<RealType, Policy>(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast<RealType>(1) / limit; // 1/eps + // for 64-bit double 1/eps = 4503599627370496 + if (df > limit) + { // Special case for really big degrees_of_freedom > 1 / eps + // - use normal distribution which is much faster and more accurate. + normal_distribution<RealType, Policy> n(0, 1); + result = pdf(n, x); + } + else + { // + RealType basem1 = x * x / df; + if(basem1 < 0.125) + { + result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2); + } + else + { + result = pow(1 / (1 + basem1), (df + 1) / 2); + } + result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy()); + } + return result; +} // pdf + +template <class RealType, class Policy> +inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x) +{ + RealType error_result; + if(false == detail::check_x( + "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy())) + return error_result; + RealType df = dist.degrees_of_freedom(); + // Error check: + + if(false == detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity. + "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy())) + return error_result; + + if (x == 0) + { // Special case with exact result. + return static_cast<RealType>(0.5); + } + if ((boost::math::isinf)(x)) + { // +infinity. + normal_distribution<RealType, Policy> n(0, 1); + RealType result = cdf(n, x); + return result; + } + RealType limit = policies::get_epsilon<RealType, Policy>(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast<RealType>(1) / limit; // 1/eps + // for 64-bit double 1/eps = 4503599627370496 + if (df > limit) + { // Special case for really big degrees_of_freedom > 1 / eps (perhaps infinite?) + // - use normal distribution which is much faster and more accurate. + normal_distribution<RealType, Policy> n(0, 1); + RealType result = cdf(n, x); + return result; + } + else + { // normal df case. + // + // Calculate probability of Student's t using the incomplete beta function. + // probability = ibeta(degrees_of_freedom / 2, 1/2, degrees_of_freedom / (degrees_of_freedom + t*t)) + // + // However when t is small compared to the degrees of freedom, that formula + // suffers from rounding error, use the identity formula to work around + // the problem: + // + // I[x](a,b) = 1 - I[1-x](b,a) + // + // and: + // + // x = df / (df + t^2) + // + // so: + // + // 1 - x = t^2 / (df + t^2) + // + RealType x2 = x * x; + RealType probability; + if(df > 2 * x2) + { + RealType z = x2 / (df + x2); + probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2; + } + else + { + RealType z = df / (df + x2); + probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2; + } + return (x > 0 ? 1 - probability : probability); + } +} // cdf + +template <class RealType, class Policy> +inline RealType quantile(const students_t_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + // + // Obtain parameters: + RealType probability = p; + + // Check for domain errors: + RealType df = dist.degrees_of_freedom(); + static const char* function = "boost::math::quantile(const students_t_distribution<%1%>&, %1%)"; + RealType error_result; + if(false == (detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity. + function, df, &error_result, Policy()) + && detail::check_probability(function, probability, &error_result, Policy()))) + return error_result; + // Special cases, regardless of degrees_of_freedom. + if (probability == 0) + return -policies::raise_overflow_error<RealType>(function, 0, Policy()); + if (probability == 1) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + if (probability == static_cast<RealType>(0.5)) + return 0; // + // +#if 0 + // This next block is disabled in favour of a faster method than + // incomplete beta inverse, but code retained for future reference: + // + // Calculate quantile of Student's t using the incomplete beta function inverse: + // + probability = (probability > 0.5) ? 1 - probability : probability; + RealType t, x, y; + x = ibeta_inv(degrees_of_freedom / 2, RealType(0.5), 2 * probability, &y); + if(degrees_of_freedom * y > tools::max_value<RealType>() * x) + t = tools::overflow_error<RealType>(function); + else + t = sqrt(degrees_of_freedom * y / x); + // + // Figure out sign based on the size of p: + // + if(p < 0.5) + t = -t; + + return t; +#endif + // + // Depending on how many digits RealType has, this may forward + // to the incomplete beta inverse as above. Otherwise uses a + // faster method that is accurate to ~15 digits everywhere + // and a couple of epsilon at double precision and in the central + // region where most use cases will occur... + // + return boost::math::detail::fast_students_t_quantile(df, probability, Policy()); +} // quantile + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c) +{ + return cdf(c.dist, -c.param); +} + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c) +{ + return -quantile(c.dist, c.param); +} + +// +// Parameter estimation follows: +// +namespace detail{ +// +// Functors for finding degrees of freedom: +// +template <class RealType, class Policy> +struct sample_size_func +{ + sample_size_func(RealType a, RealType b, RealType s, RealType d) + : alpha(a), beta(b), ratio(s*s/(d*d)) {} + + RealType operator()(const RealType& df) + { + if(df <= tools::min_value<RealType>()) + { // + return 1; + } + students_t_distribution<RealType, Policy> t(df); + RealType qa = quantile(complement(t, alpha)); + RealType qb = quantile(complement(t, beta)); + qa += qb; + qa *= qa; + qa *= ratio; + qa -= (df + 1); + return qa; + } + RealType alpha, beta, ratio; +}; + +} // namespace detail + +template <class RealType, class Policy> +RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom( + RealType difference_from_mean, + RealType alpha, + RealType beta, + RealType sd, + RealType hint) +{ + static const char* function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom"; + // + // Check for domain errors: + // + RealType error_result; + if(false == detail::check_probability( + function, alpha, &error_result, Policy()) + && detail::check_probability(function, beta, &error_result, Policy())) + return error_result; + + if(hint <= 0) + hint = 1; + + detail::sample_size_func<RealType, Policy> f(alpha, beta, sd, difference_from_mean); + tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); + boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); + std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy()); + RealType result = r.first + (r.second - r.first) / 2; + if(max_iter >= policies::get_max_root_iterations<Policy>()) + { + return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" + " either there is no answer to how many degrees of freedom are required" + " or the answer is infinite. Current best guess is %1%", result, Policy()); + } + return result; +} + +template <class RealType, class Policy> +inline RealType mode(const students_t_distribution<RealType, Policy>& /*dist*/) +{ + // Assume no checks on degrees of freedom are useful (unlike mean). + return 0; // Always zero by definition. +} + +template <class RealType, class Policy> +inline RealType median(const students_t_distribution<RealType, Policy>& /*dist*/) +{ + // Assume no checks on degrees of freedom are useful (unlike mean). + return 0; // Always zero by definition. +} + +// See section 5.1 on moments at http://en.wikipedia.org/wiki/Student%27s_t-distribution + +template <class RealType, class Policy> +inline RealType mean(const students_t_distribution<RealType, Policy>& dist) +{ // Revised for https://svn.boost.org/trac/boost/ticket/7177 + RealType df = dist.degrees_of_freedom(); + if(((boost::math::isnan)(df)) || (df <= 1) ) + { // mean is undefined for moment <= 1! + return policies::raise_domain_error<RealType>( + "boost::math::mean(students_t_distribution<%1%> const&, %1%)", + "Mean is undefined for degrees of freedom < 1 but got %1%.", df, Policy()); + return std::numeric_limits<RealType>::quiet_NaN(); + } + return 0; +} // mean + +template <class RealType, class Policy> +inline RealType variance(const students_t_distribution<RealType, Policy>& dist) +{ // http://en.wikipedia.org/wiki/Student%27s_t-distribution + // Revised for https://svn.boost.org/trac/boost/ticket/7177 + RealType df = dist.degrees_of_freedom(); + if ((boost::math::isnan)(df) || (df <= 2)) + { // NaN or undefined for <= 2. + return policies::raise_domain_error<RealType>( + "boost::math::variance(students_t_distribution<%1%> const&, %1%)", + "variance is undefined for degrees of freedom <= 2, but got %1%.", + df, Policy()); + return std::numeric_limits<RealType>::quiet_NaN(); // Undefined. + } + if ((boost::math::isinf)(df)) + { // +infinity. + return 1; + } + RealType limit = policies::get_epsilon<RealType, Policy>(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast<RealType>(1) / limit; // 1/eps + // for 64-bit double 1/eps = 4503599627370496 + if (df > limit) + { // Special case for really big degrees_of_freedom > 1 / eps. + return 1; + } + else + { + return df / (df - 2); + } +} // variance + +template <class RealType, class Policy> +inline RealType skewness(const students_t_distribution<RealType, Policy>& dist) +{ + RealType df = dist.degrees_of_freedom(); + if( ((boost::math::isnan)(df)) || (dist.degrees_of_freedom() <= 3)) + { // Undefined for moment k = 3. + return policies::raise_domain_error<RealType>( + "boost::math::skewness(students_t_distribution<%1%> const&, %1%)", + "Skewness is undefined for degrees of freedom <= 3, but got %1%.", + dist.degrees_of_freedom(), Policy()); + return std::numeric_limits<RealType>::quiet_NaN(); + } + return 0; // For all valid df, including infinity. +} // skewness + +template <class RealType, class Policy> +inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist) +{ + RealType df = dist.degrees_of_freedom(); + if(((boost::math::isnan)(df)) || (df <= 4)) + { // Undefined or infinity for moment k = 4. + return policies::raise_domain_error<RealType>( + "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)", + "Kurtosis is undefined for degrees of freedom <= 4, but got %1%.", + df, Policy()); + return std::numeric_limits<RealType>::quiet_NaN(); // Undefined. + } + if ((boost::math::isinf)(df)) + { // +infinity. + return 3; + } + RealType limit = policies::get_epsilon<RealType, Policy>(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast<RealType>(1) / limit; // 1/eps + // for 64-bit double 1/eps = 4503599627370496 + if (df > limit) + { // Special case for really big degrees_of_freedom > 1 / eps. + return 3; + } + else + { + //return 3 * (df - 2) / (df - 4); re-arranged to + return 6 / (df - 4) + 3; + } +} // kurtosis + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist) +{ + // see http://mathworld.wolfram.com/Kurtosis.html + + RealType df = dist.degrees_of_freedom(); + if(((boost::math::isnan)(df)) || (df <= 4)) + { // Undefined or infinity for moment k = 4. + return policies::raise_domain_error<RealType>( + "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)", + "Kurtosis_excess is undefined for degrees of freedom <= 4, but got %1%.", + df, Policy()); + return std::numeric_limits<RealType>::quiet_NaN(); // Undefined. + } + if ((boost::math::isinf)(df)) + { // +infinity. + return 0; + } + RealType limit = policies::get_epsilon<RealType, Policy>(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast<RealType>(1) / limit; // 1/eps + // for 64-bit double 1/eps = 4503599627370496 + if (df > limit) + { // Special case for really big degrees_of_freedom > 1 / eps. + return 0; + } + else + { + return 6 / (df - 4); + } +} + +} // namespace math +} // namespace boost + +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_STUDENTS_T_HPP diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/triangular.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/triangular.hpp new file mode 100644 index 00000000000..78ef0df7447 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/triangular.hpp @@ -0,0 +1,523 @@ +// Copyright John Maddock 2006, 2007. +// Copyright Paul A. Bristow 2006, 2007. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_TRIANGULAR_HPP +#define BOOST_STATS_TRIANGULAR_HPP + +// http://mathworld.wolfram.com/TriangularDistribution.html +// http://en.wikipedia.org/wiki/Triangular_distribution + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/expm1.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/distributions/complement.hpp> +#include <boost/math/constants/constants.hpp> + +#include <utility> + +namespace boost{ namespace math +{ + namespace detail + { + template <class RealType, class Policy> + inline bool check_triangular_lower( + const char* function, + RealType lower, + RealType* result, const Policy& pol) + { + if((boost::math::isfinite)(lower)) + { // Any finite value is OK. + return true; + } + else + { // Not finite: infinity or NaN. + *result = policies::raise_domain_error<RealType>( + function, + "Lower parameter is %1%, but must be finite!", lower, pol); + return false; + } + } // bool check_triangular_lower( + + template <class RealType, class Policy> + inline bool check_triangular_mode( + const char* function, + RealType mode, + RealType* result, const Policy& pol) + { + if((boost::math::isfinite)(mode)) + { // any finite value is OK. + return true; + } + else + { // Not finite: infinity or NaN. + *result = policies::raise_domain_error<RealType>( + function, + "Mode parameter is %1%, but must be finite!", mode, pol); + return false; + } + } // bool check_triangular_mode( + + template <class RealType, class Policy> + inline bool check_triangular_upper( + const char* function, + RealType upper, + RealType* result, const Policy& pol) + { + if((boost::math::isfinite)(upper)) + { // any finite value is OK. + return true; + } + else + { // Not finite: infinity or NaN. + *result = policies::raise_domain_error<RealType>( + function, + "Upper parameter is %1%, but must be finite!", upper, pol); + return false; + } + } // bool check_triangular_upper( + + template <class RealType, class Policy> + inline bool check_triangular_x( + const char* function, + RealType const& x, + RealType* result, const Policy& pol) + { + if((boost::math::isfinite)(x)) + { // Any finite value is OK + return true; + } + else + { // Not finite: infinity or NaN. + *result = policies::raise_domain_error<RealType>( + function, + "x parameter is %1%, but must be finite!", x, pol); + return false; + } + } // bool check_triangular_x + + template <class RealType, class Policy> + inline bool check_triangular( + const char* function, + RealType lower, + RealType mode, + RealType upper, + RealType* result, const Policy& pol) + { + if ((check_triangular_lower(function, lower, result, pol) == false) + || (check_triangular_mode(function, mode, result, pol) == false) + || (check_triangular_upper(function, upper, result, pol) == false)) + { // Some parameter not finite. + return false; + } + else if (lower >= upper) // lower == upper NOT useful. + { // lower >= upper. + *result = policies::raise_domain_error<RealType>( + function, + "lower parameter is %1%, but must be less than upper!", lower, pol); + return false; + } + else + { // Check lower <= mode <= upper. + if (mode < lower) + { + *result = policies::raise_domain_error<RealType>( + function, + "mode parameter is %1%, but must be >= than lower!", lower, pol); + return false; + } + if (mode > upper) + { + *result = policies::raise_domain_error<RealType>( + function, + "mode parameter is %1%, but must be <= than upper!", upper, pol); + return false; + } + return true; // All OK. + } + } // bool check_triangular + } // namespace detail + + template <class RealType = double, class Policy = policies::policy<> > + class triangular_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + triangular_distribution(RealType l_lower = -1, RealType l_mode = 0, RealType l_upper = 1) + : m_lower(l_lower), m_mode(l_mode), m_upper(l_upper) // Constructor. + { // Evans says 'standard triangular' is lower 0, mode 1/2, upper 1, + // has median sqrt(c/2) for c <=1/2 and 1 - sqrt(1-c)/2 for c >= 1/2 + // But this -1, 0, 1 is more useful in most applications to approximate normal distribution, + // where the central value is the most likely and deviations either side equally likely. + RealType result; + detail::check_triangular("boost::math::triangular_distribution<%1%>::triangular_distribution",l_lower, l_mode, l_upper, &result, Policy()); + } + // Accessor functions. + RealType lower()const + { + return m_lower; + } + RealType mode()const + { + return m_mode; + } + RealType upper()const + { + return m_upper; + } + private: + // Data members: + RealType m_lower; // distribution lower aka a + RealType m_mode; // distribution mode aka c + RealType m_upper; // distribution upper aka b + }; // class triangular_distribution + + typedef triangular_distribution<double> triangular; + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const triangular_distribution<RealType, Policy>& /* dist */) + { // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const triangular_distribution<RealType, Policy>& dist) + { // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + return std::pair<RealType, RealType>(dist.lower(), dist.upper()); + } + + template <class RealType, class Policy> + RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) + { + static const char* function = "boost::math::pdf(const triangular_distribution<%1%>&, %1%)"; + RealType lower = dist.lower(); + RealType mode = dist.mode(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_triangular_x(function, x, &result, Policy())) + { + return result; + } + if((x < lower) || (x > upper)) + { + return 0; + } + if (x == lower) + { // (mode - lower) == 0 which would lead to divide by zero! + return (mode == lower) ? 2 / (upper - lower) : RealType(0); + } + else if (x == upper) + { + return (mode == upper) ? 2 / (upper - lower) : RealType(0); + } + else if (x <= mode) + { + return 2 * (x - lower) / ((upper - lower) * (mode - lower)); + } + else + { // (x > mode) + return 2 * (upper - x) / ((upper - lower) * (upper - mode)); + } + } // RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) + + template <class RealType, class Policy> + inline RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) + { + static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)"; + RealType lower = dist.lower(); + RealType mode = dist.mode(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_triangular_x(function, x, &result, Policy())) + { + return result; + } + if((x <= lower)) + { + return 0; + } + if (x >= upper) + { + return 1; + } + // else lower < x < upper + if (x <= mode) + { + return ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower)); + } + else + { + return 1 - (upper - x) * (upper - x) / ((upper - lower) * (upper - mode)); + } + } // RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) + + template <class RealType, class Policy> + RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& p) + { + BOOST_MATH_STD_USING // for ADL of std functions (sqrt). + static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)"; + RealType lower = dist.lower(); + RealType mode = dist.mode(); + RealType upper = dist.upper(); + RealType result = 0; // of checks + if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_probability(function, p, &result, Policy())) + { + return result; + } + if(p == 0) + { + return lower; + } + if(p == 1) + { + return upper; + } + RealType p0 = (mode - lower) / (upper - lower); + RealType q = 1 - p; + if (p < p0) + { + result = sqrt((upper - lower) * (mode - lower) * p) + lower; + } + else if (p == p0) + { + result = mode; + } + else // p > p0 + { + result = upper - sqrt((upper - lower) * (upper - mode) * q); + } + return result; + + } // RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& q) + + template <class RealType, class Policy> + RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) + { + static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)"; + RealType lower = c.dist.lower(); + RealType mode = c.dist.mode(); + RealType upper = c.dist.upper(); + RealType x = c.param; + RealType result = 0; // of checks. + if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_triangular_x(function, x, &result, Policy())) + { + return result; + } + if (x <= lower) + { + return 1; + } + if (x >= upper) + { + return 0; + } + if (x <= mode) + { + return 1 - ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower)); + } + else + { + return (upper - x) * (upper - x) / ((upper - lower) * (upper - mode)); + } + } // RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) + + template <class RealType, class Policy> + RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) + { + BOOST_MATH_STD_USING // Aid ADL for sqrt. + static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)"; + RealType l = c.dist.lower(); + RealType m = c.dist.mode(); + RealType u = c.dist.upper(); + RealType q = c.param; // probability 0 to 1. + RealType result = 0; // of checks. + if(false == detail::check_triangular(function, l, m, u, &result, Policy())) + { + return result; + } + if(false == detail::check_probability(function, q, &result, Policy())) + { + return result; + } + if(q == 0) + { + return u; + } + if(q == 1) + { + return l; + } + RealType lower = c.dist.lower(); + RealType mode = c.dist.mode(); + RealType upper = c.dist.upper(); + + RealType p = 1 - q; + RealType p0 = (mode - lower) / (upper - lower); + if(p < p0) + { + RealType s = (upper - lower) * (mode - lower); + s *= p; + result = sqrt((upper - lower) * (mode - lower) * p) + lower; + } + else if (p == p0) + { + result = mode; + } + else // p > p0 + { + result = upper - sqrt((upper - lower) * (upper - mode) * q); + } + return result; + } // RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) + + template <class RealType, class Policy> + inline RealType mean(const triangular_distribution<RealType, Policy>& dist) + { + static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)"; + RealType lower = dist.lower(); + RealType mode = dist.mode(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) + { + return result; + } + return (lower + upper + mode) / 3; + } // RealType mean(const triangular_distribution<RealType, Policy>& dist) + + + template <class RealType, class Policy> + inline RealType variance(const triangular_distribution<RealType, Policy>& dist) + { + static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)"; + RealType lower = dist.lower(); + RealType mode = dist.mode(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) + { + return result; + } + return (lower * lower + upper * upper + mode * mode - lower * upper - lower * mode - upper * mode) / 18; + } // RealType variance(const triangular_distribution<RealType, Policy>& dist) + + template <class RealType, class Policy> + inline RealType mode(const triangular_distribution<RealType, Policy>& dist) + { + static const char* function = "boost::math::mode(const triangular_distribution<%1%>&)"; + RealType mode = dist.mode(); + RealType result = 0; // of checks. + if(false == detail::check_triangular_mode(function, mode, &result, Policy())) + { // This should never happen! + return result; + } + return mode; + } // RealType mode + + template <class RealType, class Policy> + inline RealType median(const triangular_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // ADL of std functions. + static const char* function = "boost::math::median(const triangular_distribution<%1%>&)"; + RealType mode = dist.mode(); + RealType result = 0; // of checks. + if(false == detail::check_triangular_mode(function, mode, &result, Policy())) + { // This should never happen! + return result; + } + RealType lower = dist.lower(); + RealType upper = dist.upper(); + if (mode < (upper - lower) / 2) + { + return lower + sqrt((upper - lower) * (mode - lower)) / constants::root_two<RealType>(); + } + else + { + return upper - sqrt((upper - lower) * (upper - mode)) / constants::root_two<RealType>(); + } + } // RealType mode + + template <class RealType, class Policy> + inline RealType skewness(const triangular_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // for ADL of std functions + using namespace boost::math::constants; // for root_two + static const char* function = "boost::math::skewness(const triangular_distribution<%1%>&)"; + + RealType lower = dist.lower(); + RealType mode = dist.mode(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == boost::math::detail::check_triangular(function,lower, mode, upper, &result, Policy())) + { + return result; + } + return root_two<RealType>() * (lower + upper - 2 * mode) * (2 * lower - upper - mode) * (lower - 2 * upper + mode) / + (5 * pow((lower * lower + upper + upper + mode * mode - lower * upper - lower * mode - upper * mode), RealType(3)/RealType(2))); + } // RealType skewness(const triangular_distribution<RealType, Policy>& dist) + + template <class RealType, class Policy> + inline RealType kurtosis(const triangular_distribution<RealType, Policy>& dist) + { // These checks may be belt and braces as should have been checked on construction? + static const char* function = "boost::math::kurtosis(const triangular_distribution<%1%>&)"; + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType mode = dist.mode(); + RealType result = 0; // of checks. + if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) + { + return result; + } + return static_cast<RealType>(12)/5; // 12/5 = 2.4; + } // RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist) + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist) + { // These checks may be belt and braces as should have been checked on construction? + static const char* function = "boost::math::kurtosis_excess(const triangular_distribution<%1%>&)"; + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType mode = dist.mode(); + RealType result = 0; // of checks. + if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) + { + return result; + } + return static_cast<RealType>(-3)/5; // - 3/5 = -0.6 + // Assuming mathworld really means kurtosis excess? Wikipedia now corrected to match this. + } + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_TRIANGULAR_HPP + + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/uniform.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/uniform.hpp new file mode 100644 index 00000000000..a20597a66ad --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/uniform.hpp @@ -0,0 +1,379 @@ +// Copyright John Maddock 2006. +// Copyright Paul A. Bristow 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +// TODO deal with infinity as special better - or remove. +// + +#ifndef BOOST_STATS_UNIFORM_HPP +#define BOOST_STATS_UNIFORM_HPP + +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm +// http://mathworld.wolfram.com/UniformDistribution.html +// http://documents.wolfram.com/calculationcenter/v2/Functions/ListsMatrices/Statistics/UniformDistribution.html +// http://en.wikipedia.org/wiki/Uniform_distribution_%28continuous%29 + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/distributions/complement.hpp> + +#include <utility> + +namespace boost{ namespace math +{ + namespace detail + { + template <class RealType, class Policy> + inline bool check_uniform_lower( + const char* function, + RealType lower, + RealType* result, const Policy& pol) + { + if((boost::math::isfinite)(lower)) + { // any finite value is OK. + return true; + } + else + { // Not finite. + *result = policies::raise_domain_error<RealType>( + function, + "Lower parameter is %1%, but must be finite!", lower, pol); + return false; + } + } // bool check_uniform_lower( + + template <class RealType, class Policy> + inline bool check_uniform_upper( + const char* function, + RealType upper, + RealType* result, const Policy& pol) + { + if((boost::math::isfinite)(upper)) + { // Any finite value is OK. + return true; + } + else + { // Not finite. + *result = policies::raise_domain_error<RealType>( + function, + "Upper parameter is %1%, but must be finite!", upper, pol); + return false; + } + } // bool check_uniform_upper( + + template <class RealType, class Policy> + inline bool check_uniform_x( + const char* function, + RealType const& x, + RealType* result, const Policy& pol) + { + if((boost::math::isfinite)(x)) + { // Any finite value is OK + return true; + } + else + { // Not finite.. + *result = policies::raise_domain_error<RealType>( + function, + "x parameter is %1%, but must be finite!", x, pol); + return false; + } + } // bool check_uniform_x + + template <class RealType, class Policy> + inline bool check_uniform( + const char* function, + RealType lower, + RealType upper, + RealType* result, const Policy& pol) + { + if((check_uniform_lower(function, lower, result, pol) == false) + || (check_uniform_upper(function, upper, result, pol) == false)) + { + return false; + } + else if (lower >= upper) // If lower == upper then 1 / (upper-lower) = 1/0 = +infinity! + { // upper and lower have been checked before, so must be lower >= upper. + *result = policies::raise_domain_error<RealType>( + function, + "lower parameter is %1%, but must be less than upper!", lower, pol); + return false; + } + else + { // All OK, + return true; + } + } // bool check_uniform( + + } // namespace detail + + template <class RealType = double, class Policy = policies::policy<> > + class uniform_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + uniform_distribution(RealType l_lower = 0, RealType l_upper = 1) // Constructor. + : m_lower(l_lower), m_upper(l_upper) // Default is standard uniform distribution. + { + RealType result; + detail::check_uniform("boost::math::uniform_distribution<%1%>::uniform_distribution", l_lower, l_upper, &result, Policy()); + } + // Accessor functions. + RealType lower()const + { + return m_lower; + } + + RealType upper()const + { + return m_upper; + } + private: + // Data members: + RealType m_lower; // distribution lower aka a. + RealType m_upper; // distribution upper aka b. + }; // class uniform_distribution + + typedef uniform_distribution<double> uniform; + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const uniform_distribution<RealType, Policy>& /* dist */) + { // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + 'infinity'. + // Note RealType infinity is NOT permitted, only max_value. + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const uniform_distribution<RealType, Policy>& dist) + { // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(dist.lower(), dist.upper()); + } + + template <class RealType, class Policy> + inline RealType pdf(const uniform_distribution<RealType, Policy>& dist, const RealType& x) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::pdf(const uniform_distribution<%1%>&, %1%)", lower, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_uniform_x("boost::math::pdf(const uniform_distribution<%1%>&, %1%)", x, &result, Policy())) + { + return result; + } + + if((x < lower) || (x > upper) ) + { + return 0; + } + else + { + return 1 / (upper - lower); + } + } // RealType pdf(const uniform_distribution<RealType, Policy>& dist, const RealType& x) + + template <class RealType, class Policy> + inline RealType cdf(const uniform_distribution<RealType, Policy>& dist, const RealType& x) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::cdf(const uniform_distribution<%1%>&, %1%)",lower, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_uniform_x("boost::math::cdf(const uniform_distribution<%1%>&, %1%)", x, &result, Policy())) + { + return result; + } + if (x < lower) + { + return 0; + } + if (x > upper) + { + return 1; + } + return (x - lower) / (upper - lower); // lower <= x <= upper + } // RealType cdf(const uniform_distribution<RealType, Policy>& dist, const RealType& x) + + template <class RealType, class Policy> + inline RealType quantile(const uniform_distribution<RealType, Policy>& dist, const RealType& p) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks + if(false == detail::check_uniform("boost::math::quantile(const uniform_distribution<%1%>&, %1%)",lower, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_probability("boost::math::quantile(const uniform_distribution<%1%>&, %1%)", p, &result, Policy())) + { + return result; + } + if(p == 0) + { + return lower; + } + if(p == 1) + { + return upper; + } + return p * (upper - lower) + lower; + } // RealType quantile(const uniform_distribution<RealType, Policy>& dist, const RealType& p) + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<uniform_distribution<RealType, Policy>, RealType>& c) + { + RealType lower = c.dist.lower(); + RealType upper = c.dist.upper(); + RealType x = c.param; + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::cdf(const uniform_distribution<%1%>&, %1%)", lower, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_uniform_x("boost::math::cdf(const uniform_distribution<%1%>&, %1%)", x, &result, Policy())) + { + return result; + } + if (x < lower) + { + return 1; + } + if (x > upper) + { + return 0; + } + return (upper - x) / (upper - lower); + } // RealType cdf(const complemented2_type<uniform_distribution<RealType, Policy>, RealType>& c) + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<uniform_distribution<RealType, Policy>, RealType>& c) + { + RealType lower = c.dist.lower(); + RealType upper = c.dist.upper(); + RealType q = c.param; + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::quantile(const uniform_distribution<%1%>&, %1%)", lower, upper, &result, Policy())) + { + return result; + } + if(false == detail::check_probability("boost::math::quantile(const uniform_distribution<%1%>&, %1%)", q, &result, Policy())) + if(q == 0) + { + return lower; + } + if(q == 1) + { + return upper; + } + return -q * (upper - lower) + upper; + } // RealType quantile(const complemented2_type<uniform_distribution<RealType, Policy>, RealType>& c) + + template <class RealType, class Policy> + inline RealType mean(const uniform_distribution<RealType, Policy>& dist) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::mean(const uniform_distribution<%1%>&)", lower, upper, &result, Policy())) + { + return result; + } + return (lower + upper ) / 2; + } // RealType mean(const uniform_distribution<RealType, Policy>& dist) + + template <class RealType, class Policy> + inline RealType variance(const uniform_distribution<RealType, Policy>& dist) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::variance(const uniform_distribution<%1%>&)", lower, upper, &result, Policy())) + { + return result; + } + return (upper - lower) * ( upper - lower) / 12; + // for standard uniform = 0.833333333333333333333333333333333333333333; + } // RealType variance(const uniform_distribution<RealType, Policy>& dist) + + template <class RealType, class Policy> + inline RealType mode(const uniform_distribution<RealType, Policy>& dist) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::mode(const uniform_distribution<%1%>&)", lower, upper, &result, Policy())) + { + return result; + } + result = lower; // Any value [lower, upper] but arbitrarily choose lower. + return result; + } + + template <class RealType, class Policy> + inline RealType median(const uniform_distribution<RealType, Policy>& dist) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::median(const uniform_distribution<%1%>&)", lower, upper, &result, Policy())) + { + return result; + } + return (lower + upper) / 2; // + } + template <class RealType, class Policy> + inline RealType skewness(const uniform_distribution<RealType, Policy>& dist) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::skewness(const uniform_distribution<%1%>&)",lower, upper, &result, Policy())) + { + return result; + } + return 0; + } // RealType skewness(const uniform_distribution<RealType, Policy>& dist) + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const uniform_distribution<RealType, Policy>& dist) + { + RealType lower = dist.lower(); + RealType upper = dist.upper(); + RealType result = 0; // of checks. + if(false == detail::check_uniform("boost::math::kurtosis_execess(const uniform_distribution<%1%>&)", lower, upper, &result, Policy())) + { + return result; + } + return static_cast<RealType>(-6)/5; // -6/5 = -1.2; + } // RealType kurtosis_excess(const uniform_distribution<RealType, Policy>& dist) + + template <class RealType, class Policy> + inline RealType kurtosis(const uniform_distribution<RealType, Policy>& dist) + { + return kurtosis_excess(dist) + 3; + } + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_UNIFORM_HPP + + + diff --git a/src/third_party/boost-1.56.0/boost/math/distributions/weibull.hpp b/src/third_party/boost-1.56.0/boost/math/distributions/weibull.hpp new file mode 100644 index 00000000000..da1189090cb --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/distributions/weibull.hpp @@ -0,0 +1,395 @@ +// Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_STATS_WEIBULL_HPP +#define BOOST_STATS_WEIBULL_HPP + +// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm +// http://mathworld.wolfram.com/WeibullDistribution.html + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/special_functions/gamma.hpp> +#include <boost/math/special_functions/log1p.hpp> +#include <boost/math/special_functions/expm1.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/distributions/complement.hpp> + +#include <utility> + +namespace boost{ namespace math +{ +namespace detail{ + +template <class RealType, class Policy> +inline bool check_weibull_shape( + const char* function, + RealType shape, + RealType* result, const Policy& pol) +{ + if((shape <= 0) || !(boost::math::isfinite)(shape)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Shape parameter is %1%, but must be > 0 !", shape, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_weibull_x( + const char* function, + RealType const& x, + RealType* result, const Policy& pol) +{ + if((x < 0) || !(boost::math::isfinite)(x)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Random variate is %1% but must be >= 0 !", x, pol); + return false; + } + return true; +} + +template <class RealType, class Policy> +inline bool check_weibull( + const char* function, + RealType scale, + RealType shape, + RealType* result, const Policy& pol) +{ + return check_scale(function, scale, result, pol) && check_weibull_shape(function, shape, result, pol); +} + +} // namespace detail + +template <class RealType = double, class Policy = policies::policy<> > +class weibull_distribution +{ +public: + typedef RealType value_type; + typedef Policy policy_type; + + weibull_distribution(RealType l_shape, RealType l_scale = 1) + : m_shape(l_shape), m_scale(l_scale) + { + RealType result; + detail::check_weibull("boost::math::weibull_distribution<%1%>::weibull_distribution", l_scale, l_shape, &result, Policy()); + } + + RealType shape()const + { + return m_shape; + } + + RealType scale()const + { + return m_scale; + } +private: + // + // Data members: + // + RealType m_shape; // distribution shape + RealType m_scale; // distribution scale +}; + +typedef weibull_distribution<double> weibull; + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> range(const weibull_distribution<RealType, Policy>& /*dist*/) +{ // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); +} + +template <class RealType, class Policy> +inline const std::pair<RealType, RealType> support(const weibull_distribution<RealType, Policy>& /*dist*/) +{ // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + using boost::math::tools::max_value; + using boost::math::tools::min_value; + return std::pair<RealType, RealType>(min_value<RealType>(), max_value<RealType>()); + // A discontinuity at x == 0, so only support down to min_value. +} + +template <class RealType, class Policy> +inline RealType pdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::pdf(const weibull_distribution<%1%>, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_weibull_x(function, x, &result, Policy())) + return result; + + if(x == 0) + { + if(shape == 1) + { + return 1 / scale; + } + if(shape > 1) + { + return 0; + } + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + } + result = exp(-pow(x / scale, shape)); + result *= pow(x / scale, shape - 1) * shape / scale; + + return result; +} + +template <class RealType, class Policy> +inline RealType cdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_weibull_x(function, x, &result, Policy())) + return result; + + result = -boost::math::expm1(-pow(x / scale, shape), Policy()); + + return result; +} + +template <class RealType, class Policy> +inline RealType quantile(const weibull_distribution<RealType, Policy>& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_probability(function, p, &result, Policy())) + return result; + + if(p == 1) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = scale * pow(-boost::math::log1p(-p, Policy()), 1 / shape); + + return result; +} + +template <class RealType, class Policy> +inline RealType cdf(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)"; + + RealType shape = c.dist.shape(); + RealType scale = c.dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_weibull_x(function, c.param, &result, Policy())) + return result; + + result = exp(-pow(c.param / scale, shape)); + + return result; +} + +template <class RealType, class Policy> +inline RealType quantile(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)"; + + RealType shape = c.dist.shape(); + RealType scale = c.dist.scale(); + RealType q = c.param; + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_probability(function, q, &result, Policy())) + return result; + + if(q == 0) + return policies::raise_overflow_error<RealType>(function, 0, Policy()); + + result = scale * pow(-log(q), 1 / shape); + + return result; +} + +template <class RealType, class Policy> +inline RealType mean(const weibull_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::mean(const weibull_distribution<%1%>)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + + result = scale * boost::math::tgamma(1 + 1 / shape, Policy()); + return result; +} + +template <class RealType, class Policy> +inline RealType variance(const weibull_distribution<RealType, Policy>& dist) +{ + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + static const char* function = "boost::math::variance(const weibull_distribution<%1%>)"; + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + { + return result; + } + result = boost::math::tgamma(1 + 1 / shape, Policy()); + result *= -result; + result += boost::math::tgamma(1 + 2 / shape, Policy()); + result *= scale * scale; + return result; +} + +template <class RealType, class Policy> +inline RealType mode(const weibull_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std function pow. + + static const char* function = "boost::math::mode(const weibull_distribution<%1%>)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + { + return result; + } + if(shape <= 1) + return 0; + result = scale * pow((shape - 1) / shape, 1 / shape); + return result; +} + +template <class RealType, class Policy> +inline RealType median(const weibull_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std function pow. + + static const char* function = "boost::math::median(const weibull_distribution<%1%>)"; + + RealType shape = dist.shape(); // Wikipedia k + RealType scale = dist.scale(); // Wikipedia lambda + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + { + return result; + } + using boost::math::constants::ln_two; + result = scale * pow(ln_two<RealType>(), 1 / shape); + return result; +} + +template <class RealType, class Policy> +inline RealType skewness(const weibull_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::skewness(const weibull_distribution<%1%>)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + { + return result; + } + RealType g1, g2, g3, d; + + g1 = boost::math::tgamma(1 + 1 / shape, Policy()); + g2 = boost::math::tgamma(1 + 2 / shape, Policy()); + g3 = boost::math::tgamma(1 + 3 / shape, Policy()); + d = pow(g2 - g1 * g1, RealType(1.5)); + + result = (2 * g1 * g1 * g1 - 3 * g1 * g2 + g3) / d; + return result; +} + +template <class RealType, class Policy> +inline RealType kurtosis_excess(const weibull_distribution<RealType, Policy>& dist) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::kurtosis_excess(const weibull_distribution<%1%>)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + + RealType g1, g2, g3, g4, d, g1_2, g1_4; + + g1 = boost::math::tgamma(1 + 1 / shape, Policy()); + g2 = boost::math::tgamma(1 + 2 / shape, Policy()); + g3 = boost::math::tgamma(1 + 3 / shape, Policy()); + g4 = boost::math::tgamma(1 + 4 / shape, Policy()); + g1_2 = g1 * g1; + g1_4 = g1_2 * g1_2; + d = g2 - g1_2; + d *= d; + + result = -6 * g1_4 + 12 * g1_2 * g2 - 3 * g2 * g2 - 4 * g1 * g3 + g4; + result /= d; + return result; +} + +template <class RealType, class Policy> +inline RealType kurtosis(const weibull_distribution<RealType, Policy>& dist) +{ + return kurtosis_excess(dist) + 3; +} + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_WEIBULL_HPP + + |