octave-maintainers
[Top][All Lists]
Advanced

[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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]