Fixing types for RNGs

[ChrisRackauckas]

This is related to JuliaStats/Rmath.jl#9. There are a lot of RNG functions which are returning inappropriate types, i.e. floats instead of ints. For example, rpois, rnbinom, rnbinom_mu, rgeom should all be returning ints.

[simonbyrne]

I don't think we want to change Rmath code ourselves.

[ChrisRackauckas]

This leads to obvious type problems whenever using the StatsFuns RNGs like rand(Poisson(lambda)) since one would always expect that to give an Int, and so a conversion has to occur. If Rmath is going to be continued to be used with Julia, it would make sense to fix these instabilities. However, if the plan is to drop this dependency soon, then this could be ignored.

[ChrisRackauckas]

Maybe @tkelman has a comment on whether RMath is planned to stick around with StatsFuns, or whether these RNGs are all planned to be replaced soon to get rid of the dependency on this library?

[tkelman]

That's entirely a function of the amount of work that goes into replacing it. If there's a major performance regression between what worked on 0.4 and how things are currently set up on 0.5 then we can adjust how we build and package things, but I'm the wrong person to ask any other questions here.

[ChrisRackauckas]

Got it. Wasn't sure if you were in charge of deprecating Rmath or just moving it out of Base (or what entirely is happening there). That's an issue for StatsFuns.

Is there a reason why "I don't think we want to change Rmath code ourselves"? The RNGs were already modified to use a different underlying generator of random uniform numbers. Changing the types to be the right kind would make sense as well because "everything is a float" isn't true in Julia, and in this case we have a lot of integers represented as floats.

[simonbyrne]

There are no instabilities. Rmath (the C library) returns a Cdouble (ie a Float64) and Distributions.jl converts it to an Int. I haven't timed it, but given how complicated the algorithm is, I don't think the overhead will be particularly large.

I don't think modifying Rmath code is time well spent, since it's never going to be upstreamed (due to R's inherent design), and then we have the pain of maintaining our own fork.

[simonbyrne]

The RNG changes are in a separate file (and Rmath was designed to allow for custom RNGs)

[ChrisRackauckas]

Got it. I was just tracking Poisson generation performance at this example and was wondering if more performance can be gained by this change. But it sounds like it would be a hassle to keep this change when Rmath will keep getting updated.

A related question is whether Julia's standard RNGs for these types of distributions should be using Rmath at all, and whether it should be able to switch between the different RNGs from RNG.jl (the GSoC project) as needed. But that's for StatsFuns (I'll open an issue there).