Re: Low hanging fruit - Accelerated random distribution functions

[David Bateman]

Paul Kienzle wrote:
>>
>>
>>> wblrnd
>>
>> Fast enough:
>>
>>   rnd = scale * -log(1-rand(k))) .^ (1/shape)
>
> 1-U and U are equivalent generators, so this should be
> equivalent to:
>
>    rnd = scale * rande(size) ^ (1/shape)
>
> which will be a bit faster.
>
> Note that simplifying the expression on the wikipedia page
> for the exponential distribution relating it to weibull
> distribution gives the same result.
>
>
Yes it is faster by a reasonable factor..

octave:8> function rnd = wblrnd2(scale,shape,varargin), rnd = scale .*
rande(varargin{:}) .^ (1./shape); endfunction
octave:9> tic; b1 = wblrnd2(2,3,1,1e7); toc
Elapsed time is 3.379036 seconds.
octave:10> tic; b0 = wblrnd(2,3,1,1e7); toc
Elapsed time is 5.163300 seconds.

Patch to apply this attached...

D.


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

*** ./scripts/statistics/distributions/wblrnd.m.orig45  2007-02-23 
11:32:20.435974132 +0100
--- ./scripts/statistics/distributions/wblrnd.m 2007-02-23 11:32:00.337102480 
+0100
***************
*** 79,85 ****
  
    if (isscalar (shape) && isscalar (scale))
      if ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf))
!       rnd = (scale * (- log (1 - rand (sz))) .^ (1 / shape));
      else
        rnd = NaN * ones (sz);
      endif
--- 79,85 ----
  
    if (isscalar (shape) && isscalar (scale))
      if ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf))
!       rnd = scale .* rande(sz) .^ (1./shape);
      else
        rnd = NaN * ones (sz);
      endif
***************
*** 87,94 ****
      rnd = NaN * ones (sz);
      k = find ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf));
      if (any (k))
!       rnd(k) = (scale(k)
!               .* (- log (1 - rand (size (k)))) .^ (1 ./ shape(k)));
      endif
    endif
  
--- 87,93 ----
      rnd = NaN * ones (sz);
      k = find ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf));
      if (any (k))
!       rnd(k) = scale(k) .* rande(size(k)) .^ (1./shape(k));
      endif
    endif