[Boost-users] [math] gamma distribution's quantile performance

Subject: [Boost-users] [math] gamma distribution's quantile performance
From: Christian Henning (chhenning_at_[hidden])
Date: 2009-09-04 18:16:43

Hi there, I have created a test which times the performance of
boost::math::quantile( ... ) when using a gamma distribution. I ran it
against source code we use here at work, for ages. It can be found
here:

http://people.sc.fsu.edu/~burkardt/cpp_src/dcdflib/dcdflib.html

The old source code is about 2x times faster than the boost version.

MS Visual Studio 2005:
boost: 35.4sec
att_bell:19sec

Intel 11.1:
boost: 21.4sec
att_bell: 11.2sec

Question is if there is a way to incorporate such a function into
boost::math? As far, as I can tell the results are almost identical.

Regards,
Christian

Here the code:

#include <dcdflib.h>

#include <boost/math/distributions/gamma.hpp>

#include <boost/timer.hpp>

double min_mean =      2000; //      2,000
double max_mean = 500000000; //500,000,000

double min_std =     10000;  //      10,000
double max_std = 100000000;  // 100,000,000

double min_max =  600000000; //   600,000,000
double max_max = 1000000000; // 1,000,000,000

const std::size_t max_year = 5000000; // 5,000,000

const double right = 0.999;
const double left  = 0.001;

inline
double get_rand()
{
    return static_cast< double >( std::rand() )
         / static_cast< double >( RAND_MAX    );
}

inline
void boost_( boost::math::gamma_distribution<>& d, double q )
{
    double value = boost::math::quantile( d, q );
}

inline
void att_bell( double alpha, double beta, double q )
{
        double q_Minus1 = 1 - q;
        double value = 0.0;
        double bound = 0.0;

        int which  = 2;
        int status = 0;

    cdfgam( &which
          , &q
          , &q_Minus1
          , &value
          , &alpha
          , &beta
          , &status
          , &bound
          );
}

int main()
{
    // boost
    {
        std::srand( 0 );

        boost::timer timer;
        for( std::size_t y = 0; y < max_year; ++y )
        {
            if(( y % 100000 ) == 0 )
                std::cout << y << std::endl;

            double mean = get_rand() * ( max_mean - min_mean ) + min_mean;
            double std  = get_rand() * ( max_std  - min_std  ) + min_std;

            double alpha = mean * mean / std / std; // shape parameter
            double beta  = mean / alpha; // scale parameter

            boost::math::gamma_distribution<> d( alpha, beta );
            boost_( d, right );
            boost_( d, left  );
        }

        std::cout << "Boost - Time elapsed: " << timer.elapsed() << "
sec" << std::endl;
    }

    // att bell
    {
        std::srand( 0 );

        boost::timer timer;
        for( std::size_t y = 0; y < max_year; ++y )
        {
            if(( y % 100000 ) == 0 )
                std::cout << y << std::endl;

            double mean = get_rand() * ( max_mean - min_mean ) + min_mean;
            double std  = get_rand() * ( max_std  - min_std  ) + min_std;

            double alpha = mean * mean / std / std; // shape parameter
            double beta  = mean / alpha; // scale parameter

            att_bell( alpha, beta, right );
            att_bell( alpha, beta, left  );
        }

        std::cout << "ATT Bell - Time elapsed: " << timer.elapsed() <<
" sec" << std::endl;
    }

    return 0;
}