Logpdf support

include/boost/math/distributions/arcsine.hpp

@@ -29,6 +29,7 @@
 #ifndef BOOST_MATH_DIST_ARCSINE_HPP
 #define BOOST_MATH_DIST_ARCSINE_HPP
 
+#include <cmath>
 #include <boost/math/distributions/fwd.hpp>
 #include <boost/math/distributions/complement.hpp> // complements.
 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks.

include/boost/math/distributions/chi_squared.hpp

@@ -131,6 +131,47 @@ RealType pdf(const chi_squared_distribution<RealType, Policy>& dist, const RealT
    return gamma_p_derivative(degrees_of_freedom / 2, chi_square / 2, Policy()) / 2;
 } // pdf
 
+template <class RealType, class Policy>
+RealType logpdf(const chi_squared_distribution<RealType, Policy>& 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<RealType>(
+         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<RealType>(function, 0, Policy());
+      }
+      else if(k == 2)
+      {
+         return -boost::math::constants::ln_two<RealType>();
+      }
+      else
+      {
+         return -std::numeric_limits<RealType>::infinity();
+      }
+   }
+
+   return log(pdf(dist, chi_square));
+} // logpdf
+
 template <class RealType, class Policy>
 inline RealType cdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square)
 {

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<value_type>(x));
 }
 template <class Distribution, class RealType>
+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<value_type>(x)));
+}
+template <class Distribution, class RealType>
 inline typename Distribution::value_type cdf(const Distribution& dist, const RealType& x)
 {
    typedef typename Distribution::value_type value_type;

include/boost/math/distributions/exponential.hpp

@@ -127,6 +127,24 @@ inline RealType pdf(const exponential_distribution<RealType, Policy>& dist, cons
    return result;
 } // pdf
 
+template <class RealType, class Policy>
+inline RealType logpdf(const exponential_distribution<RealType, Policy>& 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 = 0;
+   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 <class RealType, class Policy>
 inline RealType cdf(const exponential_distribution<RealType, Policy>& dist, const RealType& x)
 {

include/boost/math/distributions/extreme_value.hpp

@@ -122,6 +122,31 @@ inline RealType pdf(const extreme_value_distribution<RealType, Policy>& dist, co
    return result;
 } // pdf
 
+template <class RealType, class Policy>
+inline RealType logpdf(const extreme_value_distribution<RealType, Policy>& 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 = 0;
+   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<RealType>())
+      result = log(1/b) + e - exp(e);
+   // else.... result *must* be zero since exp(e) is infinite...
+   return result;
+} // pdf
+
 template <class RealType, class Policy>
 inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
 {

include/boost/math/distributions/gamma.hpp

@@ -17,6 +17,7 @@
 #include <boost/math/distributions/complement.hpp>
 
 #include <utility>
+#include <type_traits>
 
 namespace boost{ namespace math
 {
@@ -147,6 +148,41 @@ inline RealType pdf(const gamma_distribution<RealType, Policy>& dist, const Real
    return result;
 } // pdf
 
+template <class RealType, class Policy>
+inline RealType logpdf(const 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 gamma_distribution<%1%>&, %1%)";
+
+   RealType k = dist.shape();
+   RealType theta = dist.scale();
+
+   RealType result = 0;
+   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<RealType>::quiet_NaN();
+   }
+
+   // The following calculation does not always work with float so take the naive road out
+   BOOST_IF_CONSTEXPR(std::is_same<RealType, float>::value)
+   {
+      result = log(pdf(dist, x));
+   }
+   else
+   {
+      result = -k*log(theta) + (k-1)*log(x) - log(tgamma(k)) - (x/theta);
+   }
+
+   return result;
+} // logpdf
+
 template <class RealType, class Policy>
 inline RealType cdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
 {

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)
 {

include/boost/math/distributions/inverse_gaussian.hpp

@@ -174,6 +174,43 @@ inline RealType pdf(const inverse_gaussian_distribution<RealType, Policy>& dist,
    return result;
 } // pdf
 
+template <class RealType, class Policy>
+inline RealType logpdf(const inverse_gaussian_distribution<RealType, Policy>& 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 = 0;
+   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 0; // Convenient, even if not defined mathematically.
+   }
+
+   const RealType two_pi = boost::math::constants::two_pi<RealType>();
+
+   result = (-scale*pow(mean - x, RealType(2))/(mean*mean*x) + log(scale) - 3*log(x) - log(two_pi)) / 2;
+   return result;
+} // pdf
+
 template <class RealType, class Policy>
 inline RealType cdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
 { // Cumulative Density Function.

include/boost/math/distributions/laplace.hpp

@@ -150,6 +150,44 @@ inline RealType pdf(const laplace_distribution<RealType, Policy>& dist, const Re
    return result;
 } // pdf
 
+template <class RealType, class Policy>
+inline RealType logpdf(const laplace_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   // Checking function argument
+   RealType result = 0;
+   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 0; // pdf + and - infinity is zero.
+   }
+   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<RealType>::epsilon())
+   {
+      result = log(pdf(dist, x));
+   }
+   else
+   {
+      // General case
+      const RealType log2 = boost::math::constants::ln_two<RealType>();
+      result = -abs(x-mu)/b - log(b) - log2;
+   }
+
+   return result;
+} // logpdf
+
 template <class RealType, class Policy>
 inline RealType cdf(const laplace_distribution<RealType, Policy>& dist, const RealType& x)
 {

include/boost/math/distributions/normal.hpp

@@ -21,6 +21,7 @@
 #include <boost/math/distributions/detail/common_error_handling.hpp>
 
 #include <utility>
+#include <type_traits>
 
 namespace boost{ namespace math{
 
@@ -165,6 +166,42 @@ inline RealType pdf(const normal_distribution<RealType, Policy>& dist, const Rea
    return result;
 } // pdf
 
+template <class RealType, class Policy>
+inline RealType logpdf(const normal_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType sd = dist.standard_deviation();
+   RealType mean = dist.mean();
+
+   static const char* function = "boost::math::logpdf(const normal_distribution<%1%>&, %1%)";
+
+   RealType result = 0;
+   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 0; // pdf + and - infinity is zero.
+   }
+   if(false == detail::check_x(function, x, &result, Policy()))
+   {
+      return result;
+   }
+
+   const RealType pi = boost::math::constants::pi<RealType>();
+   const RealType half = boost::math::constants::half<RealType>();
+
+   result = -log(sd) - half*log(2*pi) - (x-mean)*(x-mean)/(2*sd*sd);
+
+   return result;
+}
+
 template <class RealType, class Policy>
 inline RealType cdf(const normal_distribution<RealType, Policy>& dist, const RealType& x)
 {