Re: betarnd function

[Paul Kienzle]

On Feb 22, 2007, at 5:32 PM, antonio palestrini wrote:

> I use the standard "inversion method" of a distrib function. That is,
>
> function p = pareto_rnd(k,a,n)
>    p = k*rand(n,1).^(-1/a);
> endfunction
>
> where "k" is the location parameter and "a" the shape parameter.
> The asymptotic properties of such distrib. requires n large. Even
> though the above function is fast since is vectorized,  do you know of
> a faster octave function?

Rewriting slightly,

   rnd = k * exp ( - log (rand(n)) / a )

The function rande returns values from the standard
exponential, distribution, which is -log(uniform),
so the following should work for pareto and be a bit
faster:

   rnd = k * exp ( rande(n) / a)

octave> tic; 5 * exp ( - log ( rand(400) ) / 20 ); toc
ans = 0.32883
octave> tic; 5 * exp ( rande(400) / 20 ); toc
ans = 0.19606

The rande() call itself is 1/4 of this.

You could probably develop your own ziggurat code for
particular pareto parameters and see similar speed
gains, but I doubt this is worth the effort --- just
summing an array this big takes as much time as rande(),
and surely your simulation code is an order of magnitude
more complex than that.

- Paul