updated code when revising the tutorial

Former-commit-id: 3dd1310 [formerly 3dd1310 [formerly 5768bd9]]
Former-commit-id: a7e0f6ed90597448b38d8fd4639ab22380d9a418
Former-commit-id: 03ec09f

Author: Johan Dahlin

r/example1-lgss.R

@@ -45,53 +45,46 @@ x    = data$x
 y    = data$y
 
 # Export plot to file
-#cairo_pdf("example1-lgss.pdf", height = 10, width = 8)
+cairo_pdf("example1-lgss.pdf", height = 10, width = 8)
 
 # Plot the latent state and observations
-layout(matrix(1:3, 3, 1, byrow = TRUE));
+layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow = TRUE));
 par(mar=c(4,5,0,0))
 grid = seq(0,T);
 
 plot(grid,y,col="#1B9E77",lwd=1.5,type="l",xlab="time",ylab="observation",bty="n",ylim=c(-6,6))
 polygon(c(grid,rev(grid)),c(y,rep(-7,T+1)),border=NA,col=rgb(t(col2rgb("#1B9E77"))/256,alpha=0.25))
 
-plot(grid,x,col="#D95F02",lwd=1.5,type="l",xlab="time",ylab="latent state",bty="n",ylim=c(-6,6))
-polygon(c(grid,rev(grid)),c(x,rep(-7,T+1)),border=NA,col=rgb(t(col2rgb("#D95F02"))/256,alpha=0.25))
-
 ###################################################################################
-# State estimation using the particle filter
+# State estimation using the particle filter and Kalman filter
 ###################################################################################
 
-grid = seq(0,T-1);
-
 # Using N = 20 particles and plot the estimate of the latent state
 N      <- 20;
 outPF  <- sm(y,phi,sigmav,sigmae,N,T,x0)
