This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

random variate from power exponential distribution: performances


First of all let me apologize for the length of the following
message. Given my deep ignorance in mathematics in general and
numerical analysis in particular, I chose to tell the whole story so
that someone wiser than me can, more easily, recognize and point out
my mistakes.


The function gsl_ran_exppow that generates pseudo-random power
exponential variates and that is actually implemented in
randist/exppow.c uses different algorithms depending on the value of the
parameter b (the so called "shape" parameter of the power exponential).
Apart from the special cases b=1 and b=2, that correspond to a Laplace
and Normal variate respectively, the function uses a transformation of
a gamma variate (Johnson M.E., 1979, Journal of Stat. Sim. and Comput.,
9,239+) if b<1 while applies a rejection method (Tadikamalla P.R., 1980,
Journal of Amer. Stat. Assoc., 75, 683+) if b>1. More precisely, here's
the table of the algo used as a function of the value of the parameter
b:

b<1   : transformation of a gamma variate [JOHNSON]
b=1   : generate a Laplace variate
1<b<2 : rejection method based on a Laplace [TADIKAMALLA]
b=2   : generate a Normal variate
b>2   : rejection method based on a Normal [TADIKAMALLA]

In the above cited paper, prof. Tadikamalla suggests that the
rejection method is more efficient for the values of b between 1 and 2
and provides evidence in terms of the time used to generate a variate
by a DEC 10 Computer.

Of course his results depended on his hardware, but also (and probably
mainly) on the specific implementations of the underlying algorithms
used to generate gamma, normal and Laplace variates.

Now I'm wondering if the implementation of Tadikamalla's suggestion
inside the present GSL library can still be considered a valid
choice. I wrote a short program to compare the various
algorithms. Here is the table of the time (in seconds) my laptop
(Linux; kernel 2.6.8-gentoo-r1, glibc2, gcc 3.3.4, gsl 1.5 ) takes
to generate 10^6 power exp. variates for different values of the
parameter b:

b value	     JOHNSON (secs)   TADIKAMALLA (secs)

0.500000	0.630000	3.810000
0.600000	1.700000	3.020000
0.700000	1.670000	2.600000
0.800000	1.600000	2.380000
0.900000	1.540000	2.260000
1.000000	0.590000	0.260000
1.100000	1.450000	2.250000
1.200000	1.450000	2.260000
1.300000	1.460000	2.250000
1.400000	1.450000	2.250000
1.500000	1.450000	2.260000
1.600000	1.440000	2.250000
1.700000	1.440000	2.260000
1.800000	1.430000	2.270000
1.900000	1.430000	2.260000
2.000000	1.420000	0.490000
2.100000	1.410000	4.610000
2.200000	1.410000	4.600000
2.300000	1.410000	4.610000
2.400000	1.400000	4.610000
2.500000	1.390000	4.620000
2.600000	1.390000	4.600000
2.700000	1.380000	4.610000
2.800000	1.380000	4.590000
2.900000	1.380000	4.590000
3.000000	1.380000	4.100000
3.200000	1.370000	4.600000
3.400000	1.360000	4.610000
3.600000	1.370000	4.600000
3.800000	1.350000	4.590000
4.000000	1.340000	4.130000

as can be seen (apart the special cases b=1 and b=2) the algo JOHNSON
constantly outperforms TADIKAMALLA. That is, the rejection method takes
longer than the transformation of gamma, irrespectively of the value of
b.

It would be nice if someone else could check my findings. To this
purpose let me explain in more details what I did:

1) cut-and-past the actual code from gsl_ran_exppow to a program,
generating two functions for the JOHNSON and TADIKAMALLA approaches

2) write a simple cli around these functions to pass values for a,b
and the number of variates

3) measure the time the two functions take using the clock() function
from libc

I attach my little program. I compiled it with NO OPTIMIZATION using

gcc -lgsl -lgslcblas testran.c -o testran

and I used a for loop from the shell to build the table reported
above.

If my analysis is confirmed, by an inspection of my code by a
numerical savvy and by results on different platforms, the actual
implementation of gsl_ran_exppow has probably to be changed.

Regards,

	Giulio.


-- 
Giulio Bottazzi <bottazzi@sssup.it>
Laboratory of Economics and Management
Sant'Anna School for Advanced Studies, Pisa, Italy
Phone: (+39)-050-883365  Fax: (+39)-050-883344
URL: http://www.sssup.it/~bottazzi/
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>

