wolffd@0: function p = dirichletpdf(x, alpha)
wolffd@0: %DIRICHLETPDF Dirichlet probability density function.
wolffd@0: %   p = dirichletpdf(x, alpha) returns the probability of vector
wolffd@0: %   x under the Dirichlet distribution with parameter vector
wolffd@0: %   alpha.
wolffd@0: %
wolffd@0: %   Author: David Ross
wolffd@0: 
wolffd@0: %-------------------------------------------------
wolffd@0: % Check the input
wolffd@0: %-------------------------------------------------
wolffd@0: error(nargchk(2,2,nargin));
wolffd@0: 
wolffd@0: % enusre alpha is a vector
wolffd@0: if min(size(alpha)) ~= 1 | ndims(alpha) > 2 | length(alpha) == 1
wolffd@0:     error('alpha must be a vector');
wolffd@0: end
wolffd@0: 
wolffd@0: % ensure x is is a vector of the same size as alpha
wolffd@0: if any(size(x) ~= size(alpha))
wolffd@0:     error('x and alpha must be the same size');
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: %-------------------------------------------------
wolffd@0: % Main
wolffd@0: %-------------------------------------------------
wolffd@0: if any(x < 0)
wolffd@0:     p = 0;
wolffd@0: elseif sum(x) ~= 1
wolffd@0:     disp(['dirichletpdf warning: sum(x)~=1, but this may be ' ...
wolffd@0:         'due to numerical issues']);
wolffd@0:     p = 0;
wolffd@0: else
wolffd@0:     z = gammaln(sum(alpha)) - sum(gammaln(alpha)); 
wolffd@0:     z = exp(z);
wolffd@0: 
wolffd@0:     p = z * prod(x.^(alpha-1));
wolffd@0: end