ENH: Ziggurat for gammas

Add ziggurat for improve generation of gamma random variables

Author: bashtage

README.md

@@ -22,22 +22,24 @@ z = rnd.randn(10,10)
 * Default random generator is identical to NumPy's RandomState (i.e., 
 same seed, same random numbers).
 * Support for random number generators that support independent streams 
-and jumping ahead so that substreams can be generated
-* Faster random number generation, especially for Normals using the 
-Ziggurat method
+and jumping ahead so that sub-streams can be generated
+* Faster random number generation, especially for normal, standard
+  exponential and standard gamma using the Ziggurat method
 
 ```python
 import randomstate as rnd
 w = rnd.standard_normal(10000, method='zig')
+x = rnd.standard_exponential(10000, method='zig')
+y = rnd.standard_gamma(5.5, 10000, method='zig')
 ```
 
 * Support for 32-bit floating randoms for core generators. 
   Currently supported:
 
     * Uniforms (`random_sample`)
-    * Exponentials (`standard_exponential`) 
+    * Exponentials (`standard_exponential`, both Inverse CDF and Ziggurat)
     * Normals (`standard_normal`, both Box-Muller and Ziggurat)
-    * Standard Gammas (via `standard_gamma`) 
+    * Standard Gammas (via `standard_gamma`, both Inverse CDF and Ziggurat)
 
   **WARNING**: The 32-bit generators are **experimental** and subject 
   to change.
@@ -64,7 +66,10 @@ The RNGs include:
 * [dSFMT](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) a 
   SSE2-aware version of the MT19937 generator that is especially fast at 
   generating doubles
