diff --git a/matlab/nonlinear-filters/traditional_resampling.m b/matlab/nonlinear-filters/traditional_resampling.m
index ab9089c8c6b92b1defcd05b177b4b64393d73ed9..fb0e2c2560ec679331525f2214ef4e1037ecd4aa 100644
--- a/matlab/nonlinear-filters/traditional_resampling.m
+++ b/matlab/nonlinear-filters/traditional_resampling.m
@@ -1,39 +1,16 @@
-function return_resample = traditional_resampling(particles,weights,noise)
-% Resamples particles.
-
-%@info:
-%! @deftypefn {Function File} {@var{indx} =} traditional_resampling (@var{weights},@var{noise})
-%! @anchor{particle/traditional_resampling}
-%! @sp 1
-%! Resamples particles (Resampling à la Kitagawa or stratified resampling).
-%! @sp 2
-%! @strong{Inputs}
-%! @sp 1
-%! @table @ @var
-%! @item weights
-%! n*1 vector of doubles, particles' weights.
-%! @item noise
-%! n*1 vector of doubles sampled from a [0,1] uniform distribution (stratified resampling) or scalar double
-%! sampled from a [0,1] uniform distribution (Kitagawa resampling).
-%! @end table
-%! @sp 2
-%! @strong{Outputs}
-%! @sp 1
-%! @table @ @var
-%! @item indx
-%! n*1 vector of intergers, indices.
-%! @end table
-%! @sp 2
-%! @strong{This function is called by:}
-%! @sp 1
-%! @ref{particle/resample}
-%! @sp 2
-%! @strong{This function calls:}
-%! @sp 2
-%! @end deftypefn
-%@eod:
-
-% Copyright © 2011-2017 Dynare Team
+function [indices, sample] = traditional_resampling(weights, noise, particles)
+
+% Resample particles.
+%
+% INPUTS:
+% - weights [double] n×1 vector of weights, sum(weights)==1.
+% - noise [double] scalar in [0,1] (kitagawa resampling) or n×1 vector of uniformly distributed number in [0,1] (stratified resampling)
+%
+% OUTPUTS:
+% - sample [double] n×1 vector, new sample of particles
+% - indices [integer] n×1 vector
+
+% Copyright © 2011-2025 Dynare Team
%
% This file is part of Dynare.
%
@@ -50,21 +27,18 @@ function return_resample = traditional_resampling(particles,weights,noise)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr
-% stephane DOT adjemian AT univ DASH lemans DOT fr
-
-% What is the number of particles?
number_of_particles = length(weights);
% Initialize the returned argument.
-indx = ones(number_of_particles,1);
+sample = NaN(number_of_particles, 1);
+indices = NaN(number_of_particles,1);
% Select method.
switch length(noise)
case 1
- kitagawa_resampling = 1;
+ kitagawa = true;
case number_of_particles
- kitagawa_resampling = 0;
+ kitagawa = false;
otherwise
error(['particle::resampling: Unknown method! The size of the second argument (' inputname(3) ') is wrong.'])
end
@@ -73,165 +47,169 @@ end
c = cumsum(weights);
% Draw a starting point.
-if kitagawa_resampling
+if kitagawa
randvec = (transpose(1:number_of_particles)-1+noise(:))/number_of_particles ;
else
randvec = fliplr(cumprod(noise.^(1./(number_of_particles:-1:1))));
end
% Start at the bottom of the CDF
-if kitagawa_resampling
+if kitagawa
j = 1;
for i=1:number_of_particles
while (randvec(i)>c(j))
j = j+1;
end
- indx(i) = j;
+ indices(i) = j;
end
else
for i=1:number_of_particles
- indx(i) = sum(randvec(i)>c);
+ indices(i) = sum(randvec(i)>c);
end
- % Matlab's indices start at 1...
- indx = indx+1;
+ indices = indices+1;
end
-if particles==0
- return_resample = indx ;
-else
- return_resample = particles(indx,:) ;
+if nargout>1
+ if nargin>2
+ sample = particles(indices,:);
+ else
+ error('Third input (a sample of particles) is missing.')
+ end
end
+
+return % --*-- Unit tests --*--
+
%@test:1
-%$ % Define the weights
-%$ weights = randn(2000,1).^2;
-%$ weights = weights/sum(weights);
-%$ % Initialize t.
-%$ t = ones(2,1);
-%$
-%$ % First, try the stratified resampling.
-%$ try
-%$ indx1 = traditional_resampling(weights,rand(2000,1));
-%$ catch
-%$ t(1) = 0;
-%$ end
-%$
-%$ % Second, try the Kitagawa resampling.
-%$ try
-%$ indx2 = traditional_resampling(weights,rand);
-%$ catch
-%$ t(2) = 0;
-%$ end
-%$
-%$ T = all(t);
+% Define the weights
+weights = randn(2000,1).^2;
+weights = weights/sum(weights);
+
+t = true(2,1);
+
+% First, try the stratified resampling.
+try
+ indx1 = traditional_resampling(weights, rand(2000,1));
+catch
+ t(1) = false;
+end
+
+% Second, try the Kitagawa resampling.
+try
+ indx2 = traditional_resampling(weights, rand);
+catch
+ t(2) = false;
+end
+
+T = all(t);
%@eof:1
%@test:2
-%$ % Define the weights
-%$ weights = exp(randn(2000,1));
-%$ weights = weights/sum(weights);
-%$ % Initialize t.
-%$ t = ones(2,1);
-%$
-%$ % First, try the stratified resampling.
-%$ try
-%$ indx1 = traditional_resampling(weights,rand(2000,1));
-%$ catch
-%$ t(1) = 0;
-%$ end
-%$
-%$ % Second, try the Kitagawa resampling.
-%$ try
-%$ indx2 = traditional_resampling(weights,rand);
-%$ catch
-%$ t(2) = 0;
-%$ end
-%$
-%$ T = all(t);
+% Define the weights
+weights = exp(randn(2000,1));
+weights = weights/sum(weights);
+
+t = true(2,1);
+
+% First, try the stratified resampling.
+try
+ indx1 = traditional_resampling(weights,rand(2000,1));
+catch
+ t(1) = false;
+end
+
+% Second, try the Kitagawa resampling.
+try
+ indx2 = traditional_resampling(weights,rand);
+catch
+ t(2) = false;
+end
+
+T = all(t);
%@eof:2
%@test:3
-%$ % Set the number of particles.
-%$ number_of_particles = 20000;
-%$
-%$ show_plot = 0;
-%$ show_time = 1;
-%$
-%$ % Define the weights
-%$ weights = randn(number_of_particles,1).^2;
-%$ weights = weights/sum(weights);
-%$
-%$ % Compute the empirical CDF
-%$ c = cumsum(weights);
-%$
-%$ % Stratified resampling.
-%$ noise = rand(number_of_particles,1);
-%$
-%$ if show_time
-%$ disp('Stratified resampling timing:')
-%$ tic
-%$ end
-%$
-%$ indx1 = traditional_resampling(weights,noise);
-%$
-%$ if show_time
-%$ toc
-%$ tic
-%$ end
-%$
-%$ indx1_ = zeros(number_of_particles,1);
-%$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles;
-%$ for i=1:number_of_particles
-%$ j = 1;
-%$ while (randvec(i)>c(j))
-%$ j = j + 1;
-%$ end
-%$ indx1_(i) = j;
-%$ end
-%$
-%$ if show_time
-%$ toc
-%$ end
-%$
-%$ % Kitagawa's resampling.
-%$ noise = rand;
-%$
-%$ if show_time
-%$ disp('Kitagawa''s resampling timing:')
-%$ tic
-%$ end
-%$
-%$ indx2 = traditional_resampling(weights,noise);
-%$
-%$ if show_time
-%$ toc
-%$ tic
-%$ end
-%$
-%$ indx2_ = zeros(number_of_particles,1);
-%$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles;
-%$ j = 1;
-%$ for i=1:number_of_particles
-%$ while (randvec(i)>c(j))
-%$ j = j + 1;
-%$ end
-%$ indx2_(i) = j;
-%$ end
-%$
-%$ if show_time
-%$ toc
-%$ end
-%$
-%$ % REMARK
-%$ % Note that the alternative code used in this test is sensibly faster than the code proposed
-%$ % in the routine for the resampling scheme à la Kitagawa...
-%$
-%$ if show_plot
-%$ plot(randvec,c,'-r'), hold on, plot([randvec(1),randvec(end)],[c(1),c(end)],'-k'), hold off, axis tight, box on
-%$ end
-%$
-%$ % Check results.
-%$ t(1) = dassert(indx1,indx1_);
-%$ t(2) = dassert(indx2,indx2_);
-%$ T = all(t);
-%@eof:3
\ No newline at end of file
+% Set the number of particles.
+number_of_particles = 20000;
+
+show_plot = false;
+show_time = true;
+
+% Define the weights
+weights = randn(number_of_particles,1).^2;
+weights = weights/sum(weights);
+
+% Compute the empirical CDF
+c = cumsum(weights);
+
+% Stratified resampling.
+noise = rand(number_of_particles,1);
+
+if show_time
+ disp('Stratified resampling timing:')
+ tic
+end
+
+indx1 = traditional_resampling(weights,noise);
+
+if show_time
+ toc
+ tic
+end
+
+indx1_ = zeros(number_of_particles,1);
+randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles;
+for i=1:number_of_particles
+ j = 1;
+ while (randvec(i)>c(j))
+ j = j + 1;
+ end
+ indx1_(i) = j;
+end
+
+if show_time
+ toc
+end
+
+% Kitagawa's resampling.
+noise = rand;
+
+if show_time
+ disp('Kitagawa''s resampling timing:')
+ tic
+end
+
+indx2 = traditional_resampling(weights,noise);
+
+if show_time
+ toc
+ tic
+end
+
+indx2_ = zeros(number_of_particles,1);
+randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles;
+j = 1;
+for i=1:number_of_particles
+ while (randvec(i)>c(j))
+ j = j + 1;
+ end
+ indx2_(i) = j;
+end
+
+if show_time
+ toc
+end
+
+% REMARK
+% Note that the alternative code used in this test is sensibly faster than the code proposed
+% in the routine for the resampling scheme à la Kitagawa...
+
+if show_plot
+ plot(randvec,c,'-r'), hold on, plot([randvec(1),randvec(end)],[c(1),c(end)],'-k'), hold off, axis tight, box on
+end
+
+% Check results.
+t(1) = dassert(indx1,indx1_);
+t(2) = dassert(indx2,indx2_);
+T = all(t);
+%@eof:3