bugGNU Scientific Library - Bugs: bug #21831, Levý random number generator for...

 
 

You are not allowed to post comments on this tracker with your current authentification level.

bug #21831: Levý random number generator for alpha < 1

Submitted by:  Brian Gough <bjg>
Submitted on:  Tue 18 Dec 2007 05:40:58 PM UTC  
 
Category: Accuracy problemSeverity: 3 - Normal
Operating System: Status: None
Assigned to: NoneOpen/Closed: Open
Release: 1.8

Tue 18 Dec 2007 05:41:30 PM UTC, comment #1:

From: -unavailable-
To: Brian Gough <bjg@network-theory.co.uk>
Cc:
Subject: Re: [Bug-gsl] Lev� random number generator
Date: Tue, 27 Mar 2007 09:35:15 -0300

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,
> -unavailable- 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.

_____________________________________________
Bug-gsl mailing list
-unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl

Brian Gough <bjg>
Project Administrator
Tue 18 Dec 2007 05:40:58 PM UTC, original submission:

From: -unavailable-
To: -unavailable-
Subject: [Bug-gsl] Levý random number generator
Date: Mon, 26 Mar 2007 19:48:01 -0300

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.

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

_____________________________________________
Bug-gsl mailing list
-unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl

Brian Gough <bjg>
Project Administrator

 

Attached Files
file #14654:  levy.c added by bjg (3kB - application/octet-stream - test program)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by bjg (Submitted the item)
  •  

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follows 1 latest change.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue 18 Dec 2007 05:42:19 PM UTCbjgAttached File-=>Added levy.c, #14654

    Back to the top


    Powered by Savane 3.1-cleanup