Hi,
I am currently working on a model which requires the use of an upper incomplete gamma function. For some parameter combinations the argument of this incomplete gamma can become negative. Unfortunately, the gamma function in Stan/boost only takes non-negative values as an input. A colleague uses Matlab and there the implemented incomplete gamma function takes negatives inputs(the documentation however makes no mention of how exactly that is computed).
I have come across this paper about an algorithm of computing the incomplete gamma with negative position and argument(Algorithm 969). It has been written in Fortran and uses series expansion if I understand it correctly(this slightly outside of my expertise and comfort zone).
I was wondering if it is straight forward to incorporate these Fortran function into my Stan code via its capability to link to functions in other languages or are there particular things to consider.
Otherwise, does anyone have any other suggestion of going about this problem? Amy help is appreciated :)
Do you know what part of the domain on a/z you need? With the special functions implementation isn’t too bad you just have to be thorough when checking the calculation over the domain you need. Part of the work is just calculating the reference values and gradients from another source or two so we can have confidence in an implementation. Then you could implement an algorithm you find as a Stan function and it would work as long as it’s accurate enough.
Great. We only consider positive a and negative values of z in the range of 300, I reckon. That would require the computation of Eq. 6 and Eq. 29 from the above algorithm reference, I think.
Would you rather explicitly write a Stan function implementing the algorithm or rather link “feed” their Fortran implementation into the Stan code? Particularly, thinking about accuracy and speed of the implementations.
You couldn’t use the Fortran implementation directly because it doesn’t calculate gradients, but if you re-wrote it in Stan code you’d get the gradients from our auto-diff so that’s probably the way to go. You might be surprised with most algorithms for how bad the (relative) accuracy gets when you put a +/-300 box around the incomplete gamma parameters. I don’t know that it currently performs well enough on the positive part of the range for parameters that large (that might be a first thing to test on your end).
Okay, I understand. The algorithm right-up states relative accuracy ∼10^-13 for down to -500 for z.
I will check about the positive part of parameters with the current implementation. What reference would you compare it to?
The algorithm for z < -50 looks like this
The fortran lines calculating the sum looks like below. How easily transferable are these lines(particularly things like the limits of computational accuracy(not sure if that is the correct terminology and if it even does what I think it does in the Fortran implementation) such as “giant” and “eps” in the lines below.
DO WHILE ((abs(t/v)>eps).AND.(m==0))
p=(1.0_r8-a+k);
t=t*p/x;
v=v+t
k=k+1
IF (t>giant) m=1
ENDDO
Sorry, if I am not clear or precise. I am new two these kind of algorithms, but I am definitely keen to get it working.
I’ve used Mathematica, Matlab, octave, and boost::math to get reference values. You calculate and store a grid of values and then compare your calculations to that.
eps is a relative tolerance, it just cuts off the calculation once it’s close enough. In Stan you can usually use something like 10^-8 but experiment with larger values so you can see when it’s good enough. Giant probably just detects if the summation has diverged. You can look at the incomplete gamma is calculate in Stan (or the f32 special function) and you’ll see a similar pattern applied.
I’m currently grappling with this problem myself. No breakthroughs just yet but I have a quadrature implementation I can share if you want [needs to be slightly modified]. Not very classy, I know, but should get the job done.
I’m still looking for some symmetry argument or obscure property of \Gamma(s, x) that can be exploited to avoid computing the quadrature.
In the file for the hypergeometric function (3,2 a.k.a. F32) there’s an example of what this sort of algorithm looks like in Stan (it’s a short function). In the file for the the gradient of the lower incomplete gamma there’s another example that’s more complicated mostly because it has two implementations that cover two parts of the parameter space. In the test file for the gradient of the lower incomplete gamma is how I used the reference values (calculated from Mathematica I believe) to test my implementation. @maxbiostat@jmac I’m sharing this because the actual C++ code involved here is relatively short and it’s the math and the testing that are most time consuming. For the lower incomplete gamma gradients I tried out a half dozen implementations (a few from papers + variations) to come up with the two that did a decent job of covering the parameter space (papers are not always clear or right about which part of parameter space their algorithm fails in). If you can come up with a table of reference values and an implementation in R that we can run to demonstrate satisfactory accuracy over a known parameter range then I’m happy to turn the R into C++.
I saw on the threads @maxbiostat linked that there is a range of solutions and some distributions that need the calculation so I think that’s likely enough reason for @syclik to ok including the functionality.
One risk here that might double the amount of work on your end (and mine) is that the auto-diff through these algorithms doesn’t always work all that well. I did the work on the algorithm specifically for the gradients of the lower incomplete gamma because for a distribution I derived the values were good enough accuracy but the gradients were less much accurate in parts of the parameter space and it leads the HMC acceptance rate to crash (therefore stepsize and effective sample size too).
I realised my problem is distinct from @jmac’s, in that I needed to compute the incomplete gamma function for negative first argument. I think the discussion could benefit from precise definitions so we don’t get confused about what’s being computed. Following Wikipedia, what I needed to compute was
which I find a bit weird, but OK.
I needed \Gamma(a, z) with a <0 and got it done following the implementation in one of the links. I have no idea how to use Stan-only functions to compute \gamma^\star(a, z) for z < 0, but if you want we can try and implement the quadrature for use with Stan. I know how to use it with Rstan and maybe Pystan.
@sakrejda, I got the impression you’d like to see this implemented in the math library, which is of course more work. But if the distributions you mention do not involve negative arguments, the function in (*) can be easily obtained from tgamma and gamma_q, right?
The request from @jmac was about negative function arguments, that’s the only use case that isn’t implemented*. Showing that your function does what you think it does on the math side is the bulk of the work, that’s why I shared the math lib code (it’s relatively straightfoward). If somebody doesn’t care whether their model works that leaves a lot more leeway in terms of implementation.
It’s just normalized which helps (can help) with computational accuracy. If you want to blow up the scale you’re working on you can always multiply by tgamma. It’s really common in these papers to skip the normalized/non-normalized, lower/upper, complete/incomplete language and just let the reader figure it out. Oh well.
