Replace code depending on Rmath.jl

[dmbates]

There has been discussion on discourse and on slack regarding replacing the last dependencies of StatsFuns on Rmath.jl. In the discourse thread the VBA code by Ian J. Smith was mentioned. See https://github.com/dmbates/IanSmithdistmath.jl

@andreasnoack commented in the slack discussion that

The gamma and the beta cover a lot of cases. One thing missing for those two is the inverse log cdf. I might be possible to adjust the non-log versions from SpecialFunctions. The most important completely unimplemented distribution is probably the hypergeometric. Maybe it can be composed from things in SpecialFunctions and HypergeometricFunctions. I’d have to look closer.
The all the non-central distributions are not yet covered and unfortunately I’m not sure that they’ll be anytime soon since 1. few people care about them and 2. they will probably be demanding to implement. Maybe we should split StatsFuns and Distributions into packages with the current names and then StatsFunsExtra and DistributionsExtra where we keep the exotic stuff that would then still depend on Rmath while StatsFuns and Distributions could be made GPL free.

It appears that some of the noncentral distributions are covered in the VBA code.

[dmbates]

I have added some automatic translation to https://github.com/dmbates/IanSmithdistmath.jl but am reaching the point where I need advice from someone who knows VBA. @ararslan ?
I am assuming that the Dim declarations, at least for scalars, can safely be omitted. Is that correct?

Also, I am not sure what a return value of [#VALUE!] means.

I assumed that the Sqr(x) function is abs2(x) but it seems to occur more than sqrt(x) and I see other places where statements like x2 = x * x occur. Should I have translated it as sqrt(x)?

[palday]

#VALUE! is the sentinel return value for a domain or similar error in the inputs in Excel.

[ararslan]

Should I have translated it as sqrt(x)?

Yes, Sqr is indeed sqrt.

[ararslan]

I am assuming that the Dim declarations, at least for scalars, can safely be omitted. Is that correct?

Yep. Dim X will define X as a scalar (but otherwise basically Any) without a value, Dim X(i_1, ..., i_n) defines X as an n-dimensional array. Dim X As T or Dim X(...) As T will define X as a value with type or array with eltype T, i.e. X::T or X::Array{T,N}.

[dmbates]

I will start work on the hypergeometric distribution functions in a branch of StatsFuns.jl using translations of the VBA code from Ian Smith.

[dmbates]

I will start work on the hypergeometric distribution functions in a branch of StatsFuns.jl using translations of the VBA code from Ian Smith.

This ended up being a "Heart of Darkness" descent into madness where each function depends on several other functions that probably have equivalents in SpecialFunctions.jl either under gamma or beta sources but the documentation is so sparse that it is difficult to tell. There are just enough VBA idiosyncracies that I can't be sure I am doing the right thing and I haven't found tests so that I can check.

I am afraid I am going to have to abandon this approach.

[palday]

@dmbates is there a reason you started with hypergeometric? In other words, if I give parts of this a shot, is that a good place for me to start as well?

[dmbates]

@dmbates is there a reason you started with hypergeometric? In other words, if I give parts of this a shot, is that a good place for me to start as well?

The reason I started with the hypergeometric is because it is one of the few distributions that needs a full suite of functions to replace the Rmath code. But, as I said, it led into a descent into madness trying to decipher all the Stirling formula approximations for the loggamma, etc.

According to the R documentation, the hypergeometric pdf in Rmath now uses Loader's saddle-point approximation, which is interesting in that it is related to deviance calculations for binomial glms. Her paper (https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf) contains C code to implement this method for the binomial (log)pdf. However, there is no copyright statement so I am not sure whether MIT-licensed code can be derived from it. The R documentation also says that the cdf calculations are based on methods from Ian Smith so this code may be relevant. I just wish it was possible to step back from the implementation details to get an overview of the methodology.

[dmbates]

I just noticed SpecialFunctions.stirling_error which is one piece of Loader's saddlepoint approximation that is already implemented.

[dmbates]

As I have said, I got bogged down trying to translate Ian Smith's VBA code. It is a single 12,600 line source code file. I find it hard to pull apart the pieces associated with a particular function that we would use.

As I mentioned, the R-project provides a PDF of Catherine Loader's paper that gives a description of her algorithm and includes both equations and code. She describes the implementation of (log)pdf and (log)cdf for binomial and Poisson.

I am beginning to think that the most effective way of writing unencumbered Julia implementations of functions related to the hypergeometric distribution would be to have one person, possibly me, read the code and write a description of the algorithm, then have someone else code from the description. I don't think it would be appropriate for me to be the one writing the code because of my close ties to the R Project.

Does this seem like a workable approach or are there other avenues we haven't tried yet?

[palday]

Clean room implementations are a lot of work, but I would be willing to take on the Julia side of things.