function [prob] = DirichletDistribution2(alpha, Q, P)
%[prob] = DirichletDistribution2(alpha, Q, P)
%   Evaluate the density function of the Dirichlet distribution
%   with parameters alpha (scalar) and Q (vector) at the columns of
%   P (each column of P must be a vector of the same size as Q) 
%   and return the result as a row matrix with number of columns of
%   P entries. 
%   ($RCSfile: DirichletDistribution2.m,v $ $Revision: 1.1.1.1 $)

%%% Check arguments:
Q=Q(:); alpha=alpha(1,1);
if ( sum([ P(:); Q; alpha] <0) )
  error('DirichletDistribution2: arguments must be nonnegative.')
end
% sum(P) summiert P spaltenweise und ergibt Zeilenvektor mit den
% Summen. sum(P) ~= 1 ist dann ein boolscher Vektor und die Summe
% dar"uber ist 0 genau dann wenn alle Spaltensummen von P eins
% sind. 
if ( sum((sum(P) ~= 1)) || (sum(Q) ~= 1) )
  error('DirichletDistribution2: Sum of Q and of each column of P must be equal to 1.')
end
% size(P,1) gibt die Zeilenzahl von P zur"uck
if ( size(P,1)~=length(Q) )
  error('DirichletDistribution2: Number of rows of P must equal length of Q.')
end
%%% Do computation:
% Exponenten hat gleiche Gr"o"se wie P, denn (alpha*Q-1) ist ein
% Spaltenvektor der L"ange wie die Zeilenzahl von P und
% ones(1,size(P,2)) ergibt eine Matrix voller Einsen mit einer
% Zeilen und sovie Spalten wie P hat. 
Exponenten=(alpha*Q-1)*ones(1,size(P,2));
% prod() auf eine Matrix angewendet gibt einen Zeilenvektor mit den
% Spaltenprodukten zur"uck, genau wie hier n"otig.
prob=gamma(alpha)/prod(gamma(alpha*Q))*prod(P.^Exponenten);
return;