Re: [Boost-users] [random] present library status and doc

Subject: Re: [Boost-users] [random] present library status and doc
From: Marco Guazzone (marco.guazzone_at_[hidden])
Date: 2010-02-02 04:36:37

On Tue, Feb 2, 2010 at 4:40 AM, er <erwann.rogard_at_[hidden]> wrote:
> Glad to hear a doc is coming for random.
>
>>
>> In principle this should be possible by means of the "scaling property":
>> Gamma(shape,scale) ~ scale*Gamma(shape,1)
>>
>
>>
>
> FYI, you  might recall that I started a mapping from math::distribution to
> random which I mention again because I've been refining it a bit. A
> Kolmogorov Smirnov accumulator is used for each distribution to verify that
> the mapping is correct for a particular parameter, for example,
>
> boost::mt19937 urng;
> typedef kolmovorov_smirnov::check_convergence<double> check_;
>
> math::gamma_distribution<double> d( 2, 3 );
>
> check_()(6,10,10,d,make_random_generator(urng,d),os)
>
> prints
>
> gamma(2,3)
> (10,0.195231)
> (110,0.0720463)
> (1110,0.0190203)
> (11110,0.00706733)
> (111110,0.00400383)
> (1111110,0.000770538)
>
> https://svn.boost.org/svn/boost/sandbox/statistics/distribution_toolkit/
>

So does your code solve this problem?

I've looked at your code but I fail to create a test case. I would
like to compare the result of your code with other tools like R,
octave.
For instance using this
  https://svn.r-project.org/R/trunk/src/nmath/rgamma.c
adapted accordingly in order to remove R dependencies.

--- [code] ---
#include <boost/math/distributions/gamma.hpp>
#include <boost/random/gamma_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/statistics/detail/distribution_common/meta/random/generator.hpp>
#include <boost/typeof/typeof.hpp>
#include <iostream>

template <typename T, typename GeneratorT>
T rgamma_rnd(T shape, T scale, GeneratorT rng)
{
    typedef boost::gamma_distribution<T> dist_type;

    dist_type d(shape);

    boost::variate_generator<GeneratorT&, dist_type> rvg(rng,d);

    return scale*rvg();
}

template <typename T, typename GeneratorT>
T rgamma_math(T shape, T scale, GeneratorT rng)
{
    typedef boost::math::gamma_distribution<T> dist_type;
    dist_type d(shape, scale);

    BOOST_AUTO(
        rvg,
        boost::statistics::detail::distribution::make_random_generator(rng,d)
    );

    return rvg();
}

int main()
{
    typedef double value_type;
    typedef boost::mt19937 urng_type;

    value_type shape = 1.0;
    value_type scale = 3.0;
    unsigned int seed = 5489u;

    std::cout << "Using Boost.Random ..." << ::std::endl;
    urng_type urng(seed);
    std::cout << rgamma_rnd(shape, scale, urng) << std::endl;
    std::cout << rgamma_rnd(shape, scale, urng) << std::endl;
    std::cout << rgamma_rnd(shape, scale, urng) << std::endl;

    std::cout << "Using Boost.Random + Boost.Statistic ..." << ::std::endl;
    urng = urng_type(seed);
    std::cout << rgamma_math(shape, scale, urng) << std::endl;
    std::cout << rgamma_math(shape, scale, urng) << std::endl;
    std::cout << rgamma_math(shape, scale, urng) << std::endl;
}
--- [/code] ---

Compiling with GCC 4.4.2 + Boost 1.39 + boost-statistics_rev_59420:

$ g++ -ansi -Wall -pedantic -I ./distribution_common -I
./distribution_toolkit -I /usr/include/boost -o test_gamma
test_gamma.cpp

I get these errors:

--- [error] ---
test_gamma.cpp:27:   instantiated from ‘T rgamma_math(T, T,
GeneratorT) [with T = double, GeneratorT =
boost::random::mersenne_twister<unsigned int, 32, 624, 397, 31,
2567483615u, 11, 7, 2636928640u, 15, 4022730752u, 18, 3346425566u>]’
test_gamma.cpp:52:   instantiated from here
./distribution_common/boost/statistics/detail/distribution_common/meta/random/generator.hpp:24:
error: no type named ‘type’ in ‘struct
boost::statistics::detail::distribution::meta::random_distribution<boost::math::gamma_distribution<double,
boost::math::policies::policy<boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy> > >’
./distribution_common/boost/statistics/detail/distribution_common/meta/random/generator.hpp:25:
error: no type named ‘type’ in ‘struct
boost::statistics::detail::distribution::meta::random_distribution<boost::math::gamma_distribution<double,
boost::math::policies::policy<boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy> > >’
./distribution_common/boost/statistics/detail/distribution_common/meta/random/generator.hpp:27:
error: no type named ‘type’ in ‘struct
boost::statistics::detail::distribution::meta::random_distribution<boost::math::gamma_distribution<double,
boost::math::policies::policy<boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy,
boost::math::policies::default_policy> > >’
test_gamma.cpp: In function ‘T rgamma_math(T, T, GeneratorT) [with T =
double, GeneratorT = boost::random::mersenne_twister<unsigned int, 32,
624, 397, 31, 2567483615u, 11, 7, 2636928640u, 15, 4022730752u, 18,
3346425566u>]’:
test_gamma.cpp:52:   instantiated from here
test_gamma.cpp:27: error: no matching function for call to
‘make_random_generator(boost::random::mersenne_twister<unsigned int,
32, 624, 397, 31, 2567483615u, 11, 7, 2636928640u, 15, 4022730752u,
18, 3346425566u>&, rgamma_math(T, T, GeneratorT) [with T = double,
GeneratorT = boost::random::mersenne_twister<unsigned int, 32, 624,
397, 31, 2567483615u, 11, 7, 2636928640u, 15, 4022730752u, 18,
3346425566u>]::dist_type&)’
test_gamma.cpp:27: error: no matching function for call to
‘make_random_generator(boost::random::mersenne_twister<unsigned int,
32, 624, 397, 31, 2567483615u, 11, 7, 2636928640u, 15, 4022730752u,
18, 3346425566u>&, rgamma_math(T, T, GeneratorT) [with T = double,
GeneratorT = boost::random::mersenne_twister<unsigned int, 32, 624,
397, 31, 2567483615u, 11, 7, 2636928640u, 15, 4022730752u, 18,
3346425566u>]::dist_type&)’
--- [/error] ---

Where I am wrong?

> I've only dealt with what I needed but am open to suggestions.
>
>> However, when shape == 1 the implementation uses the fact that
>> Gamma(1) ~ Exp(1).
>> This is rigth when scale==1 but not when scale != 1 since the exact
>> relation is
>> Gamma(1,scale) ~ Exp(1/scale)
>>
>
> Some textbooks define Gamma(x|a,b) as prop to x^(a-1) exp(-b x) where b is
> called the "inverse scale" such as BDA by Andrew Gelman. Others, use x^(a-1)
> exp(- x/b) such as Wikipedia. I did not think through the implications here.
> It's just a mention in passing.
>

Ok! Once this is clarified this is not a problem (it's like the use of
"inverse rate" rather than "rate" parameter for the Exponential
distribution).
The most important thing is to make 2 parameters Gamma generation
working right ;)

Best,

-- Marco