[boost] More distributions

From: John Maddock (john_at_[hidden])
Date: 2007-12-16 05:02:28

Andrew Sutton wrote:
>> So you can call it whatever you like and still be right :-)
>>
>> BTW I believe the min and max cases are basically just mirror
>> images of each
>> other about the location parameter?
>
> So names aside, it's all the extreme value distribution, but there
> are minimum and maximum variations of it.
>
> After doing more math than I've done in some time, I'll agree to
> that :)

:-)

> Actually, the min/max case are so similar that it might be
> worth considering generalizing the extreme value distribution to
> handle both cases. It comes down to adding a third parameter, c which
> is either 1 or -1. Naturally, 1 for max, -1 for min.

I think this is along the right lines, but I don't much like the +-1 
parameter: especially as it's a floating point type argument.  How about 
using an enum here:

enum extreme_value_type
{
   type_1_maximum = 1,
   type_1_minimum = -1
};

Then the constructor looks like:

extreme_value_distribution(RealType a = RealType(0), RealType b = 
RealType(1), extreme_value_type s = type_1_maximum)

And the user can only pass valid values for the distribution type.  We also 
leave the door open for other extreme value types should the need arise.

> c(x) = (c * (a - x)) / b
> pdf(x) = exp(c(x) - exp(cx))
>
> if c == 1, then it's the maximum case. If c == -1, then it's the
> minimum because (a - x) becomes (x - a). Unfortunately, the cdf's are
> a little different...
>
> cdf(x) = {
> exp(-exp(c(x)), for c == 1
> 1 - exp(-exp(c(x)), for c== -1
> }
>
> Which is kind of fun because the cdf for the minimum is the
> complement of the cdf for the maximum (and vice versa for respective
> complements). The characteristics (outside the mean) are all the
> same. The mean can be redefined given as:
>
> mean = a + c * b * euler
>
> I don't know what the effect on the quantile functions would be,
> since I'm not sure how to compute those.

Just reverse the CDF:

Let c = sign == 1 ? log(-log(p)) ? log(-boost::math::log1p(-q, Policy()));

Then x = a - (sign/b)* c;

I hope :-)

All these need tests and sanity checks: preferably some sanity checks from 
an independent source if we can find one.  Mathematica would be a good 
choice, but I don't have access to that here :-(

> I spent the last hour rewriting the extreme_value_distribution to
> accept a sign argument, allowing it to act as either the min or max
> case of the distribution. It defaults to the max. I also figured out
> how to tweak the random number generator to generate numbers for both
> cases (also taking the additional parameter). The code is here and
> here:
>
> http://warhol.sdml.cs.kent.edu/trac/miniboost/browser/trunk/boost/
> math/distributions/extreme_value.hpp
> http://warhol.sdml.cs.kent.edu/trac/miniboost/browser/trunk/boost/
> random/extreme_value_distribution.hpp
>
> Thoughts?

Looks good so far.  Some further comments:

Median looks like it needs adjusting to reflect the "sign".
Skewness for min case is presumably the negative of the max case (needs 
checking)?

Many thanks for taking this on, I notice you've started on the geometric 
distribution as well, that's also on my list of TODO's that I've just added 
to the end of here: 
http://svn.boost.org/trac/boost/browser/sandbox/math_toolkit/libs/math/doc/sf_and_dist/issues.qbk 
so more power to your elbow I say :-)

BTW the sandbox/math_toolkit directory is open for breaking changes again 
now (Boost.Math in 1.35 is effectively frozen bar trivial fixes), so when 
you have these worked up feel free to commit to the sandbox.  Also let me 
know if you get stumped for test cases.

Cheers, John.