Resamples particles.
0001 function indx = traditional_resampling(weights,noise) 0002 % Resamples particles. 0003 0004 %@info: 0005 %! @deftypefn {Function File} {@var{indx} =} traditional_resampling (@var{weights},@var{noise}) 0006 %! @anchor{particle/traditional_resampling} 0007 %! @sp 1 0008 %! Resamples particles (Resampling à la Kitagawa or stratified resampling). 0009 %! @sp 2 0010 %! @strong{Inputs} 0011 %! @sp 1 0012 %! @table @ @var 0013 %! @item weights 0014 %! n*1 vector of doubles, particles' weights. 0015 %! @item noise 0016 %! n*1 vector of doubles sampled from a [0,1] uniform distribution (stratified resampling) or scalar double 0017 %! sampled from a [0,1] uniform distribution (Kitagawa resampling). 0018 %! @end table 0019 %! @sp 2 0020 %! @strong{Outputs} 0021 %! @sp 1 0022 %! @table @ @var 0023 %! @item indx 0024 %! n*1 vector of intergers, indices. 0025 %! @end table 0026 %! @sp 2 0027 %! @strong{This function is called by:} 0028 %! @sp 1 0029 %! @ref{particle/resample} 0030 %! @sp 2 0031 %! @strong{This function calls:} 0032 %! @sp 2 0033 %! @end deftypefn 0034 %@eod: 0035 0036 % Copyright (C) 2011, 2012 Dynare Team 0037 % 0038 % This file is part of Dynare. 0039 % 0040 % Dynare is free software: you can redistribute it and/or modify 0041 % it under the terms of the GNU General Public License as published by 0042 % the Free Software Foundation, either version 3 of the License, or 0043 % (at your option) any later version. 0044 % 0045 % Dynare is distributed in the hope that it will be useful, 0046 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0047 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0048 % GNU General Public License for more details. 0049 % 0050 % You should have received a copy of the GNU General Public License 0051 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0052 0053 % AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr 0054 % stephane DOT adjemian AT univ DASH lemans DOT fr 0055 0056 % What is the number of particles? 0057 number_of_particles = length(weights); 0058 0059 % Initialize the returned argument. 0060 indx = ones(number_of_particles,1); 0061 0062 % Select method. 0063 switch length(noise) 0064 case 1 0065 kitagawa_resampling = 1; 0066 case number_of_particles 0067 kitagawa_resampling = 0; 0068 otherwise 0069 error(['particle::resampling: Unknown method! The size of the second argument (' inputname(noise) ') is wrong.']) 0070 end 0071 0072 % Get the empirical CDF. 0073 c = cumsum(weights); 0074 0075 % Draw a starting point. 0076 randvec = (transpose(1:number_of_particles)-1+noise(:))/number_of_particles ; 0077 0078 % Start at the bottom of the CDF 0079 if kitagawa_resampling 0080 j = 1; 0081 for i=1:number_of_particles 0082 while (randvec(i)>c(j)) 0083 j = j+1; 0084 end 0085 indx(i) = j; 0086 end 0087 else 0088 for i=1:number_of_particles 0089 indx(i) = sum(randvec(i)>c); 0090 end 0091 % Matlab's indices start at 1... 0092 indx = indx+1; 0093 end 0094 0095 %@test:1 0096 %$ % Define the weights 0097 %$ weights = randn(2000,1).^2; 0098 %$ weights = weights/sum(weights); 0099 %$ % Initialize t. 0100 %$ t = ones(2,1); 0101 %$ 0102 %$ % First, try the stratified resampling. 0103 %$ try 0104 %$ indx1 = traditional_resampling(weights,rand(2000,1)); 0105 %$ catch 0106 %$ t(1) = 0; 0107 %$ end 0108 %$ 0109 %$ % Second, try the Kitagawa resampling. 0110 %$ try 0111 %$ indx2 = traditional_resampling(weights,rand); 0112 %$ catch 0113 %$ t(2) = 0; 0114 %$ end 0115 %$ 0116 %$ T = all(t); 0117 %@eof:1 0118 0119 %@test:2 0120 %$ % Define the weights 0121 %$ weights = exp(randn(2000,1)); 0122 %$ weights = weights/sum(weights); 0123 %$ % Initialize t. 0124 %$ t = ones(2,1); 0125 %$ 0126 %$ % First, try the stratified resampling. 0127 %$ try 0128 %$ indx1 = traditional_resampling(weights,rand(2000,1)); 0129 %$ catch 0130 %$ t(1) = 0; 0131 %$ end 0132 %$ 0133 %$ % Second, try the Kitagawa resampling. 0134 %$ try 0135 %$ indx2 = traditional_resampling(weights,rand); 0136 %$ catch 0137 %$ t(2) = 0; 0138 %$ end 0139 %$ 0140 %$ T = all(t); 0141 %@eof:2 0142 0143 %@test:3 0144 %$ % Set the number of particles. 0145 %$ number_of_particles = 20000; 0146 %$ 0147 %$ show_plot = 0; 0148 %$ show_time = 1; 0149 %$ 0150 %$ % Define the weights 0151 %$ weights = randn(number_of_particles,1).^2; 0152 %$ weights = weights/sum(weights); 0153 %$ 0154 %$ % Compute the empirical CDF 0155 %$ c = cumsum(weights); 0156 %$ 0157 %$ % Stratified resampling. 0158 %$ noise = rand(number_of_particles,1); 0159 %$ 0160 %$ if show_time 0161 %$ disp('Stratified resampling timing:') 0162 %$ tic 0163 %$ end 0164 %$ 0165 %$ indx1 = traditional_resampling(weights,noise); 0166 %$ 0167 %$ if show_time 0168 %$ toc 0169 %$ tic 0170 %$ end 0171 %$ 0172 %$ indx1_ = zeros(number_of_particles,1); 0173 %$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles; 0174 %$ for i=1:number_of_particles 0175 %$ j = 1; 0176 %$ while (randvec(i)>c(j)) 0177 %$ j = j + 1; 0178 %$ end 0179 %$ indx1_(i) = j; 0180 %$ end 0181 %$ 0182 %$ if show_time 0183 %$ toc 0184 %$ end 0185 %$ 0186 %$ % Kitagawa's resampling. 0187 %$ noise = rand; 0188 %$ 0189 %$ if show_time 0190 %$ disp('Kitagawa''s resampling timing:') 0191 %$ tic 0192 %$ end 0193 %$ 0194 %$ indx2 = traditional_resampling(weights,noise); 0195 %$ 0196 %$ if show_time 0197 %$ toc 0198 %$ tic 0199 %$ end 0200 %$ 0201 %$ indx2_ = zeros(number_of_particles,1); 0202 %$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles; 0203 %$ j = 1; 0204 %$ for i=1:number_of_particles 0205 %$ while (randvec(i)>c(j)) 0206 %$ j = j + 1; 0207 %$ end 0208 %$ indx2_(i) = j; 0209 %$ end 0210 %$ 0211 %$ if show_time 0212 %$ toc 0213 %$ end 0214 %$ 0215 %$ % REMARK 0216 %$ % Note that the alternative code used in this test is sensibly faster than the code proposed 0217 %$ % in the routine for the resampling scheme à la Kitagawa... 0218 %$ 0219 %$ if show_plot 0220 %$ plot(randvec,c,'-r'), hold on, plot([randvec(1),randvec(end)],[c(1),c(end)],'-k'), hold off, axis tight, box on 0221 %$ end 0222 %$ 0223 %$ % Check results. 0224 %$ t(1) = dyn_assert(indx1,indx1_); 0225 %$ t(2) = dyn_assert(indx2,indx2_); 0226 %$ T = all(t); 0227 %@eof:3