diff --git a/.gitignore b/.gitignore index 901283a3f7..5ff21efb64 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,6 @@ Makefile **CTestTestfile.cmake DartConfiguration.tcl cmake-build-debug/* +.cmake/* +build.ninja +.ninja* diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 5cb2c05f9b..a8fcbbc05f 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -29,6 +29,7 @@ #ifndef BOOST_MATH_DIST_ARCSINE_HPP #define BOOST_MATH_DIST_ARCSINE_HPP +#include #include #include // complements. #include // error checks. diff --git a/include/boost/math/distributions/chi_squared.hpp b/include/boost/math/distributions/chi_squared.hpp index e97bee7ce7..7b65a0da68 100644 --- a/include/boost/math/distributions/chi_squared.hpp +++ b/include/boost/math/distributions/chi_squared.hpp @@ -23,10 +23,10 @@ template > class chi_squared_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - chi_squared_distribution(RealType i) : m_df(i) + explicit chi_squared_distribution(RealType i) : m_df(i) { RealType result; detail::check_df( @@ -53,7 +53,7 @@ class chi_squared_distribution RealType m_df; // degrees of freedom is a positive real number. }; // class chi_squared_distribution -typedef chi_squared_distribution chi_squared; +using chi_squared = chi_squared_distribution; #ifdef __cpp_deduction_guides template @@ -66,7 +66,7 @@ chi_squared_distribution(RealType)->chi_squared_distribution -inline const std::pair range(const chi_squared_distribution& /*dist*/) +inline std::pair range(const chi_squared_distribution& /*dist*/) { // Range of permissible values for random variable x. if (std::numeric_limits::has_infinity) { @@ -84,7 +84,7 @@ inline const std::pair range(const chi_squared_distribution< #endif template -inline const std::pair support(const chi_squared_distribution& /*dist*/) +inline std::pair support(const chi_squared_distribution& /*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(static_cast(0), tools::max_value()); // 0 to + infinity. @@ -224,17 +224,6 @@ inline RealType mode(const chi_squared_distribution& 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( - // 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( @@ -244,25 +233,12 @@ inline RealType mode(const chi_squared_distribution& dist) return df - 2; } -//template -//inline RealType median(const chi_squared_distribution& dist) -//{ // Median is given by Quantile[dist, 1/2] -// RealType df = dist.degrees_of_freedom(); -// if(df <= 1) -// return tools::domain_error( -// 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 inline RealType skewness(const chi_squared_distribution& dist) { BOOST_MATH_STD_USING // For ADL RealType df = dist.degrees_of_freedom(); - return sqrt (8 / df); // == 2 * sqrt(2 / df); + return sqrt (8 / df); } template diff --git a/include/boost/math/distributions/detail/derived_accessors.hpp b/include/boost/math/distributions/detail/derived_accessors.hpp index 9101c8ae1f..e2eca511bf 100644 --- a/include/boost/math/distributions/detail/derived_accessors.hpp +++ b/include/boost/math/distributions/detail/derived_accessors.hpp @@ -110,6 +110,13 @@ inline typename Distribution::value_type pdf(const Distribution& dist, const Rea return pdf(dist, static_cast(x)); } template +inline typename Distribution::value_type logpdf(const Distribution& dist, const RealType& x) +{ + using std::log; + typedef typename Distribution::value_type value_type; + return log(pdf(dist, static_cast(x))); +} +template inline typename Distribution::value_type cdf(const Distribution& dist, const RealType& x) { typedef typename Distribution::value_type value_type; diff --git a/include/boost/math/distributions/exponential.hpp b/include/boost/math/distributions/exponential.hpp index ba7eae927f..5214575a64 100644 --- a/include/boost/math/distributions/exponential.hpp +++ b/include/boost/math/distributions/exponential.hpp @@ -60,10 +60,10 @@ template > class exponential_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - exponential_distribution(RealType l_lambda = 1) + explicit exponential_distribution(RealType l_lambda = 1) : m_lambda(l_lambda) { RealType err; @@ -76,7 +76,7 @@ class exponential_distribution RealType m_lambda; }; -typedef exponential_distribution exponential; +using exponential = exponential_distribution; #ifdef __cpp_deduction_guides template @@ -84,7 +84,7 @@ exponential_distribution(RealType)->exponential_distribution -inline const std::pair range(const exponential_distribution& /*dist*/) +inline std::pair range(const exponential_distribution& /*dist*/) { // Range of permissible values for random variable x. if (std::numeric_limits::has_infinity) { @@ -98,7 +98,7 @@ inline const std::pair range(const exponential_distribution< } template -inline const std::pair support(const exponential_distribution& /*dist*/) +inline std::pair support(const exponential_distribution& /*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; @@ -127,6 +127,24 @@ inline RealType pdf(const exponential_distribution& dist, cons return result; } // pdf +template +inline RealType logpdf(const exponential_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::logpdf(const exponential_distribution<%1%>&, %1%)"; + + RealType lambda = dist.lambda(); + RealType result = -std::numeric_limits::infinity(); + if(0 == detail::verify_lambda(function, lambda, &result, Policy())) + return result; + if(0 == detail::verify_exp_x(function, x, &result, Policy())) + return result; + + result = log(lambda) - lambda * x; + return result; +} // logpdf + template inline RealType cdf(const exponential_distribution& dist, const RealType& x) { diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp index d503b31bc9..b2e8bea966 100644 --- a/include/boost/math/distributions/extreme_value.hpp +++ b/include/boost/math/distributions/extreme_value.hpp @@ -53,10 +53,10 @@ template > class extreme_value_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - extreme_value_distribution(RealType a = 0, RealType b = 1) + explicit extreme_value_distribution(RealType a = 0, RealType b = 1) : m_a(a), m_b(b) { RealType err; @@ -68,10 +68,11 @@ class extreme_value_distribution RealType scale()const { return m_b; } private: - RealType m_a, m_b; + RealType m_a; + RealType m_b; }; -typedef extreme_value_distribution extreme_value; +using extreme_value = extreme_value_distribution; #ifdef __cpp_deduction_guides template @@ -81,7 +82,7 @@ extreme_value_distribution(RealType,RealType)->extreme_value_distribution -inline const std::pair range(const extreme_value_distribution& /*dist*/) +inline std::pair range(const extreme_value_distribution& /*dist*/) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair( @@ -90,7 +91,7 @@ inline const std::pair range(const extreme_value_distributio } template -inline const std::pair support(const extreme_value_distribution& /*dist*/) +inline std::pair support(const extreme_value_distribution& /*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; @@ -122,6 +123,31 @@ inline RealType pdf(const extreme_value_distribution& dist, co return result; } // pdf +template +inline RealType logpdf(const extreme_value_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::logpdf(const extreme_value_distribution<%1%>&, %1%)"; + + RealType a = dist.location(); + RealType b = dist.scale(); + RealType result = -std::numeric_limits::infinity(); + if(0 == detail::verify_scale_b(function, b, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if((boost::math::isinf)(x)) + return 0.0f; + if(0 == detail::check_x(function, x, &result, Policy())) + return result; + RealType e = (a - x) / b; + if(e < tools::log_max_value()) + result = log(1/b) + e - exp(e); + // else.... result *must* be zero since exp(e) is infinite... + return result; +} // logpdf + template inline RealType cdf(const extreme_value_distribution& dist, const RealType& x) { diff --git a/include/boost/math/distributions/gamma.hpp b/include/boost/math/distributions/gamma.hpp index 8a3414d312..28b7c55b0b 100644 --- a/include/boost/math/distributions/gamma.hpp +++ b/include/boost/math/distributions/gamma.hpp @@ -17,6 +17,7 @@ #include #include +#include namespace boost{ namespace math { @@ -71,10 +72,10 @@ template > class gamma_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - gamma_distribution(RealType l_shape, RealType l_scale = 1) + explicit gamma_distribution(RealType l_shape, RealType l_scale = 1) : m_shape(l_shape), m_scale(l_scale) { RealType result; @@ -108,14 +109,14 @@ gamma_distribution(RealType,RealType)->gamma_distribution -inline const std::pair range(const gamma_distribution& /* dist */) +inline std::pair range(const gamma_distribution& /* dist */) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); } template -inline const std::pair support(const gamma_distribution& /* dist */) +inline std::pair support(const gamma_distribution& /* 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; @@ -147,6 +148,33 @@ inline RealType pdf(const gamma_distribution& dist, const Real return result; } // pdf +template +inline RealType logpdf(const gamma_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + using boost::math::lgamma; + + static const char* function = "boost::math::logpdf(const gamma_distribution<%1%>&, %1%)"; + + RealType k = dist.shape(); + RealType theta = dist.scale(); + + RealType result = -std::numeric_limits::infinity(); + if(false == detail::check_gamma(function, theta, k, &result, Policy())) + return result; + if(false == detail::check_gamma_x(function, x, &result, Policy())) + return result; + + if(x == 0) + { + return std::numeric_limits::quiet_NaN(); + } + + result = -k*log(theta) + (k-1)*log(x) - lgamma(k) - (x/theta); + + return result; +} // logpdf + template inline RealType cdf(const gamma_distribution& dist, const RealType& x) { diff --git a/include/boost/math/distributions/inverse_gamma.hpp b/include/boost/math/distributions/inverse_gamma.hpp index 9266bc22f6..8c9e4763d5 100644 --- a/include/boost/math/distributions/inverse_gamma.hpp +++ b/include/boost/math/distributions/inverse_gamma.hpp @@ -28,6 +28,7 @@ #include #include +#include namespace boost{ namespace math { @@ -88,10 +89,10 @@ template > class inverse_gamma_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1) + explicit inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1) : m_shape(l_shape), m_scale(l_scale) { RealType result; @@ -117,10 +118,9 @@ class inverse_gamma_distribution RealType m_scale; // distribution scale }; -typedef inverse_gamma_distribution inverse_gamma; +using inverse_gamma = inverse_gamma_distribution; // typedef - but potential clash with name of inverse gamma *function*. -// but there is a typedef for gamma -// typedef boost::math::gamma_distribution gamma; +// but there is a typedef for the gamma distribution (gamma) #ifdef __cpp_deduction_guides template @@ -132,14 +132,14 @@ inverse_gamma_distribution(RealType,RealType)->inverse_gamma_distribution -inline const std::pair range(const inverse_gamma_distribution& /* dist */) +inline std::pair range(const inverse_gamma_distribution& /* dist */) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); } template -inline const std::pair support(const inverse_gamma_distribution& /* dist */) +inline std::pair support(const inverse_gamma_distribution& /* 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; @@ -190,11 +190,47 @@ inline RealType pdf(const inverse_gamma_distribution& dist, co } result /= (x * x); } - // better than naive - // result = (pow(scale, shape) * pow(x, (-shape -1)) * exp(-scale/x) ) / tgamma(shape); + return result; } // pdf +template +inline RealType logpdf(const inverse_gamma_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + using boost::math::lgamma; + + static const char* function = "boost::math::logpdf(const inverse_gamma_distribution<%1%>&, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = -std::numeric_limits::infinity(); + if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) + { // distribution parameters bad. + return result; + } + if(x == 0) + { // Treat random variate zero as a special case. + return result; + } + else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy())) + { // x bad. + return result; + } + result = scale / x; + if(result < tools::min_value()) + return result; // random variable is infinite or so close as to make no difference. + + // x * x may under or overflow, likewise our result + if (!(boost::math::isfinite)(x*x)) + { + return policies::raise_overflow_error(function, "PDF is infinite.", Policy()); + } + + return shape * log(scale) + (-shape-1)*log(x) - lgamma(shape) - (scale/x); +} // pdf + template inline RealType cdf(const inverse_gamma_distribution& dist, const RealType& x) { @@ -269,7 +305,6 @@ inline RealType cdf(const complemented2_type #include #include // for gamma function -// using boost::math::gamma_p; #include -//using std::tr1::tuple; -//using std::tr1::make_tuple; #include -//using boost::math::tools::newton_raphson_iterate; #include @@ -71,10 +67,10 @@ template > class inverse_gaussian_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - inverse_gaussian_distribution(RealType l_mean = 1, RealType l_scale = 1) + explicit inverse_gaussian_distribution(RealType l_mean = 1, RealType l_scale = 1) : m_mean(l_mean), m_scale(l_scale) { // Default is a 1,1 inverse_gaussian distribution. static const char* function = "boost::math::inverse_gaussian_distribution<%1%>::inverse_gaussian_distribution"; @@ -113,7 +109,7 @@ class inverse_gaussian_distribution RealType m_scale; // distribution standard deviation or scale, aka lambda. }; // class normal_distribution -typedef inverse_gaussian_distribution inverse_gaussian; +using inverse_gaussian = inverse_gaussian_distribution; #ifdef __cpp_deduction_guides template @@ -123,14 +119,14 @@ inverse_gaussian_distribution(RealType,RealType)->inverse_gaussian_distribution< #endif template -inline const std::pair range(const inverse_gaussian_distribution& /*dist*/) +inline std::pair range(const inverse_gaussian_distribution& /*dist*/) { // Range of permissible values for random variable x, zero to max. using boost::math::tools::max_value; return std::pair(static_cast(0.), max_value()); // - to + max value. } template -inline const std::pair support(const inverse_gaussian_distribution& /*dist*/) +inline std::pair support(const inverse_gaussian_distribution& /*dist*/) { // Range of supported values for random variable x, zero to max. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -174,6 +170,43 @@ inline RealType pdf(const inverse_gaussian_distribution& dist, return result; } // pdf +template +inline RealType logpdf(const inverse_gaussian_distribution& dist, const RealType& x) +{ // Probability Density Function + BOOST_MATH_STD_USING // for ADL of std functions + + RealType scale = dist.scale(); + RealType mean = dist.mean(); + RealType result = -std::numeric_limits::infinity(); + static const char* function = "boost::math::logpdf(const inverse_gaussian_distribution<%1%>&, %1%)"; + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, mean, &result, Policy())) + { + return result; + } + if(false == detail::check_x_gt0(function, mean, &result, Policy())) + { + return result; + } + if(false == detail::check_positive_x(function, x, &result, Policy())) + { + return result; + } + + if (x == 0) + { + return std::numeric_limits::quiet_NaN(); // Convenient, even if not defined mathematically. log(0) + } + + const RealType two_pi = boost::math::constants::two_pi(); + + result = (-scale*pow(mean - x, RealType(2))/(mean*mean*x) + log(scale) - 3*log(x) - log(two_pi)) / 2; + return result; +} // pdf + template inline RealType cdf(const inverse_gaussian_distribution& dist, const RealType& x) { // Cumulative Density Function. @@ -203,10 +236,7 @@ inline RealType cdf(const inverse_gaussian_distribution& dist, { return 0; // Convenient, even if not defined mathematically. } - // Problem with this formula for large scale > 1000 or small x, - //result = 0.5 * (erf(sqrt(scale / x) * ((x / mean) - 1) / constants::root_two(), Policy()) + 1) - // + exp(2 * scale / mean) / 2 - // * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two(), Policy())); + // Problem with this formula for large scale > 1000 or small x // so use normal distribution version: // Wikipedia CDF equation http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution. @@ -270,22 +300,20 @@ namespace detail template inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1) { // guess at random variate value x for inverse gaussian quantile. - BOOST_MATH_STD_USING - using boost::math::policies::policy; - // Error type. - using boost::math::policies::overflow_error; - // Action. - using boost::math::policies::ignore_error; + BOOST_MATH_STD_USING + using boost::math::policies::policy; + // Error type. + using boost::math::policies::overflow_error; + // Action. + using boost::math::policies::ignore_error; - typedef policy< - overflow_error // Ignore overflow (return infinity) - > no_overthrow_policy; + using no_overthrow_policy = policy>; RealType x; // result is guess at random variate value x. RealType phi = lambda / mu; if (phi > 2.) { // Big phi, so starting to look like normal Gaussian distribution. - // x=(qnorm(p,0,1,true,false) - 0.5 * sqrt(mu/lambda)) / sqrt(lambda/mu); + // // Whitmore, G.A. and Yalovsky, M. // A normalising logarithmic transformation for inverse Gaussian random variables, // Technometrics 20-2, 207-208 (1978), but using expression from @@ -300,14 +328,12 @@ namespace detail using boost::math::gamma_distribution; // Define the distribution, using gamma_nooverflow: - typedef gamma_distribution gamma_nooverflow; + using gamma_nooverflow = gamma_distribution; gamma_nooverflow g(static_cast(0.5), static_cast(1.)); - // gamma_nooverflow g(static_cast(0.5), static_cast(1.)); - // R qgamma(0.2, 0.5, 1) 0.0320923 + // R qgamma(0.2, 0.5, 1) = 0.0320923 RealType qg = quantile(complement(g, p)); - //RealType qg1 = qgamma(1.- p, 0.5, 1.0, true, false); x = lambda / (qg * 2); // if (x > mu/2) // x > mu /2? @@ -350,7 +376,7 @@ inline RealType quantile(const inverse_gaussian_distribution& { // overflow result = policies::raise_overflow_error(function, "probability parameter is 1, but must be < 1!", Policy()); - return result; // std::numeric_limits::infinity(); + return result; // infinity; } RealType guess = detail::guess_ig(p, dist.mean(), dist.scale()); @@ -378,21 +404,7 @@ inline RealType cdf(const complemented2_type::has_infinity && x == std::numeric_limits::infinity()) - //{ // cdf complement +infinity is zero. - // return 0; - //} - //if(std::numeric_limits::has_infinity && x == -std::numeric_limits::infinity()) - //{ // cdf complement -infinity is unity. - // return 1; - //} + RealType result = 0; if(false == detail::check_scale(function, scale, &result, Policy())) return result; diff --git a/include/boost/math/distributions/laplace.hpp b/include/boost/math/distributions/laplace.hpp index 9b268f3f33..cc922a879a 100644 --- a/include/boost/math/distributions/laplace.hpp +++ b/include/boost/math/distributions/laplace.hpp @@ -36,13 +36,13 @@ class laplace_distribution // ---------------------------------- // public Types // ---------------------------------- - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; // ---------------------------------- // Constructor(s) // ---------------------------------- - laplace_distribution(RealType l_location = 0, RealType l_scale = 1) + explicit laplace_distribution(RealType l_location = 0, RealType l_scale = 1) : m_location(l_location), m_scale(l_scale) { RealType result; @@ -78,7 +78,7 @@ class laplace_distribution // // Convenient type synonym for double. -typedef laplace_distribution laplace; +using laplace = laplace_distribution; #ifdef __cpp_deduction_guides template @@ -90,9 +90,9 @@ laplace_distribution(RealType,RealType)->laplace_distribution -inline const std::pair range(const laplace_distribution&) +inline std::pair range(const laplace_distribution&) { - if (std::numeric_limits::has_infinity) + if (std::numeric_limits::has_infinity) { // Can use infinity. return std::pair(-std::numeric_limits::infinity(), std::numeric_limits::infinity()); // - to + infinity. } @@ -105,7 +105,7 @@ inline const std::pair range(const laplace_distribution -inline const std::pair support(const laplace_distribution&) +inline std::pair support(const laplace_distribution&) { if (std::numeric_limits::has_infinity) { // Can Use infinity. @@ -150,6 +150,48 @@ inline RealType pdf(const laplace_distribution& dist, const Re return result; } // pdf +template +inline RealType logpdf(const laplace_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + // Checking function argument + RealType result = -std::numeric_limits::infinity(); + const char* function = "boost::math::logpdf(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 result; // pdf + and - infinity is zero so logpdf is -INF + } + if (false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType mu = dist.scale(); + const RealType b = dist.location(); + + // if b is 0 avoid divde by 0 error + if(abs(b) < std::numeric_limits::epsilon()) + { + result = log(pdf(dist, x)); + } + else + { + // General case + const RealType log2 = boost::math::constants::ln_two(); + result = -abs(x-mu)/b - log(b) - log2; + } + + return result; +} // logpdf + template inline RealType cdf(const laplace_distribution& dist, const RealType& x) { @@ -201,14 +243,14 @@ inline RealType quantile(const laplace_distribution& dist, con { result = policies::raise_overflow_error(function, "probability parameter is 0, but must be > 0!", Policy()); - return -result; // -std::numeric_limits::infinity(); + return -result; // -inf } if(p == 1) { result = policies::raise_overflow_error(function, "probability parameter is 1, but must be < 1!", Policy()); - return result; // std::numeric_limits::infinity(); + return result; // inf } // Calculate Quantile RealType scale( dist.scale() ); @@ -238,8 +280,6 @@ inline RealType cdf(const complemented2_type #include +#include namespace boost{ namespace math{ @@ -28,10 +29,10 @@ template > class normal_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - normal_distribution(RealType l_mean = 0, RealType sd = 1) + explicit 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"; @@ -69,7 +70,7 @@ class normal_distribution RealType m_sd; // distribution standard deviation or scale. }; // class normal_distribution -typedef normal_distribution normal; +using normal = normal_distribution; // // Deduction guides, note we don't check the @@ -91,7 +92,7 @@ normal_distribution(RealType)->normal_distribution -inline const std::pair range(const normal_distribution& /*dist*/) +inline std::pair range(const normal_distribution& /*dist*/) { // Range of permissible values for random variable x. if (std::numeric_limits::has_infinity) { @@ -105,7 +106,7 @@ inline const std::pair range(const normal_distribution -inline const std::pair support(const normal_distribution& /*dist*/) +inline std::pair support(const normal_distribution& /*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::has_infinity) { @@ -145,11 +146,6 @@ inline RealType pdf(const normal_distribution& dist, const Rea { return 0; // pdf + and - infinity is zero. } - // Below produces MSVC 4127 warnings, so the above used instead. - //if(std::numeric_limits::has_infinity && abs(x) == std::numeric_limits::infinity()) - //{ // pdf + and - infinity is zero. - // return 0; - //} if(false == detail::check_x(function, x, &result, Policy())) { return result; @@ -165,6 +161,42 @@ inline RealType pdf(const normal_distribution& dist, const Rea return result; } // pdf +template +inline RealType logpdf(const normal_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + const RealType sd = dist.standard_deviation(); + const RealType mean = dist.mean(); + + static const char* function = "boost::math::logpdf(const normal_distribution<%1%>&, %1%)"; + + RealType result = -std::numeric_limits::infinity(); + 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 result; // pdf + and - infinity is zero so logpdf is -inf + } + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType pi = boost::math::constants::pi(); + const RealType half = boost::math::constants::half(); + + result = -log(sd) - half*log(2*pi) - (x-mean)*(x-mean)/(2*sd*sd); + + return result; +} + template inline RealType cdf(const normal_distribution& dist, const RealType& x) { @@ -187,15 +219,6 @@ inline RealType cdf(const normal_distribution& dist, const Rea if(x < 0) return 0; // -infinity return 1; // + infinity } - // These produce MSVC 4127 warnings, so the above used instead. - //if(std::numeric_limits::has_infinity && x == std::numeric_limits::infinity()) - //{ // cdf +infinity is unity. - // return 1; - //} - //if(std::numeric_limits::has_infinity && x == -std::numeric_limits::infinity()) - //{ // cdf -infinity is zero. - // return 0; - //} if(false == detail::check_x(function, x, &result, Policy())) { return result; @@ -249,15 +272,6 @@ inline RealType cdf(const complemented2_type::has_infinity && x == std::numeric_limits::infinity()) - //{ // cdf complement +infinity is zero. - // return 0; - //} - //if(std::numeric_limits::has_infinity && x == -std::numeric_limits::infinity()) - //{ // cdf complement -infinity is unity. - // return 1; - //} if(false == detail::check_x(function, x, &result, Policy())) return result; diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp index 533c31a80d..5507360e82 100644 --- a/include/boost/math/distributions/poisson.hpp +++ b/include/boost/math/distributions/poisson.hpp @@ -47,6 +47,7 @@ #include #include +#include namespace boost { @@ -144,10 +145,10 @@ namespace boost class poisson_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - poisson_distribution(RealType l_mean = 1) : m_l(l_mean) // mean (lambda). + explicit 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( @@ -165,7 +166,7 @@ namespace boost RealType m_l; // mean number of occurrences. }; // template class poisson_distribution - typedef poisson_distribution poisson; // Reserved name of type double. + using poisson = poisson_distribution; // Reserved name of type double. #ifdef __cpp_deduction_guides template @@ -175,14 +176,14 @@ namespace boost // Non-member functions to give properties of the distribution. template - inline const std::pair range(const poisson_distribution& /* dist */) + inline std::pair range(const poisson_distribution& /* dist */) { // Range of permissible values for random variable k. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); // Max integer? } template - inline const std::pair support(const poisson_distribution& /* dist */) + inline std::pair support(const poisson_distribution& /* 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; @@ -202,15 +203,7 @@ namespace boost return floor(dist.mean()); } - //template - //inline RealType median(const poisson_distribution& dist) - //{ // median = approximately lambda + 1/3 - 0.2/lambda - // RealType l = dist.mean(); - // return dist.mean() + static_cast(0.3333333333333333333333333333333333333333333333) - // - static_cast(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. + // Median now implemented via quantile(half) in derived accessors. template inline RealType variance(const poisson_distribution& dist) @@ -218,7 +211,6 @@ namespace boost return dist.mean(); } - // RealType standard_deviation(const poisson_distribution& dist) // standard_deviation provided by derived accessors. template @@ -281,6 +273,42 @@ namespace boost return boost::math::gamma_p_derivative(k+1, mean, Policy()); } // pdf + template + RealType logpdf(const poisson_distribution& dist, const RealType& k) + { + BOOST_FPU_EXCEPTION_GUARD + + BOOST_MATH_STD_USING // for ADL of std functions. + using boost::math::lgamma; + + RealType mean = dist.mean(); + // Error check: + RealType result = -std::numeric_limits::infinity(); + 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 std::numeric_limits::quiet_NaN(); + } + + // Special case where k and lambda are both positive + if(k > 0 && mean > 0) + { + return -lgamma(k+1) + k*log(mean) - mean; + } + + result = log(pdf(dist, k)); + return result; + } + template RealType cdf(const poisson_distribution& dist, const RealType& k) { // Cumulative Distribution Function Poisson. @@ -320,8 +348,8 @@ namespace boost return 0; } if (k == 0) - { // return pdf(dist, static_cast(0)); - // but mean (and k) have already been checked, + { + // mean (and k) have already been checked, // so this avoids unnecessary repeated checks. return exp(-mean); } @@ -414,9 +442,10 @@ namespace boost { return policies::raise_overflow_error(function, 0, Policy()); } - typedef typename Policy::discrete_quantile_type discrete_type; + using discrete_type = typename Policy::discrete_quantile_type; std::uintmax_t max_iter = policies::get_max_root_iterations(); - RealType guess, factor = 8; + RealType guess; + RealType factor = 8; RealType z = dist.mean(); if(z < 1) guess = z; @@ -484,9 +513,10 @@ namespace boost { return 0; // Exact result regardless of discrete-quantile Policy } - typedef typename Policy::discrete_quantile_type discrete_type; + using discrete_type = typename Policy::discrete_quantile_type; std::uintmax_t max_iter = policies::get_max_root_iterations(); - RealType guess, factor = 8; + RealType guess; + RealType factor = 8; RealType z = dist.mean(); if(z < 1) guess = z; diff --git a/include/boost/math/distributions/rayleigh.hpp b/include/boost/math/distributions/rayleigh.hpp index cbe934471f..cbbf39876d 100644 --- a/include/boost/math/distributions/rayleigh.hpp +++ b/include/boost/math/distributions/rayleigh.hpp @@ -56,10 +56,10 @@ template > class rayleigh_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - rayleigh_distribution(RealType l_sigma = 1) + explicit rayleigh_distribution(RealType l_sigma = 1) : m_sigma(l_sigma) { RealType err; @@ -75,7 +75,7 @@ class rayleigh_distribution RealType m_sigma; }; // class rayleigh_distribution -typedef rayleigh_distribution rayleigh; +using rayleigh = rayleigh_distribution; #ifdef __cpp_deduction_guides template @@ -83,14 +83,14 @@ rayleigh_distribution(RealType)->rayleigh_distribution -inline const std::pair range(const rayleigh_distribution& /*dist*/) +inline std::pair range(const rayleigh_distribution& /*dist*/) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : max_value()); } template -inline const std::pair support(const rayleigh_distribution& /*dist*/) +inline std::pair support(const rayleigh_distribution& /*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; @@ -122,6 +122,32 @@ inline RealType pdf(const rayleigh_distribution& dist, const R return result; } // pdf +template +inline RealType logpdf(const rayleigh_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std function exp. + + const RealType sigma = dist.sigma(); + RealType result = -std::numeric_limits::infinity(); + static const char* function = "boost::math::logpdf(const rayleigh_distribution<%1%>&, %1%)"; + + if(false == detail::verify_sigma(function, sigma, &result, Policy())) + { + return result; + } + if(false == detail::verify_rayleigh_x(function, x, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + return result; + } + + result = -(x*x)/(2*sigma*sigma) - 2*log(sigma) + log(x); + return result; +} // logpdf + template inline RealType cdf(const rayleigh_distribution& dist, const RealType& x) { @@ -265,31 +291,26 @@ inline RealType median(const rayleigh_distribution& dist) template inline RealType skewness(const rayleigh_distribution& /*dist*/) { - // using namespace boost::math::constants; return static_cast(0.63111065781893713819189935154422777984404221106391L); // Computed using NTL at 150 bit, about 50 decimal digits. - // return 2 * root_pi() * pi_minus_three() / pow23_four_minus_pi(); + // 2 * sqrt(pi) * (pi-3) / pow(4, 2/3) - pi } template inline RealType kurtosis(const rayleigh_distribution& /*dist*/) { - // using namespace boost::math::constants; return static_cast(3.2450893006876380628486604106197544154170667057995L); // Computed using NTL at 150 bit, about 50 decimal digits. - // return 3 - (6 * pi() * pi() - 24 * pi() + 16) / - // (four_minus_pi() * four_minus_pi()); + // 3 - (6*pi*pi - 24*pi + 16) / pow(4-pi, 2) } template inline RealType kurtosis_excess(const rayleigh_distribution& /*dist*/) { - //using namespace boost::math::constants; - // Computed using NTL at 150 bit, about 50 decimal digits. return static_cast(0.2450893006876380628486604106197544154170667057995L); - // return -(6 * pi() * pi() - 24 * pi() + 16) / - // (four_minus_pi() * four_minus_pi()); -} // kurtosis + // Computed using NTL at 150 bit, about 50 decimal digits. + // -(6*pi*pi - 24*pi + 16) / pow(4-pi,2) +} // kurtosis_excess template inline RealType entropy(const rayleigh_distribution& dist) diff --git a/include/boost/math/distributions/weibull.hpp b/include/boost/math/distributions/weibull.hpp index 76246e4ff2..724cce6ed8 100644 --- a/include/boost/math/distributions/weibull.hpp +++ b/include/boost/math/distributions/weibull.hpp @@ -70,10 +70,10 @@ template > class weibull_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - weibull_distribution(RealType l_shape, RealType l_scale = 1) + explicit weibull_distribution(RealType l_shape, RealType l_scale = 1) : m_shape(l_shape), m_scale(l_scale) { RealType result; @@ -97,7 +97,7 @@ class weibull_distribution RealType m_scale; // distribution scale }; -typedef weibull_distribution weibull; +using weibull = weibull_distribution; #ifdef __cpp_deduction_guides template @@ -107,14 +107,14 @@ weibull_distribution(RealType,RealType)->weibull_distribution -inline const std::pair range(const weibull_distribution& /*dist*/) +inline std::pair range(const weibull_distribution& /*dist*/) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); } template -inline const std::pair support(const weibull_distribution& /*dist*/) +inline std::pair support(const weibull_distribution& /*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; @@ -157,6 +157,40 @@ inline RealType pdf(const weibull_distribution& dist, const Re return result; } +template +inline RealType logpdf(const weibull_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::logpdf(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 log(1 / scale); + } + if(shape > 1) + { + return 0; + } + return policies::raise_overflow_error(function, 0, Policy()); + } + + result = log(shape) - shape * log(scale) + log(x) * (shape - 1) - pow(x / scale, shape); + + return result; +} + template inline RealType cdf(const weibull_distribution& dist, const RealType& x) { @@ -342,12 +376,11 @@ inline RealType skewness(const weibull_distribution& dist) { 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)); + RealType g1 = boost::math::tgamma(1 + 1 / shape, Policy()); + RealType g2 = boost::math::tgamma(1 + 2 / shape, Policy()); + RealType g3 = boost::math::tgamma(1 + 3 / shape, Policy()); + RealType d = pow(g2 - g1 * g1, RealType(1.5)); result = (2 * g1 * g1 * g1 - 3 * g1 * g2 + g3) / d; return result; @@ -367,15 +400,13 @@ inline RealType kurtosis_excess(const weibull_distribution& di 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; + RealType g1 = boost::math::tgamma(1 + 1 / shape, Policy()); + RealType g2 = boost::math::tgamma(1 + 2 / shape, Policy()); + RealType g3 = boost::math::tgamma(1 + 3 / shape, Policy()); + RealType g4 = boost::math::tgamma(1 + 4 / shape, Policy()); + RealType g1_2 = g1 * g1; + RealType g1_4 = g1_2 * g1_2; + RealType d = g2 - g1_2; d *= d; result = -6 * g1_4 + 12 * g1_2 * g2 - 3 * g2 * g2 - 4 * g1 * g3 + g4; diff --git a/test/compile_test/test_compile_result.hpp b/test/compile_test/test_compile_result.hpp index a86d390fac..0afdced913 100644 --- a/test/compile_test/test_compile_result.hpp +++ b/test/compile_test/test_compile_result.hpp @@ -9,14 +9,12 @@ // detect missing includes). // -static const float f = 0; -static const double d = 0; -static const long double l = 0; -static const unsigned u = 0; -static const int i = 0; - -//template -//inline void check_result_imp(T, T){} +static constexpr float f = 0; +static constexpr double d = 0; +static constexpr long double l = 0; +static constexpr unsigned u = 0; +static constexpr int i = 0; +static constexpr unsigned long li = 1; inline void check_result_imp(float, float){} inline void check_result_imp(double, double){} @@ -38,12 +36,12 @@ inline void check_result_imp(bool, bool){} template struct local_is_same { - enum{ value = false }; + static constexpr bool value = false; }; template struct local_is_same { - enum{ value = true }; + static constexpr bool value = true; }; template @@ -54,7 +52,7 @@ inline void check_result_imp(T1, T2) #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-local-typedefs" #endif - typedef BOOST_MATH_ASSERT_UNUSED_ATTRIBUTE int static_assertion[local_is_same::value ? 1 : 0]; + using static_assertion = int[local_is_same::value ? 1 : 0]; #if defined(__GNUC__) && ((__GNUC__ > 4) || ((__GNUC__ == 4) && (__GNUC_MINOR__ >= 8))) #pragma GCC diagnostic pop #endif @@ -68,23 +66,12 @@ inline void check_result(T2) return check_result_imp(a, b); } -union max_align_type -{ - char c; - short s; - int i; - long l; - double d; - long double ld; - long long ll; -}; - template struct DistributionConcept { static void constraints() { - typedef typename Distribution::value_type value_type; + using value_type = typename Distribution::value_type; const Distribution& dist = DistributionConcept::get_object(); @@ -93,6 +80,7 @@ struct DistributionConcept check_result(cdf(dist, x)); check_result(cdf(complement(dist, x))); check_result(pdf(dist, x)); + check_result(logpdf(dist, x)); check_result(quantile(dist, x)); check_result(quantile(complement(dist, x))); check_result(mean(dist)); @@ -116,6 +104,7 @@ struct DistributionConcept check_result(cdf(dist, f)); check_result(cdf(complement(dist, f))); check_result(pdf(dist, f)); + check_result(logpdf(dist, f)); check_result(quantile(dist, f)); check_result(quantile(complement(dist, f))); check_result(hazard(dist, f)); @@ -123,6 +112,7 @@ struct DistributionConcept check_result(cdf(dist, d)); check_result(cdf(complement(dist, d))); check_result(pdf(dist, d)); + check_result(logpdf(dist, d)); check_result(quantile(dist, d)); check_result(quantile(complement(dist, d))); check_result(hazard(dist, d)); @@ -130,6 +120,7 @@ struct DistributionConcept check_result(cdf(dist, l)); check_result(cdf(complement(dist, l))); check_result(pdf(dist, l)); + check_result(logpdf(dist, l)); check_result(quantile(dist, l)); check_result(quantile(complement(dist, l))); check_result(hazard(dist, l)); @@ -137,20 +128,32 @@ struct DistributionConcept check_result(cdf(dist, i)); check_result(cdf(complement(dist, i))); check_result(pdf(dist, i)); + check_result(logpdf(dist, i)); check_result(quantile(dist, i)); check_result(quantile(complement(dist, i))); check_result(hazard(dist, i)); check_result(chf(dist, i)); - unsigned long li = 1; check_result(cdf(dist, li)); check_result(cdf(complement(dist, li))); check_result(pdf(dist, li)); + check_result(logpdf(dist, li)); check_result(quantile(dist, li)); check_result(quantile(complement(dist, li))); check_result(hazard(dist, li)); check_result(chf(dist, li)); } private: + union max_align_type + { + char c; + short s; + int i; + long l; + double d; + long double ld; + long long ll; + }; + static void* storage() { static max_align_type storage[sizeof(Distribution)]; diff --git a/test/test_arcsine.cpp b/test/test_arcsine.cpp index 68036686e5..0c2d847a96 100644 --- a/test/test_arcsine.cpp +++ b/test/test_arcsine.cpp @@ -105,6 +105,14 @@ void test_ignore_policy(RealType) BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + // Logpdf + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + // Mean BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(-nan, 0)))); BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(+nan, 0)))); @@ -239,6 +247,7 @@ void test_spots(RealType) using boost::math::arcsine_distribution; using ::boost::math::cdf; using ::boost::math::pdf; + using ::boost::math::logpdf; using ::boost::math::complement; using ::boost::math::quantile; @@ -292,6 +301,16 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(pdf(arcsine_01, static_cast(1) - tolerance), 1 /(sqrt(tolerance) * boost::math::constants::pi()), 2 * tolerance); // + // Log PDF + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.000001), static_cast(5.7630258931329868780772138043668005779060097243996L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.000005), static_cast(4.9583089369219367114435788047327747268154560240604L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.05), static_cast(0.37878289812137058928728250884555529541061717942415L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.5), static_cast(-0.45158270528945486472619522989488214357179467855506L), tolerance); + // Note loss of significance when x is near x_max. + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.95), static_cast(0.37878289812137058928728250884555529541061717942415L), 8 * tolerance); // Less accurate. + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.999995), static_cast(4.9583089369219367114435788047327747268154560240604L), 50000 * tolerance); // Much less accurate. + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.999999), static_cast(5.7630258931329868780772138043668005779060097243996L), 100000 * tolerance);// Even less accurate. + // CDF BOOST_CHECK_CLOSE_FRACTION(cdf(arcsine_01, 0.000001), static_cast(0.00063661987847092448418377367957384866092127786060574L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(arcsine_01, 0.000005), static_cast(0.0014235262731079289297302426454125318201831474507326L), tolerance); @@ -353,6 +372,10 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.5), static_cast(0.36755259694786136634088433220864629426492432024443L), tolerance); BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.95), static_cast(1.0194074882503562519812229448639426942621591013381L), 2 * tolerance); // Less accurate. + BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.05), static_cast(-1.1434783207403409089630164813372974217316704642782L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.5), static_cast(-1.0008888496235097104238178483561449958955399574664L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.95), static_cast(0.019221564639767605567429885545559909302927558782238L), 100 * tolerance); // Less accurate. + BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.05), static_cast(0.51592213323666034437274347433261364289389772737836L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.5), static_cast(0.66666666666666666666666666666666666666666666666667L), 2 * tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.95), static_cast(0.89891737589574013042121018491729701360300248368629L), tolerance); // Not less accurate. @@ -445,6 +468,31 @@ void test_spots(RealType) arcsine_distribution(static_cast(0), static_cast(1)), // bad x > 1. static_cast(999)), std::domain_error); + BOOST_MATH_CHECK_THROW( // For various bad arguments. + logpdf( + arcsine_distribution(static_cast(+1), static_cast(-1)), // min_x > max_x + static_cast(1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(1), static_cast(0)), // bad constructor parameters. + static_cast(1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(1), static_cast(-1)), // bad constructor parameters. + static_cast(1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(1), static_cast(1)), // equal constructor parameters. + static_cast(-1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(0), static_cast(1)), // bad x > 1. + static_cast(999)), std::domain_error); + // Checks on things that are errors. // Construction with 'bad' parameters. @@ -453,6 +501,7 @@ void test_spots(RealType) arcsine_distribution<> dist; BOOST_MATH_CHECK_THROW(pdf(dist, -1), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(dist, -1), std::domain_error); BOOST_MATH_CHECK_THROW(cdf(dist, -1), std::domain_error); BOOST_MATH_CHECK_THROW(cdf(complement(dist, -1)), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error); @@ -463,6 +512,8 @@ void test_spots(RealType) // Various combinations of bad constructor and member function parameters. BOOST_MATH_CHECK_THROW(pdf(boost::math::arcsine_distribution(0, 1), -1), std::domain_error); BOOST_MATH_CHECK_THROW(pdf(boost::math::arcsine_distribution(-1, 1), +2), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(boost::math::arcsine_distribution(0, 1), -1), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(boost::math::arcsine_distribution(-1, 1), +2), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(boost::math::arcsine_distribution(1, 1), -1), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(boost::math::arcsine_distribution(1, 1), 2), std::domain_error); @@ -484,6 +535,7 @@ void test_spots(RealType) arcsine_distribution w(RealType(-1), RealType(+1)); // NaN parameters to member functions should throw. BOOST_MATH_CHECK_THROW(pdf(w, +nan), std::domain_error); // x = NaN + BOOST_MATH_CHECK_THROW(logpdf(w, +nan), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(w, +nan), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(complement(w, +nan)), std::domain_error); // x = + nan BOOST_MATH_CHECK_THROW(quantile(w, +nan), std::domain_error); // p = + nan @@ -511,6 +563,7 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(arcsine_distribution(1, inf), std::domain_error); #endif BOOST_MATH_CHECK_THROW(pdf(w, +inf), std::domain_error); // x = inf + BOOST_MATH_CHECK_THROW(logpdf(w, +inf), std::domain_error); // x = inf BOOST_MATH_CHECK_THROW(cdf(w, +inf), std::domain_error); // x = inf BOOST_MATH_CHECK_THROW(cdf(complement(w, +inf)), std::domain_error); // x = + inf BOOST_MATH_CHECK_THROW(quantile(w, +inf), std::domain_error); // p = + inf diff --git a/test/test_chi_squared.cpp b/test/test_chi_squared.cpp index 41e4dd5e50..cc7747a6c0 100644 --- a/test/test_chi_squared.cpp +++ b/test/test_chi_squared.cpp @@ -36,6 +36,8 @@ using std::cout; using std::endl; #include using std::numeric_limits; +#include +using std::log; template RealType naive_pdf(RealType df, RealType x) @@ -59,6 +61,8 @@ void test_spot( cdf(dist, cs), P, tol); BOOST_CHECK_CLOSE( pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), cs), tol); + BOOST_CHECK_CLOSE( + logpdf(dist, cs), log(pdf(dist, cs)), tol); if((P < 0.99) && (Q < 0.99)) { // diff --git a/test/test_exponential_dist.cpp b/test/test_exponential_dist.cpp index c0423bf0de..12d9fcad19 100644 --- a/test/test_exponential_dist.cpp +++ b/test/test_exponential_dist.cpp @@ -189,6 +189,31 @@ void test_spots(RealType T) static_cast(9.0799859524969703071183031121101e-5L), // probability. tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(0.5), + static_cast(0.125)), // x + log(static_cast(0.46970653140673789305985541231115L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(0.5), + static_cast(5)), // x + log(static_cast(0.04104249931194939758476433723358L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(2), + static_cast(0.125)), // x + log(static_cast(1.5576015661428097364903405339566L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(2), + static_cast(5)), // x + log(static_cast(9.0799859524969703071183031121101e-5L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( ::boost::math::mean( exponential_distribution(2)), diff --git a/test/test_extreme_value.cpp b/test/test_extreme_value.cpp index f61731647e..2493f466df 100644 --- a/test/test_extreme_value.cpp +++ b/test/test_extreme_value.cpp @@ -120,6 +120,25 @@ void test_spots(RealType) static_cast(0.11522236828583456431277265757312L), // probability. tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + extreme_value_distribution(0.5, 2), + static_cast(0.125)), // x + log(static_cast(0.18052654830890205978204427757846L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + extreme_value_distribution(1, 3), + static_cast(5)), // x + log(static_cast(0.0675057324099851209129017326286L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + extreme_value_distribution(1, 3), + static_cast(0)), // x + log(static_cast(0.11522236828583456431277265757312L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( ::boost::math::mean( extreme_value_distribution(2, 3)), diff --git a/test/test_gamma_dist.cpp b/test/test_gamma_dist.cpp index 682535c283..b7776c79cb 100644 --- a/test/test_gamma_dist.cpp +++ b/test/test_gamma_dist.cpp @@ -88,6 +88,14 @@ void check_gamma(RealType shape, RealType scale, RealType x, RealType p, RealTyp x), // random variable. NaivePDF(shape, scale, x), // PDF tol); // %tolerance. + + // LOGPDF: + BOOST_CHECK_CLOSE( + boost::math::logpdf( + gamma_distribution(shape, scale), // distribution. + x), // random variable. + log(boost::math::pdf(gamma_distribution(shape, scale), x)), // PDF + tol); // %tolerance. } template diff --git a/test/test_inverse_gamma_distribution.cpp b/test/test_inverse_gamma_distribution.cpp index bad1bb158b..68b238fbc8 100644 --- a/test/test_inverse_gamma_distribution.cpp +++ b/test/test_inverse_gamma_distribution.cpp @@ -72,6 +72,9 @@ void test_spot( BOOST_CHECK_CLOSE_FRACTION( // Compare to naive formula (might be less accurate). pdf(dist, x), naive_pdf(dist.shape(), dist.scale(), x), tol); + BOOST_CHECK_CLOSE_FRACTION( // Compare direct logpdf to naive log(pdf()) + logpdf(dist, x), log(pdf(dist,x)), tol); + BOOST_CHECK_CLOSE_FRACTION( // Compare to expected CDF. cdf(dist, x), P, tol); diff --git a/test/test_inverse_gaussian.cpp b/test/test_inverse_gaussian.cpp index 2bb5303122..68012d48a2 100644 --- a/test/test_inverse_gaussian.cpp +++ b/test/test_inverse_gaussian.cpp @@ -36,6 +36,8 @@ using std::endl; using std::setprecision; #include using std::numeric_limits; +#include +using std::log; template void check_inverse_gaussian(RealType mean, RealType scale, RealType x, RealType p, RealType q, RealType tol) @@ -207,21 +209,29 @@ BOOST_AUTO_TEST_CASE( test_main ) // formatC(SuppDists::dinverse_gaussian(1, 1, 1), digits=17) ... BOOST_CHECK_CLOSE_FRACTION( // x = 1 pdf(w11, 1.), static_cast(0.3989422804014327), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( // x = 1 + logpdf(w11, 1.), static_cast(log(0.3989422804014327)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 1.), static_cast(0.66810200122317065), 10 * tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 0.1), static_cast(0.21979480031862672), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 0.1), static_cast(log(0.21979480031862672)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.1), static_cast(0.0040761113207110162), 10 * tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( // small x pdf(w11, 0.01), static_cast(2.0811768202028392e-19), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( // small x + logpdf(w11, 0.01), static_cast(log(2.0811768202028392e-19)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.01), static_cast(4.122313403318778e-23), 10 * tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( // smaller x pdf(w11, 0.001), static_cast(2.4420044378793562e-213), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( // smaller x + logpdf(w11, 0.001), static_cast(log(2.4420044378793562e-213)), tolfeweps); // pdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.001), static_cast(4.8791443010851493e-219), 1000 * tolfeweps); // cdf // 4.8791443010859224e-219 versus 4.8791443010851493e-219 so still 14 decimal digits. @@ -240,25 +250,35 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 0.5), static_cast(0.87878257893544476), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 0.5), static_cast(log(0.87878257893544476)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.5), static_cast(0.3649755481729598), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 2), static_cast(0.10984782236693059), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 2), static_cast(log(0.10984782236693059)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 2), static_cast(.88547542598600637), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 10), static_cast(0.00021979480031862676), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 10), static_cast(log(0.00021979480031862676)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 10), static_cast(0.99964958546279115), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 100), static_cast(2.0811768202028246e-25), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 100), static_cast(log(2.0811768202028246e-25)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 100), static_cast(1), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 1000), static_cast(2.4420044378793564e-222), 10 * tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 1000), static_cast(log(2.4420044378793564e-222)), 10 * tolfeweps); // pdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 1000), static_cast(1.), tolfeweps); // cdf @@ -295,26 +315,37 @@ BOOST_AUTO_TEST_CASE( test_main ) // =================== BOOST_CHECK_CLOSE_FRACTION( // formatC(SuppDists::dinvGauss(1, 2, 3), digits=17) "0.47490884963330904" pdf(w23, 1.), static_cast(0.47490884963330904), tolfeweps ); // pdf - + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 1.), static_cast(log(0.47490884963330904)), tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( pdf(w23, 0.1), static_cast(2.8854207087665401e-05), tolfeweps * 2); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 0.1), static_cast(log(2.8854207087665401e-05)), tolfeweps * 2); // logpdf //2.8854207087665452e-005 2.8854207087665401e-005 BOOST_CHECK_CLOSE_FRACTION( pdf(w23, 10.), static_cast(0.0019822751498574636), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 10.), static_cast(log(0.0019822751498574636)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( pdf(w23, 10.), static_cast(0.0019822751498574636), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 10.), static_cast(log(0.0019822751498574636)), tolfeweps); // logpdf // Bigger changes in mean and scale. inverse_gaussian w012(0.1, 2); BOOST_CHECK_CLOSE_FRACTION( pdf(w012, 1.), static_cast(3.7460367141230404e-36), tolfeweps ); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w012, 1.), static_cast(log(3.7460367141230404e-36)), tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w012, 1.), static_cast(1), tolfeweps ); // pdf inverse_gaussian w0110(0.1, 10); BOOST_CHECK_CLOSE_FRACTION( pdf(w0110, 1.), static_cast(1.6279643678071011e-176), 100 * tolfeweps ); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w0110, 1.), static_cast(log(1.6279643678071011e-176)), 100 * tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w0110, 1.), static_cast(1), tolfeweps ); // cdf BOOST_CHECK_CLOSE_FRACTION( @@ -323,6 +354,8 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_CLOSE_FRACTION( pdf(w0110, 0.1), static_cast(39.894228040143268), tolfeweps ); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w0110, 0.1), static_cast(log(39.894228040143268)), tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w0110, 0.1), static_cast(0.51989761564832704), 10 * tolfeweps ); // cdf diff --git a/test/test_laplace.cpp b/test/test_laplace.cpp index 9167547f91..716c738736 100644 --- a/test/test_laplace.cpp +++ b/test/test_laplace.cpp @@ -67,6 +67,8 @@ Test 8: test_extreme_function_arguments() #include #include "test_out_of_range.hpp" using boost::math::laplace_distribution; +#include +using std::log; /* #include @@ -211,6 +213,12 @@ void test_pdf_cdf_ocatave() static_cast(0.067667641618306345946999747486242201703815773119812L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(-2.L)), + // static_cast(0.06766764161831L), + log(static_cast(0.067667641618306345946999747486242201703815773119812L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(-2.L)), //static_cast(0.06766764161831L), @@ -223,6 +231,12 @@ void test_pdf_cdf_ocatave() static_cast(0.18393972058572116079776188508073043372290556554506L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(-1.L)), + //static_cast(0.18393972058572L), + log(static_cast(0.18393972058572116079776188508073043372290556554506L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(-1.L)), static_cast(0.18393972058572L), @@ -245,6 +259,11 @@ void test_pdf_cdf_ocatave() static_cast(0.5L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(0.0L)), + log(static_cast(0.5L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(0.0L)), static_cast(0.5L), @@ -256,6 +275,11 @@ void test_pdf_cdf_ocatave() static_cast(0.30326532985631671180189976749559022672095906778368L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(0.5L)), + log(static_cast(0.30326532985631671180189976749559022672095906778368L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(0.5L)), // static_cast(0.69673467014368L), @@ -268,6 +292,12 @@ void test_pdf_cdf_ocatave() static_cast(0.18393972058572116079776188508073043372290556554506L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(1.0L)), + // static_cast(0.18393972058572L), + log(static_cast(0.18393972058572116079776188508073043372290556554506L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(1.00000000000000L)), // static_cast(0.81606027941428L), @@ -280,6 +310,12 @@ void test_pdf_cdf_ocatave() static_cast(0.067667641618306345946999747486242201703815772944649L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(2.0L)), + // static_cast(0.06766764161831L), + log(static_cast(0.067667641618306345946999747486242201703815772944649L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(2.0L)), // static_cast(0.93233235838169L), diff --git a/test/test_normal.cpp b/test/test_normal.cpp index b8d0024872..ef984d5e63 100644 --- a/test/test_normal.cpp +++ b/test/test_normal.cpp @@ -41,6 +41,8 @@ using std::setprecision; #include using std::numeric_limits; +#include + using std::log; template RealType NaivePDF(RealType mean, RealType sd, RealType x) @@ -126,6 +128,7 @@ void test_spots(RealType) { // No longer allow x to be NaN, then these tests should throw. BOOST_MATH_CHECK_THROW(pdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN + BOOST_MATH_CHECK_THROW(logpdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(complement(N01, +std::numeric_limits::quiet_NaN())), std::domain_error); // x = + infinity BOOST_MATH_CHECK_THROW(quantile(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // p = + infinity @@ -215,6 +218,31 @@ void test_spots(RealType) static_cast(0.3989422804014326779399460599343818684759L / 5), tolerance); + // + // Tests for logpdf + // + RealType temp_tol = tolerance; + + BOOST_IF_CONSTEXPR (std::is_same::value) + { + tolerance *= 100; + } + + BOOST_CHECK_CLOSE( + logpdf(normal_distribution(), static_cast(0)), + log(static_cast(0.3989422804014326779399460599343818684759L)), // 1/sqrt(2*pi) + tolerance); + BOOST_CHECK_CLOSE( + logpdf(normal_distribution(3), static_cast(3)), + log(static_cast(0.3989422804014326779399460599343818684759L)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(normal_distribution(3, 5), static_cast(3)), + log(static_cast(0.3989422804014326779399460599343818684759L / 5)), + tolerance); + + tolerance = temp_tol; + // // Spot checks for mean = -5, sd = 6: // @@ -307,6 +335,8 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(pdf(normal_distribution(0, 0), 0), std::domain_error); BOOST_MATH_CHECK_THROW(pdf(normal_distribution(0, -1), 0), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(normal_distribution(0, 0), 0), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(normal_distribution(0, -1), 0), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(normal_distribution(0, 1), -1), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(normal_distribution(0, 1), 2), std::domain_error); } // template void test_spots(RealType) diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 934395fe8e..9b75ce162f 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -208,6 +208,31 @@ void test_spots(RealType) static_cast(20)), // K>> mean static_cast(8.277463646553730E-009), // probability. tolerance); + + // LOGPDF + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(4)), // mean 4. + static_cast(0)), + log(static_cast(1.831563888873410E-002)), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(4)), // mean 4. + static_cast(2)), + log(static_cast(1.465251111098740E-001)), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(20)), // mean big. + static_cast(1)), // k small + log(static_cast(4.122307244877130E-008)), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(4)), // mean 4. + static_cast(20)), // K>> mean + log(static_cast(8.277463646553730E-009)), // probability. + tolerance); // CDF BOOST_CHECK_CLOSE( diff --git a/test/test_rayleigh.cpp b/test/test_rayleigh.cpp index 388adc68fc..6aa2d2fb83 100644 --- a/test/test_rayleigh.cpp +++ b/test/test_rayleigh.cpp @@ -26,6 +26,8 @@ using std::cout; using std::endl; using std::setprecision; +#include + using std::log; template void test_spot(RealType s, RealType x, RealType p, RealType q, RealType tolerance) @@ -159,6 +161,25 @@ void test_spots(RealType T) static_cast(exp_minus_half() /2), // probability. tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + rayleigh_distribution(1.L), + static_cast(1.L)), // x + log(static_cast(exp_minus_half())), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + rayleigh_distribution(0.5L), + static_cast(0.5L)), // x + log(static_cast(2 * exp_minus_half())), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + rayleigh_distribution(2.L), + static_cast(2.L)), // x + log(static_cast(exp_minus_half() /2)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( ::boost::math::mean( rayleigh_distribution(1.L)), @@ -255,6 +276,27 @@ BOOST_AUTO_TEST_CASE( test_main ) static_cast(exp_minus_half() /2 ), // p 1e-15); // % + BOOST_CHECK_CLOSE_FRACTION( + ::boost::math::logpdf( + rayleigh_distribution(1.), + static_cast(1)), // x + log(static_cast(exp_minus_half())), // p + 1e-15); // % + + BOOST_CHECK_CLOSE_FRACTION( + ::boost::math::logpdf( + rayleigh_distribution(0.5), + static_cast(0.5)), // x + log(static_cast(2 * exp_minus_half())), // p + 1e-15); // % + + BOOST_CHECK_CLOSE_FRACTION( + ::boost::math::logpdf( + rayleigh_distribution(2.), + static_cast(2)), // x + log(static_cast(exp_minus_half() /2 )), // p + 1e-15); // % + BOOST_CHECK_CLOSE_FRACTION( ::boost::math::cdf( rayleigh_distribution(1.), diff --git a/test/test_weibull.cpp b/test/test_weibull.cpp index 17e2bffb04..6cacf5402a 100644 --- a/test/test_weibull.cpp +++ b/test/test_weibull.cpp @@ -29,6 +29,8 @@ using std::setprecision; #include using std::numeric_limits; +#include + using std::log; template void check_weibull(RealType shape, RealType scale, RealType x, RealType p, RealType q, RealType tol) @@ -244,6 +246,50 @@ void test_spots(RealType) static_cast(0.551819), tolerance); + // + // Tests for logpdf + // + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.25, 0.5), static_cast(0.1)), + log(static_cast(0.856579)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.25, 0.5), static_cast(0.5)), + log(static_cast(0.183940)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.25, 0.5), static_cast(5)), + log(static_cast(0.015020)), + tolerance * 10); // fewer digits in test value + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.5, 2), static_cast(0.1)), + log(static_cast(0.894013)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.5, 2), static_cast(0.5)), + log(static_cast(0.303265)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.5, 2), static_cast(1)), + log(static_cast(0.174326)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(2, 0.25), static_cast(0.1)), + log(static_cast(2.726860)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(2, 0.25), static_cast(0.5)), + log(static_cast(0.293050)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(3, 2), static_cast(1)), + log(static_cast(0.330936)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(3, 2), static_cast(2)), + log(static_cast(0.551819)), + tolerance); + // // These test values were obtained using the formulas at // http://en.wikipedia.org/wiki/Weibull_distribution