Gamma quantile function

Developers
1

Am I right that we don’t have the quantile function for the gamma density? I was looking at putting priors on the median rather than the scale parameter for generalized gamma models to make them more comparable with gamma models for survival time so I’m looking into adding it.

2

We don’t have any of the inverse CDFs. They would be
useful if we could code them with derivatives. But
they don’t have simple closed forms.

Looks like MATLAB uses Newton’s method to calculate them:

https://fr.mathworks.com/help/stats/gaminv.html

So maybe if we come up with a general differentiable optimizer
then some of these functions will be easier to implement.

  • Bob
3

Thanks for confirming, sometimes what’s in stan::math and what’s not is still hard for me to figure out. It looks like R’s qgamma uses a chi-squared approximation (that does have a closed form) and does Newton from there. AFAICT it’s all going to work with autodiff. Would a one-off for this function be acceptable? (yes with tests and everything).

4

Absolutely. But if it involves differentiating a Newton solver,
it’d be great if that piece were modular.

We just have to come up with the naming convention for inverse
cdfs. I suppose these aren’t going to be on the log-odds or log
scale, so it’ll just be something like

gamma_inv_cdf().

My own preference is for inv_cdf over _quantile. And here I
don’t think we want any of the vertical bar notation.

  • Bob
5

I like foo_icdf. But we need to decide what we want to do with inverse CDFs generally. There are many distributions that don’t have them in closed form. And there are some where the inverse CDF exists but the PDF does not. And when the PDF does exist in Stan, there are some distributions where you would want a different parameterization for the inverse CDF.

6

For known shape, you don’t need to autodiff. The derivative of the scale parameter would be the quantile with scale = 1 and the derivative with respect to the quantile is the reciprocal of the PDF (including constants). But I don’t know what to do about the shape parameter.

7

Also, Boost has numerical approximations to the inverse of the incomplete Gamma function in

boost/math/special_functions/gamma.hpp

so you wouldn’t have to do your own Newton just for this.

8

Thanks! I looked at GSL and R but didn’t think of Boost!

9

Looks like the Boost versions invert on ‘a’ so no dice, we need the inverse of the incomplete gamma w.r.t. z to get the quantile function.

10

Boost has all 4 kinds. The one you are looking for is gamma_p_inv:

http://www.boost.org/doc/libs/1_62_0/libs/math/doc/html/math_toolkit/dist_ref/dists/gamma_dist.html#math_toolkit.dist_ref.dists.gamma_dist.implementation

although remember to divide instead of multiply if you are using the shape-rate parameterization.

11

Found it this time, it was in a different header, thanks.

12

Or you could try incorporating

13

I’d much prefer Boost to another dependency if it’s workable.
If not, that linked quantus package is at least MIT licensed.
It’s a very focused package—just the inverse gamma distribution!

  • Bob
14

Boost is certainly workable, but its algorithm works slowly.

15

The quantus code is really small—I think we (I) could just rewrite the
algorithm for the stan::math library.

16

We can’t just “rewrite” someone else’s code. It’s copyrighted.
If there’s a published algorithm somewhere you could work from
to come up with a different implementation, that’d be different.

If you really think it’s that much better we can include it, but it
makes all of our notes and everything more complicated. Our current
dependencies are minimal:

Boost
Eigen
Sundials

plus googltest for testing.

Before adding it, I’d like to see a speed and accuracy comparison
vs. Boost.

  • Bob
17

Bob_Carpenter http://discourse.mc-stan.org/users/bob_carpenter Developer
December 4

We can’t just “rewrite” someone else’s code. It’s copyrighted.

Sure, it’s MIT license so you can merge/modify the code if you keep the
license header with the header. This is not something like boost where we
have to worry about tracking version so I don’t think it’s much of a
dependency.

If there’s a published algorithm somewhere you could work from

18

Gotta say adding the boost algorithm looks pretty trivial in comparison so
it’s worth seeing if it’s adequate.

19

That’s what I’m looking for now. Maybe somebody’s done something usable but I haven’t run into it yet.

20

My suggestion: do the easy thing first.

If there’s an implementation in a library that we already have a dependency on, please start there. Don’t worry about speed until you have something working. (I can always give you something that’s really fast and wrong.) Get the easiest thing working, write enough tests to be sure it works, then write a pull request for that. Then go and try to implement a fast version using the existing tests to help make sure it actually works. And who knows… maybe the slow version will be fast enough for whatever use case you have.