random variate from power exponential distribution: performances

Mon Sep 27 00:32:00 GMT 2004

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/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testran.c
Type: text/x-csrc
Size: 3823 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20040927/0807530b/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20040927/0807530b/attachment.sig>