Fix Von Mises random number generation algorithm to match that of Python and R.

1 files changed, 15 insertions, 14 deletions

diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c
index 528d8f332..bbf02a061 100644
--- a/numpy/random/mtrand/distributions.c
+++ b/numpy/random/mtrand/distributions.c
@@ -549,6 +549,12 @@ double rk_standard_t(rk_state *state, double df)
     return X;
 }
 
+/* Uses the rejection algorithm compared against the wrapped Cauchy
+   distribution suggested by Best and Fisher and documented in 
+   Chapter 9 of Luc's Non-Uniform Random Variate Generation.
+   http://cg.scs.carleton.ca/~luc/rnbookindex.html
+   (but corrected to match the algorithm in R and Python)
+*/
 double rk_vonmises(rk_state *state, double mu, double kappa)
 {
     double r, rho, s;
@@ -567,32 +573,27 @@ double rk_vonmises(rk_state *state, double mu, double kappa)
 
         while (1)
         {
-            U = 2*rk_double(state) - 1;
-            V = 2*rk_double(state) - 1;
+	    U = rk_double(state);
             Z = cos(M_PI*U);
             W = (1 + s*Z)/(s + Z);
             Y = kappa * (s - W);
+            V = rk_double(state);
             if ((Y*(2-Y) - V >= 0) || (log(Y/V)+1 - Y >= 0))
             {
                 break;
             }
         }
 
-        if (U < 0)
-        {
-            result = acos(W);
-        }
-        else
+	U = rk_double(state);
+
+	result = acos(W);
+        if (U < 0.5)
         {
-            result = -acos(W);
+	    result = -result;
         }
-        result += mu + M_PI;
+        result += mu;
         mod = fmod(result, 2*M_PI);
-        if (mod && (mod < 0))
-        {
-            mod += 2*M_PI;
-        }
-        return mod - M_PI;
+        return mod;
     }
 }