From 570b40e14aae9091da5aebd454ccce41c9bc210b Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Thu, 24 Feb 2022 13:08:50 +0100
Subject: [PATCH] Add logpdf to inverse gamma distribution

[ci skip]
---
 .../math/distributions/inverse_gamma.hpp      | 38 +++++++++++++++++++
 test/test_inverse_gamma_distribution.cpp      |  3 ++
 2 files changed, 41 insertions(+)

diff --git a/include/boost/math/distributions/inverse_gamma.hpp b/include/boost/math/distributions/inverse_gamma.hpp
index 9266bc22f..139e8982a 100644
--- a/include/boost/math/distributions/inverse_gamma.hpp
+++ b/include/boost/math/distributions/inverse_gamma.hpp
@@ -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)
 {
diff --git a/test/test_inverse_gamma_distribution.cpp b/test/test_inverse_gamma_distribution.cpp
index bad1bb158..68b238fbc 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);