-plot(grid,outPF$xh,col="#7570B3",lwd=1.5,type="l",xlab="time",ylab="state estimate",bty="n",ylim=c(-6,6))
-polygon(c(grid,rev(grid)),c(outPF$xh,rep(-7,T)),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
+outKF  <- kf(y,phi,sigmav,sigmae,x0,0.01);
+diff   <- outPF$xh-outKF$xh[-(T+1)];
 
-###################################################################################
-# State estimation using the Kalman filter
-###################################################################################
+grid = seq(0,T-1);
+plot(grid,diff,col="#7570B3",lwd=1.5,type="l",xlab="time",ylab="difference in state estimate",bty="n",ylim=c(-0.1,0.1))
+polygon(c(grid,rev(grid)),c(diff,rep(-0.10,T)),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
 
-outKF  <- kf(y,phi,sigmav,sigmae,x0,0.01);
-lines( grid, outKF$xh[-(T+1)], col="black", lty="dashed",lwd=1.5)
-
-#dev.off()
-
-# Compute bias
-# log( mean( abs(sm(y,phi,sigmav,sigmae,10,T,x0)$xh-outKF$xh[-(T+1)]) ) )
-# log( mean( abs(sm(y,phi,sigmav,sigmae,20,T,x0)$xh-outKF$xh[-(T+1)]) ) )
-# log( mean( abs(sm(y,phi,sigmav,sigmae,50,T,x0)$xh-outKF$xh[-(T+1)]) ) )
-# log( mean( abs(sm(y,phi,sigmav,sigmae,100,T,x0)$xh-outKF$xh[-(T+1)]) ) )
-# log( mean( abs(sm(y,phi,sigmav,sigmae,200,T,x0)$xh-outKF$xh[-(T+1)]) ) )
-
-# Compute MSE
-# log( mean( (sm(y,phi,sigmav,sigmae,10,T,x0)$xh-outKF$xh[-(T+1)])^2 ) )
-# log( mean( (sm(y,phi,sigmav,sigmae,20,T,x0)$xh-outKF$xh[-(T+1)])^2 ) )
-# log( mean( (sm(y,phi,sigmav,sigmae,50,T,x0)$xh-outKF$xh[-(T+1)])^2 ) )
-# log( mean( (sm(y,phi,sigmav,sigmae,100,T,x0)$xh-outKF$xh[-(T+1)])^2 ) )
-# log( mean( (sm(y,phi,sigmav,sigmae,200,T,x0)$xh-outKF$xh[-(T+1)])^2 ) )
+# Compute bias and MSE
+logBiasMSE = matrix( 0, nrow=7, ncol=2);
+gridN      = c( 10, 20, 50, 100, 200, 500, 1000);
+
+for ( ii in 1:length(gridN) ) {
+  logBiasMSE[ ii, 1 ] = log( mean( abs(sm(y,phi,sigmav,sigmae,gridN[ii],T,x0)$xh-outKF$xh[-(T+1)]) ) )  
+  logBiasMSE[ ii, 2 ] = log( mean( (sm(y,phi,sigmav,sigmae,gridN[ii],T,x0)$xh-outKF$xh[-(T+1)])^2 ) )
+}
+
+# Plot the bias and MSE for comparison
+plot(   gridN,logBiasMSE[ , 1 ],col="#E7298A",lwd=1.5,type="l",xlab="no. particles (N)",ylab="log-bias",bty="n",ylim=c(-7,-3))
+points( gridN,logBiasMSE[ , 1 ],col="#E7298A",pch=19)
+plot(   gridN,logBiasMSE[ , 2 ],col="#66A61E",lwd=1.5,type="l",xlab="no. particles (N)",ylab="log-MSE",bty="n",ylim=c(-12,-6))
+points( gridN,logBiasMSE[ , 2 ],col="#66A61E",pch=19)
+
+dev.off()
 
 ###################################################################################
 # End of file

r/example2-lgss.R

@@ -58,11 +58,10 @@ nPart    = 100;
 nBurnIn  = 1000;
 nIter    = 5000;
 
-# The standard deviation in the random walk proposal
-stepSize = 0.10;
-
 # Run the PMH algorithm
-res = pmh(data$y,initPar,sigmav,sigmae,nPart,T,x0,nIter,stepSize)
+res1 = pmh(data$y,initPar,sigmav,sigmae,nPart,T,x0,nIter,stepSize = 0.01)
+res2 = pmh(data$y,initPar,sigmav,sigmae,nPart,T,x0,nIter,stepSize = 0.10)
+res3 = pmh(data$y,initPar,sigmav,sigmae,nPart,T,x0,nIter,stepSize = 0.50)
 
 ##############################################################################
 # Plot the results
@@ -71,32 +70,66 @@ res = pmh(data$y,initPar,sigmav,sigmae,nPart,T,x0,nIter,stepSize)
 # Export plot to file
 #cairo_pdf("example2-lgss.pdf", height = 10, width = 8)
 
+layout(matrix(1:9, 3, 3, byrow = TRUE)); par(mar=c(4,5,0,0))
+
 # Plot the parameter posterior estimate
-# Solid black line indicate posterior mean
-layout(matrix(1:3, 3, 1, byrow = TRUE)); par(mar=c(4,5,0,0))
-hist(res[nBurnIn:nIter], breaks=floor(sqrt(nIter-nBurnIn)), col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25),border=NA,xlab=expression(phi),ylab="posterior estimate",main="",xlim=c(0.5,0.8),ylim=c(0,12),freq=F)
-lines(density(res[nBurnIn:nIter],kernel="e",from=0.5,to=0.8),lwd=2,col="#7570B3")
+hist(res1[nBurnIn:nIter], breaks=floor(sqrt(nIter-nBurnIn)), col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25),border=NA,xlab=expression(phi),ylab="posterior estimate",main="",xlim=c(0.5,0.8),ylim=c(0,12),freq=F)
+lines(density(res1[nBurnIn:nIter],kernel="e",from=0.5,to=0.8),lwd=2,col="#7570B3")
+lines(seq(0.5,1.0,0.01),dnorm(seq(0.5,1.0,0.01),0,1),lwd=1,col="darkgrey")
+abline(v=mean(res1[nBurnIn:nIter]),lwd=1,lty="dotted")
+
+hist(res2[nBurnIn:nIter], breaks=floor(sqrt(nIter-nBurnIn)), col=rgb(t(col2rgb("#E7298A"))/256,alpha=0.25),border=NA,xlab=expression(phi),ylab="posterior estimate",main="",xlim=c(0.5,0.8),ylim=c(0,12),freq=F)
+lines(density(res2[nBurnIn:nIter],kernel="e",from=0.5,to=0.8),lwd=2,col="#E7298A")
 lines(seq(0.5,1.0,0.01),dnorm(seq(0.5,1.0,0.01),0,1),lwd=1,col="darkgrey")
