# Gamma quantile function

**URL:** <https://discourse.mc-stan.org/t/gamma-quantile-function/152>\
**Category:** Developers\
**Tags:** features\
**Created:** [December 2, 2016, 5:46pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152 "2016-12-02T17:46:32Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 2, 2016, 5:46pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/1 "2016-12-02T17:46:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 2, 2016, 6:39pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/2 "2016-12-02T18:39:28Z")

</div>

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](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

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 2, 2016, 6:53pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/3 "2016-12-02T18:53:53Z")

</div>

> [@Bob\_Carpenter](#):
>
> So maybe if we come up with a general differentiable optimizer then some of these functions will be easier to implement.

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).

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 2, 2016, 7:21pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/4 "2016-12-02T19:21:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [December 2, 2016, 7:30pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/5 "2016-12-02T19:30:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [December 2, 2016, 8:02pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/6 "2016-12-02T20:02:52Z")

</div>

> [@sakrejda](#):
>
> AFAICT it’s all going to work with autodiff

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.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [December 2, 2016, 8:13pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/7 "2016-12-02T20:13:15Z")

</div>

> [@sakrejda](#):
>
> Would a one-off for this function be acceptable?

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.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 2, 2016, 8:20pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/8 "2016-12-02T20:20:02Z")

</div>

> [@bgoodri](#):
>
> boost/math/special\_functions/gamma.hpp

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

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 2, 2016, 9:07pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/9 "2016-12-02T21:07:03Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [December 3, 2016, 12:44am UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/10 "2016-12-03T00:44:45Z")

</div>

> [@sakrejda](#):
>
> 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.

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](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.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 3, 2016, 12:56am UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/11 "2016-12-03T00:56:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [December 3, 2016, 6:35am UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/12 "2016-12-03T06:35:16Z")

</div>

Or you could try incorporating

> **[thomasluu/quantus](https://github.com/thomasluu/quantus)**
>
> Efficient and accurate parallel inversion of the gamma distribution - thomasluu/quantus

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 3, 2016, 6:36pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/13 "2016-12-03T18:36:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [December 3, 2016, 7:49pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/14 "2016-12-03T19:49:27Z")

</div>

Boost is certainly workable, but its algorithm works slowly.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 3, 2016, 10:22pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/15 "2016-12-03T22:22:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 4, 2016, 8:41pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/16 "2016-12-04T20:41:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 4, 2016, 9:10pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/17 "2016-12-04T21:10:28Z")

</div>

> Bob\_Carpenter [http://discourse.mc-stan.org/users/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

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 4, 2016, 9:22pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/18 "2016-12-04T21:22:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 5, 2016, 4:09pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/19 "2016-12-05T16:09:19Z")

</div>

> [@bgoodri](#):
>
> But I don’t know what to do about the shape parameter.

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

---

<div class="post-metadata">

**Author:** ![syclik](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/syclik/32/6_2.png) [@syclik](https://discourse.mc-stan.org/u/syclik)\
**Post date:** [December 5, 2016, 8:16pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/20 "2016-12-05T20:16:03Z")

</div>

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.

[Next page](https://discourse.mc-stan.org/t/gamma-quantile-function/152.md?page=2)