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

include/boost/math/distributions/poisson.hpp

@@ -47,6 +47,7 @@
 #include <boost/math/distributions/detail/inv_discrete_quantile.hpp>
 
 #include <utility>
+#include <limits>
 
 namespace boost
 {
@@ -281,6 +282,42 @@ namespace boost
       return boost::math::gamma_p_derivative(k+1, mean, Policy());
     } // pdf
 
+    template <class RealType, class Policy>
+    RealType logpdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k)
+    {
+      BOOST_FPU_EXCEPTION_GUARD
+
+      BOOST_MATH_STD_USING // for ADL of std functions.
+      using boost::math::tgamma;
+
+      RealType mean = dist.mean();
+      // Error check:
+      RealType result = 0;
+      if(false == poisson_detail::check_dist_and_k(
+        "boost::math::pdf(const poisson_distribution<%1%>&, %1%)",
+        mean,
+        k,
+        &result, Policy()))
+      {
+        return result;
+      }
+
+      // Special case of mean zero, regardless of the number of events k.
+      if (mean == 0)
+      { // Probability for any k is zero.
+        return std::numeric_limits<RealType>::quiet_NaN();
+      }
+
+      // Special case where k and lambda are both positive
+      if(k > 0 && mean > 0)
+      {
+        return -log(tgamma(k+1)) + k*log(mean) - mean;
+      }
+
+      result = log(pdf(dist, k));
+      return result;
+    }
+
     template <class RealType, class Policy>
     RealType cdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k)
     { // Cumulative Distribution Function Poisson.

test/test_arcsine.cpp

@@ -105,6 +105,14 @@ void test_ignore_policy(RealType)
       BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(0, 1), static_cast <RealType>(+2))));  // x > x_max
       BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(-1, 1), static_cast <RealType>(+2)))); // x > x_max
 
+      // Logpdf
+      BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), std::numeric_limits<RealType>::infinity()))); // x == infinity
+      BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), std::numeric_limits<RealType>::infinity()))); // x == infinity
+      BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast <RealType>(-2))));  // x < xmin
+      BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast <RealType>(-2))));  // x < xmin
+      BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast <RealType>(+2))));  // x > x_max
+      BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast <RealType>(+2)))); // x > x_max
+
       // Mean
       BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(-nan, 0))));
       BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(+nan, 0))));
@@ -239,6 +247,7 @@ void test_spots(RealType)
     using boost::math::arcsine_distribution;
     using  ::boost::math::cdf;
     using  ::boost::math::pdf;
+    using  ::boost::math::logpdf;
     using  ::boost::math::complement;
     using  ::boost::math::quantile;
 
@@ -292,6 +301,16 @@ void test_spots(RealType)
     BOOST_CHECK_CLOSE_FRACTION(pdf(arcsine_01, static_cast<RealType>(1) - tolerance),
       1 /(sqrt(tolerance) * boost::math::constants::pi<RealType>()), 2 * tolerance); //
 
+    // Log PDF
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.000001), static_cast<RealType>(5.7630258931329868780772138043668005779060097243996L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.000005), static_cast<RealType>(4.9583089369219367114435788047327747268154560240604L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.05), static_cast<RealType>(0.37878289812137058928728250884555529541061717942415L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.5), static_cast<RealType>(-0.45158270528945486472619522989488214357179467855506L), tolerance);
+    // Note loss of significance when x is near x_max.
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.95), static_cast<RealType>(0.37878289812137058928728250884555529541061717942415L), 8 * tolerance); // Less accurate.
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.999995), static_cast<RealType>(4.9583089369219367114435788047327747268154560240604L), 50000 * tolerance); // Much less accurate.
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.999999), static_cast<RealType>(5.7630258931329868780772138043668005779060097243996L), 100000 * tolerance);// Even less accurate.
+
     // CDF
     BOOST_CHECK_CLOSE_FRACTION(cdf(arcsine_01, 0.000001), static_cast<RealType>(0.00063661987847092448418377367957384866092127786060574L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(cdf(arcsine_01, 0.000005), static_cast<RealType>(0.0014235262731079289297302426454125318201831474507326L), tolerance);
@@ -353,6 +372,10 @@ void test_spots(RealType)
     BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.5), static_cast<RealType>(0.36755259694786136634088433220864629426492432024443L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.95), static_cast<RealType>(1.0194074882503562519812229448639426942621591013381L), 2 * tolerance); // Less accurate.
 
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.05), static_cast<RealType>(-1.1434783207403409089630164813372974217316704642782L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.5), static_cast<RealType>(-1.0008888496235097104238178483561449958955399574664L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.95), static_cast<RealType>(0.019221564639767605567429885545559909302927558782238L), 100 * tolerance); // Less accurate.
+
     BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.05), static_cast<RealType>(0.51592213323666034437274347433261364289389772737836L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.5), static_cast<RealType>(0.66666666666666666666666666666666666666666666666667L), 2 * tolerance);
     BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.95), static_cast<RealType>(0.89891737589574013042121018491729701360300248368629L), tolerance); //  Not less accurate.
@@ -445,6 +468,31 @@ void test_spots(RealType)
       arcsine_distribution<RealType>(static_cast<RealType>(0), static_cast<RealType>(1)), // bad x > 1.
       static_cast<RealType>(999)), std::domain_error);
 
