[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Low hanging fruit - Accelerated random distribution functions
From: |
David Bateman |
Subject: |
Low hanging fruit - Accelerated random distribution functions |
Date: |
Thu, 22 Feb 2007 17:46:57 +0100 |
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
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.
- Low hanging fruit - Accelerated random distribution functions,
David Bateman <=