# [Numpy-discussion] Quick Question about Optimization

Charles R Harris charlesr.harris@gmail....
Mon May 19 18:39:59 CDT 2008

```On Mon, May 19, 2008 at 4:36 PM, Robert Kern <robert.kern@gmail.com> wrote:

> On Mon, May 19, 2008 at 5:27 PM, Charles R Harris
> <charlesr.harris@gmail.com> wrote:
> > The latest versions of Matlab use the ziggurat method to generate random
> > normals and it is faster than the method used in numpy. I have ziggurat
> code
> > at hand, but IIRC, Robert doesn't trust the method ;)
>
> Well, I outlined the tests that would satisfy me, but I don't think
> you ever responded.
>
>  http://projects.scipy.org/pipermail/scipy-dev/2005-December/004405.html
>  which references
>  http://projects.scipy.org/pipermail/scipy-dev/2005-December/004400.html
>

It's been tested in the literature. It's basically just sampling a Gaussian
but done so most of the tests for being under the curve are trivial and
there are few misses, i.e., the Gaussian is covered with a stack (ziggurat)
of slices of equal areas, each slice is randomly chosen, then the position
along the slice is randomly chosen. Most of those last points will be under
the curve except at the ends, and it is those last that require computation.
However, like all sampling it has to be carefully implemented and the
samples are discretized differently than for the current way. Floats are
strange that way because they are on a log scale. The tails will be fine,
the real question is how much precision you want when doubles are returned,
i.e., how fine the discetization of the resulting samples should be. The
same method also works for the exponential distribution.

I don't feel this is a pressing issue, when I need fast normals I use my own
code. But if we are in competition with Matlab maybe we should give it a
shot.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/numpy-discussion/attachments/20080519/70ba5541/attachment.html
```