diff --git a/include/boost/math/distributions/weibull.hpp b/include/boost/math/distributions/weibull.hpp index 76246e4ff..572c96716 100644 --- a/include/boost/math/distributions/weibull.hpp +++ b/include/boost/math/distributions/weibull.hpp @@ -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 * pow(scale, -shape) * pow(x, shape - 1)) - pow(x / scale, shape); + + return result; +} + template inline RealType cdf(const weibull_distribution& dist, const RealType& x) { diff --git a/test/test_weibull.cpp b/test/test_weibull.cpp index 17e2bffb0..6cacf5402 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