function s = resampstr(p);
%RESAMPSTR Stratified resampling
% S = RESAMPSTR(P) returns a new set of indices according to
% the probabilities P. P is array of probabilities, which are
% not necessarily normalized, though they must be non-negative,
% and not all zero. The size of S is the size of P.
%
% Default is to use no-sort resampling. For sorted resampling use
% [PS,PI]=SORT(P);
% S=PI(RESAMPSTR(PS));
% Sorted re-sampling is slower but has slightly smaller variance.
% Stratified resampling is unbiased, almost as fast as
% deterministic resampling (RESAMPDET), and has only slightly
% larger variance.
%
% In stratified resampling indices are sampled using random
% numbers u_j~U[(j-1)/n,j/n], where n is length of P. Compare
% this to simple random resampling where u_j~U[0,1]. See,
% Kitagawa, G., Monte Carlo Filter and Smoother for
% Non-Gaussian Nonlinear State Space Models, Journal of
% Computational and Graphical Statistics, 5(1):1-25, 1996.
%
% See also RESAMPSIM, RESAMPRES, RESAMPDET
% Copyright (c) 2003-2004 Aki Vehtari
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.
[m,n]=size(p);
mn=m.*n;
pn=p./sum(p(:)).*mn;
fpn=floor(pn);
s=zeros(m,n);
r=rand(1,mn);
k=0;
c=0;
for i=1:mn
c=c+pn(i);
if c>=1
a=floor(c);
c=c-a;
s(k+[1:a])=i;
k=k+a;
end
if k<mn && c>=r(k+1)
c=c-1;
k=k+1;
s(k)=i;
end
end