root/applications/dual/SIDP/smc29/resample_residual.m @ 1287

Revision 1256, 1.3 kB (checked in by smidl, 14 years ago)

resampling procedures

Line 
1function outIndex = resample_residual(inIndex,q);
2% PURPOSE : Performs the resampling stage of the SIR
3%           in order(number of samples) steps. It uses Liu's
4%           residual resampling algorithm and Niclas' magic line.
5% INPUTS  : - inIndex = Input particle indices.
6%           - q = Normalised importance ratios.
7% OUTPUTS : - outIndex = Resampled indices.
8
9% AUTHORS  : Arnaud Doucet, Nando de Freitas and Neil Gordon
10% DATE     : 08-09-98
11
12if nargin < 2, error('Not enough input arguments.'); end
13
14[S,arb] = size(q);  % S = Number of particles.
15
16% RESIDUAL RESAMPLING:
17% ====================
18N_babies= zeros(1,S);
19% first integer part
20q_res = S.*q';
21N_babies = fix(q_res);
22% residual number of particles to sample
23N_res=S-sum(N_babies);
24if (N_res~=0)
25  q_res=(q_res-N_babies)/N_res;
26  cumDist= cumsum(q_res);   
27  % generate N_res ordered random variables uniformly distributed in [0,1]
28  u = fliplr(cumprod(rand(1,N_res).^(1./(N_res:-1:1))));
29  j=1;
30  for i=1:N_res
31    while (u(1,i)>cumDist(1,j))
32      j=j+1;
33    end
34    N_babies(1,j)=N_babies(1,j)+1;
35  end;
36end;
37
38% COPY RESAMPLED TRAJECTORIES: 
39% ============================
40index=1;
41for i=1:S
42  if (N_babies(1,i)>0)
43    for j=index:index+N_babies(1,i)-1
44      outIndex(j) = inIndex(i);
45    end;
46  end;   
47  index= index+N_babies(1,i);   
48end
49
50
51
52
53
54
55
56
57
58
59
60
Note: See TracBrowser for help on using the browser.