-abline(v=mean(res[nBurnIn:nIter]),lwd=1,lty="dotted")
+abline(v=mean(res2[nBurnIn:nIter]),lwd=1,lty="dotted")
+
+hist(res3[nBurnIn:nIter], breaks=floor(sqrt(nIter-nBurnIn)), col=rgb(t(col2rgb("#66A61E"))/256,alpha=0.25),border=NA,xlab=expression(phi),ylab="posterior estimate",main="",xlim=c(0.5,0.8),ylim=c(0,12),freq=F)
+lines(density(res3[nBurnIn:nIter],kernel="e",from=0.5,to=0.8),lwd=2,col="#66A61E")
+lines(seq(0.5,1.0,0.01),dnorm(seq(0.5,1.0,0.01),0,1),lwd=1,col="darkgrey")
+abline(v=mean(res3[nBurnIn:nIter]),lwd=1,lty="dotted")
+
+# Plot the trace of the Markov chain during 1000 iterations after the burn-in
+plot(seq(nBurnIn,nBurnIn+1000,1),res1[nBurnIn:(nBurnIn+1000)],col = '#7570B3', type="l",xlab="iteration",ylab=expression(phi),bty="n",ylim=c(0.5,0.8))
+grid=seq(nBurnIn,nBurnIn+1000,1)
+polygon(c(grid,rev(grid)),c(res1[nBurnIn:(nBurnIn+1000)],rep(0,length(res1[nBurnIn:(nBurnIn+1000)]))),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
+abline(h=mean(res1[nBurnIn:nIter]),lwd=1,lty="dotted")
+
+plot(seq(nBurnIn,nBurnIn+1000,1),res2[nBurnIn:(nBurnIn+1000)],col = '#E7298A', type="l",xlab="iteration",ylab=expression(phi),bty="n",ylim=c(0.5,0.8))
+grid=seq(nBurnIn,nBurnIn+1000,1)
+polygon(c(grid,rev(grid)),c(res2[nBurnIn:(nBurnIn+1000)],rep(0,length(res2[nBurnIn:(nBurnIn+1000)]))),border=NA,col=rgb(t(col2rgb("#E7298A"))/256,alpha=0.25))
+abline(h=mean(res2[nBurnIn:nIter]),lwd=1,lty="dotted")
 
-# Plot the trace of the Markov chain after burn-in
-# Solid black line indicate posterior mean
-plot(seq(nBurnIn,nIter,1),res[nBurnIn:nIter],col = '#7570B3', type="l",xlab="iteration",ylab=expression(phi),bty="n",ylim=c(0.5,0.8))
-grid=seq(nBurnIn,nIter,1)
-polygon(c(grid,rev(grid)),c(res[nBurnIn:nIter],rep(0,length(res[nBurnIn:nIter]))),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
-abline(h=mean(res[nBurnIn:nIter]),lwd=1,lty="dotted")
+plot(seq(nBurnIn,nBurnIn+1000,1),res3[nBurnIn:(nBurnIn+1000)],col = '#66A61E', type="l",xlab="iteration",ylab=expression(phi),bty="n",ylim=c(0.5,0.8))
+grid=seq(nBurnIn,nBurnIn+1000,1)
+polygon(c(grid,rev(grid)),c(res3[nBurnIn:(nBurnIn+1000)],rep(0,length(res3[nBurnIn:(nBurnIn+1000)]))),border=NA,col=rgb(t(col2rgb("#66A61E"))/256,alpha=0.25))
+abline(h=mean(res3[nBurnIn:nIter]),lwd=1,lty="dotted")
 
 # Plot the ACF of the Markov chain
-foo=acf( res[nBurnIn:nIter], plot=F)
-plot(foo$lag,foo$acf,col = '#7570B3', type="l",xlab="iteration",ylab="ACF",bty="n",lwd=2)
+foo=acf( res1[nBurnIn:nIter], plot=F, lag.max=60)
+plot(foo$lag,foo$acf,col = '#7570B3', type="l",xlab="iteration",ylab="ACF",bty="n",lwd=2,ylim=c(-0.2,1))
 polygon(c(foo$lag,rev(foo$lag)),c(foo$acf,rep(0,length(foo$acf))),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
 abline(h=1.96/sqrt(nIter-nBurnIn),lty="dotted")
 abline(h=-1.96/sqrt(nIter-nBurnIn),lty="dotted")
 
+foo=acf( res2[nBurnIn:nIter], plot=F, lag.max=60)
+plot(foo$lag,foo$acf,col = '#E7298A', type="l",xlab="iteration",ylab="ACF",bty="n",lwd=2,ylim=c(-0.2,1))
+polygon(c(foo$lag,rev(foo$lag)),c(foo$acf,rep(0,length(foo$acf))),border=NA,col=rgb(t(col2rgb("#E7298A"))/256,alpha=0.25))
+abline(h=1.96/sqrt(nIter-nBurnIn),lty="dotted")
+abline(h=-1.96/sqrt(nIter-nBurnIn),lty="dotted")
+
+
+foo=acf( res3[nBurnIn:nIter], plot=F, lag.max=60)
+plot(foo$lag,foo$acf,col = '#66A61E', type="l",xlab="iteration",ylab="ACF",bty="n",lwd=2,ylim=c(-0.2,1))
+polygon(c(foo$lag,rev(foo$lag)),c(foo$acf,rep(0,length(foo$acf))),border=NA,col=rgb(t(col2rgb("#66A61E"))/256,alpha=0.25))
+abline(h=1.96/sqrt(nIter-nBurnIn),lty="dotted")
+abline(h=-1.96/sqrt(nIter-nBurnIn),lty="dotted")
+
 #dev.off()
 
 # Estimate the parameter posterior mean
-mean( res[nBurnIn:nIter] )
+mean( res1[nBurnIn:nIter] )
+mean( res2[nBurnIn:nIter] )
+mean( res3[nBurnIn:nIter] )
 
 ##############################################################################
 # End of file