On this site you will find different implementations of probability distributions in C, used in the SILAS-model. Some of them build up on each other and the original sources and links are provided. A seed for the random-number generation can be set for all distributions.
Standard normal/Gaußian distribution
The C code for the standard-normal distribution uses the polar method by Marsaglia and Bray (1964)
double rnormal(int seed){
double u1;
double u2;
double w;
double mult;
double x1;
srand((unsigned)seed);
do {
u1 = -1 + ((double) rand()/ (double) RAND_MAX)*2;
u2 = -1 + ((double) rand()/ (double) RAND_MAX)*2;
w = pow(u1, 2) + pow(u2, 2);
}
while(w>=1 || w==0);
mult = sqrt((-2*log(w))/w);
x1 = u1*mult;
return((double)x1);
}
Gamma distribution
The receiving a gamma-distribted random variable builds up on the standard normal and is based on the article of Marsaglia and Tsang (2000).
double rgamma(int seed, double a){
double rnormal(int seed);
double d;
double c;
double x;
double v;
double u;
double bound;
d = a-1/3;
c = 1 / sqrt(9*d);
bound = -( sqrt(9*a-3) );
for(;;){
do{
x = rnormal((int)seed);
seed++;
}
while(x <= bound);
v = 1+c*x;
v = v*v*v;
u = ((double) rand()/ (double) (RAND_MAX) );
if(u < 1-0.0331*(x*x)*(x*x)){
return (d*v);}
if(log(u) < 0.5*x*x+d*(1-v+log(v))){
return (d*v);}
}
}
Beta distribution
Generating a beta-distributed random variable is using two gammas, which we draw from Georgii but can also be found on the wikipedia-page of the beta-distribution. The parametrization is identical to the gamlss-package.
double rbeta(int seed, double a, double b){
double rgamma();
double x;
double y;
double z;
if(a < 1){
a = (1 + a);
x = rgamma(seed, a);
x = x*pow(((double)rand()/(double)(RAND_MAX)),1/a);
}
else{x = rgamma(seed, a);}
if(b < 1){
b = (1 + b);
y = rgamma(seed*seed, b);
y = y*pow(((double)rand()/(double)(RAND_MAX)),1/a);
}
else {y = rgamma(seed*seed, b);}
z = x/(x+y);
x = 0;
y = 0;
return z;
}
Beta-inflated distribution
As we needed not only a regular bet-distribution with support from (0,1), but also a 0 and 1 inflated beta distribution, we applied the following code according to Ospina and Ferrari (2010) (where c in our code is P1 in the paper and d is P0). This parametrization is also used in the gamlss-package.
double rbetainf(int seed, double a, double b, double c, double d){
double rbeta();
double beta;
double bern;
double u;
beta = rbeta(seed, a, b);
srand((unsigned)seed);
u = ((double) rand()/ (double) (RAND_MAX) );
if(u (1-d)){bern = 0;}
else{bern = (1-c-d)*beta;}
}
return bern;
}
Weibull distribution
For all survival-type like analyses (e.g. duration of relationships) we used a Weibull-Distribution (http://www.taygeta.com/random/weibull.html). Again this parametrization is also used in the gamlss-package.
double rweibull(int seed, double wscale, double wshape){
srand((unsigned)seed);
double zufall;
zufall = ((double) rand()/ (double) (RAND_MAX) );
return(wscale*pow((-log(zufall)),(1/wshape)));
}