[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] gsl_ran_gamma wrong for certain parameters
From: |
Tino Kluge |
Subject: |
[Bug-gsl] gsl_ran_gamma wrong for certain parameters |
Date: |
Thu, 6 Apr 2006 12:25:26 +0100 |
User-agent: |
Mutt/1.4.2.1i |
gsl_ran_gamma doesn't cope with parameters a,b producing a mean of 1
and very small standard deviations. See the programme below. I'm
using fedora core precompiled packages (gsl-1.6-2 and
gsl-devel-1.6-2).
---------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
int main(char** argv, int args) {
int i;
double a,b;
gsl_rng * rgen = gsl_rng_alloc(gsl_rng_taus);
a=5E9;
b=1.0/a;
printf("a=%e, b=%e, mean=%f, std dev=%f\n",a,b,a*b,sqrt(a)*b);
for(i=1;i<10;i++){
printf("%f\n",gsl_ran_gamma(rgen,a,b));
}
return 1;
}
---------------------------------------------------------------------
It produces seemingly correct samples till about a=4E9.
$ gcc gsl_gamma_test.c -lgsl -lgslcblas -lm -o gsl_gamma_test
$ ./gsl_gamma_test
a=4.000000e+09, b=2.500000e-10, mean=1.000000, std dev=0.000016
1.000020
1.000007
1.000000
1.000017
1.000023
0.999999
0.999995
0.999998
0.999994
---------------------------------------------------------------------
But completely changes behaviour at about a=5E9.
$ gcc gsl_gamma_test.c -lgsl -lgslcblas -lm -o gsl_gamma_test
$ ./gsl_gamma_test
a=5.000000e+09, b=2.000000e-10, mean=1.000000, std dev=0.000014
0.141013
0.141001
0.141012
0.141002
0.141006
0.141009
0.141012
0.141005
0.141002
- [Bug-gsl] gsl_ran_gamma wrong for certain parameters,
Tino Kluge <=