重采样主要是为了解决经典蒙特卡洛方法中出现的粒子匮乏现象。其主要思想是对粒子和其相应的权值表示的概率密度函数重新进行采样。通过增加权值较大粒子和减少权值较小粒子来实现。重采样虽然可以改善粒子匮乏现象,但也降低了粒子的多样性。因此,重采样过程中一般选取一些准则来判断有效粒子的个数,通过这个个数来判断是否进行重采样。一般的判断准则为:
其中Neff为有效粒子个数,表示粒子权值,N为粒子个数。通过将Neff与预先设定的个数进行比较来决定是否重采样。一般Neff<2/3*N时候进行重采样。
重采样具体过程伪代码形式如下:
按照如此方法,权值比较大的粒子会被多次复制,那些权值小对计算后验概率函数贡献非常小的粒子就会被剔除。上面的过程其实就是相当于一个转转盘的过程,如下图所示:
这个转盘被分成了N分,即从w1到wn,和为1。每个区间就是一个粒子的权值,权值越大,区间的弧度就越大。每次都从w1左边的竖线开始转,显而易见,转到w1的概率为0-w1之间的随机数,转到w2的概率为w1到w1+w2之间的随机数,转到w3的概率为w1+w2到w1+w2+w3之间的随机数,后面的依次类推。当某个粒子的权值较大的时候,产生的随机数落在相应区间的概率就会增大,从而实现了对权值大粒子的多次复制,权值小的粒子的剔除。这样重采样的过程中不是全部复制大权值粒子,也有可能对小权值粒子进行了复制,因为区间虽小,但也有可能转到这个区间。不过这样在一定程度上也保证了粒子的多样性。
MATLAB代码如下:
function w_new=resample_particles1(w)
w_new=w;
Neff=1/sum(w.*w);
N=length(w);
if Neff<75 %75为预先设定阈值
for i = 1 : N
u = rand;
qtempsum = 0;
for j = 1 : N
qtempsum = qtempsum + w(j);
if qtempsum >= u
w_new(i)=w(j);
break;
end
end
end
end
今天就整理到这,如有问题,请指教。