Random number generation : Von Mises distribution

Alexandre Campo Alexandre.Campo@ulb.ac.be
Thu Mar 9 14:42:00 GMT 2006

```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;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20060309/9405f353/attachment.sig>
```