Add logcdf to distributions

include/boost/math/distributions/detail/derived_accessors.hpp

@@ -122,6 +122,13 @@ inline typename Distribution::value_type cdf(const Distribution& dist, const Rea
    typedef typename Distribution::value_type value_type;
    return cdf(dist, static_cast<value_type>(x));
 }
+template <class Distribution, class Realtype>
+inline typename Distribution::value_type logcdf(const Distribution& dist, const Realtype& x)
+{
+   using std::log;
+   using value_type = typename Distribution::value_type;
+   return log(cdf(dist, static_cast<value_type>(x)));
+}
 template <class Distribution, class RealType>
 inline typename Distribution::value_type quantile(const Distribution& dist, const RealType& x)
 {

include/boost/math/distributions/extreme_value.hpp

@@ -174,6 +174,32 @@ inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, co
    return result;
 } // cdf
 
+template <class RealType, class Policy>
+inline RealType logcdf(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::logcdf(const extreme_value_distribution<%1%>&, %1%)";
+
+   if((boost::math::isinf)(x))
+      return x < 0 ? 0.0f : 1.0f;
+   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(0 == detail::check_finite(function, a, &result, Policy()))
+      return result;
+   if(0 == detail::check_x("boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))
+      return result;
+
+   result = -exp((a-x)/b);
+
+   return result;
+} // logcdf
+
 template <class RealType, class Policy>
 RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p)
 {

include/boost/math/distributions/laplace.hpp

@@ -226,6 +226,50 @@ inline RealType cdf(const laplace_distribution<RealType, Policy>& dist, const Re
    return result;
 } // cdf
 
+template <class RealType, class Policy>
+inline RealType logcdf(const laplace_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING  // For ADL of std functions.
+
+   RealType result = 0;
+   // Checking function argument.
+   const char* function = "boost::math::logcdf(const laplace_distribution<%1%>&, %1%)";
+   // Check scale and location.
+   if (false == dist.check_parameters(function, &result)) 
+   {
+      return result;
+   }
+
+   // Special cdf values:
+   if((boost::math::isinf)(x))
+   {
+      if(x < 0) 
+      {
+         return 0; // -infinity.
+      }
+      return 1; // + infinity.
+   }
+
+   if (false == detail::check_x(function, x, &result, Policy())) 
+   {
+      return result;
+   }
+
+   // General cdf  values
+   RealType scale( dist.scale() );
+   RealType location( dist.location() );
+
+   if (x < location)
+   {
+      result = ((x - location) / scale) - boost::math::constants::ln_two<RealType>();
+   }
+   else
+   {
+      result = log(2 - exp((location - x) / scale)) - boost::math::constants::ln_two<RealType>();
+   }
+
+   return result;
+} // logcdf
 
 template <class RealType, class Policy>
 inline RealType quantile(const laplace_distribution<RealType, Policy>& dist, const RealType& p)

include/boost/math/distributions/logistic.hpp

@@ -146,6 +146,49 @@ namespace boost { namespace math {
        if(power < -tools::log_max_value<RealType>())
           return 1;
        return 1 / (1 + exp(power)); 
+    }
+
+    template <class RealType, class Policy>
+    inline RealType logcdf(const logistic_distribution<RealType, Policy>& dist, const RealType& x)
+    {
+       RealType scale = dist.scale();
+       RealType location = dist.location();
+       RealType result = 0; // of checks.
+       static const char* function = "boost::math::cdf(const logistic_distribution<%1%>&, %1%)";
+       if(false == detail::check_scale(function, scale, &result, Policy()))
+       {
+          return result;
+       }
+       if(false == detail::check_location(function, location, &result, Policy()))
+       {
+          return result;
+       }
+
+       if((boost::math::isinf)(x))
+       {
+          if(x < 0) 
+          {
+            return 0; // -infinity
+          }
+          return 1; // + infinity
+       }
+
+       if(false == detail::check_x(function, x, &result, Policy()))
+       {
+          return result;
+       }
+       BOOST_MATH_STD_USING
+       RealType power = (location - x) / scale;
+       if(power > tools::log_max_value<RealType>())
+       {
+          return 0;
+       }
+       if(power < -tools::log_max_value<RealType>())
+       {
+          return 1;
+       }
+
+       return -log(1 + exp(power));
     } 
 
     template <class RealType, class Policy>

test/test_extreme_value.cpp

@@ -23,7 +23,7 @@
    using std::setprecision;
 
 template <class RealType>
-void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealType tolerance)
+void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealType logp, RealType tolerance)
 {
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
@@ -37,6 +37,12 @@ void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealT
          x)),
          q,
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         extreme_value_distribution<RealType>(a, b),      
+         x),
+         logp,
+         tolerance); // %
    if((p < 0.999) && (p > 0))
    {
       BOOST_CHECK_CLOSE(
@@ -78,27 +84,31 @@ void test_spots(RealType)
       static_cast<RealType>(0.125), // x
       static_cast<RealType>(0.27692033409990891617007608217222L), // p
       static_cast<RealType>(0.72307966590009108382992391782778L), //q
+      static_cast<RealType>(-1.2840254166877414840734205680624364583362808652814L), // Log(p)
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
       static_cast<RealType>(2), // b
       static_cast<RealType>(-5), // x
       static_cast<RealType>(1.6087601139887776413169427645933e-7L), // p
       static_cast<RealType>(0.99999983912398860112223586830572L), //q
+      static_cast<RealType>(-15.6426318841881716102126980461566588450380350341076L), // Log(p)
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
       static_cast<RealType>(0.25), // b
       static_cast<RealType>(0.75), // x
       static_cast<RealType>(0.69220062755534635386542199718279L), // p
       static_cast<RealType>(0.30779937244465364613457800281721), //q
+      static_cast<RealType>(-0.36787944117144232159552377016146086744581113103177L), // Log(p)
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
       static_cast<RealType>(0.25), // b
       static_cast<RealType>(5), // x
       static_cast<RealType>(0.99999998477002037126351248727041L), // p
       static_cast<RealType>(1.5229979628736487512729586276294e-8L), //q
+      static_cast<RealType>(-1.52299797447126284361366292335174318621748e-8L), // Log(p)
       tolerance);
 
    BOOST_CHECK_CLOSE(
@@ -195,6 +205,9 @@ void test_spots(RealType)
    BOOST_MATH_CHECK_THROW(
        cdf(extreme_value_distribution<RealType>(0, -1), RealType(1)),
        std::domain_error);
+   BOOST_MATH_CHECK_THROW(
+       logcdf(extreme_value_distribution<RealType>(0, -1), RealType(1)),
+       std::domain_error);
    BOOST_MATH_CHECK_THROW(
        quantile(dist, RealType(-1)),
        std::domain_error);
@@ -211,6 +224,8 @@ void test_spots(RealType)
       BOOST_CHECK_EQUAL(cdf(extreme_value_distribution<RealType>(), inf), 1);
       BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution<RealType>(), -inf)), 1);
       BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution<RealType>(), inf)), 0);
