wolffd@0: function r = multirnd(theta,k)
wolffd@0: %MULTIRND - Random vector from multinomial distribution.
wolffd@0: %   r = multirnd(theta,k) returns a vector randomly selected
wolffd@0: %   from the multinomial distribution with parameter vector
wolffd@0: %   theta, and count k (i.e. sum(r) = k).
wolffd@0: %
wolffd@0: %   Note: if k is unspecified, then it is assumed k=1.
wolffd@0: %
wolffd@0: %   Author: David Ross
wolffd@0: %
wolffd@0: 
wolffd@0: %--------------------------------------------------------
wolffd@0: % Check the arguments.
wolffd@0: %--------------------------------------------------------
wolffd@0: error(nargchk(1,2,nargin));
wolffd@0: 
wolffd@0: % make sure theta is a vector
wolffd@0: if ndims(theta) > 2 | all(size(theta) > 1)
wolffd@0:     error('theta must be a vector');
wolffd@0: end
wolffd@0: 
wolffd@0: % if theta is a row vector, convert it to a column vector
wolffd@0: if size(theta,1) == 1
wolffd@0:     theta = theta';
wolffd@0: end
wolffd@0: 
wolffd@0: % make sure k is a scalar?
wolffd@0: 
wolffd@0: % if the number of samples has not been provided, set
wolffd@0: % it to one
wolffd@0: if nargin == 1
wolffd@0:     k = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: %--------------------------------------------------------
wolffd@0: % Main...
wolffd@0: %--------------------------------------------------------
wolffd@0: n = length(theta);
wolffd@0: theta_cdf = cumsum(theta);
wolffd@0: 
wolffd@0: r = zeros(n,1);
wolffd@0: random_vals = rand(k,1);
wolffd@0: 
wolffd@0: for j = 1:k
wolffd@0:     index = min(find(random_vals(j) <= theta_cdf));
wolffd@0:     r(index) = r(index) + 1;
wolffd@0: end