2
0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00

Add logpdf to inverse gamma distribution

[ci skip]
This commit is contained in:
Matt Borland
2022-02-24 13:08:50 +01:00
parent 208a1bef46
commit 570b40e14a
2 changed files with 41 additions and 0 deletions

View File

@@ -28,6 +28,7 @@
#include <boost/math/distributions/complement.hpp>
#include <utility>
#include <cfenv>
namespace boost{ namespace math
{
@@ -195,6 +196,43 @@ inline RealType pdf(const inverse_gamma_distribution<RealType, Policy>& dist, co
return result;
} // pdf
template <class RealType, class Policy>
inline RealType logpdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
{
BOOST_MATH_STD_USING // for ADL of std functions
using boost::math::tgamma;
static const char* function = "boost::math::logpdf(const inverse_gamma_distribution<%1%>&, %1%)";
RealType shape = dist.shape();
RealType scale = dist.scale();
RealType result = 0;
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 0;
}
else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
{ // x bad.
return result;
}
result = scale / x;
if(result < tools::min_value<RealType>())
return 0; // 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<RealType, Policy>(function, "PDF is infinite.", Policy());
}
return shape * log(scale) + (-shape-1)*log(x) - log(tgamma(shape)) - (scale/x);
} // pdf
template <class RealType, class Policy>
inline RealType cdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
{

View File

@@ -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);