# Implementation of Dirichlet cdf?

I need to compute the Dirichlet CDF, but I can only find implementations of the PDF.

Do you guys know of any library (preferably in R) implementing it?

• Not directly aware of any. But there may be something that can be done. What do you need to do with it? Apr 25, 2013 at 22:58
• I need to take the complementary of the CDF and consider it as my p value. Apr 26, 2013 at 8:40
• Hmm. So yeah, if you need $1-P(X_1\leq x_1,X_2\leq x_2, ...,X_k\leq x_k)$, you kind of do need the cdf. Zen's idea of simulation is certainly a way to do it (and the higher the number of dimensions, the better it starts to look), but if you do that, use one of the packages with built-in implementations of rdirichlet. If it's only 3-variate or possibly 4-variate (the last component, of course, being redundant) it may be worth trying numerical quadrature. Apr 27, 2013 at 0:50

Remember that, if $Y_i$ are independent $\mathrm{Gamma}(a_i,b)$, for $i=1,\dots,k$, then $$(X_1,\dots,X_k) = \left(\frac{Y_1}{\sum_{j=1}^k Y_j}, \dots, \frac{Y_k}{\sum_{j=1}^k Y_j} \right) \sim \mathrm{Dirichlet}(a_1,\dots,a_k) \, .$$

The proof can be found on page 594 of Luc Devroye's book.

Therefore, one possibility is to compute a Monte Carlo approximation of $$F_{X_1,\dots,X_k}(t_1,\dots,t_k)=P\left\{X_1\leq t_1,\dots, X_k\leq t_k\right\} \, ,$$ starting with gammas. In R, try this:

pdirichlet <- function(a, t) {
N <- 10000
rdirichlet <- function(a) { y <- rgamma(length(a), a, 1); y / sum(y) }
x <- replicate(N, rdirichlet(a), simplify = FALSE)
sum(sapply(x, function(x) prod(x <= t))) / N
}


I didn't check the code. Use it carefully. If you find any bugs, please tell us.

• There's a vectorized rdirichlet function in R already - in fact several of them (in gtools, MCMCpack and dirmult for example). Apr 27, 2013 at 0:43
• @Zen I have a <- c(6, 20,2) how to obtain the Drichelt cdf? is t 2 by 2 matrix? Oct 4, 2018 at 13:33
• I used the above code, but its throwing an error. Oct 4, 2018 at 19:08
• @score324 A Dirichlet vector usually has each element in $[0,1]$ and their sum being $1$, so $(6, 20,2)$ is outside the support Mar 10, 2020 at 15:28

Any library? Mathematica has it. Here's the code for an example plot of a Dirichlet CDF from the documentation:

Plot3D[CDF[DirichletDistribution[{1, 3, 2}], {x, y}], {x, 0, 1}, {y, 0, 1}]

• How to get expression for CDF[DirichletDistribution[{1, 3, 2}]] ?
– AIB
May 20, 2013 at 10:30