# 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:** 10\
**Page:** 2

<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, 8:25pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/21 "2016-12-05T20:25:46Z")

</div>

I wasn’t too concerned with the speed of the inverse calculation but with auto-diff to get the gradients since the inverse gamma calculation is an approximation to begin with. I don’t think there’s any way to get around that (no published algorithm that I can find that doesn’t invoke something awkward like Meijer’s G function) so I’m going to take your advice. I can test code I write against values from R’s qgamma and gradients against finite diff? If there are other suggestions on testing I’d welcome them.

---

<div class="post-metadata">

**Author:** ![thel](https://avatars.discourse-cdn.com/v4/letter/t/97f17d/32.png) [@thel](https://discourse.mc-stan.org/u/thel)\
**Post date:** [December 8, 2016, 3:43am UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/22 "2016-12-08T03:43:09Z")

</div>

Regarding the derivatives, for any distribution we have icdf’(y) = 1 / pdf(icdf(y)) – this is just the chain rule. So there’s no extra work to get the derivatives, assuming we have a pdf and a good approximation of icdf. If we were willing to do newton searches we could do a generic implementation inverse cdfs with derivatives for all the univariate distributions at once, by implementing inverses for univariate functions. If you need the derivative with respect to the shape parameters, there’s a similar generic computation using the multivariate chain rule. But I do expect a lot of overflows out in the tails.

---

<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 8, 2016, 6:39am UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/23 "2016-12-08T06:39:34Z")

</div>

> thel Developer  
> December 8  
> Regarding the derivatives, for any distribution we have icdf’(y) = 1 / pdf(icdf(y)) – this is just the chain rule. So there’s no extra work to get the derivatives, assuming we have a pdf and a good approximation of icdf. If we were willing to do newton searches we could do a generic implementation inverse cdfs with derivatives for all the univariate distributions at once, by implementing inverses for univariate functions.

I think this’d be worth exploring. It’d be nice to have  
the inverse CDFs.

It’d be nice to expose the actual optimizer and its derivatives,  
I think.

> If you need the derivative with respect to the shape parameters,

Yes, we need derivatives with respect to everything. At least  
with the way functions are architected and documented now.  
They can be used anywhere the Stan types match.

- Bob

---

<div class="post-metadata">

**Author:** ![thel](https://avatars.discourse-cdn.com/v4/letter/t/97f17d/32.png) [@thel](https://discourse.mc-stan.org/u/thel)\
**Post date:** [December 8, 2016, 7:26am UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/24 "2016-12-08T07:26:33Z")

</div>

In this approach the story with the derivative with respect to a parameter q is to just to write

cdf(q, icdf(q, x)) = x

differentiating with respect to q and writing partial\_q for the derivative we find

partial\_q cdf(q, icdf(q,x)) + pdf(q, icdf(q,x)) \* partial\_q icdf(q,x) = 0

so

partial\_q icdf(q,x) = - partial\_q cdf(q, icdf(q,x)) / pdf(q, icdf(q,x))

again I worry about overflows but I think they’re “true” overflows in the sense that icdf is genuinely badly behaved in the tails (swap horizontal for vertical). But the pattern is the same as before: if I can evaluate the function, its inverse, and the derivatives of the inverse, I can evaluate all derivatives of the function with generic formulae.

The downside of this is of course that these formulas depend on the inverse function relationship being exact, so the “derivatives” would not actually be derivatives of the approximation, nor are they obviously the best approximation to the derivative. I’m starting to doubt it. Looking at those pdfs in the denominator I’m quite worried.

Anyway, a generic enough approach but I’m suspicious of the accuracy. Might be worth trying if we’ve got time to do a lot of testing.

-thel

---

<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 8, 2016, 3:38pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/25 "2016-12-08T15:38:38Z")

</div>

This all makes sense, I like the idea of getting a newton root solver in and having an implementation for all the inverse cdf’s. I would be interested in working on that sometime in January but won’t have any time for it before that really. In the next week or so I will get Boost’s inverse gamma and derivatives done at least on a branch and I’ll use this idea for the shape derivative so maybe that can be a test case.

---

<div class="post-metadata">

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

</div>

Maybe there are some ideas to take from this blog post:  
[http://www.shapeoperator.com/2014/02/22/bisecting-floats/](http://www.shapeoperator.com/2014/02/22/bisecting-floats/)  
where some floating point traps in the context of bisection are described

Julia implementation of some methods are there:

> **[JuliaMath/Roots.jl](https://github.com/JuliaMath/Roots.jl)**
>
> Root finding functions for Julia. Contribute to JuliaMath/Roots.jl development by creating an account on GitHub.

This derivative free method looks especially interesting:

> **[New Eighth-Order Derivative-Free Methods for Solving Nonlinear Equations](https://www.hindawi.com/journals/ijmms/2012/493456/)**
>
> A new family of eighth-order derivative-free methods for solving nonlinear equations is presented. It is proved that these methods have the convergence order of eight. These new methods are derivative-free and only use four evaluations of the...

---

<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 13, 2016, 1:55pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/27 "2016-12-13T13:55:19Z")

</div>

> [@vincent-picaud](#):
>
> some floating point traps in the context of bisection are described

That’s fun stuff! I hadn’t run into such a clear example of the discreteness of floating point numbers yet.

---

<div class="post-metadata">

**Author:** ![vincent-picaud](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/vincent-picaud/32/24_2.png) [@vincent-picaud](https://discourse.mc-stan.org/u/vincent-picaud)\
**Post date:** [December 13, 2016, 6:04pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/28 "2016-12-13T18:04:21Z")

</div>

> That’s fun stuff! I hadn’t run into such a clear example of the discreteness of floating point numbers yet.

That is also the case for me, the example is quite illuminating.

I have not implemented [New Eighth-Order Derivative-Free Methods for Solving Nonlinear Equations](http://www.hindawi.com/journals/ijmms/2012/493456/) yet. To be honest I fear a little bit about numerical stability. However if it works in a stable way, the 8th order of CV must be impresive.

---

<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 13, 2016, 7:52pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/29 "2016-12-13T19:52:35Z")

</div>

The usual place we run into it is that 1 - epsilon evaluates  
to 1 for any epsilon with -1e-15 \<\< epsilon \<\< 1e-15 or so.  
There’s much better resolution around 0, hence the need for  
complementary CDFs.

---

<div class="post-metadata">

**Author:** ![jstag](https://avatars.discourse-cdn.com/v4/letter/j/9dc877/32.png) [@jstag](https://discourse.mc-stan.org/u/jstag)\
**Post date:** [July 1, 2020, 7:42pm UTC](https://discourse.mc-stan.org/t/gamma-quantile-function/152/30 "2020-07-01T19:42:30Z")

</div>

Were you ever able to incorporate Boost’s inverse gamma as you mention here? I have been struggling to do the same.

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