function x=multinomrnd(n,p)
% generates a multinomial (n,p)
x=[]; m=length(p); tot=0;
for i=1:m
    q=p(i)/sum(p(i:m));
    x=[x binornd(n-tot,q,1,1)];
    tot=tot+x(i);
end