[Boost-bugs] [Boost C++ Libraries] #3408: gamma distribution's quantile performance

Subject: [Boost-bugs] [Boost C++ Libraries] #3408: gamma distribution's quantile performance
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2009-09-06 18:50:32

#3408: gamma distribution's quantile performance
---------------------------------+------------------------------------------
 Reporter:  chhenning@…          |       Owner:  johnmaddock 
     Type:  Tasks                |      Status:  new         
Milestone:  Boost 1.41.0         |   Component:  math        
  Version:  Boost 1.40.0         |    Severity:  Optimization
 Keywords:                       |  
---------------------------------+------------------------------------------
 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.

 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;
 }

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/3408>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.