+    BOOST_MATH_CHECK_THROW( // For various bad arguments.
+      logpdf(
+      arcsine_distribution<RealType>(static_cast<RealType>(+1), static_cast<RealType>(-1)), // min_x > max_x
+      static_cast<RealType>(1)), std::domain_error);
+
+    BOOST_MATH_CHECK_THROW(
+      logpdf(
+      arcsine_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(0)), // bad constructor parameters.
+      static_cast<RealType>(1)), std::domain_error);
+
+    BOOST_MATH_CHECK_THROW(
+      logpdf(
+      arcsine_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(-1)), // bad constructor parameters.
+      static_cast<RealType>(1)), std::domain_error);
+
+    BOOST_MATH_CHECK_THROW(
+      logpdf(
+      arcsine_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)), // equal constructor parameters.
+      static_cast<RealType>(-1)), std::domain_error);
+
+    BOOST_MATH_CHECK_THROW(
+      logpdf(
+      arcsine_distribution<RealType>(static_cast<RealType>(0), static_cast<RealType>(1)), // bad x > 1.
+      static_cast<RealType>(999)), std::domain_error);
+
     // Checks on things that are errors.
 
     // Construction with 'bad' parameters.
@@ -453,6 +501,7 @@ void test_spots(RealType)
 
     arcsine_distribution<> dist;
     BOOST_MATH_CHECK_THROW(pdf(dist, -1), std::domain_error);
+    BOOST_MATH_CHECK_THROW(logpdf(dist, -1), std::domain_error);
     BOOST_MATH_CHECK_THROW(cdf(dist, -1), std::domain_error);
     BOOST_MATH_CHECK_THROW(cdf(complement(dist, -1)), std::domain_error);
     BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error);
@@ -463,6 +512,8 @@ void test_spots(RealType)
     // Various combinations of bad constructor and member function parameters.
     BOOST_MATH_CHECK_THROW(pdf(boost::math::arcsine_distribution<RealType>(0, 1), -1), std::domain_error);
     BOOST_MATH_CHECK_THROW(pdf(boost::math::arcsine_distribution<RealType>(-1, 1), +2), std::domain_error);
+    BOOST_MATH_CHECK_THROW(logpdf(boost::math::arcsine_distribution<RealType>(0, 1), -1), std::domain_error);
+    BOOST_MATH_CHECK_THROW(logpdf(boost::math::arcsine_distribution<RealType>(-1, 1), +2), std::domain_error);
     BOOST_MATH_CHECK_THROW(quantile(boost::math::arcsine_distribution<RealType>(1, 1), -1), std::domain_error);
     BOOST_MATH_CHECK_THROW(quantile(boost::math::arcsine_distribution<RealType>(1, 1), 2), std::domain_error);
 
@@ -484,6 +535,7 @@ void test_spots(RealType)
       arcsine_distribution<RealType> w(RealType(-1), RealType(+1));
       // NaN parameters to member functions should throw.
       BOOST_MATH_CHECK_THROW(pdf(w, +nan), std::domain_error); // x = NaN
+      BOOST_MATH_CHECK_THROW(logpdf(w, +nan), std::domain_error); // x = NaN
       BOOST_MATH_CHECK_THROW(cdf(w, +nan), std::domain_error); // x = NaN
       BOOST_MATH_CHECK_THROW(cdf(complement(w, +nan)), std::domain_error); // x = + nan
       BOOST_MATH_CHECK_THROW(quantile(w, +nan), std::domain_error); // p = + nan
@@ -511,6 +563,7 @@ void test_spots(RealType)
       BOOST_MATH_CHECK_THROW(arcsine_distribution<RealType>(1, inf), std::domain_error);
 #endif
       BOOST_MATH_CHECK_THROW(pdf(w, +inf), std::domain_error); // x = inf
+      BOOST_MATH_CHECK_THROW(logpdf(w, +inf), std::domain_error); // x = inf
       BOOST_MATH_CHECK_THROW(cdf(w, +inf), std::domain_error); // x = inf
       BOOST_MATH_CHECK_THROW(cdf(complement(w, +inf)), std::domain_error); // x = + inf
       BOOST_MATH_CHECK_THROW(quantile(w, +inf), std::domain_error); // p = + inf

test/test_exponential_dist.cpp

@@ -189,6 +189,31 @@ void test_spots(RealType T)
          static_cast<RealType>(9.0799859524969703071183031121101e-5L),                // probability.
          tolerance); // %
 
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logpdf(
+         exponential_distribution<RealType>(0.5),
+         static_cast<RealType>(0.125)),              // x
+         log(static_cast<RealType>(0.46970653140673789305985541231115L)),                // probability.
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logpdf(
+         exponential_distribution<RealType>(0.5),
+         static_cast<RealType>(5)),              // x
+         log(static_cast<RealType>(0.04104249931194939758476433723358L)),                // probability.
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logpdf(
+         exponential_distribution<RealType>(2),
+         static_cast<RealType>(0.125)),              // x
+         log(static_cast<RealType>(1.5576015661428097364903405339566L)),                // probability.
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logpdf(
+         exponential_distribution<RealType>(2),
+         static_cast<RealType>(5)),              // x
+         log(static_cast<RealType>(9.0799859524969703071183031121101e-5L)),                // probability.
+         tolerance); // %
+
    BOOST_CHECK_CLOSE(
       ::boost::math::mean(
          exponential_distribution<RealType>(2)),

test/test_gamma_dist.cpp

@@ -88,6 +88,14 @@ void check_gamma(RealType shape, RealType scale, RealType x, RealType p, RealTyp
          x),                                            // random variable.
          NaivePDF(shape, scale, x),                     // PDF
          tol);                                          // %tolerance.
+
+   // LOGPDF:
+   BOOST_CHECK_CLOSE(
+      boost::math::logpdf(
+         gamma_distribution<RealType>(shape, scale),                             // distribution.
+         x),                                                                     // random variable.
+         log(boost::math::pdf(gamma_distribution<RealType>(shape, scale), x)),   // PDF
+         tol);                                                                   // %tolerance.
 }
 
 template <class RealType>