function i = resampling(q)

% Resample indeces using Multinomial sampling.
% q is a vector of probabilities.
% Returns a vector of indeces "i" (positive integer values).
% The function samples length(q) times with replacement from the discrete 
% distribution taking values on {1,2,...,length(q)} with associated probabilities {q(1),q(2),...}.
% In the context of sequential Monte Carlo this implies producing indeces corresponding 
% to the particles we will sample at next time.

% See https://en.wikipedia.org/wiki/Multinomial_distribution#Sampling_from_a_multinomial_distribution

N = length(q);
% Multinomial resampling
u = rand(N,1);
qc = cumsum(q);
qc = qc(:);
qc=qc/qc(N);
[~,ind1]=sort([u;qc]);
ind2=find(ind1<=N);
i=ind2'-(0:N-1);
end