Re: Low hanging fruit - Accelerated random distribution functions

[Daniel J Sebald]

David Bateman wrote:
> Paul Kienzle wrote:
>
> > > 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)));
>
> This will only be a win in the case where mu or sigma are a matrix.

It's still a win even for scalar mu and sigma, isn't it?  A pair of exp and 
multiply loses to an add and exp.  The extra exp is significant, probably the 
overhead for the octave interpretation masks that.  Anyway...

Why the restriction of mu > 0?  The underlying exponentially distributed r.v. 
could have a negative mean with no consequence.  Naturally, sigma needs to be greater 
than 0 from a definition standpoint.  (And I wonder about the > 0 part because if 
variance is zero, that simply means the variable is a constant, exp(mu)).  And then 
even if the user specifies a sigma less than zero, it is effectively still a log 
normal random variable.  I.e., flip the sign of

sigma(k) .* randn (1, length(k))

and it is an r.v. with the same distribution.  (OK, may still want to take out the 
sigma < 0 case.)

Also, I wonder if the special conditioning is necessary.  If mu = Inf, exp(Inf) 
= Inf.  If mu = -Inf, exp(-Inf) = 0.  How about avoiding the use of the index?  
In most cases I'd think that the user will have no pathological cases.  I.e.,

    rnd = exp(mu + sigma .* randn (1, length(sigma))); 
    k = find ((sigma < 0) | (sigma == Inf));
    if (any (k))
        rnd(k) = NaN
    endif

instead of

    k = find ((mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf));
    if (any (k))
      rnd(k) = exp(mu(k) + sigma(k) .* randn (1, length(k))); 
    endif

Does that help any for when k is the full array?

Dan