Boost logo

Boost :

Subject: Re: [boost] [math distributions]
From: Thijs van den Berg (thijs_at_[hidden])
Date: 2008-11-29 09:15:10


Paul A. Bristow wrote:
>
>> -----Original Message-----
>> From: boost-bounces_at_[hidden] [mailto:boost-bounces_at_[hidden]] On
>> Behalf Of Thijs van den Berg
>> Sent: 29 November 2008 09:59
>> To: boost_at_[hidden]
>> Subject: Re: [boost] [math distributions]
>>
>>>>> PS, I see no LambertW function in math::special_function. I'm sure
>>>>> Knuth is going to be very upset! :)
>>>>>
>
> I'm entirely ignorant of other uses this has apart from this
>
> http://www.apmaths.uwo.ca/~rcorless/frames/PAPERS/LambertW/knuth&me.jpg
>
> (And some very pretty 3D colorplots ;-)
>
> http://en.wikipedia.org/wiki/Lambert_W_function
>
> gives some clues:
>
> "The Lambert W function cannot be expressed in terms of elementary functions. It is useful in combinatorics, for instance in the enumeration of trees. It can be used to solve various equations involving exponentials and also occurs in the solution of time-delayed differential equations, such as y'(t) = a y(t − 1)."
>
> But I'm sure you are right - Knuth will be very, very disappointed.
>
> Anyone with nothing better to do fancy implementing this? :-))
>
> Paul
>
well I was thinking about implementing it using this ref
http://www.whim.org/nebula/math/lambertw.html

but I'm not sure about the complex number version (which is very
important, because it's needed for the colorplots!)

Another thing is that the multivariate Gaussian has higher priority on
my list. It's something I actually need to do quite soon (in- or
out-side boost), plus the fact that adding a new function is *way* more
work than just implementing a simple formula.

Essentially lapalace is nothing more than

T function pdf_lapace(const T& x, const T& a, const T& b) return
exp(-abs(x-a)/b)/b/2;

so I thought, I can do that in 10 minutes!
..but I ended up writing an additional 800 lines. aagghh. A little think
I could try to do is decompose the unit test into a distribution
invariant part (which is reusable) and a laplace specific part. How
about common_test.hpp, as an analogy to common_error_handling.hpp? That
would make adding a test other as simple as adding a single line. Your
idea on testing all non-member functions for the way they handle
infinity could e.g. be done like this:

// The centralized part. Testing various values, is a bit configurable
with a couple of boolean

template <class DistributionType>
template distribution_check_throw_on_inf(const DistributionType dist,
bool x_neg_inf, bool x_pos_inf, ...)
{
   typedef DistributionType::RealType realtype;

   if (x_neg_inf) //check if an error gets thrown on x == -inf
   {
      BOOST_CHECK_THROW(pdf(dist,
-std::numeric_limits<RealType>::infinity()), std::domain_error);
      BOOST_CHECK_THROW(cdf(dist,
-std::numeric_limits<RealType>::infinity()), std::domain_error);

      ... other non-member functions...

   }

   if (x_pos_inf)
   {
      ...
   }

   ... other tests ..., e.g. p=0, p=1 in quantile
}

And then have a test that uses it like this:

BOOST_AUTO_TEST_CASE( dist_test )
{
   distribution_check_throw_on_inf(normal<float>(0,1), true, true, ..
   distribution_check_throw_on_inf(normal<double>(0,1), true, true, ..
   distribution_check_throw_on_inf(normal<long double>(0,1), true, true, ..
}

we could even make a single centralized test for *all* distributions! A
author of a new distribution will have to add a couple of lines to this

BOOST_AUTO_TEST_CASE( dist_throw_test )
{
   distribution_check_throw_on_inf(normal<float>(0,1), true, true, ..
   distribution_check_throw_on_inf(normal<double>(0,1), true, true, ..
   ...

   distribution_check_throw_on_inf(laplace<float>(0,1), true, true, ..
   ...
}

something like this would make me much more confident that all tests are
performed on all distribution (nothing forgotten). What do you think? It
also makes the burden of thinking about all possible cases much less, as
the test cases are centralized...

-- 
SITMO Quantitative Financial Consultancy - Software Development
M.A. (Thijs) van den Berg
Tel.+31 (0)6 2411 0061
Fax.+31 (0)15 285 1984
thijs_at_[hidden] <mailto:thijs_at_[hidden]> - www.sitmo.com 
<http://www.sitmo.com>


Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk