function pdf = dirichlet_pdf ( x, n, a ) %*****************************************************************************80 % %% DIRICHLET_PDF evaluates the Dirichlet PDF. % % Discussion: % % PDF(X)(N,A) = Product ( 1 <= I <= N ) X(I)**( A(I) - 1 ) % * Gamma ( A_SUM ) / A_PROD % % where % % 0 < A(I) for all I; % 0 <= X(I) for all I; % Sum ( 1 <= I <= N ) X(I) = 1; % A_SUM = Sum ( 1 <= I <= N ) A(I). % A_PROD = Product ( 1 <= I <= N ) Gamma ( A(I) ) % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 30 October 2004 % % Author: % % John Burkardt % % Parameters: % % Input, real X(N), the argument of the PDF. Each X(I) should % be greater than 0.0, and the X(I)'s must add up to 1.0. % % Input, integer N, the number of components. % % Input, real A(N), the probabilities for each component. % Each A(I) should be positive. % % Output, real PDF, the value of the PDF. % tol = 0.0001; for i = 1 : n if ( x(i) <= 0.0 ) fprintf ( 1, '\n' ); fprintf ( 1, 'DIRICHLET_PDF - Fatal error!\n' ); fprintf ( 1, ' X(I) <= 0.\n' ); end end x_sum = sum ( x(1:n) ); if ( tol < abs ( x_sum - 1.0 ) ) fprintf ( 1, '\n' ); fprintf ( 1, 'DIRICHLET_PDF - Fatal error!\n' ); fprintf ( 1, ' SUM X(I) =/= 1.\n' ); end a_sum = sum ( a(1:n) ); a_prod = 1.0; for i = 1 : n a_prod = a_prod * gamma ( a(i) ); end pdf = gamma ( a_sum ) / a_prod; for i = 1 : n pdf = pdf * x(i)^( a(i) - 1.0 ); end return end