%?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