function X=multsamp(N,theta)
%MULTSAMP   Samples from a multinomial, parameterized by theta.
%   This is probably a naive implementation as it is O(N).
%       (Probably could be O(M), where M=length(theta).)
%
%   TODO: error checking
%
%   Arguments:
%       N      Number of times to sample the multinomial.
%       theta  Multinomial parameters.
%
%   Returns:
%       X      samples
%
%   Authors:
%       Joshua V. Dillon; Feb. 21, 2007; jvdillon AAT purdue DDOT edu

% jvd: always fix a broken theta
if sum(theta)~=1
	% jvd: only warn if its really bad
	if abs(1-sum(theta))>2*eps
		warning('normalizing theta');
	end
	theta = theta./sum(theta);
end

M = length(theta);

% replicated for each of N (naive?)
foo   = (repmat(rand(1,N),M,1)<repmat(cumsum(theta),1,N))';
[i j] = find( xor([foo ones(N,1)],[zeros(N,1) foo]) );
X = histc(j,1:M);