wolffd@0: function x = sample(p, n)
wolffd@0: % SAMPLE    Sample from categorical distribution.
wolffd@0: % Returns a row vector of integers, sampled according to the probability
wolffd@0: % distribution p.
wolffd@0: % Uses the stick-breaking algorithm.
wolffd@0: % Much faster algorithms are also possible.
wolffd@0: 
wolffd@0: if nargin < 2
wolffd@0:   n = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: cdf = cumsum(p(:));
wolffd@0: for i = 1:n
wolffd@0:   x(i) = sum(cdf < rand) + 1;
wolffd@0: end