[NumPy-Tickets] [NumPy] #1450: Patch with Ziggurat method for Normal distribution
NumPy Trac
numpy-tickets@scipy....
Sat Apr 10 11:55:23 CDT 2010
#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
More information about the NumPy-Tickets
mailing list