Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

dirmultnom_pdf(n, a, b, c)
function [pdf] = dirmultnom_pdf(n, a, b, c)

%*******************************************************************************
%
%! DIRICHLET_MULTINOMIAL_PDF evaluates a Dirichlet Multinomial PDF.
%
%
%  Formula:
%
%    PDF(X)(A,B,C) = Comb(A,B,X) * ( Gamma(C_Sum) / Gamma(C_Sum+A) )
%      Product ( 1 <= I <= B ) Gamma(C(I)+X(I)) / Gamma(C(I))
%
%    where:
%
%      Comb(A,B,X) is the multinomial coefficient C( A; X(1), X(2), ..., X(B) ),
%      C_Sum = Sum ( 1 <= I <= B ) C(I)
%
%  Reference:
%
%    Kenneth Lange,
%    Mathematical and Statistical Methods for Genetic Analysis,
%    Springer, 1997, page 45.
%
%  Modified:
%
%    17 December 1999 (pdf evaluation for single vector, x), 
%    30 March 2011 (pdf evaluation for all permutation vectors of 'b' outcomes 
%    in 'a' trials)
%
%  Authors:
%
%    John Burkardt (Fortran 90); Jim Huntley (Matlab)
%
%  Parameters:
%
%    Input, integer X(B); X(I) counts the number of occurrences of
%    outcome I, out of the total of A trials.
%
%    Input, integer A, the total number of trials.
%
%    Input, integer B, the number of different possible outcomes on
%    one trial.
%
%    Input, integer C(B); C(I) is the Dirichlet parameter associated
%    with outcome I.
%
%    Output, real PDF, the value of the Dirichlet multinomial PDF.
%
%   NOTE: SINCE THE ORDER OF THE PERMUTATION VECTORS IS ARBITRARY,
%         THE PDF VALUES ARE SORTED INTO DESCENDING ORDER TO PRODUCE 
%         THE SMOOTHEST PDF AND CDF!!!
%
%    Calls: 'dist_nex.m'
%
  
  [d,jmax] = dist_nex(b,a);
  c_sum = sum( c(1:b) );
  pdf_log0 = - gammaln( c_sum + real( a ) ) + gammaln( c_sum ) ...
                + gammaln( real( a + 1 ) );

  for j = 1:jmax
      pdf_log = pdf_log0;
      x = d(j,:);      
      for i = 1: b
          pdf_log = pdf_log + gammaln( c(i) + real ( x(i) ) ) ...
                    - gammaln( c(i) ) - gammaln( real ( x(i) + 1 ) );
      end   
      pdf(j) = exp( pdf_log );
  end
  pdf = sort(pdf,'descend');
 
  return

Contact us