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