Commit 195ad9f7 authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Added Routines for resampling (particle filter).

parent 11134746
function indx = resample(weights,method)
% Resamples particles.
%@info:
%! @deftypefn {Function File} {@var{indx} =} resample (@var{weights}, @var{method})
%! @anchor{particle/resample}
%! @sp 1
%! Resamples particles.
%!
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item weights
%! n*1 vector of doubles, particles' weights.
%! @item method
%! string equal to 'residual' or 'traditional'.
%! @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:}
%! @ref{particle/sequantial_importance_particle_filter}
%! @sp 2
%! @strong{This function calls:}
%! @ref{residual_resampling}, @ref{traditional_resampling}
%! @end deftypefn
%@eod:
% Copyright (C) 2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr
% stephane DOT adjemian AT univ DASH lemans DOT fr
switch method
case 'residual'
indx = residual_resampling(weights);
case 'traditional'
indx = traditional_resampling(weights);
otherwise
error('particle::resample: Unknown method!')
end
\ No newline at end of file
function indx = residual_resampling(weights)
%@info:
%! @deftypefn {Function File} {@var{indx} =} resample (@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:}
%! @ref{particle/resample}
%! @sp 2
%! @strong{This function calls:}
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://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);
% Set vectors of indices.
jndx = 1:number_of_particles;
indx = zeros(1,number_of_particles);
% Multiply the weights by the number of particles.
WEIGHTS = number_of_particles*weights;
% Compute the integer part of the normalized weights.
iWEIGHTS = fix(WEIGHTS);
% Compute the number of resample
number_of_trials = number_of_particles-sum(iWEIGHTS);
if number_of_trials
WEIGHTS = (WEIGHTS-iWEIGHTS)/number_of_trials;
EmpiricalCDF = cumsum(WEIGHTS);
u = fliplr(cumprod(rand(1,number_of_trials).^(1./(number_of_trials:-1:1))));
j=1;
for i=1:number_of_trials
while (u(i)>EmpiricalCDF(j))
j=j+1;
end
iWEIGHTS(j)=iWEIGHTS(j)+1;
end
end
k=1;
for i=1:number_of_particles
if (iWEIGHTS(i)>0)
for j=k:k+iWEIGHTS(i)-1
indx(j) = jndx(i);
end
end
k = k + iWEIGHTS(i);
end
\ No newline at end of file
function indx = traditional_resampling(weights,noise)
% Resamples particles.
%@info:
%! @deftypefn {Function File} {@var{indx} =} resample (@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:}
%! @ref{particle/resample}
%! @sp 2
%! @strong{This function calls:}
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://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);
% Select method.
switch length(noise)
case 1
kitagawa_resampling = 1;
case number_of_particles
kitagawa_resampling = 0;
otherwise
error(['particle::resampling: Unknown method! The size of the second argument (' inputname(noise) ') is wrong.'])
end
% Get the empirical CDF.
c = cumsum(weights);
% Draw a starting point.
randvec = (transpose(1:number_of_particles)-1+noise(:))/number_of_particles ;
% Start at the bottom of the CDF
if kitagawa_resampling
j = 1;
for i=1:number_of_particles
while (randvec(i)>c(j))
j = j+1;
end
indx(i) = j;
end
else
for i=1:number_of_particles
indx(i) = sum(randvec(i)>c);
end
% Matlab's indices start at 1...
indx = indx+1;
end
%@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);
%@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);
%@eof:2
%@test:3
%$ % Set the number of particles.
%$ number_of_particles = 20000;
%$
%$ show_plot = 0;
%$ show_time = 0;
%$
%$ % 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
%$
%$ 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) = dyn_assert(indx1,indx1_);
%$ t(2) = dyn_assert(indx2,indx2_);
%$ T = all(t);
%@eof:3
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment