Re: [Bug-gsl] Lev� random number generator

[rafael]

Thanks for your quick answer, and sorry about my poor english, it is  
not my natural language.

The code below generates 10^6 random numbers, and makes a normalized  
histogram wich is compared to the levy pdf. To get the levy pdf, it  
numericaly integrates the characteristic function for levy process  
(the function f in the code). The n parameter just adds a series of  
levy numbers to get better precision. The code saves 2 files:  
lhist-$alpha-$n (the normalized histogram) and lpdf-$alpha (the pdf  
for the levy process). It also prints to the stdout the absolute error  
(square of the difference) between the histogram and the pdf.

The function levy skew shows problems for alpha<1.

With this code, you can also check the problems related to the  
gsl_ran_levy function (just change gsl_ran_levy_skew by gsl_ran_levy,  
cutting the last paramenter).

I am using the pre-compiled gsl that comes with debian etch (gsl  
version 1.8.2).

If you are interested, I also encoded a routine to generate levy skew  
random numbers, it is not fully tested, but it works for beta=0 and  
alpha<1 (it suffers from the same precision problem as gsl_ran_levy  
function for small alpha)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_integration.h>

double f (double x, void *params)
{
        double alpha = *(double *) params;
        double f = exp(-pow(x,alpha))/M_PI;
        return f;
}

double *levy_pdf(double alpha,int hist_size,double a,double b)
{
        double abserr,*lpdf,dx;
        gsl_function F;
        int i;
        gsl_integration_workspace *w=gsl_integration_workspace_alloc(1000);
        gsl_integration_workspace *cw=gsl_integration_workspace_alloc(1000);
	gsl_integration_qawo_table  
*wf=gsl_integration_qawo_table_alloc(0.0,1.0,GSL_INTEG_COSINE,200);
        lpdf=(double*)calloc(hist_size,sizeof(double));
        F.function=&f;
        F.params=&alpha;
        dx=(double)(b-a)/(double)hist_size;
        for (i=0;i<hist_size;i++)
        {
                gsl_integration_qawo_table_set(wf,i*dx+a,1.0,GSL_INTEG_COSINE);
                gsl_integration_qawf 
(&F,0.0,1e-10,1000,w,cw,wf,&lpdf[i],&abserr);
        }
        gsl_integration_qawo_table_free(wf);
        gsl_integration_workspace_free(w);
        gsl_integration_workspace_free(cw);
        return (lpdf);
}

int main (int argc,char *argv[])
{
        double *l,*lpdf,a=-20,b=20,alpha,dx,n,errabs=0.0;
        unsigned long int rnd_seed;
        int i,j,rand_numbers=1e6,hist_size=400;
        gsl_histogram *h;
        gsl_rng *r;
        struct timeval *tv;
        struct timezone *tz;
        char filename[50];
        FILE *f1,*f2;
        if(argc!=3)
        {
                printf("\nThe program must be called with 2 parameters: alpha and 
n\n \n");
                exit(1);
        }
        dx=(double)(b-a)/(double)hist_size;
        h=gsl_histogram_alloc(hist_size);
        alpha=atof(argv[1]);
        n=atof(argv[2]);
        strcpy(filename,"lhist-");
        strcat(filename,argv[1]);
        strcat(filename,"-");
        strcat(filename,argv[2]);
        f1=fopen(filename,"w+");
        strcpy(filename,"lpdf-");
        strcat(filename,argv[1]);
        f2=fopen(filename,"w+");
        l=(double*)calloc(rand_numbers,sizeof(double));
        lpdf=(double*)calloc(hist_size,sizeof(double));
        r=gsl_rng_alloc (gsl_rng_mt19937);
        gettimeofday(tv,tz);
        rnd_seed=(unsigned long int)tv->tv_usec;
        gsl_rng_set(r,rnd_seed);
        i=0;
        do
        {
                l[i]=0.0;
                for (j=1;j<n;j++) l[i]+=gsl_ran_levy_skew(r,1.0,alpha,0.0);
                l[i]=l[i]/pow(n,1.0/alpha);
		if (abs(l[i])<=20) i++;//picks only random numbers in the interval  
[a,b] to get good precision in the histogram
        }while(i<rand_numbers);
        gsl_histogram_set_ranges_uniform(h,a,b);
        for(i=0;i<rand_numbers;i++) gsl_histogram_increment(h,l[i]);
        gsl_histogram_scale(h,(double)hist_size/((b-a)*gsl_histogram_sum(h)));
        gsl_histogram_fprintf(f1,h,"%g","%g");
        lpdf=levy_pdf(alpha,hist_size,a,b);
        for (i=0;i<hist_size;i++) fprintf(f2,"%e\t%e\n",i*dx+a,lpdf[i]);
        for (i=0;i<hist_size;i++) 
errabs+=pow((gsl_histogram_get(h,i)-lpdf[i]),2.0);
        printf("%e\n",errabs/(double)hist_size);
        gsl_histogram_free(h);
        gsl_rng_free(r);
        fclose(f1);
        exit (0);
}


> At Mon, 26 Mar 2007 19:48:01 -0300,
> address@hidden wrote:
> > The Lev� skew random number generator (gsl_ran_levy_skew) does not
> > procuce a Lev� random number when beta=0 (symmetric case), and the
> > gsl_ran_levy function does not work as stated in the docs. I made some
> > histograms from 10^6 samples to check the accuracy of the algorithms,
> > by comparison agaisnt the numerical integration of the equation of
> > Lev�'s PDF. For the gsl_ran_levy function there is a good precison for
> > alpha [1,2], for alpha (0.3,1) you must sum a series of random numbers
> > to get the same precision (tipicaly 100 or more gsl_ran_levy numbers).
> > For alpha<=0.3 the algorithm does not work properly, even worse, the
> > error increases as you add more random numbers. This contradicts the
> > manual that says "the algoritm only works for alpha (0,2]". The
> > function gsl_ran_levy_skew does not produce levy random numbers when
> > beta=0, instead the pdf of the random numbers is a linear (?!?!) one.
>
> Thanks for your email.  Please can you send a small example program
> which demonstrates the problem.
>
> Note that the Levy skew generator is tested in the GSL test suite for
> several cases where beta=0 -- if you have not done so, can you run
> "make check" and confirm that it works for these cases.
>
> --
> Brian Gough
> (GSL Maintainer)
>
> Network Theory Ltd,
> Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/



----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.