Random number generation : Von Mises distribution

[Alexandre Campo <Alexandre dot Campo at ulb dot ac dot be>]

Hi,

I have been looking for a function to generate random numbers
out of a von mises distribution (circular stats).

I did not manage to find such a thing in GSL therefore 
I adapted the R version. 

The function is provided below (GNU GPL License of course).

Cheers,
Alexandre Campo.

// int main(int argc, char** argv)
// {
//     init_rng(&rng);

//     // run 1000 times rvm ...
//     for (int i = 0; i < 100000; i++)
//     {
// 	double x = rvm(M_PI, 1.0);
// 	printf("%lf\n", x);
//     }

//     del_rng(rng);
// }


double rvm (double mean, double k) 
{
    double result = 0.0;

    double a = 1.0 + sqrt(1 + 4.0 * (k * k));
    double b = (a - sqrt(2.0 * a))/(2.0 * k);
    double r = (1.0 + b * b)/(2.0 * b);

    while (1)
    {
	double U1 = gsl_ran_flat(rng, 0.0, 1.0);
	double z = cos(M_PI * U1);
	double f = (1.0 + r * z)/(r + z);
	double c = k * (r - f);
	double U2 = gsl_ran_flat(rng, 0.0, 1.0);
	
	if (c * (2.0 - c) - U2 > 0.0) 
	{
	    double U3 = gsl_ran_flat(rng, 0.0, 1.0);
	    double sign = 0.0;
	    if (U3 - 0.5 < 0.0)
		sign = -1.0;
	    if (U3 - 0.5 > 0.0)
		sign = 1.0;
	    result = sign * acos(f) + mean;
	    while (result >= 2.0 * M_PI)
		result -= 2.0 * M_PI;
	    break;
	}
	else 
	{
	    if(log(c/U2) + 1.0 - c >= 0.0) 
	    {
		double U3 = gsl_ran_flat(rng, 0.0, 1.0);
		double sign = 0.0;
		if (U3 - 0.5 < 0.0)
		    sign = -1.0;
		if (U3 - 0.5 > 0.0)
		    sign = 1.0;
		result = sign * acos(f) + mean;
		while (result >= 2.0 * M_PI)
		    result -= 2.0 * M_PI;
		break;
	    }
	}
    }
    return result;
}

Attachment: pgp00000.pgp
Description: PGP signature