Resamples particles.
0001 function indx = residual_resampling(weights) 0002 % Resamples particles. 0003 0004 %@info: 0005 %! @deftypefn {Function File} {@var{indx} =} residual_resampling (@var{weights}) 0006 %! @anchor{particle/residual_resampling} 0007 %! @sp 1 0008 %! Resamples particles. 0009 %! @sp 2 0010 %! @strong{Inputs} 0011 %! @sp 1 0012 %! @table @ @var 0013 %! @item weights 0014 %! n*1 vector of doubles, particles' weights. 0015 %! @end table 0016 %! @sp 2 0017 %! @strong{Outputs} 0018 %! @sp 1 0019 %! @table @ @var 0020 %! @item indx 0021 %! n*1 vector of intergers, indices. 0022 %! @end table 0023 %! @sp 2 0024 %! @strong{This function is called by:} 0025 %! @sp 1 0026 %! @ref{particle/resample} 0027 %! @sp 2 0028 %! @strong{This function calls:} 0029 %! @sp 2 0030 %! @end deftypefn 0031 %@eod: 0032 0033 % Copyright (C) 2011, 2012 Dynare Team 0034 % 0035 % This file is part of Dynare. 0036 % 0037 % Dynare is free software: you can redistribute it and/or modify 0038 % it under the terms of the GNU General Public License as published by 0039 % the Free Software Foundation, either version 3 of the License, or 0040 % (at your option) any later version. 0041 % 0042 % Dynare is distributed in the hope that it will be useful, 0043 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0044 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0045 % GNU General Public License for more details. 0046 % 0047 % You should have received a copy of the GNU General Public License 0048 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0049 0050 % AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr 0051 % stephane DOT adjemian AT univ DASH lemans DOT fr 0052 0053 % What is the number of particles? 0054 number_of_particles = length(weights); 0055 0056 % Set vectors of indices. 0057 jndx = 1:number_of_particles; 0058 indx = zeros(1,number_of_particles); 0059 0060 % Multiply the weights by the number of particles. 0061 WEIGHTS = number_of_particles*weights; 0062 0063 % Compute the integer part of the normalized weights. 0064 iWEIGHTS = fix(WEIGHTS); 0065 0066 % Compute the number of resample 0067 number_of_trials = number_of_particles-sum(iWEIGHTS); 0068 0069 if number_of_trials 0070 WEIGHTS = (WEIGHTS-iWEIGHTS)/number_of_trials; 0071 EmpiricalCDF = cumsum(WEIGHTS); 0072 u = fliplr(cumprod(rand(1,number_of_trials).^(1./(number_of_trials:-1:1)))); 0073 j=1; 0074 for i=1:number_of_trials 0075 while (u(i)>EmpiricalCDF(j)) 0076 j=j+1; 0077 end 0078 iWEIGHTS(j)=iWEIGHTS(j)+1; 0079 end 0080 end 0081 0082 k=1; 0083 for i=1:number_of_particles 0084 if (iWEIGHTS(i)>0) 0085 for j=k:k+iWEIGHTS(i)-1 0086 indx(j) = jndx(i); 0087 end 0088 end 0089 k = k + iWEIGHTS(i); 0090 end 0091 0092 %@test:1 0093 %$ % Define the weights 0094 %$ weights = randn(2000,1).^2; 0095 %$ weights = weights/sum(weights); 0096 %$ % Initialize t. 0097 %$ t = ones(1,1); 0098 %$ 0099 %$ try 0100 %$ indx1 = residual_resampling(weights); 0101 %$ catch 0102 %$ t(1) = 0; 0103 %$ end 0104 %$ 0105 %$ T = all(t); 0106 %@eof:1 0107 0108 %@test:2 0109 %$ % Define the weights 0110 %$ weights = exp(randn(2000,1)); 0111 %$ weights = weights/sum(weights); 0112 %$ % Initialize t. 0113 %$ t = ones(1,1); 0114 %$ 0115 %$ try 0116 %$ indx1 = residual_resampling(weights); 0117 %$ catch 0118 %$ t(1) = 0; 0119 %$ end 0120 %$ 0121 %$ T = all(t); 0122 %@eof:2