function zzz=SampleFromMultinomial(p)
% SampleFromMultinomial - Samples from a multinomial distribution
% If the elements of p do not sum to one, p gets normalised
%
% INPUT
% p(1:K) - array of probabilities
%
% OUTPUT
% zzz - Element in (1,2,...,K) = sample from the multinomial
%
% USAGE
% zzz=SampleFromMultinomial(p)

% -- Error checking --
if sum(p<0) > 0
	error('p contains negative elements');
end
Z=sum(p);
if Z==0
	error('All elements in p are zero.')
end
if Z~=1
	p=p/Z;
end

K= length(p);

cumulative_dist=zeros(1,K);
for i=1:K-1
	cumulative_dist(i+1)=cumulative_dist(i)+p(i);
end
% This is an array with treshold(1)=0, cumulative_dist(K)=1-p(K)

x=rand;

zzz=sum(x>=cumulative_dist);