Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

 

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