This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Other format: | [Raw text] |
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
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |