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