10

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?

Ricky Robinson
  • 489
  • 5
  • 16
  • 1
    Not directly aware of any. But there may be something that can be done. What do you need to do with it? – Glen_b Apr 25 '13 at 22:58
  • 2
    I need to take the complementary of the CDF and consider it as my p value. – Ricky Robinson Apr 26 '13 at 08:40
  • 3
    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. – Glen_b Apr 27 '13 at 00:50

2 Answers2

10

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

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 = 10^5) {
    rdirichlet <- function(a) { y <- rgamma(length(a), a, 1); y / sum(y) }
    x <- replicate(N, rdirichlet(a), simplify = FALSE)
    mean(sapply(x, function(x) prod(x <= t)))
}
Zen
  • 24,121
  • 3
    There's a vectorized rdirichlet function in R already - in fact several of them (in gtools, MCMCpack and dirmult for example). – Glen_b Apr 27 '13 at 00:43
  • @Zen I have a <- c(6, 20,2) how to obtain the Drichelt cdf? is t 2 by 2 matrix? – score324 Oct 04 '18 at 13:33
  • I used the above code, but its throwing an error. – score324 Oct 04 '18 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 – Henry Mar 10 '20 at 15:28
2

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}]
Mike Z.
  • 209