6
$\begingroup$

First, let's assume that we want to generated from a Dirichlet(1,1,1,1) distribution. Would the following method be correct?:

  • generate three variates from a Uniform(0,1). Call them $x_1$, $x_2$, $x_3$.
  • then, order these such that $0 \leq x_{(1)} \leq x_{(2)} \leq x_{(3)} \leq 1$
  • then, return the differences as our Dirichlet variate: $(x_{(1)}, x_{(2)}-x_{(1)}, x_{(3)}-x_{(2)}, 1-x_{(3)})$

Is this correct? I have a feeling it is correct, but I'm not sure and this does not seem to be the same as either method described on Wikipedia or any other search I did. Maybe it is slow or has other problems, but I'm curious if it is correct.

Assuming this is correct, can it be extended to non-uniform Dirichlets, such as Dirichlet(a,b,c,d)?

Extra note: I am not simply asking how to generate a Dirichlet; there is plenty of information about that already. I'm just curious to see if the method for uniforms can be extended. Is there a more general method that involves drawing from a distribution, then ordering those numbers, then using the gaps?

Improve this question
$\endgroup$
6
  • $\begingroup$ I think this is exactly the "method described elsewhere", just in disguise. Hint: What do you know about the Poisson process? $\endgroup$ – cardinal Feb 23 '12 at 23:06
  • $\begingroup$ possible duplicate of Generate uniformly distributed weights that sum to unity? $\endgroup$ – cardinal Feb 23 '12 at 23:18
  • $\begingroup$ @cardinal, I can't see the relevance of Poisson. The probability of the event happening in a Poisson process in independent of how long it has been since the last event. One of the methods on Wikipedia is based on Gamma, which I know is the conjugate prior for Poisson - but I can't really see any direct connection. I can see some vague connections, but nothing that helps with Dirichlets other than Dirichlet(1,1,...) $\endgroup$ – Aaron McDaid Feb 23 '12 at 23:19
  • 1
    $\begingroup$ Specifically, there is a property of the homogeneous Poisson process that often goes by the name of the order-statistic property. I suspect you should immediately see its relevance. In fact it can be viewed as the intrinsic property that leads to the definition of a Poisson process on more general spaces than the nonnegative half-line. $\endgroup$ – cardinal Feb 23 '12 at 23:25
  • $\begingroup$ The previous comments were not intended to address your last question, only the first. $\endgroup$ – cardinal Feb 23 '12 at 23:26
6
$\begingroup$

If $Y_i$ are independent $\mathrm{Gamma}(\alpha_i,\beta)$, 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}(\alpha_1,\dots,\alpha_k) \, .$$ So, in R just do something like

rdirichlet <- function(a) {
    y <- rgamma(length(a), a, 1)
    return(y / sum(y))
}

And use it uniformily

> rdirichlet(c(1, 1, 1, 1))
[1] 0.40186737 0.03924152 0.37070316 0.18818796

or non-uniformily

> rdirichlet(c(3, 2.5, 9, 7))
[1] 0.1377426 0.1043081 0.4701179 0.2878314

The proof is given on page 594 of Luc Devroye's beautiful book:

http://luc.devroye.org/rnbookindex.html

P.S. Thanks to @cardinal for the R hacking tips.

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ You shouldn't need the sapply call, I don't believe; rgamma should allow a vector to be supplied to the parameter argument. In the special case described in the question, using rexp would be preferable. $\endgroup$ – cardinal Jun 23 '12 at 1:06
  • $\begingroup$ rgamma doesn't seem to vectorize automatically. $\endgroup$ – Zen Jun 23 '12 at 1:22
  • 1
    $\begingroup$ Instead of passing 1 as the first argument, did you try passing length(a)? $\endgroup$ – cardinal Jun 23 '12 at 14:09
  • 1
    $\begingroup$ That works! Cool! Edited. $\endgroup$ – Zen Jun 23 '12 at 17:44
  • $\begingroup$ I'm not simply looking for methods to generate a Dirichlet, it is not difficult to find information on that. I'm specifically asking if there is a method which involves drawing a set of numbers, ordering them, and using the gaps. Such a method works for the uniform case, I wonder if there is a more general version of this. This is just a curiousity - in hindsight I think I shouldn't be wasting all our time with this :-) $\endgroup$ – Aaron McDaid Jun 25 '12 at 10:37
1
$\begingroup$

Let $U_{(1)},U_{(2)}, \ldots, U_{(n)}$ be the order statistics from $U(0,1)$ distribution. Let $W_0 = U_{(1)},W_i = U_{(i+1)}-U{(i)}$, $1 \leq i \leq n-1$.

Then $(W_0, W_1, \ldots, W_{n-1}) \sim \cal{D}(\alpha)$, a dirichlet distribution with parameter $\alpha_{n+1,1} = (1,1,\ldots,1)$. That is, $(W_0, W_1, \ldots, W_{n-1})$ is uniformly distributed in the n-dimensional simplex.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ Something's not quite right here, because $U_{(n)}\ne1$ almost surely and therefore $W_0+W_1+\cdots+W_{n-1}=U_{(n)}$ does not even lie on the unit simplex. $\endgroup$ – whuber♦ Feb 23 '12 at 18:16
  • 3
    $\begingroup$ @VitalStatistix, I think your answer is just a rephrasing of my question. Do you have a reference or proof to back up your answer? $\endgroup$ – Aaron McDaid Feb 23 '12 at 21:24
  • 1
    $\begingroup$ @Aaron. This is a theorem from Probability for Statistics and Machine Learning by Anirban Dasgupta, Theorem 6.6. Consider the joint density of order statistics of uniform $f(U_{(1)}, \ldots, U_{(n)})= n!$. Now make the transformations $W_0=U{(1)}, W_{(i)}=U_{(i+1)}-U_{(i)}, 1\leq i \leq n-1$, the Jacobian of the transformation is 1. We can easily see from the Change of Variable theorem that the joint density of the $W_i$'s is the following: $f(w_0,w_2, ..., w_{n-1})= n! \{w_i\geq 0 \forall i; w_o+\ldots+w_{n-1} \leq 1\}$ - which we recognize as the Dirichlet density function. $\endgroup$ – VitalStatistix Feb 24 '12 at 14:06
  • $\begingroup$ @Whuber: Sorry for not being clear about what I mean by uniformly distributed on unit simplex. Let us consider the simplest Dirichlet density. $f(x_1, x_2) \propto {x_1}^{\alpha_1-1} {x_2}^{\alpha_2-1} x_3^{\alpha_3-1}$ where, $0 \leq x_1,x_2,x_3 \leq 1; x_1+x_2+x_3 =1 $. Here, $x_1,x_2,x_3$ lie in the unit simplex but the density is on $x_1,x_2$. Let me know if I could answer your question. $\endgroup$ – VitalStatistix Feb 24 '12 at 14:17
  • $\begingroup$ Can we generalize this to Dirichlet distribution with other $\alpha$ parameters? $\endgroup$ – Blade Aug 15 '20 at 15:40

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.