GNU Scientific Library - Bugs: bug #21831, Levý random number generator for...
You are not allowed to post comments on this tracker with your current authentication level.
bug #21831: Levý random number generator for alpha < 1
Submitter: | -Deleted Account- <bjg> | ||
Submitted: | Tue 18 Dec 2007 05:40:58 PM UTC | ||
Category: | Accuracy problem | Severity: | 3 - Normal |
Operating System: | Status: | None | |
Assigned to: | None | Open/Closed: | Open |
Release: | 1.8 |
Tue 18 Dec 2007 05:41:30 PM UTC, comment #1: |
-Deleted Account- <bjg> |
Tue 18 Dec 2007 05:40:58 PM UTC, original submission:
From: -email is unavailable- |
-Deleted Account- <bjg> |
Depends on the following items: None found
Items that depend on this one: None found
Follows 1 latest change.
Date | Changed by | Updated Field | Previous Value | => | Replaced by |
---|---|---|---|---|---|
2007-12-18 | bjg | Attached File | - | Added levy.c, #14654 |
Powered by Savane 3.13-f8d8.
Corresponding source code
From: -email is 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=α
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);
}
----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.
_____________________________________________
Bug-gsl mailing list
-email is unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl