Low hanging fruit - Accelerated random distribution functions

[David Bateman]

Dear All,

I took a look at the random distribution functions as they are quite
slow with the thought of trying to use rande, randp and randg as much as
possible to accelerate them. The attached patch is the result and
converts the ones that Paul discussed in the randg function's help plus
a couple of other simple one (exprnd, poissrnd).  In particular it
converts betarnd, chi2rnd, exprnd, frnd, gamrnd, poissrnd and trnd.
Furthermore the functions binornd,  cauchy_rnd, discrete_rnd,
empirical_rnd, normrnd and unifrnd were already reasonably good
implementations.  I haven't tried to convert the functions geornd,
hygernd, laplace_rnd, logistic_rnd, lognrnd, pascal_rnd, stdnormal_rnd,
wblrnd and wienrnd yet, so if someone has any suggests then please
propose them.

This patch however raises a question in that it will alter the sequences
of the converted functions. This was already effectively done at the
time we introduced the Mersenne Twister and Ziggurat into octave and no
one complained then. However I believe that before this patch is applied
the fact that it will alter the sequences for a particular state or seed
must be known.. Note that typical speed ups with this patch are a factor
of 100...

Regards
David

-- 
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/frnd.m.orig43    2006-10-10 
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/frnd.m   2007-02-22 17:06:19.785403761 
+0100
***************
*** 80,86 ****
  
    if (isscalar (m) && isscalar (n))
      if ((m > 0) && (m < Inf) && (n > 0) && (n < Inf))
!       rnd =  finv (rand (sz), m, n);
      else
        rnd = NaN * ones (sz);
      endif
--- 80,86 ----
  
    if (isscalar (m) && isscalar (n))
      if ((m > 0) && (m < Inf) && (n > 0) && (n < Inf))
!       rnd = n ./ m .* randg(m/2,sz) ./ randg(n/2,sz);
      else
        rnd = NaN * ones (sz);
      endif
***************
*** 96,102 ****
      k = find ((m > 0) & (m < Inf) &
                (n > 0) & (n < Inf));
      if (any (k))
!       rnd(k) = finv (rand (size (k)), m(k), n(k));
      endif
    endif
  
--- 96,102 ----
      k = find ((m > 0) & (m < Inf) &
                (n > 0) & (n < Inf));
      if (any (k))
!       rnd(k) = n(k) ./ m(k) .* randg(m(k)./2,size(k)) ./ 
randg(n(k)./2,size(k));
      endif
    endif
  
*** ./scripts/statistics/distributions/exprnd.m.orig43  2006-10-10 
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/exprnd.m 2007-02-22 17:06:38.812439582 
+0100
***************
*** 69,75 ****
  
    if (isscalar (l))
      if ((l > 0) && (l < Inf))
!       rnd = - log (1 - rand (sz)) ./ l;
      else
        rnd = NaN * ones (sz);
      endif
--- 69,75 ----
  
    if (isscalar (l))
      if ((l > 0) && (l < Inf))
!       rnd = rande(sz) / l;
      else
        rnd = NaN * ones (sz);
      endif
***************
*** 81,87 ****
      endif
      k = find ((l > 0) & (l < Inf));
      if (any (k))
!       rnd(k) = - log (1 - rand (size (k))) ./ l(k);
      endif
    endif
  
--- 81,87 ----
      endif
      k = find ((l > 0) & (l < Inf));
      if (any (k))
!       rnd(k) = rande(size(k)) / l(k);
      endif
    endif
  
*** ./scripts/statistics/distributions/gamrnd.m.orig43  2006-10-10 
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/gamrnd.m 2007-02-22 17:07:33.587664666 
+0100
***************
*** 82,88 ****
      if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
        rnd = NaN * ones (sz);
      else
!       rnd =  gaminv (rand (sz), a, b);
      endif
    else 
      k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
--- 82,88 ----
      if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
        rnd = NaN * ones (sz);
      else
!       rnd = randg(a,sz)/b;
      endif
    else 
      k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
***************
*** 91,97 ****
      endif
      k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
      if (any (k))
!       rnd(k) = gaminv (rand (size (k)), a(k), b(k));
      endif
    endif
  
--- 91,97 ----
      endif
      k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
      if (any (k))
!       rnd(k) = randg(a(k),size(k))/b(k);
      endif
    endif
  
*** ./scripts/statistics/distributions/trnd.m.orig43    2006-10-10 
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/trnd.m   2007-02-22 17:23:34.021050514 
+0100
***************
*** 70,76 ****
      if (!(n > 0) || !(n < Inf))
        rnd = NaN * ones (sz);
      elseif ((n > 0) && (n < Inf))
!       rnd = tinv (rand (sz), n);
      else
        rnd = zeros (size (n));
      endif
--- 70,76 ----
      if (!(n > 0) || !(n < Inf))
        rnd = NaN * ones (sz);
      elseif ((n > 0) && (n < Inf))
!       rnd = randn(sz) ./ sqrt(2*randg(n/2,sz)./n); 
      else
        rnd = zeros (size (n));
      endif
***************
*** 84,90 ****
  
      k = find ((n > 0) & (n < Inf));
      if (any (k))
!       rnd(k) = tinv (rand (size (k)), n(k));
      endif
    endif
  
--- 84,90 ----
  
      k = find ((n > 0) & (n < Inf));
      if (any (k))
!       rnd(k) = randn(size(k)) ./ sqrt(2*randg(n(k)/2,size(k))./n(k)); 
      endif
    endif
  
*** ./scripts/statistics/distributions/poissrnd.m.orig43        2006-10-10 
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/poissrnd.m       2007-02-22 
17:16:03.800856596 +0100
***************
*** 69,86 ****
      if (!(l >= 0) | !(l < Inf))
        rnd = NaN * ones (sz);
      elseif ((l > 0) & (l < Inf))
!       num = zeros (sz);
!       sum = - log (1 - rand (sz)) ./ l;
!       while (1)
!       ind = find (sum < 1);
!       if (any (ind))
!           sum(ind) = (sum(ind) - log (1 - rand (size (ind))) / l);
!           num(ind) = num(ind) + 1;
!       else
!           break;
!       endif
!       endwhile
!       rnd = num;
      else
        rnd = zeros (sz);
      endif
--- 69,75 ----
      if (!(l >= 0) | !(l < Inf))
        rnd = NaN * ones (sz);
      elseif ((l > 0) & (l < Inf))
!       rnd = randp(l, sz);
      else
        rnd = zeros (sz);
      endif
***************
*** 94,113 ****
  
      k = find ((l > 0) & (l < Inf));
      if (any (k))
!       l = l(k);
!       num = zeros (size (k));
!       sum = - log (1 - rand (size (k))) ./ l;
!       while (1)
!       ind = find (sum < 1);
!       if (any (ind))
!           sum(ind) = (sum(ind)
!                       - log (1 - rand (size (ind))) ./ l(ind));
!           num(ind) = num(ind) + 1;
!       else
!           break;
!       endif
!       endwhile
!       rnd(k) = num;
      endif
    endif
  
--- 83,89 ----
  
      k = find ((l > 0) & (l < Inf));
      if (any (k))
!       rnd(k) = randp(l(k), size(k));
      endif
    endif
  
*** ./scripts/statistics/distributions/chi2rnd.m.orig43 2006-10-10 
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/chi2rnd.m        2007-02-22 
13:48:59.311691304 +0100
***************
*** 69,75 ****
       if (find (!(n > 0) | !(n < Inf)))
         rnd = NaN * ones (sz);
       else
!        rnd =  chi2inv (rand (sz), n);
       endif
    else
      [retval, n, dummy] = common_size (n, ones (sz));
--- 69,75 ----
       if (find (!(n > 0) | !(n < Inf)))
         rnd = NaN * ones (sz);
       else
!        rnd = 2 * randg(n/2, sz)
       endif
    else
      [retval, n, dummy] = common_size (n, ones (sz));
***************
*** 85,91 ****
  
      k = find ((n > 0) & (n < Inf));
      if (any (k))
!       rnd(k) = chi2inv (rand (size (k)), n(k));
      endif
    endif
  
--- 85,91 ----
  
      k = find ((n > 0) & (n < Inf));
      if (any (k))
!       rnd(k) = 2 * randg(n(k)/2, size(k))
      endif
    endif
  
*** ./scripts/statistics/distributions/betarnd.m.orig43 2006-10-10 
18:10:29.000000000 +0200
--- ./scripts/statistics/distributions/betarnd.m        2007-02-22 
13:35:32.725157769 +0100
***************
*** 79,85 ****
      if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
        rnd = NaN * ones (sz);
      else
!       rnd = betainv (rand(sz), a, b);
      endif
    else
      rnd = zeros (sz);
--- 79,86 ----
      if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
        rnd = NaN * ones (sz);
      else
!       r1 = randg(a,sz); 
!       rnd = r1 ./ (r1 + randg(b,sz));
      endif
    else
      rnd = zeros (sz);
***************
*** 91,97 ****
  
      k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
      if (any (k))
!       rnd(k) = betainv (rand (size (k)), a(k), b(k));
      endif
    endif
  
--- 92,99 ----
  
      k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
      if (any (k))
!       r1 = randg(a(k),size(k)); 
!       rnd(k) = r1 ./ (r1 + randg(b(k),size(k)));
      endif
    endif
  

2007-02-22  David Bateman  <address@hidden>

        * frnd.m, exprnd.m, gamrnd.m, trnd.m, poissrnd.m, chi2rnd.m,
        betarnd.m: Convert to use randg, rande and randp to accelerate.