Opened 9 years ago

Closed 9 years ago

Last modified 4 years ago

#9672 closed Feature Requests (fixed)

PDF and CDF of a laplace distribution throwing domain_error

Reported by: HS <tan@…> Owned by: Paul A. Bristow
Milestone: Boost 1.56.0 Component: math
Version: Boost Development Trunk Severity: Showstopper
Keywords: laplace Cc:

Description

Hello,

I would expect the following to behave just like a normal distribution (or a cauchy distribution) at x=+-infinity:

// #define BOOST_MATH_DOMAIN_ERROR_POLICY ignore_error

#include <cstdlib>
#include <cstdio>

#include <limits>

#include <boost/math/distributions/cauchy.hpp>
#include <boost/math/distributions/laplace.hpp>
#include <boost/math/distributions/normal.hpp>

static const double INFTY = std::numeric_limits< double >::infinity();

template< typename T >
void TestBoundary( const char* distName, const T& t )
{
   std::printf(
         "%15s ---------------------\n"
         "%15s %'*g %'*g\n",
         distName,
         "pdf",
         boost::math::pdf( t, -INFTY ), 10,
         boost::math::pdf( t, INFTY ), 10 );

   std::printf(
         "%15s %'*g %'*g\n",
         "cdf",
         boost::math::cdf( t, -INFTY ), 10,
         boost::math::cdf( t, INFTY ), 10 );

   std::printf(
         "%15s %'*g %'*g\n",
         "co-cdf",
         boost::math::cdf( boost::math::complement( t, -INFTY ) ), 10,
         boost::math::cdf( boost::math::complement( t, INFTY ) ), 10 );
};

int main()
{
   std::printf(
         "%15s %10s %10s\n",
         "Distribution",
         "x=-infty",
         "x=+infty" );

   double _location = 3.2;
   double _scale = 0.7;

   /*
   TestBoundary(
      "Normal",
      boost::math::normal_distribution<>( _location, _scale ) );

   TestBoundary(
      "Cauchy",
      boost::math::cauchy_distribution<>( _location, _scale ) );
   */

   // To NOT throw domain_error exception at infinity?
   TestBoundary(
      "Laplace",
      boost::math::laplace_distribution<>( _location, _scale ) );

   return EXIT_SUCCESS;
}

Simply compile with "g++ filename.cpp" and obtain the following when the resulting executable is ran:

terminate called after throwing an instance of 'boost::exception_detail::clone_impl<boost::exception_detail::error_info_injector<std::domain_error> >'
what():  Error in function boost::math::pdf(const laplace_distribution<d>&, d)): Random variate x is inf, but must be finite!
Abort

If a random variable is documented to have a domain [-infty, +infty]---notice the closed interval, does it mean that the implementation will honor it (in that it would return the valid PDF or CDF at any value of the domain)---see the line above the graph of the PDF on http://www.boost.org/doc/libs/1_55_0/libs/math/doc/html/math_toolkit/dist_ref/dists/laplace_dist.html; I assumed it meant to say the "domain of the random variable" instead of "the range of the random variable" (the last expression under Non-member Accessors paragraph has it right).

Thank you.
HS

P.S. The following is from "git log -2" inside libs/math:

commit ec438ff39d855f810e236dbaf14b8755f7157869
Author: Beman <bdawes@acm.org>
Date:   Sun Dec 1 09:14:12 2013 -0500

    End-of-line normalization. Most normalization was taken care of by .gitattributes, but a few files originally committed to svn with incorrect eol tags need explicit normalization. See .gitattributes man page and extensive list discussion.

commit fb52d2de42ff7e1cbc61103a25fbc82c5e034a90
Author: John Maddock <john@johnmaddock.co.uk>
Date:   Thu Oct 24 08:35:41 2013 +0000

    Fix initialization of power series so that we don't get a spurious overflow error from tgamma when the res
ult is actually zero.
    
    [SVN r86415]

Attachments (1)

laplace.hpp (9.7 KB ) - added by Paul A. Bristow 9 years ago.

Download all attachments as: .zip

Change History (14)

comment:1 by John Maddock, 9 years ago

Owner: changed from John Maddock to Paul A. Bristow

Originally we did not allow non-finite arguments anywhere, however, in the last release that should have been - mostly - fixed to allow infinite random variables when that made sense. Paul, I think you made those changes, were these distributions included?

comment:2 by Paul A. Bristow, 9 years ago

I have reordered the checks in laplace.hpp to now allow + and - infinity.

Please can I email you a revised version for you to check that it meets your needs? Please email me at pbristow@…

comment:3 by anonymous, 9 years ago

Paul, you can attach files direct to this ticket.

by Paul A. Bristow, 9 years ago

Attachment: laplace.hpp added

in reply to:  1 comment:4 by HS <tan@…>, 9 years ago

Replying to johnmaddock:

