Re: Low hanging fruit - Accelerated random distribution functions

[Paul Kienzle]

On Feb 22, 2007, at 11:46 AM, David Bateman wrote:

> I haven't tried to convert the functions geornd,

Octave uses the following; it would be hard to beat:

   rnd(k) = floor (log (rand (size (k))) ./ log (1 - p(k)));


> hygernd,

This needs a fast discrete_rnd, which could be implemented
using lookup as [untested!]

rnd = discrete_rnd(p,v,n)
        edges = cumsum(p(1:end-1))/sum(p);
   rnd = v(lookup(edges,rand(n,1))+1);

The current code is approximately:

  edges = cumsum(p)/sum(p);
  u = rand(1,n);
  rnd = v (1 + sum( edges*ones(1,n) <= ones(m,1)*u));

Given n numbers generated from m possibilities:
 current code takes O(m*n) space and O(m*n) time.
 lookup takes O(m+n) space and O((m+n)*log(m+n)) time.
so lookup saves both space and time.

Note that the documentation says p is a probability
whereas it is actually a proportion.

> laplace_rnd,

The heart of this will be fast enough:

  tmp = rand (sz);
  rnd = (tmp<1/2) .* log(2*tmp) - (tmp>1/2) .* log (2*(1-tmp));

A C implementation would run in 1/2 the time with
1/4? the memory but I don't think we care.

> logistic_rnd,

Fast enough:

  rnd = - log (1 ./ rand (sz) - 1);


> lognrnd,

Current code can be improved a bit from:

 rnd(k) = exp(mu(k)) .* exp(sigma(k).*randn(1,length(k)));

to:

 rnd(k) = exp(mu(k) + sigma(k).*randn(1,length(k)));


> pascal_rnd,

Needs work.

> stdnormal_rnd,

Fast enough:

  rnd = randn(sz)

> wblrnd

Fast enough:

  rnd = scale * -log(1-rand(k))) .^ (1/shape)

>  and wienrnd

Fast enough:

  rnd = cumsum(randn(n*t,d))/sqrt(n)


- Paul