Random Number Generation Performance Issues

[ChrisRackauckas]

There are some issues with how the random number generating procedures are wrapped. To explain, look at how rpois is wrapped:

    ## Distributions with 1 parameter and no default
    macro libRmath_1par_0d(base)
        dd = Symbol("d", base)
        pp = Symbol("p", base)
        qq = Symbol("q", base)
        rr = Symbol("r", base)
        quote
            global $dd, $pp, $qq, $rr
            $dd(x::Number, p1::Number, give_log::Bool) = 
                ccall(($(string(dd)),libRmath), Float64,
                      (Float64,Float64,Int32), x, p1, give_log)
            $pp(q::Number, p1::Number, lower_tail::Bool, log_p::Bool) =
                ccall(($(string(pp)),libRmath), Float64,
                      (Float64,Float64,Int32,Int32), q, p1, lower_tail, log_p)
            $qq(p::Number, p1::Number, lower_tail::Bool, log_p::Bool) =
                ccall(($(string(qq)),libRmath), Float64,
                      (Float64,Float64,Int32,Int32), p, p1, lower_tail, log_p)
            $rr(nn::Integer, p1::Number) =
                [ccall(($(string(rr)),libRmath), Float64, (Float64,), p1) for i=1:nn]
            @libRmath_1par_0d_aliases $base
        end
    end
@libRmath_1par_0d pois

For reference, using rpois is like this: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Poisson.html

Two things to point out. One is that using a comprehensions gets rid of a lot of potential speedups. For example, it can re-use values if it makes more than 1 random number at a time, and so it should actually pass nn in and do something more as mentioned here. Also, Poisson random numbers should return ints, not Floats. Lastly, you should be able to pass it a p1 which is a vector. Also, to work like other random number algorithms in Julia, it should return a number and not a vector when nn is 1.

I suspect the other random number generating algorithms have similar problems.

[tkelman]

sorry my mistake, nevermind that comment. but note that the julia code here is very old and almost entirely replaced elsewhere, e.g. StatsFuns.jl

[simonbyrne]

Perhaps we should get rid of most of the Julia code, if it's not being used?

[tkelman]

that may make sense

[ChrisRackauckas]

That would be helpful, since I first tracked down rand(Poisson(lambda)) to here accidentally, instead of to StatsFuns.RFunctions.poisrand(lambda) defined here.