Originally we did not allow non-finite arguments anywhere, however, in the last release that should have been - mostly - fixed to allow infinite random variables when that made sense. Paul, I think you made those changes, were these distributions included?

Dear John,

Thank you for your comment.

I actually would expect every distribution function to be defined on the whole real line. And since each distribution function is often for us a "terminal function"---eg. determining P-values or critical values of statistical tests, I would also expect its implementation to support infinities (due to overflow when calculating a test statistic) and NaN (from 0/0 say). However, I do expect the distribution parameters to take on valid finite values only---conceptually, we cannot use a "normal distribution with infinite mean"; here, I would expect an implementation to throw the domain_error like your team is doing, or---my preference---to return NaN.

But that is just my opinion as a statistician.

With best regards,
HS

in reply to:  2 comment:5 by HS <tan@…>, 9 years ago

Replying to pbristow:

I have reordered the checks in laplace.hpp to now allow + and - infinity.

Please can I email you a revised version for you to check that it meets your needs? Please email me at pbristow@…

Dear Paul,

Thank you very much for the fixes; it has cleared my test deck.

With best regards,
HS

comment:6 by Paul A. Bristow, 9 years ago

Good :-)

Useful to have a professional statistician's view on when infinity should always be supported.

In a sense, we already provide a way of achieving this view policies. If you implement a policy to ignore domain and or overflow errors, then you will get a NaN. You can achieve this now but studying the examples of using policies.

When Boost.Math was first developed, the implementation of infinity and NaN was poorly and/or inconsistently supported, so our instinct was to duck the issues.

But then John devised the 'policies' control of what happens when overflow or domain error is encountered. (This allows us to mimic what happens when you have C view of error handling and error return codes).

Also other parts of Boost.Math have been extended to make it appear that infinity and NaN are fully supported. So we have started to permit infinity in some cases, but not all because of the work involved to changes things (and the many tests) in all the many functions. It has become a bit labyrinthine as a consequence. I'll put a fuller review on my TODO list.

comment:7 by Paul A. Bristow, 9 years ago

Keywords: laplace added
Milestone: To Be DeterminedBoost 1.56.0
Resolution: fixed
Status: newclosed

Fixed and pushed to develop.

SHA-1: bf40c602967e3f87a375b3dc9b684245e270c2c1

  • Random variate can be infinite.

in reply to:  6 comment:8 by HS <tan@…>, 9 years ago

Replying to pbristow:

Good :-)

Useful to have a professional statistician's view on when infinity should always be supported.

In a sense, we already provide a way of achieving this view policies. If you implement a policy to ignore domain and or overflow errors, then you will get a NaN. You can achieve this now but studying the examples of using policies.

When Boost.Math was first developed, the implementation of infinity and NaN was poorly and/or inconsistently supported, so our instinct was to duck the issues.

But then John devised the 'policies' control of what happens when overflow or domain error is encountered. (This allows us to mimic what happens when you have C view of error handling and error return codes).

Also other parts of Boost.Math have been extended to make it appear that infinity and NaN are fully supported. So we have started to permit infinity in some cases, but not all because of the work involved to changes things (and the many tests) in all the many functions. It has become a bit labyrinthine as a consequence. I'll put a fuller review on my TODO list.

Dear Paul,

Sorry for the late response---trying to get use to modular-boost and "git submodule".

Thanks for the reminder about policies---it is a beautiful pattern and I feel it is very well implemented in boost::math.

Yes, I do use my own policies. While exceptions and/or error numbers are great tools in programming, I (as a statistician) do not find them useful here for me with distributions: for example, if any of the distribution parameters are invalid, then any caller (programmer or actual user) should get a "non-answer" (implemented by silent NaN which is wonderful since it "propagates")---no need to "leak" an exception; and each "proper" distribution (meaning all its parameters are valid) should also returns NaN when the input value is NaN.

Thanks for including +-infinity as possible input values ... they do represent the sets of "large" (and "small") numbers ... and we do want to return 0 as the P-value when a test statistic is "infinite" (or "large" in mathematical sense).

Let me know if you would like an extra pair of eyes to look over your proposed changes: I am not very good in reading code but I think I am a better reader of documents. :)

With best regards,
HS

in reply to:  2 comment:9 by anonymous, 5 years ago

Replying to Paul A. Bristow:

I have reordered the checks in laplace.hpp to now allow + and - infinity.

Please can I email you a revised version for you to check that it meets your needs? Please email me at pbristow@…

comment:10 by merciu, 5 years ago

Last edited 4 years ago by Paul A. Bristow (previous) (diff)

comment:11 by anonymous, 5 years ago

Last edited 4 years ago by Paul A. Bristow (previous) (diff)

comment:12 by anonymous, 4 years ago

Last edited 4 years ago by Paul A. Bristow (previous) (diff)

comment:13 by anonymous, 4 years ago

Last edited 4 years ago by Paul A. Bristow (previous) (diff)
Note: See TracTickets for help on using tickets.