HP PPL and Statistical distributions

[compaqdrew]

I did some digging on this. The root of the issue is a typical problem of mathematical calculation. Expressions which are nice for math, are not always nice numerically. Even my Mathematica struggles on it. You need either a simpler expression or a fast approximate method. Of course the nature of either is it will be less general than the pure form, but how to handle that depends a lot on your application.

On the expression side, Wikipedia has nice closed-form solutions for special values of k, if you can be satisfied with a less general solution. If those don't work it might be possible to get some bounds for your application and ask a big boy CAS like Mathematica if it can find a nice solution. I did ask it to simplify the general formula on the offchance it can prove some new math for us, it appears it cannot.

The simplest and typical implementation is to use tables and interpolate. For example, this is the approach given here. I also saw an approach like this recommended in one of the older papers on this function. Basically, load up a matrix with the table, do lookups and average nearby values. One of the papers I read recommended harmonic interpolation between table entries but that may be overkill.

Another approach is to use a numerical method to slice open the nasty integrals, such as some kind of quadrature rule. Glancing at the formula, I'd guess you use some Laguerre method on one integral and a Hermite method on the other. Fortunately we don't have to guess, there are known algorithms for this function. Here is an older method (indeed, seems to be Laguerrish) and here is a modern one from R. Some of these methods have limitations on the parameters. There's also this survey of similar methods which has a large appendix with implementations. Unfortunately none of these are HP PPL, but it's probably straightforward to translate.