# Differentiation of incomplete Beta function

**URL:** <https://discourse.julialang.org/t/differentiation-of-incomplete-beta-function/49675>\
**Category:** Numerics\
**Tags:** question, differentiation, forwarddiff\
**Created:** [November 6, 2020, 10:42am UTC](https://discourse.julialang.org/t/differentiation-of-incomplete-beta-function/49675 "2020-11-06T10:42:41Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![arzwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arzwa/32/12886_2.png) [@arzwa](https://discourse.julialang.org/u/arzwa)\
**Post date:** [November 6, 2020, 10:42am UTC](https://discourse.julialang.org/t/differentiation-of-incomplete-beta-function/49675/1 "2020-11-06T10:42:41Z")

</div>

First, to give some context: my use case is equivalent to the following contrived example, in which I would like to define a (`Turing.jl`) probabilistic program that uses a discretized Beta mixture model.

```julia
@model themodel(X, K) = begin
    a ~ Exponential()
    b ~ Exponential()
    p = discretize(Beta(a, b), K)
    X ~ MixtureModel([Bernouilli(p[i]) for i=1:K])
end

function discretize(d, K)
    qstart = 1.0/2K
    qend = 1. - 1.0/2K
    xs = quantile.(d, qstart:(1/K):qend)
    xs *= mean(d)*K/sum(xs) # rescale by factor mean(d)/mean(xs)
end

```

Now, this does not work, because currently the quantile function for the Beta distribution (which is the inverse of the incomplete Beta function, or `beta_inc_inv` in `SpecialFunctions.jl`) does not work with AD. (This brings me to a first question, why do many methods in `SpecialFunctions` only take `Float64` type args?). I experimented (lazily) with adapting the arg types for the relevant `SpecialFunctions` methods and using AD to get my gradient of interest, but this gets too slow.

Then I’ve implemented an algorithm for the gradient of the CDF of the beta distribution (i.e. `beta_inc`) with respect to the `a` and `b` parameters (which I found [here](https://www.jstatsoft.org/article/view/v003i01)). This gives me the gradient of `beta_inc(a, b, x)` w.r.t. `a` and `b` fast. I would expect that this would enable faster computations of the `beta_inc_inv` gradient (which mainly calls `beta_inc` in a numerical inversion routine if I understand correctly).

However, now I’m not sure how to proceed (I am no analysis guru, and all the stuff in `JuliaDiff` is a bit daunting). It seems to me I should implement some diffrule for `beta_inc`? Say my function is `beta_inc_grad(a,b,x)` and returns a tuple with `beta_inc(a,b,x)` and its partial derivatives w.r.t. `a` and `b`, is

```julia
@define_diffrule SpecialFunctions.beta_inc(a, b, x) =
    :( beta_inc_grad($a, $b, $x)[2] ), :( beta_inc_grad($a, $b, $x)[3] ), :NaN

```

what I should implement to make this work with AD (it doesn’t seem to be enough to make `beta_inc` play with ForwardDiff)? (Also, would this call `beta_inc_grad` twice?)

---

<div class="post-metadata">

**Author:** ![bdeonovic](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bdeonovic/32/3928_2.png) [@bdeonovic](https://discourse.julialang.org/u/bdeonovic)\
**Post date:** [June 9, 2021, 4:40pm UTC](https://discourse.julialang.org/t/differentiation-of-incomplete-beta-function/49675/2 "2021-06-09T16:40:59Z")

</div>

Would you be willing to share your `beta_inc_grad` code?

---

<div class="post-metadata">

**Author:** ![arzwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arzwa/32/12886_2.png) [@arzwa](https://discourse.julialang.org/u/arzwa)\
**Post date:** [June 9, 2021, 7:22pm UTC](https://discourse.julialang.org/t/differentiation-of-incomplete-beta-function/49675/3 "2021-06-09T19:22:06Z")

</div>

Hi Benjamin, I haven’t followed up on this, but I put what I had implemented back when I asked this question here [https://github.com/arzwa/IncBetaDer](https://github.com/arzwa/IncBetaDer). The `beta_inc_grad` code is [here](https://github.com/arzwa/IncBetaDer/blob/main/src/beta_inc_grad.jl), I hope it can be useful. I’m not at all a numerical algorithms guy, so it might not be very good or correct currently, although I do think I had some tests. (Needless to say, I haven’t prepared this code for release or so, so comments etc. should be read as notes to myself).

---

<div class="post-metadata">

**Author:** ![bdeonovic](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bdeonovic/32/3928_2.png) [@bdeonovic](https://discourse.julialang.org/u/bdeonovic)\
**Post date:** [June 9, 2021, 7:31pm UTC](https://discourse.julialang.org/t/differentiation-of-incomplete-beta-function/49675/4 "2021-06-09T19:31:18Z")

</div>

I’m not a numerical algo guy either, but I also need these derivatives (in my case because I want to use TDistribution for some Bayesian stuff in Turing, but it doesn’t work with AD because of the `beta_inc`).

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [October 5, 2025, 12:55pm UTC](https://discourse.julialang.org/t/differentiation-of-incomplete-beta-function/49675/5 "2025-10-05T12:55:22Z")

</div>

@arzwa I have ported your code for integration into SpecialFunctions.jl in [this PR](https://github.com/JuliaMath/SpecialFunctions.jl/pull/506) but we require your aggrement for licensing it under MIT of course.