The GSL library has a large number of random number generators. In addition there are a large set of functions for turning uniform random numbers into discrete or continuous probabilty distributions. The random number generator algorithms vary in terms of quality numbers output, speed of execution and maximium number output. For more information see the GSL documentation. The algorithm and seed are set via environment variables, these are picked up by the ncap2
code.
Setup
The number algorithm is set by the environment variable GSL_RNG_TYPE
. If this variable isn't set then the default rng algorithm is gsl_rng_19937. The seed is set with the environment variable GSL_RNG_SEED
. The following wrapper functions in ncap2 provide information about the chosen algorithm.
gsl_rng_min()
gsl_rng_max()
Uniformly Distributed Random Numbers
gsl_rng_get(var_in)
gsl_rng_uniform_int(var_in)
gsl_rng_uniform(var_in)
gsl_rng_uniform_pos(var_in)
Below are examples of gsl_rng_get()
and gsl_rng_uniform_int()
in action.
export GSL_RNG_TYPE=ranlux export GSL_RNG_SEED=10 ncap2 -v -O -s 'a1[time]=0;a2=gsl_rng_get(a1);' in.nc foo.nc // 10 random numbers from the range 0 - 16777215 // a2=9056646, 12776696, 1011656, 13354708, 5139066, 1388751, 11163902, 7730127, 15531355, 10387694 ; ncap2 -v -O -s 'a1[time]=21;a2=gsl_rng_uniform_int(a1).sort();' in.nc foo.nc // 10 random numbers from the range 0 - 20 a2 = 1, 1, 6, 9, 11, 13, 13, 15, 16, 19 ;
The following example produces an ncap2
runtime error. This is because the chose rng algorithm has a maximium value greater than NC_MAX_INT=2147483647
; the wrapper functions to gsl_rng_get()
and gsl_rng_uniform_int()
return variable of type NC_INT
. Please be aware of this when using random number distribution functions functions from the GSL library which return unsigned int
. Examples of these are gsl_ran_geometric()
and gsl_ran_pascal()
.
export GSL_RNG_TYPE=mt19937 ncap2 -v -O -s 'a1[time]=0;a2=gsl_rng_get(a1);' in.nc foo.nc
To find the maximium value of the chosen rng algorithm use the following code snippet.
ncap2 -v -O -s 'rng_max=gsl_rng_max();print(rng_max)' in.nc foo.nc
Random Number Distributions
The GSL library has a rich set of random number disribution functions. The library also provides cumulative distribution functions and inverse cumulative distribution functions sometimes referred to a quantile functions. To see whats available on your build use the shell command ncap2 -f|grep -e _ran -e _cdf
.
The following examples all return variables of type NC_INT
defdim("out",15); a1[$out]=0.5; a2=gsl_ran_binomial(a1,30).sort(); //a2 = 10, 11, 12, 12, 13, 14, 14, 15, 15, 16, 16, 16, 16, 17, 22 ; a3=gsl_ran_geometric(a2).sort(); //a2 = 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 4, 5 ; a4=gsl_ran_pascal(a2,50); //a5 = 37, 40, 40, 42, 43, 45, 46, 49, 52, 58, 60, 62, 62, 65, 67 ;
The following all return variables of type NC_DOUBLE
;
defdim("b1",1000); b1[$b1]=0.8; b2=gsl_ran_exponential(b1); b2_avg=b2.avg(); print(b2_avg); // b2_avg = 0.756047976787 b3=gsl_ran_gaussian(b1); b3_avg=b3.avg(); b3_rms=b3.rms(); print(b3_avg); // b3_avg= -0.00903446534258 ; print(b3_rms); // b3_rms= 0.81162979889 ; b4[$b1]=10.0; b5[$b1]=20.0; b6=gsl_ran_flat(b4,b5); b6_avg=b6.avg(); print(b6_avg); // b6_avg=15.0588129413