[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: randn benchmarks
From: |
David Bateman |
Subject: |
Re: randn benchmarks |
Date: |
Fri, 23 Jan 2004 13:36:03 +0100 |
User-agent: |
Mutt/1.4.1i |
According to Paul Kienzle <address@hidden> (on 01/23/04):
> Yeah, I saw it there, but I didn't see any indication of a
> license in the ziggurat code when I looked. This is a
> user contribution rather than a base package, so it
> hasn't necessarily gone through the same scrutiny as
> others. That said, I didn't ask the package author if he
> obtained permission before using it.
>
> The license on the jstatsoft page says that the paper
> is freely distributed. It does not say that it is freely
> redistributable.
>
> It is probably okay for us to use it, but I prefer to get
> permission first.
There is also an implementation as a C++ class
http://www.codeproject.com/cpp/zigurat.asp
I don't think Marsaglia is going to reply, give that he was writing
books in the 60's there is a strong chance he has retired... You'd
probably be better of contacting teh co-author
http://www.csis.hku.hk/~tsang/
and getting permission from him. In any case, the ziggurat code in
this paper is for 2^32 bit seeds, which is a bit short for my liking.
I'd prefer to target a 64 bit seed as does matlab, so a possible
solution is just to rewrite the code from the paper for a 64 bit
seed..
In any case, before that the randn in octave-forge is really terrible.
I sent another mail last night, but forgot reply-all (so only you got
it), which is partially quoted here
>
> Ok, randn in octave-forge is slow but rand is pretty good. I tried a simple
> code to generate randn based on some fortran code on netlib to compare it
> with what is currently in octave-forge... The code snippet is
>
> {
> double u,v;
> while(1)
> {
> u = randu.randExc();
> v = randu.randExc();
> v = 1.7156 * (v - 0.5);
> double x = u - 0.449871;
> double y = abs(v) + 0.386595;
> double q = x*x + y*(0.196*y - 0.25472*x);
> if (q < 0.27597)
> break;
> if ((q <= 0.27846) && (v*v < - 4. * log(u) *u*u))
> break;
> }
> X(r,c) = v / u;
> }
>
> I've done a simple test with "a=randn(1,1e6);hist(a,[-5:0.2:5])" to see
> if its normally distributed, and to my eye it appears to be. The interesting
> thing is this
>
> before
> -----
> octave:2> tic; a = randn(1500, 1500); toc
>
> after
> -----
> octave:2> tic; a = randn(1500, 1500); toc
> ans = 0.45562
>
> matlab R12
> ----------
> >> tic; a = rand(1500, 1500); toc
>
> elapsed_time =
>
> 0.1762
>
> So, this simplistic algorithm is not as fast as matlab, but its still a good
> 4 times faster than what is currently in octave-forge.
>
> In fact even, the most basic way of generating a normal distribution
>
> {
> double sq, u, v;
> do {
> u = randu.randExc();
> v = randu.randExc();
> sq = u*u + v*v;
> } while (sq > 1.);
> X(r,c) = sqrt(-2*log(sq)/sq)*u;
> }
>
> is faster than what is in octave-forge. Is there any reason to
> use the code for randn that is already in octave-forge, over one
> of these?
>
> octave:2> tic; a = randn(1500, 1500); toc
> ans = 0.86871
>
Couldn't we at least replace the randn code in octave-forge with one of
these to generate the distribution from the Mersenne Twister code? This
would be a stop-gap measure until we can deal with the Ziggurat code...
Cheers
David
--
David Bateman address@hidden
Motorola CRM +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE
The information contained in this communication has been classified as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
- Re: benchmarks, (continued)
- Re: benchmarks - sort, Schloegl Alois, 2004/01/06
- benchmarks, John W. Eaton, 2004/01/06
- packages (was: benchmarks), John W. Eaton, 2004/01/06
- randn benchmarks, David Bateman, 2004/01/22
- Re: randn benchmarks, Paul Kienzle, 2004/01/23
- Re: randn benchmarks, Dirk Eddelbuettel, 2004/01/23
- Re: randn benchmarks, Dmitri A. Sergatskov, 2004/01/23
- Re: randn benchmarks, Dirk Eddelbuettel, 2004/01/24
- Re: randn benchmarks, Paul Kienzle, 2004/01/24
- Re: randn benchmarks, Paul Kienzle, 2004/01/24
- Re: randn benchmarks, Dmitri A. Sergatskov, 2004/01/24
- Re: randn benchmarks, David Bateman, 2004/01/25
eigenvalues 3 times speedup patch [Was: benchmarks], David Bateman, 2004/01/23