A3.3. inverse Gaussian random variate

The details of this procedure are contained in Chapter 4, appendix d.

#include <math.h>

double drand48();
void srand48();
 
/* returns a random variate from the standard normal */
/* distribution. Taken from Numerical Recipes in C. */ 
double normal()
{
        static int iset = 0;
        static double gset;
        double fac, r, v1, v2;
 
        if (iset == 0){
                do {
                        v1 = 2.0 * drand48() - 1.0;
                        v2 = 2.0 * drand48() - 1.0;
                        r = v1*v1 + v2*v2;
                } while (r >= 1.0);
                fac = sqrt(-2.0 * log((float)r)/(float)r);
                gset = v1 * fac;
                iset = 1;
                return v2*fac;
        } else {
                iset = 0;
                return gset;
        }
}

/* returns a random variate from the Inverse Gaussian */
/* distribution. Details of the algorithm are contained in */
/* appendix 4.d. */
float travel(mu,lam)
        float mu, lam; /* model parameters */
{
        float n, v, w, c, x;
        float p;

        n = normal();
        v = n*n;
        w = mu*v;
        c = mu/(2.0*lam);

        x = mu + c*(w - sqrt(w*(4.0*lam + w)));
        p = mu/(mu + x);

        if (p > drand48()) return(x);
        else return(mu*mu/x);
}