function [ x, seed ] = multinomial_sample ( a, b, c, seed ) %*****************************************************************************80 % %% MULTINOMIAL_SAMPLE samples the Multinomial PDF. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 08 October 2004 % % Author: % % John Burkardt % % Reference: % % Luc Devroye, % Non-Uniform Random Variate Generation, % Springer-Verlag, New York, 1986, page 559. % % Parameters: % % Input, integer A, the total number of trials. % 0 <= A. % % Input, integer B, the number of outcomes possible on one trial. % 1 <= B. % % Input, real C(B). C(I) is the probability of outcome I on % any trial. % 0.0 <= C(I) <= 1.0, % Sum ( 1 <= I <= B) C(I) = 1.0. % % Input, integer SEED, a seed for the random number generator. % % Output, integer X(B); X(I) is the number of % occurrences of event I during the N trials. % % Output, integer SEED, an updated seed for the random number generator. % ntot = a; sum2 = 1.0; x(1:b) = 0; for ifactor = 1 : b - 1 prob = c(ifactor) / sum2; % % Generate a binomial random deviate for NTOT trials with % single trial success probability PROB. % [ x(ifactor), seed ] = binomial_sample ( ntot, prob, seed ); ntot = ntot - x(ifactor); if ( ntot <= 0 ) x = x'; return end sum2 = sum2 - c(ifactor); end % % The last factor gets what's left. % x(b) = ntot; x = x'; return end