function ix = genmul ( n, p, ncat ) %*****************************************************************************80 % %% GENMUL generates a multinomial random deviate. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 31 March 2013 % % Author: % % Original FORTRAN77 version by Barry Brown, James Lovato. % MATLAB version by John Burkardt. % % Reference: % % Luc Devroye, % Non-Uniform Random Variate Generation, % Springer, 1986, % ISBN: 0387963057, % LC: QA274.D48. % % Parameters: % % Input, integer N, the number of events, which will be % classified into one of the NCAT categories. % % Input, real P(NCAT-1). P(I) is the probability that an event % will be classified into category I. Thus, each P(I) must be between % 0.0 and 1.0. Only the first NCAT-1 values of P must be defined since % P(NCAT) would be 1.0 minus the sum of the first NCAT-1 P's. % % Input, integer NCAT, the number of categories. % % Output, integer IX(NCAT), a random observation from % the multinomial distribution. All IX(i) will be nonnegative and their % sum will be N. % if ( n < 0 ) fprintf ( 1, ' \n' ); fprintf ( 1, 'GENMUL - Fatal error!\n' ); fprintf ( 1, ' N < 0\n' ); error ( 'GENMUL - Fatal error!\n' ); end if ( ncat <= 1 ) fprintf ( 1, ' \n' ); fprintf ( 1, 'GENMUL - Fatal error!\n' ); fprintf ( 1, ' NCAT <= 1\n' ); error ( 'GENMUL - Fatal error!\n' ); end do i = 1, ncat - 1 if ( p(i) < 0.0 ) fprintf ( 1, ' \n' ); fprintf ( 1, 'GENMUL - Fatal error!\n' ); fprintf ( 1, ' Some P(i) < 0.\n' ); error ( 'GENMUL - Fatal error!\n' ); end if ( 1.0 < p(i) ) fprintf ( 1, ' \n' ); fprintf ( 1, 'GENMUL - Fatal error!\n' ); fprintf ( 1, ' Some 1 < P(i).\n' ); error ( 'GENMUL - Fatal error!\n' ); end end ptot = 0.0; for i = 1 : ncat - 1 ptot = ptot + p(i); end if ( 0.99999 < ptot ) fprintf ( 1, ' \n' ); fprintf ( 1, 'GENMUL - Fatal error!\n' ); fprintf ( 1, ' 1 < Sum of P().\n' ); error ( 'GENMUL - Fatal error!\n' ); end % % Initialize variables. % ntot = n; ptot = 1.0; ix = zeros ( ncat, 1 ); % % Generate the observation. % for icat = 1 : ncat - 1 prob = p(icat) / ptot; ix(icat) = ignbin ( ntot, prob ); ntot = ntot - ix(icat); if ( ntot <= 0 ) return end if ptot = ptot - p(icat); end ix(ncat) = ntot; return end