# (approximate) PDFs/CDFs suitable for automatic differentiation?

**URL:** <https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027>\
**Category:** Statistics\
**Created:** [April 3, 2017, 5:10pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027 "2017-04-03T17:10:09Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![evanfields](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evanfields/32/1744_2.png) [@evanfields](https://discourse.julialang.org/u/evanfields)\
**Post date:** [April 3, 2017, 5:10pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027/1 "2017-04-03T17:10:09Z")

</div>

I’m working on a project where, given some observations of transformed random variables, we would like to estimate the parameters of the distribution those observations came from.

As a simple example, suppose `Y = min(X_1, X_2)` where `X_i ~ Gamma(a,b)`. We see `Y`, and we’d like to estimate `a` and `b` via maximum likelihood estimation.

We can write down the [log]pdf of `Y`, but it involves the CDF of a gamma random variable. The gamma CDF available in the JuliaStats ecosystem is a Rmath call, and thus won’t accept DualNumbers as input. Therefore ForwardDiff cannot compute a derivative, and I expect the same would be true for any other differentiation package since the CDF is not computed in Julia.

Is there a pure Julia implementation of either CDFs of common distributions or the special functions needed which would support automatic differentiation? Or am I better off just fitting the maximum likelihood parameters via a derivative-free method?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 3, 2017, 5:37pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027/2 "2017-04-03T17:37:57Z")

</div>

Given that the derivative of the CDF is the PDF which is also available, you can implement a method that works on dual numbers.

---

<div class="post-metadata">

**Author:** ![evanfields](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evanfields/32/1744_2.png) [@evanfields](https://discourse.julialang.org/u/evanfields)\
**Post date:** [April 3, 2017, 5:52pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027/3 "2017-04-03T17:52:55Z")

</div>

That’s a great point for the example given. But unless I’m mistaken (quite possible), sometimes you do need the CDF. Modify the example so that `Y = min(X_1, X_2, X_3)` where the `X_i` are again iid. Now the pdf of `Y` has a `F^2` (`F` the CDF of `X`) term in it, and the derivative of the likelihood still has a `F` term, right?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 3, 2017, 6:00pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027/4 "2017-04-03T18:00:41Z")

</div>

But if I understand correctly, you have both the CDF and the PDF, so you can still define a method for it. Also, if you just define one for `F(::ForwardDiff.Dual)`, it will take care of `F(...)^2` etc. Look at the manual of ForwardDiff.

---

<div class="post-metadata">

**Author:** ![evanfields](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evanfields/32/1744_2.png) [@evanfields](https://discourse.julialang.org/u/evanfields)\
**Post date:** [April 3, 2017, 6:26pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027/5 "2017-04-03T18:26:10Z")

</div>

On further reflection I must not understand what you mean here. Denote by `f(a,x)` the PDF with parameters `a` at point `x`, and likewise for `F` the CDF. Then `f(a,x)` is the derivative of `F(a,x)` with respect to `x`. But for this MLE problem, we need the derivative with respect to `a`, the parameters.

Am I missing something obvious?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 3, 2017, 6:36pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027/6 "2017-04-03T18:36:32Z")

</div>

Sure, but if I understand your problem correctly, the only tricky part of the CDF is the (incomplete) gamma function, which has recursively defined derivatives, so a method can be defined for `ForwardDiff.Dual` using existing building blocks.

---

<div class="post-metadata">

**Author:** ![piever](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/piever/32/1815_2.png) [@piever](https://discourse.julialang.org/u/piever)\
**Post date:** [April 3, 2017, 9:56pm UTC](https://discourse.julialang.org/t/approximate-pdfs-cdfs-suitable-for-automatic-differentiation/3027/7 "2017-04-03T21:56:42Z")

</div>

I needed to differentiate a definite integral of a pdf of a gamma for my (in progress) package survival and I think I have a reasonable solution. First one needs to realize that:  
`gradient integral pdf(Gamma(a,b),x)dx = integral pdf(Gamma(a,b),x) gradlog(pdf(Gamma(a,b),x)dx`

Then I’ve implemented an efficient way to compute `integral f(x)pdf(dist,x)dx` (for one dimensional functions, so you should maybe try one component of the graldog at a time). First you get a vector of coefficients you’ll need for the actual function and then you compute it. Let’s say `dist` is you distribution object and `f` your function (gradlog in your case):

Pkg.clone([https://github.com/piever/Survival.jl.git](https://github.com/piever/Survival.jl.git))

```julia
using Survival
const coefs = Survival.clenshaw_coefs(dist, f)
# value of indefinite integral at t is
indefinite_integral_at_t = Survival.clenshaw_asin(cdf(dist,t),coefs)

```

Once you have computed the coefficients, evaluating the indefinite\_integral should be pretty fast (ideally not much slower than computing the cdf). There is a bit of a speed accuracy trade-off, handled by an optional argument. You cal also call

```julia
const coefs = Survival.clenshaw_coefs(dist, f, Val{N}())

```

where larger N gives more accuracy but takes longer. Default is 50. The code is [here](https://github.com/piever/Survival.jl/blob/master/src/fast_integral.jl). Sorry for the scarce documentation, hope it helps!