% DIRPDF   Dirichlet distribution density
%
% P should be n*d matrix, with n rows of d-dimensional points
% A should be 1*d row vector of parameters.
%
function f=dirpdf(P, A)

[n,d] = size(P);

% compute normalization coefficient
c = exp(sum(gammaln(A)) - gammaln(sum(A)));

% compute the density
f = prod(P.^(A-1), 2) / c;

% set to zero outside the triangular area
f(any(P<0, 2)) = 0;
f(sum(P,2) > 1) = 0;