Next: , Previous: GSL statistics, Up: ncap2 netCDF Arithmetic Processor


4.1.23 GSL random number generation

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()
the minimium value returned by the rng algorithm.
gsl_rng_max()
the maximium value returned by the rng algorithm.

Uniformly Distributed Random Numbers

gsl_rng_get(var_in)
This function returns var_in with integers from the chosen rng algorithm. The min and max values depend uoon the chosen rng algorthm.
gsl_rng_uniform_int(var_in)
This function returns var_in with random integers from 0 to n-1. The value n must be less than or equal to the maximium value of the chosen rng algorithm.
gsl_rng_uniform(var_in)
This function returns var_in with double precision numbers in the range [0.0,1). The range include 0.0 but excludes 1.0.
gsl_rng_uniform_pos(var_in)
This function returns var_in with double precision numbers in the range (0.0,1) - excluding both 0.0 and 1.0.

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