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?

$\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$ 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.

$\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$ 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.

$\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$ 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$ 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$ 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

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .