[Numpy-svn] r4027 - in trunk/numpy: core random/mtrand

numpy-svn@scip... numpy-svn@scip...
Sat Sep 8 23:44:16 CDT 2007


Author: oliphant
Date: 2007-09-08 23:44:11 -0500 (Sat, 08 Sep 2007)
New Revision: 4027

Modified:
   trunk/numpy/core/numerictypes.py
   trunk/numpy/random/mtrand/distributions.c
Log:
Fix Von Mises random number generation algorithm to match that of Python and R.

Modified: trunk/numpy/core/numerictypes.py
===================================================================
--- trunk/numpy/core/numerictypes.py	2007-09-04 16:28:04 UTC (rev 4026)
+++ trunk/numpy/core/numerictypes.py	2007-09-09 04:44:11 UTC (rev 4027)
@@ -250,8 +250,8 @@
                   ('int0', 'intp'),
                   ('uint0', 'uintp'),
                   ('single', 'float'),
+                  ('csingle', 'cfloat'),
                   ('singlecomplex', 'cfloat'),
-                  ('csingle', 'cfloat'),
                   ('float_', 'double'),
                   ('intc', 'int'),
                   ('uintc', 'uint'),

Modified: trunk/numpy/random/mtrand/distributions.c
===================================================================
--- trunk/numpy/random/mtrand/distributions.c	2007-09-04 16:28:04 UTC (rev 4026)
+++ trunk/numpy/random/mtrand/distributions.c	2007-09-09 04:44:11 UTC (rev 4027)
@@ -549,6 +549,12 @@
     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 @@
 
         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)
+	U = rk_double(state);
+
+	result = acos(W);
+        if (U < 0.5)
         {
-            result = acos(W);
+	    result = -result;
         }
-        else
-        {
-            result = -acos(W);
-        }
-        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;
     }
 }
 



More information about the Numpy-svn mailing list