double
ran_gamma(const gsl_rng * r, const double a, const double b)
{
  
  double u = gsl_rng_uniform (r) ;
  double v = gsl_ran_gamma (r, 1/b, 1.0) ;
  double z = a * pow(v, 1/b) ;
  
  if (u > 0.5) 
    {
      return z ;
    } 
  else 
    {
      return -z ;
    }
  
}

double
ran_reject(const gsl_rng * r, const double a, const double b)
{

  if (b == 1) 
    {
      /* Laplace distribution */
      return gsl_ran_laplace (r, a) ;
    }
  else if (b < 2) 
    {
      /* Use laplace distribution for rejection method */

      double x, y, h, ratio, u ;
      
      /* Scale factor chosen by upper bound on ratio at b = 2 */

      double s = 1.4489 ; 
      do 
        {
          x = gsl_ran_laplace (r, a) ;
          y = gsl_ran_laplace_pdf (x,a) ;
          h = gsl_ran_exppow_pdf (x,a,b) ;
          ratio = h/(s * y) ;
          u = gsl_rng_uniform (r) ;
        } 
      while (u > ratio) ;
      
      return x ;
    }
  else if (b == 2)
    {
      /* Gaussian distribution */
      return gsl_ran_gaussian (r, a/sqrt(2.0)) ;
    }
  else
    {
      /* Use gaussian for rejection method */

      double x, y, h, ratio, u ;
      const double sigma = a / sqrt(2.0) ;

      /* Scale factor chosen by upper bound on ratio at b = infinity.
         This could be improved by using a rational function
         approximation to the bounding curve. */

      double s = 2.4091 ;  /* this is sqrt(pi) e / 2 */

      do 
        {
          x = gsl_ran_gaussian (r, sigma) ;
          y = gsl_ran_gaussian_pdf (x, sigma) ;
          h = gsl_ran_exppow_pdf (x, a, b) ;
          ratio = h/(s * y) ;
          u = gsl_rng_uniform (r) ;
        } 
      while (u > ratio) ;

      return x;
    }

}


int
main(int argc,char* argv[])
{

  unsigned i;
  
  /* parameters */
  double a=1;
  double b=1;
  unsigned int n=1000000;

  /* GSL RNG */
  gsl_rng * rng;
  long seed=0; /* the seed of the RNG */

  /* time */
  clock_t start, end;
  double gamma_time_used,reject_time_used;
 
  
  /* COMMAND LINE PROCESSING */
  
  /* variables for reading command line options */
  /* ------------------------------------------ */
  char opt;
  extern char *optarg;
/*   extern int optind, opterr, optopt; */
  extern int optopt;
  /* ------------------------------------------ */
    
    
  /* read the command line */    
  while((opt=getopt(argc,argv,"a:b:N:R:"))!=EOF){
	if(opt=='?'){
	    fprintf(stderr,"option %c not recognized\n",optopt);
	    return(-1);
	}
	else if(opt=='N'){
	  /*number of rv to generate*/
	  n= (unsigned) atoi(optarg);
	}
	else if(opt=='a'){
	  a=atof(optarg);
	}
	else if(opt=='b'){
	    b=atof(optarg);
	}
	else if(opt=='R'){
	  seed= (long) atoi(optarg);
	}
	else if(opt=='h'){
	    /*print short help*/
	  printf("Compare method to generate rv\n\n");
	  printf("Usage: %s [options]\n",argv[0]);
	  printf(" Options:\n");
	  printf("-N\tnumber of points [100]\n");
	  printf("-a\tset the scale [1.0]\n");
	  printf("-b\tset the exponent [2.0]\n");
	  printf("-R\tset the seed [0]\n");
	  exit(0);
	}
    }

  /* initialize RNG */
  rng = gsl_rng_alloc (gsl_rng_mt19937);
  gsl_rng_set (rng,seed);


  /* benchmarking gamma */
  start = clock();
  for (i=0;i<n;i++){
    ran_gamma(rng,a,b);
  }
  end = clock();
  gamma_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

  /* benchmarking reject */
  start = clock();
  for (i=0;i<n;i++){
    ran_reject(rng,a,b);
  }
  end = clock();
  reject_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

  /* print results */
/*   printf("#a\t\tb\t\tgamma\t\treject\n",a,b,gamma_time_used,reject_time_used); */
  printf("%f\t%f\t%f\t%f\n",a,b,gamma_time_used,reject_time_used);
  

}

Attachment: pgp00000.pgp
Description: PGP signature


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]