-* [xorshift128+](http://xorshift.di.unimi.it/), 
+* [SFMT](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) a
+  SSE2-aware version of the MT19937 generator that is optimized for
+  integer values
+* [xorshift128+](http://xorshift.di.unimi.it/),
   [xoroshiro128+](http://xoroshiro.di.unimi.it/) and
   [xorshift1024*](http://xorshift.di.unimi.it/)
 * [PCG32](http://www.pcg-random.org/) and 
@@ -76,14 +81,18 @@ The RNGs include:
 
 ### New Features
 * `standard_normal`, `normal`, `randn` and `multivariate_normal` all 
-support an additional `method` keyword argument which can be `bm` or 
-`zig` where `bm` corresponds to the current method using the Box-Muller 
-transformation and `zig` uses the much faster (100%+) Ziggurat method.
+  support an additional `method` keyword argument which can be `bm` or
+  `zig` where `bm` corresponds to the current method using the Box-Muller
+  transformation and `zig` uses the much faster (100%+) Ziggurat method.
+* `standard_exponential` and `standard_gamma` both support an additional
+  `method` keyword argument which can be `inv` or
+  `zig` where `inv` corresponds to the current method using the inverse
+  CDF and `zig` uses the much faster (100%+) Ziggurat method.
 * Core random number generators can produce either single precision
-(`np.float32`) or double precision (`np.float64`, the default) using
-an the optional keyword argument `dtype`
+  (`np.float32`) or double precision (`np.float64`, the default) using
+  an the optional keyword argument `dtype`
 * Core random number generators can fill existing arrays using the
-`out` keyword argument
+  `out` keyword argument
 
 
 ### New Functions
@@ -110,8 +119,8 @@ will produce an identical sequence of random numbers for a given seed.
   * Linux 32/64 bit, Python 2.7, 3.4, 3.5, 3.6 (probably works on 2.6 and 3.3)
   * PC-BSD (FreeBSD) 64-bit, Python 2.7
   * OSX 64-bit, Python 2.7
-  * Windows 32/64 bit (only tested on Python 2.7 and 3.5, but should 
-    work on 3.3/3.4)
+  * Windows 32/64 bit (only tested on Python 2.7, 3.5 and 3.6, but
+    should work on 3.3/3.4)
 
 ## Version
 The version matched the latest version of NumPy where 
@@ -138,7 +147,7 @@ Building requires:
   * Cython (0.22, 0.23, 0.24, 0.25)
   * tempita (0.5+), if not provided by Cython
 
-Testing requires nose (1.3+).
+Testing requires pytest (3.0+).
 
 **Note:** it might work with other versions but only tested with these 
 versions. 

README.rst

@@ -22,22 +22,26 @@ Features
 -  Default random generator is identical to NumPy's RandomState (i.e.,
    same seed, same random numbers).
 -  Support for random number generators that support independent streams
-   and jumping ahead so that substreams can be generated
--  Faster random number generation, especially for Normals using the
-   Ziggurat method
+   and jumping ahead so that sub-streams can be generated
+-  Faster random number generation, especially for normal, standard
+   exponential and standard gamma using the Ziggurat method
 
 .. code:: python
 
     import randomstate as rnd
     w = rnd.standard_normal(10000, method='zig')
+    x = rnd.standard_exponential(10000, method='zig')
+    y = rnd.standard_gamma(5.5, 10000, method='zig')
 
 -  Support for 32-bit floating randoms for core generators. Currently
    supported:
 
    -  Uniforms (``random_sample``)
-   -  Exponentials (``standard_exponential``)
+   -  Exponentials (``standard_exponential``, both Inverse CDF and
+      Ziggurat)
    -  Normals (``standard_normal``, both Box-Muller and Ziggurat)
-   -  Standard Gammas (via ``standard_gamma``)
+   -  Standard Gammas (via ``standard_gamma``, both Inverse CDF and
+      Ziggurat)
 
 **WARNING**: The 32-bit generators are **experimental** and subject to
 change.
@@ -64,6 +68,9 @@ in addition to the MT19937 that is included in NumPy. The RNGs include:
 -  `dSFMT <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/>`__ a
    SSE2-aware version of the MT19937 generator that is especially fast
    at generating doubles
+-  `SFMT <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/>`__ a
+   SSE2-aware version of the MT19937 generator that is optimized for
+   integer values
 -  `xorshift128+ <http://xorshift.di.unimi.it/>`__,
    `xoroshiro128+ <http://xoroshiro.di.unimi.it/>`__ and
    `xorshift1024\* <http://xorshift.di.unimi.it/>`__
@@ -83,6 +90,10 @@ New Features
    argument which can be ``bm`` or ``zig`` where ``bm`` corresponds to
    the current method using the Box-Muller transformation and ``zig``
    uses the much faster (100%+) Ziggurat method.
+-  ``standard_exponential`` and ``standard_gamma`` both support an
+   additional ``method`` keyword argument which can be ``inv`` or
+   ``zig`` where ``inv`` corresponds to the current method using the
+   inverse CDF and ``zig`` uses the much faster (100%+) Ziggurat method.
 -  Core random number generators can produce either single precision
    (``np.float32``) or double precision (``np.float64``, the default)
    using an the optional keyword argument ``dtype``
@@ -117,8 +128,8 @@ Status
    3.3)
 -  PC-BSD (FreeBSD) 64-bit, Python 2.7
 -  OSX 64-bit, Python 2.7
--  Windows 32/64 bit (only tested on Python 2.7 and 3.5, but should work
-   on 3.3/3.4)
+-  Windows 32/64 bit (only tested on Python 2.7, 3.5 and 3.6, but should
+   work on 3.3/3.4)
 
 Version
 -------
@@ -152,7 +163,7 @@ Building requires:
 -  Cython (0.22, 0.23, 0.24, 0.25)
 -  tempita (0.5+), if not provided by Cython
 
-Testing requires nose (1.3+).
+Testing requires pytest (3.0+).
 
 **Note:** it might work with other versions but only tested with these
 versions.

doc/source/change-log.rst

@@ -5,6 +5,20 @@ Change Log
 
 Since Version 1.13
 ------------------
+* Add Ziggurat generation for standard gamma
+  (:meth:`~randomstate.prng.mt19937.standard_gamma`) for both floats and
+  doubles.  The gamma generator uses a rejection sampler that
+  depends on random double, normal and/or exponential values.
+
+.. ipython:: python
+
+   import numpy as np
+   import randomstate as rs
+   rs.seed(23456)
+   rs.standard_gamma(3, method='zig') # New method
+
+   rs.standard_gamma(3, method='inv') # Old method
+
 * Add Ziggurat generator for standard exponential
   (:meth:`~randomstate.prng.mt19937.standard_exponential`) for both floats and
   doubles
@@ -18,7 +32,6 @@ Since Version 1.13
 
    rs.standard_exponential(3, method='inv') # Old method
 
-
 * Add SIMD-oriented Fast Mersenne Twister
   (`SFMT <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/>`_) generator.
 * Add complex normal (:meth:`~randomstate.prng.mt19937.complex_normal`)
@@ -83,19 +96,16 @@ Version 1.11.1
 
 Version 1.11
 ------------
-
 * Update to recent changes in NumPy's RandomState
 * Expose system entropy through :meth:`randomstate.entropy.random_entropy`
 * Add vector initialization for all PRNGs
 
 Version 1.10.1
 --------------
-
 * Added support for jumping the MRG32K3A generator
 * Added support for jumping the dSFMT generator
 * Update to recent changes in NumPy's RandomState
 
 Version 1.10
 ------------
-
 * This is the initial release with compatibility with NumPy 1.10

doc/source/index.rst

@@ -10,16 +10,18 @@ What's New or Different
   ``/dev/urandom`` on Unix).
 * Simulate from the complex normal distribution
   (:meth:`~randomstate.prng.mt19937.complex_normal`)
-* The both the normal and exponential generators support 256-step Ziggurat
-  methods which are 2-6 times faster than NumPy's default implementation in
-  ``standard_normal`` or ``standard_exponential``.  This Ziggurat generator can
-  be accessed by passing the keyword argument ``method='zig'``.
+* The normal, exponential and gamma generators support 256-step Ziggurat
+  methods which are 2-10 times faster than NumPy's default implementation in
+  ``standard_normal``, ``standard_exponential`` or ``standard_gamma``.
+  The Ziggurat generator can be accessed by passing the keyword
+  argument ``method='zig'``.
 
 .. ipython:: python
 
   from randomstate import standard_normal, standard_exponential
-  standard_normal(100000, method='zig')
-  standard_exponential(100000, method='zig')
+  x = standard_normal(100000, method='zig')
+  y = standard_exponential(100000, method='zig')
+  z = standard_gamma(100000, method='zig')
 
 * Optional ``dtype`` argument that accepts ``np.float32`` or ``np.float64``
   to produce either single or double prevision uniform random variables for

randomstate/distributions.c

@@ -9,6 +9,10 @@
  *
  */
 
+/* Internal Prototypes */
+static inline double standard_exponential_zig_double(aug_state* state);
+static inline float standard_exponential_zig_float(aug_state* state);
+
 static inline float random_float(aug_state* state)
 {
     return (random_uint32(state) >> 9) * (1.0f / 8388608.0f);
@@ -1558,8 +1562,6 @@ void random_bounded_bool_fill(aug_state *state, npy_bool off, npy_bool rng, npy_
     }
 }
 
-static inline double standard_exponential_zig_double(aug_state* state);
-
 static double standard_exponential_zig_double_unlikely(aug_state* state, uint8_t idx, double x)
 {
     if (idx == 0)
@@ -1605,9 +1607,7 @@ void random_standard_exponential_zig_double_fill(aug_state* state, npy_intp coun
     }
 }
 
-static inline float standard_exponential_zig_float(aug_state* state);
-
-static double standard_exponential_zig_float_unlikely(aug_state* state, uint8_t idx, float x)
+static float standard_exponential_zig_float_unlikely(aug_state* state, uint8_t idx, float x)
 {
     if (idx == 0)
     {
@@ -1653,3 +1653,121 @@ void random_standard_exponential_zig_float_fill(aug_state* state, npy_intp count
         out[i] = standard_exponential_zig_float(state);
     }
 }
+
+static inline double standard_gamma_zig_double(aug_state* state, double shape)
+{
+    double b, c;
+    double U, V, X, Y;
+
+    if (shape == 1.0)
+    {
+        return standard_exponential_zig_double(state);
+    }
+    else if (shape < 1.0)
+    {
+        for (;;)
+        {
+            U = random_double(state);
+            V = standard_exponential_zig_double(state);
+            if (U <= 1.0 - shape)
+            {
+                X = pow(U, 1./shape);
+                if (X <= V)
+                {
+                    return X;
+                }
+            }
+            else
+            {
+                Y = -log((1-U)/shape);
+                X = pow(1.0 - shape + shape*Y, 1./shape);
+                if (X <= (V + Y))
+                {
+                    return X;
+                }
+            }
+        }
+    }
+    else
+    {
+        b = shape - 1./3.;
+        c = 1./sqrt(9*b);
+        for (;;)
+        {
+            do
+            {
+                X = gauss_zig_double(state);
+                V = 1.0 + c*X;
+            } while (V <= 0.0);
+
+            V = V*V*V;
+            U = random_double(state);
+            if (U < 1.0 - 0.0331*(X*X)*(X*X)) return (b*V);
+            if (log(U) < 0.5*X*X + b*(1. - V + log(V))) return (b*V);
+        }
+    }
+}
+
+static inline float standard_gamma_zig_float(aug_state* state, float shape)
+{
+    float b, c;
+    float U, V, X, Y;
+
+    if (shape == 1.0f)
+    {
+        return standard_exponential_zig_float(state);
+    }
+    else if (shape < 1.0f)
+    {
+        for (;;)
+        {
+            U = random_float(state);
+            V = standard_exponential_zig_float(state);
+            if (U <= 1.0f - shape)
+            {
+                X = powf(U, 1.0f/shape);
+                if (X <= V)
+                {
+                    return X;
+                }
+            }
+            else
+            {
+                Y = -logf((1.0f-U)/shape);
+                X = powf(1.0f - shape + shape*Y, 1.0f/shape);
+                if (X <= (V + Y))
+                {
+                    return X;
+                }
+            }
+        }
+    }
+    else
+    {
+        b = shape - 1.0f/3.0f;
+        c = 1.0f / sqrtf(9.0f*b);
+        for (;;)
+        {
+            do
+            {
+                X = gauss_zig_float(state);
+                V = 1.0f + c*X;
+            } while (V <= 0.0f);
+
+            V = V*V*V;
+            U = random_float(state);
+            if (U < 1.0f - 0.0331f * (X*X)*(X*X)) return (b*V);
+            if (logf(U) < 0.5f * X*X + b*(1.0f - V + logf(V))) return (b*V);
+        }
+    }
+}
+
+double random_standard_gamma_zig_double(aug_state* state, double shape)
+{
+    return standard_gamma_zig_double(state, shape);
+}
+
+float random_standard_gamma_zig_float(aug_state* state, float shape)
+{
+    return standard_gamma_zig_float(state, shape);
+}

randomstate/distributions.h

@@ -170,4 +170,8 @@ extern void random_standard_exponential_zig_double_fill(aug_state* state, npy_in
 
 extern float random_standard_exponential_zig_float(aug_state* state);
 
-extern void random_standard_exponential_zig_float_fill(aug_state* state, npy_intp count, float *out);
+extern void random_standard_exponential_zig_float_fill(aug_state* state, npy_intp count, float *out);
+
+extern double random_standard_gamma_zig_double(aug_state* state, double shape);
+
+extern float random_standard_gamma_zig_float(aug_state* state, float shape);