[NumPy-Tickets] [NumPy] #1450: Patch with Ziggurat method for Normal distribution

[NumPy Trac]

#1450: Patch with Ziggurat method for Normal distribution
---------------------------------+------------------------------------------
 Reporter:  ilan                 |       Owner:  somebody
     Type:  enhancement          |      Status:  new     
 Priority:  normal               |   Milestone:          
Component:  numpy.random         |     Version:  1.4.0   
 Keywords:  normal distribution  |  
---------------------------------+------------------------------------------

Comment(by charris):

 A few additional comments.

 1) The random variable

 {{{
 u = 2.0 * rk_double(state) - 1.0
 }}}

 isn't evenly distributed between positive and negative. I don't think that
 makes any difference in practice given the other roundoff errors, but is
 worth noting because the steps in the ziggurat are symmetric about zero.
 Accepting the roundoff errors, it can be fixed by shifting and scaling the
 ziggurat layer to [0, 1-eps], using the positive double, and shifting and
 scaling the return value. This also avoids the call to fabs and saves one
 floating multiplication.

 2)
 {{{
        /* first try the rectangular boxes */
        if (fabs(u) < s_adZigR[i])
            return u * s_adZigX[i];
 }}}
 uses a double in the comparison. One of the virtues of the original
 ziggurat was using integers here, which is most likely faster even on
 modern hardware.

 3) It would be nice to move the initialization to load time so that there
 is no need to check it for every call. That isn't likely to make much
 difference given the overhead of the call to the rng but it looks a bit
 out of place.

 4) The function zig_NormalTail can be removed if an improved ziggurat is
 used.

 5) You can avoid the occasional extra call to rk_random if you handle the
 double conversion yourself. There are enough bits left over from the
 double to deal with the ziggurat layers. Having the pre-conversion integer
 form will also allow you to use integer comparisons in 2).

 Chuck

-- 
Ticket URL: <http://projects.scipy.org/numpy/ticket/1450#comment:4>
NumPy <http://projects.scipy.org/numpy>
My example project