% 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, on logarithmic scale
logc = sum(gammaln(A)) - gammaln(sum(A));

% compute the logarithmic density
logf = sum(log(P).*(A-1), 2) - logc;

% convert to actual density
f = exp(logf);

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