blob: f11074db15cb61bb422924220f78c9ba0db7d72a (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
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: dirichletrnd.m,v 1.1.1.1 2005/05/22 23:32:12 yozhik Exp $
%-------------------------------------------------
% 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;
|
