[random] present library status and doc
Hi, What is the present library status of the random library,and in how far are the documentation and the code base out-of-sync ? I just tried a variate_generator with mt19937 as engine and a normal_distribution as distribution. IIRC, then some years ago this did not work as the normal_distribution requires an input in the [0,1] range, but now it seems to work. Looking into the header files I concluded under the hood there has been quite some development going on. Also for uniform_01 things apparently changed quite a lot, there is now also for example a new_uniform_01 etc. My main questions: Is the library presently under active development ? Shall I rely on the docs, or not ? For example, is above normal distribution now deliberately working, or am I relying on a deliberately undocumented feature which might be removed later ? Because many more distributions now exist I think the docs are just totally out of date. Any plans of updating them ? cheers, Thomas
AMDG Thomas Mang wrote:
What is the present library status of the random library,and in how far are the documentation and the code base out-of-sync ?
I just tried a variate_generator with mt19937 as engine and a normal_distribution as distribution. IIRC, then some years ago this did not work as the normal_distribution requires an input in the [0,1] range, but now it seems to work. Looking into the header files I concluded under the hood there has been quite some development going on.
variate_generator is supposed to work with any generator and any distribution.
Also for uniform_01 things apparently changed quite a lot, there is now also for example a new_uniform_01 etc.
It hasn't changed much. The mess there is to make uniform_01 consistent with the other distributions without breaking backwards compatibility.
My main questions: Is the library presently under active development ?
Mostly bug fixes. I'd like to sync the library with the standard at some point, but that's likely to be a while.
Shall I rely on the docs, or not ? For example, is above normal distribution now deliberately working, or am I relying on a deliberately undocumented feature which might be removed later ?
If I recall correctly, the docs are mostly correct although incomplete.
Because many more distributions now exist I think the docs are just totally out of date. Any plans of updating them ?
Yes. There's a draft at http://www.cs.hmc.edu/~swatanabe/random_doc/libs/random/doc/html/ In Christ, Steven Watanabe
On Mon, Feb 1, 2010 at 4:46 PM, Steven Watanabe <watanabesj@gmail.com> wrote:
AMDG
Thomas Mang wrote:
[cut]
Yes. There's a draft at http://www.cs.hmc.edu/~swatanabe/random_doc/libs/random/doc/html/
Good news! But for what regards Gamma distribution I would to point out that actually the code used to generate random Gamma variates cannot be used to generate variates for two-parameters Gamma distributions (e.g., http://en.wikipedia.org/wiki/Gamma_distribution or http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm). In principle this should be possible by means of the "scaling property": Gamma(shape,scale) ~ scale*Gamma(shape,1) This would simply translate into boost::mt19937 rng; boost::gamma_distribution<> gamma(shape); // 1-parameter Gamma distribution boost::variate_generator<boost::mt19937&, boost::gamma_distribution<>
rvg(rng, gamma); double r = scale*rvg();
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) It is possibile to fix this problem? Also, it seems that the C++0x is rethinking the way to model the Gamma distribution as a two parameter distribution (see Sec. 26.5.8.3.3 [rand.dist.pois.gamma] at http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2009/n3000.pdf). Maybe it is the case that also Boost think about it, do you? Best, -- Marco
In Christ, Steven Watanabe
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
AMDG Marco Guazzone wrote:
But for what regards Gamma distribution I would to point out that actually the code used to generate random Gamma variates cannot be used to generate variates for two-parameters Gamma distributions (e.g., http://en.wikipedia.org/wiki/Gamma_distribution or http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm).
In principle this should be possible by means of the "scaling property": Gamma(shape,scale) ~ scale*Gamma(shape,1)
This would simply translate into boost::mt19937 rng; boost::gamma_distribution<> gamma(shape); // 1-parameter Gamma distribution boost::variate_generator<boost::mt19937&, boost::gamma_distribution<>
rvg(rng, gamma);
double r = scale*rvg();
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)
It is possibile to fix this problem? Also, it seems that the C++0x is rethinking the way to model the Gamma distribution as a two parameter distribution (see Sec. 26.5.8.3.3 [rand.dist.pois.gamma] at http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2009/n3000.pdf). Maybe it is the case that also Boost think about it, do you?
Please file a trac ticket. In Christ, Steven Watanabe
On Mon, Feb 1, 2010 at 5:59 PM, Steven Watanabe <watanabesj@gmail.com> wrote:
AMDG
Marco Guazzone wrote:
But for what regards Gamma distribution I would to point out that actually the code used to generate random Gamma variates cannot be used to generate variates for two-parameters Gamma distributions (e.g., http://en.wikipedia.org/wiki/Gamma_distribution or http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm).
In principle this should be possible by means of the "scaling property": Gamma(shape,scale) ~ scale*Gamma(shape,1)
This would simply translate into boost::mt19937 rng; boost::gamma_distribution<> gamma(shape); // 1-parameter Gamma distribution boost::variate_generator<boost::mt19937&, boost::gamma_distribution<>
rvg(rng, gamma);
double r = scale*rvg();
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)
It is possibile to fix this problem? Also, it seems that the C++0x is rethinking the way to model the Gamma distribution as a two parameter distribution (see Sec. 26.5.8.3.3 [rand.dist.pois.gamma] at http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2009/n3000.pdf). Maybe it is the case that also Boost think about it, do you?
Please file a trac ticket.
Just done! See https://svn.boost.org/trac/boost/ticket/3888 Thanks!!! Best, -- Marco
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/ 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.
On Tue, Feb 2, 2010 at 4:40 AM, er <erwann.rogard@gmail.com> 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
So does your code solve this problem?
The current mapping applies a scaling factor to boost::gamma_distribution and apparently the resulting empirical cdf converges to the true cdf...
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,
Sure, double checking is a good idea. Among your included directories you should have: #include <boost/statistics/detail/distribution_toolkit/distributions/gamma/include.hpp> OR #include <boost/statistics/detail/distribution_toolkit/distributions/gamma/random.hpp> These have to be in the search path. /sandbox/statistics/non_parametric/ /sandbox/statistics/distribution_common/ /sandbox/statistics/distribution_toolkit/ /sandbox/statistics/random/ /usr/local/boost_1_41_0/ but I have just removed a deprecated dependency so please check out again. To be certain, here's a small subset that I just compiled: #include <boost/typeof/typeof.hpp> #include <boost/random/mersenne_twister.hpp> #include <boost/statistics/detail/distribution_toolkit/distributions/gamma/include.hpp> #include <boost/statistics/detail/distribution_common/meta/random/generator.hpp> void example_gamma(){ namespace stat = boost::statistics::detail; namespace dist = stat::distribution; typedef double val_; typedef boost::mt19937 urng_; typedef boost::math::gamma_distribution<val_> dist_; const val_ shape = 2.0; const val_ scale = 3.0; dist_ d( shape, scale ); urng_ urng; BOOST_AUTO( vg, dist::make_random_generator(urng,d) ); } Please let me know if you have any other issues.
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.
The code is a bit cryptic, but is it possible that the 2 different algorithms for a<1 and a>1 are for efficiency reasons? This is probably obvious to you but I'll say it anway : I wouldn't expect, nor is it required, that the generated values agree but their sample distribution, however, should be the same.
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; }
I ran these and they agree (PS: I think you meant to pass the generator by reference, not value).
In principle this should be possible by means of the "scaling property": Gamma(shape,scale) ~ scale*Gamma(shape,1)
This would simply translate into boost::mt19937 rng; boost::gamma_distribution<> gamma(shape); // 1-parameter Gamma distribution boost::variate_generator<boost::mt19937&, boost::gamma_distribution<>
rvg(rng, gamma); double r = scale*rvg();
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)
An example of this would be gamma(1,0.5) or gamma(1,2.0)? FYI I ran these kolmogorov-smirnov statistics which seem to converge to zero: gamma(1,0.5) (10,0.0974325) (110,0.0464457) (1110,0.0173321) (11110,0.0108157) (111110,0.00259073) (1111110,0.000503333) gamma(1,2) (10,0.0974325) (110,0.0464457) (1110,0.0173321) (11110,0.0108157) (111110,0.00259073) (1111110,0.000503333) To rule out a fluke I've compared make_random_generator(urng,d1) to d where d and d1 only differ by 0.01 in their shape xor scale, and the statistic is larger by a factor of 10 and 100 respectively, at the last iteration... HTH
participants (4)
-
er -
Marco Guazzone -
Steven Watanabe -
Thomas Mang