+      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType>(), -inf), 0);
+      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType>(), inf), 1);
    }
    //
    // Bug reports:

test/test_laplace.cpp

@@ -225,6 +225,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.067667641618306345946999747486242201703815773119812L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(-2.L)),
+      //static_cast<RealType>(0.06766764161831L),
+      static_cast<RealType>(-2.6931471805599453094172321214581765680755001319215L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(-1.L)),
       //static_cast<RealType>(0.18393972058572L),
@@ -239,7 +245,14 @@ void test_pdf_cdf_ocatave()
 
    BOOST_CHECK_CLOSE(
       cdf(laplace_distribution<RealType>(), static_cast<RealType>(-1.L)),
-      static_cast<RealType>(0.18393972058572L),
+      //static_cast<RealType>(0.18393972058572L),
+      static_cast<RealType>(0.18393972058572116079776188508073043372290556554506L),
+   tolerance);
+
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(-1.L)),
+      //log(static_cast<RealType>(0.18393972058572L)),
+      static_cast<RealType>(-1.6931471805599453094172321214581765680755001342016L),
    tolerance);
 
    BOOST_CHECK_CLOSE(
@@ -254,6 +267,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.30326532985631671180189976749559022672095906778368L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(-0.5L)),
+      //static_cast<RealType>(0.30326532985632L),
+      static_cast<RealType>(-1.1931471805599453094172321214581765680755001342281L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(0.0L)),
       static_cast<RealType>(0.5L),
@@ -269,6 +288,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.5L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(0.0L)),
+      //log(static_cast<RealType>(0.5L)),
+      static_cast<RealType>(-0.69314718055994530941723212145817656807550013436026L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(0.5L)),
       //static_cast<RealType>(0.30326532985632L),
@@ -286,6 +311,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.69673467014368328819810023250440977327904093221632L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(0.5L)),
+      // static_cast<RealType>(0.69673467014368L),
+      static_cast<RealType>(-0.36135061480875908299180460835979393734587166765723L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(1.0L)),
     //  static_cast<RealType>(0.18393972058572L),
@@ -304,6 +335,14 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.81606027941427883920223811491926956627709443427977L),
    tolerance);
 
+   // Tolerance of naive implementation requires tol * 1e6
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(1.00000000000000L)),
+     // static_cast<RealType>(0.81606027941428L),
+      static_cast<RealType>(-0.2032670549151953327040710243903843274153282048153L),
+      //log(static_cast<RealType>(0.81606027941427883920223811491926956627709443427977L)),
+   tolerance*1e5);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(2.0L)),
     //  static_cast<RealType>(0.06766764161831L),
@@ -322,6 +361,13 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.93233235838169365405300025251375779829618422688019L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(2.0L)),
+   //   static_cast<RealType>(0.93233235838169L),
+      static_cast<RealType>(-0.07006592016028141016812825150355188141186366978811L),
+      //log(static_cast<RealType>(0.93233235838169365405300025251375779829618422688019L)),
+   tolerance*1e6);
+
    check_out_of_range<laplace_distribution<RealType> >(0, 1);
    BOOST_MATH_CHECK_THROW(laplace_distribution<RealType>(0, 0), std::domain_error);
    BOOST_MATH_CHECK_THROW(laplace_distribution<RealType>(0, -1), std::domain_error);