From 93593fb8f2856a4e5eb9d69a894d131f138451a8 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Thu, 24 Feb 2022 12:19:15 +0100 Subject: [PATCH] Add logpdf to chi squared distribution [ci skip] --- .../boost/math/distributions/chi_squared.hpp | 41 +++++++++++++++++++ test/test_chi_squared.cpp | 4 ++ 2 files changed, 45 insertions(+) diff --git a/include/boost/math/distributions/chi_squared.hpp b/include/boost/math/distributions/chi_squared.hpp index e97bee7ce..a81197a85 100644 --- a/include/boost/math/distributions/chi_squared.hpp +++ b/include/boost/math/distributions/chi_squared.hpp @@ -131,6 +131,47 @@ RealType pdf(const chi_squared_distribution& dist, const RealT return gamma_p_derivative(degrees_of_freedom / 2, chi_square / 2, Policy()) / 2; } // pdf +template +RealType logpdf(const chi_squared_distribution& dist, const RealType& chi_square) +{ + BOOST_MATH_STD_USING // for ADL of std functions + RealType k = dist.degrees_of_freedom(); + // Error check: + RealType error_result; + + static const char* function = "boost::math::logpdf(const chi_squared_distribution<%1%>&, %1%)"; + + if(false == detail::check_df(function, k, &error_result, Policy())) + { + return error_result; + } + + if((chi_square < 0) || !(boost::math::isfinite)(chi_square)) + { + return policies::raise_domain_error( + function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy()); + } + + if(chi_square == 0) + { + // Handle special cases: + if(k < 2) + { + return policies::raise_overflow_error(function, 0, Policy()); + } + else if(k == 2) + { + return -boost::math::constants::ln_two(); + } + else + { + return -std::numeric_limits::infinity(); + } + } + + return log(pdf(dist, chi_square)); +} // logpdf + template inline RealType cdf(const chi_squared_distribution& dist, const RealType& chi_square) { diff --git a/test/test_chi_squared.cpp b/test/test_chi_squared.cpp index 41e4dd5e5..cc7747a6c 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)) { //