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