10
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ Not directly aware of any. But there may be something that can be done. What do you need to do with it? $\endgroup$
    – Glen_b
    Apr 25, 2013 at 22:58
  • 2
    $\begingroup$ I need to take the complementary of the CDF and consider it as my p value. $\endgroup$ Apr 26, 2013 at 8:40
  • 2
    $\begingroup$ 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. $\endgroup$
    – Glen_b
    Apr 27, 2013 at 0:50

2 Answers 2

9
$\begingroup$

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.

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

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}]
$\endgroup$
1
  • 1
    $\begingroup$ How to get expression for CDF[DirichletDistribution[{1, 3, 2}]] ? $\endgroup$
    – AIB
    May 20, 2013 at 10:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.