Distributions in C

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)));
}

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: