序贯蒙特卡洛

可以参考的sis(序贯蒙特卡洛),大家相互学习，一起加油！

% Demo Sequential Importance Sampling for linear gaussian model

T=100;
sw=2;
alpha=0.5;

% simulate data
x=zeros（1T）;
y=zeros（1T）;
x=randn（1）;
for k=2:T
    x（1k）=alpha*x（1k-1）+randn（1）;
end
y=x+sw*randn（1T）;

N=1000;

% prior
xs1=zeros（TN）;
lw1=zeros（TN）;
w1=zeros（TN）;
wnorm1=zeros（TN）;
ess1=zeros（1T）;

% optimal
xs2=zeros（TN）;
lw2=zeros（TN）;
w2=zeros（TN）;
wnorm2=zeros（TN）;
ess2=zeros（1T）;

% SIS using prior
for k=1:T
    if （k==1）
        xs1（k:）=randn（1N）;
        lw1（k:）=-0.5*log（2*pi*sw^2）.*ones（1N）-0.5*（y（1k）*ones（1N）-xs1（k:））.^2./（2*sw^2）;
    else
        xs1（k:）=alpha.*xs1（k-1:）+randn（1N）;
        lw1（k:）=lw1（k-1:）-0.5*log（2*pi*sw^2）.*ones（1N）-0.5*（y（1k）*ones（1N）-xs1（k:））.^2./（2*sw^2）;
    end
    lmax=max（lw1（k:））;
    w1（k:）=exp（lw1（k:）-lmax）;  % correct only up to a multiplicative factor for unnormalized weights
    wnorm1（k:）=w1（k:）./sum（w1（k:））;
    ess1（1k）=1/sum（wnorm1（k:）.^2）;
end

% SIS using optimal
ss2=sw^2/（sw^2+1）;
ss=sqrt（ss2）;
for k=1:T
    if （k==1）
        xs2（k:）=ss2*y（1k）*ones（1N）/sw^2+ss.*randn（1N）;
        lw2（k:）=-0.5*log（2*pi*（sw^2+1））.*ones（1N）-0.5*（y（1k）*ones（1N））.^2./（2*（sw^2+1））;
    else
        xs2（k:）=ss2.*（alpha.*xs2（k-1:）+y（1k）/sw^2）+ss.*randn（1N）;
        lw2（k:）=lw2（k-1:）-0.5*log（2*pi*（sw^2+1））.*ones（1N）-0.5*（y（1k）*ones（1N）-alpha.*xs2（k-1:））.^2./（2*sw^2）;
    end
    lmax=max（lw2（k:））;
    w2（k:）=exp（lw2（k:）-lmax）;  % correct only up to a multiplicative factor for unnormalized weights
    wnorm2（k:）=w2（k:）./sum（w2（k:））;
    ess2（1k）=1/sum（wnorm2（k:）.^2）;
end

std（lw1（T:））
std（lw2（T:））

figure（1）
plot（ess1‘b‘）
hold on
plot （ess2‘r‘）
hold on
plot（sum（xs1.*wnorm1）‘k‘）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       1776  2010-03-30 11:10  sis.m

----------- ---------  ---------- -----  ----

                 1776                    1