Language Reference

RANDDIRICHLET Function

RANDDIRICHLET (N, Shape);

The RANDDIRICHLET function is part of the IMLMLIB library . The RANDDIRICHLET function generates a random sample from a Dirichlet distribution, which is a multivariate generalization of the beta distribution.

The input parameters are as follows:

N

is the number of observations to sample.

Shape

is a $1 \times (p+1)$ vector of shape parameters for the distribution, $\mbox{Shape}[i]>0$.

The RANDDIRICHLET function returns an $N \times p$ matrix that contains N random draws from the Dirichlet distribution.

If $X=\{ X_1 \, X_2 \, \ldots \, X_ p\} $ with $\sum _{i=1}^ p {X_ i} < 1$ and $X_ i>0$ follows a Dirichlet distribution with shape parameter $\alpha =\{ \alpha _1 \, \alpha _2 \, \ldots \, \alpha _{p+1}\} $, then

  • the probability density function for x is

    \[ f(x; \alpha ) = \frac{\Gamma (\sum _{i=1}^{p+1}{\alpha _ i})}{\prod _{i=1}^{p+1} \Gamma (\alpha _ i) } \prod _{i=1}^ p { {x_ i}^{\alpha _ i -1}(1-x_1-x_2- \ldots -x_ p)^{\alpha _{p+1}-1} } \]

  • if $p=1$, the probability distribution is a beta distribution.

  • if $\alpha _0 = \Sigma _{i=1}^{p+1} \alpha _ i$, then

    • the expected value of $X_ i$ is $\alpha _ i / \alpha _0$.

    • the variance of $X_ i$ is $\alpha _ i(\alpha _0-\alpha _ i) / (\alpha _0^2 (\alpha _0+1))$.

    • the covariance of $X_ i$ and $X_ j$ is $-\alpha _ i \alpha _ j / (\alpha _0^2 (\alpha _0+1))$.

The following example generates 1,000 samples from a two-dimensional Dirichlet distribution. Each row of the returned matrix x is a row vector sampled from the Dirichlet distribution. The following example computes the sample mean and covariance and compares them with the expected values: 

call randseed(1);
n = 1000;
Shape = {2, 1, 1};
x = RandDirichlet(n,Shape);
d = nrow(Shape)-1;
s = Shape[1:d];
Shape0 = sum(Shape);
Mean = s`/Shape0;
Cov = -s*s` / (Shape0##2*(Shape0+1));
/* replace diagonal elements with variance */
Variance = s#(Shape0-s) / (Shape0##2*(Shape0+1));
do i = 1 to d;
   Cov[i,i] = Variance[i];
end;

SampleMean = mean(x);
SampleCov = cov(x);
print SampleMean Mean, SampleCov Cov;

Figure 25.301: Estimated Mean and Covariance Matrix

SampleMean   Mean  
0.4992449 0.2485677 0.5 0.25

SampleCov   Cov  
0.0502652 -0.026085 0.05 -0.025
-0.026085 0.0393922 -0.025 0.0375 




            

            


            
For further details about sampling from the Dirichlet distribution, see Kotz, Balakrishnan, and Johnson (2000); Gentle (2003); or Devroye (1986).