#--COMPARISON OF DIFFERENT FILTERS(RANDOM WALK)--
#two cases: (sigma, sigma.meas)=(3,1) or (1,3)
# all methods resample every time
#load library
source("filter_mod.R")
#set model variables
x0<-0
tau<-100
sigma<-1 #1
sigma.meas<-3 #3
N<-500
M<-100
#generate states and measurements
set.seed(123)
x<-rand.walk.1D(tau=tau,x0=x0,sigma=sigma)
#generate simulated data y
y.mat<-matrix(NA,nrow=tau,ncol=M)
for(k in 1:M){
y.mat[,k]<-rand.y.1D(x,sigma.meas=sigma.meas)
}
#------------------------------------BOOTSTRAP---------------------------------------
# resamples at every step
# run filter M times
MSE.boot.mat<-matrix(NA,tau,M)
N.eff.boot<-matrix(NA,tau,M)
for(k in 1:M){
obj.boot<-particle.filter(N=N,x=x,y=y.mat[,k],x0=5,sigma=sigma,sigma.meas=sigma.meas,resample.type="standard",N.thr.per=1)
MSE.boot.mat[,k]<-obj.boot$MSE.k.out
N.eff.boot[,k]<-obj.boot$N.eff.out
}
#compute MSE and AMSE
MSE.boot<-colMeans(MSE.boot.mat)
AMSE.boot<-mean(MSE.boot)
MSE.boot<-rowMeans(MSE.boot.mat)
N.eff.boot<-rowMeans(N.eff.boot)
#assign mean
m.boot<-obj.boot$m.out
#------------------------------------Exact---------------------------------------
# no resampling for this filter
# run filter M times
MSE.exact.mat<-matrix(NA,tau,M)
N.eff.exact<-matrix(NA,tau,M)
for(k in 1:M){
obj.exact<-exact.filter(N=N,x=x,y=y.mat[,k],x0=5,sigma=sigma,sigma.meas=sigma.meas,resample.type="standard")
MSE.exact.mat[,k]<-obj.exact$MSE.k.out
N.eff.exact[,k]<-obj.exact$N.eff.out
}
#compute MSE and AMSE
MSE.exact<-colMeans(MSE.exact.mat)
AMSE.exact<-mean(MSE.exact)
MSE.exact<-rowMeans(MSE.exact.mat)
N.eff.exact<-rowMeans(N.eff.exact)
#assign mean
m.exact<-obj.exact$m.out
#------------------------------------Optimum---------------------------------------
# resample 50% of the time
#compute MSE
MSE.opt.mat<-matrix(NA,tau,M)
N.eff.opt<-matrix(NA,tau,M)
for(k in 1:M){
obj.opt<-optimum.filter(N=N,x=x,y=y.mat[,k],x0=5,sigma=sigma,sigma.meas=sigma.meas,resample.type="standard",N.thr=1)
MSE.opt.mat[,k]<-obj.opt$MSE.k.out
N.eff.opt[,k]<-obj.opt$N.eff.out
}
#compute MSE and AMSE
MSE.opt<-colMeans(MSE.opt.mat)
AMSE.opt<-mean(MSE.opt)
MSE.opt<-rowMeans(MSE.opt.mat)
N.eff.opt<-rowMeans(N.eff.opt)
#assign mean
m.opt<-obj.opt$m.out
#------------------------------------Likelihood---------------------------------------
# resample 50% of the time
#compute MSE
MSE.like.mat<-matrix(NA,tau,M)
N.eff.like<-matrix(NA,tau,M)
for(k in 1:M){
obj.like<-likelihood.filter(N=N,x=x,y=y.mat[,k],x0=5,sigma=sigma,sigma.meas=sigma.meas,resample.type="standard",N.thr.per=1)
MSE.like.mat[,k]<-obj.like$MSE.k.out
N.eff.like[,k]<-obj.like$N.eff.out
}
#compute MSE and AMSE
MSE.like<-colMeans(MSE.like.mat)
AMSE.like<-mean(MSE.like)
AMSE.like
MSE.like<-rowMeans(MSE.like.mat)
N.eff.like<-rowMeans(N.eff.like)
#assign mean
m.like<-obj.like$m.out
#------------------------------------RESULTS---------------------------------------
#------(sigma, sigma.meas)=(3,1)-------
#display AMSE
AMSE.like #2.79889(1,3) 0.8977392(3,1)
AMSE.opt #2.611359 (1,3) 0.8954476 (3,1)
AMSE.boot #2.625013(1,3) 0.8964912 (3,1)
AMSE.exact #2.626761 (1,3) 0.894724 (3,1)
#plot means(There is hardly any difference)
plot(x,main=paste("Random Walk (",N, "particles)"),xlab="t")
lines(m.boot)
lines(m.opt,col="blue")
lines(m.like, col="red")
lines(m.exact, col="orange")
legend(20, 25, legend=c("Bootstrap", "PF Optimum","Likelihood","Exact"),
col=c("black", "blue","red","orange"),lty=c(1,1,1,1), cex=0.8)
#plot effective sample size
#obj.boot$N.eff.out
plot(N.eff.boot,ylab="N.eff",xlab="Time,t", main=" Effective Sample Size",type="l",ylim=c(1,500))
lines(N.eff.like,col="red")
lines(N.eff.opt,col="blue")
legend(70, 100, legend=c("Bootstrap", "PF Optimum","Likelihood"),
col=c("black", "blue","red"),lty=c(1,1,1,1), cex=0.8)
#plot MSE(There is hardly any difference)
plot(MSE.boot,main="MSE",xlab="Time,t",type="l",ylim=c(0,5))
lines(MSE.like,col="red")
lines(MSE.opt,col="blue")
lines(MSE.exact,col="orange")
#------(sigma, sigma.meas)=(1,3)-------
#display AMSE
AMSE.like #2.745584 (1,3) 0.8984577(3,1)
AMSE.opt #2.629193 (1,3) 0.8950922(3,1)
AMSE.boot #2.625013(1,3) 0.8964912(3,1)
AMSE.exact #2.626761(1,3) 0.894724(3,1)
#old results(keep just in case)
#AMSE.like #2.79889(1,3) 0.8977392(3,1) 0.6573997(1,1)
#AMSE.opt #2.611359 (1,3) 0.8954476 (3,1) 0.623635(1,1)
#AMSE.boot #2.625013(1,3) 0.8964912 (3,1) 0.6381266(1,1)
#AMSE.exact #2.626761 (1,3) 0.894724 (3,1) 0.6228132(1,1)
#plot means(There is hardly any difference)
plot(x,main=paste("Random Walk (",N, "particles)"),xlab="t")
lines(m.boot)
lines(m.opt,col="blue")
lines(m.like, col="red")
lines(m.exact, col="orange")
#plot effective sample size
plot(N.eff.boot,ylab="N.eff",xlab="Time,t", main=" Effective Sample Size",type="l",ylim=c(1,500))
lines(N.eff.like,col="red")
lines(N.eff.opt,col="blue")
legend(70, 100, legend=c("Bootstrap", "PF Optimum","Likelihood"),
col=c("black", "blue","red"),lty=c(1,1,1), cex=0.8)
#plot MSE(There is hardly any difference)
plot(MSE.boot,main="MSE",xlab="Time,t",type="l",ylim=c(0,5))
lines(MSE.like,col="red")
lines(MSE.opt,col="blue")
lines(MSE.exact,col="orange")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.