www.pudn.com > HMM1.zip > dirichletrnd.m
function x = dirichletrnd(alpha)
>DIRICHLETRND Random vector from a dirichlet distribution.
> x = dirichletrnd(alpha) returns a vector randomly selected
> from the Dirichlet distribution with parameter vector alpha.
>
> The algorithm used is the following:
> For each alpha(i), generate a value s(i) with distribution
> Gamma(alpha(i),1). Now x(i) = s(i) / sum_j s(j).
>
> The above algorithm was recounted to me by Radford Neal, but
> a reference would be appreciated...
> Do the gamma parameters always have to be 1?
>
> Author: David Ross
> $Id$
>-------------------------------------------------
> Check the input
>-------------------------------------------------
error(nargchk(1,1,nargin));
if min(size(alpha)) ~= 1 | length(alpha) < 2
error('alpha must be a vector of length at least 2');
end
>-------------------------------------------------
> Main
>-------------------------------------------------
gamma_vals = gamrnd(alpha, ones(size(alpha)), size(alpha));
denom = sum(gamma_vals);
x = gamma_vals / denom;