Re: Low hanging fruit - Accelerated random distribution functions

[David Bateman]

Daniel J Sebald wrote:
>> But its an exp of a scalar plus a scalar multiply, which is
>> insignificant.. Which is why I couldn't get an advantage
> 
> Principle though.  Exp is a function call (probably recursive
> approximation) while +/-/* are simple one instruction machine cycles.

Ok, given :-)


>> The restriction are reproduced from what was in the original code for
>> geometric_rnd.
> 
> Oh, well geometric and log normal are both positive support (one's
> discrete the other continuous).  Quite different otherwise.  That
> explains things.

All of the restrictions in the code come from Kurt Hornik's original
code. I meant geometric_rnd for geornd and lognormal_rnd for lognrnd, etc..

>> Looking at what the matlab code does it restricts sigma >
>> 1 and has no restrictions on mu.. So yes I think you are right and we
>> should remove the restriction.
> 
> Sigma > 1?  There's no reason sigma can't be in (0,1].

Sorry typo, I meant Sigma in [0,Inf], and is what I used in the code...


> Rather than tic/toc, use cputime().  tic/toc is too susceptible to other
> system processes.  Also, be sure to run the routine once before taking
> the measurement.  (First use is slower.)

I wanted a faster gross answer as I was looking for factors of about 100
originally, and couldn't be bothered writing a full benchmarking script.
Don't worry, when I'm looking for improvements of a few 10's of % I also
use cputime and multiple runs..

> Let's look at indexing.  Changing the mu/sigma to matrices... oops,
> getting a bug here with the CVS version:
> 
> t=cputime(); y=lognrnd(mu,sigma,400,400); cputime()-t
> error: product: nonconformant arguments (op1 is 160000x1, op2 is 1x160000)
> error: evaluating binary operator `.*' near line 97, column 45
> 
> The documentation says "Both MU and SIGMA must be scalar or of size R by
> C." so I assume what I wrote above is OK.

Yes I fixed a bug here in passing. In the original code randn (1, length
(k)) should have read randn(sz)... I can also get my code to run with

octave:1> mu = ones(400,400);
octave:2> sigma = ones(400,400);
octave:3> t=cputime(); y=lognrnd(mu,sigma,400,400); cputime()-t
ans =  0.094986

> In summary, indexing seems to have most effect (55 percent?), the mu
> tests had rough 5-6 percent, the math structuring 8 to 10 percent.  Of
> course, the most common use by the user will be scalar mu and sigma.

If I understood correctly, then you think the version I sent is ok? If
not what would you change?

D.