RANDDIRICHLET Function :: SAS/IML(R) 13.1 User's Guide

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 24.296: 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).