rmvn <- function(mn,var) { tmp <- chol(var) mn + t(tmp)%*%rnorm(length(mn)) } ellip <-function(mn,var) { x <- cos(seq(from=0,to=2*pi,length=41)) y <- sin(seq(from=0,to=2*pi,length=41)) tmp <- chol(var)%*%rbind(x,y) + matrix(mn, 2, 41) points(tmp[1,], tmp[2,], type="l") } TYP<- 2 set.seed(13) ### `parameters' F <- diag(2); G <- matrix(c(1,.75,0,.5),2,2) V <- .2*diag(2); W <- 0.2*diag(2) ### initial values t.hat <- c(0,0); sigma <- diag(2) theta <- rmvn(t.hat, sigma) ### set-up plot postscript("kalman.ps") plot(0,0,xlim=c(-.5,1.5),ylim=c(-3,2), type="n") for (tm in 1:10) { ### evolve theta and y theta <- G %*% theta + rmvn(c(0,0),W) y <- rmvn(theta, V) text(theta[1],theta[2],paste("",as.character(tm),sep="")) ### filter R <- G %*% sigma %*% t(G) + W t.hat <- G %*% t.hat + R %*% t(F) %*% solve(F%*%R%*%t(F)+V) %*% (y-F%*%G%*%t.hat) sigma <- R - R%*%t(F)%*%solve(F%*%R%*%t(F)+V)%*%F%*%R if (TYP==1) { points(t.hat[1],t.hat[2], pch=1) points(y[1],y[2], pch=2) lines(c(theta[1],t.hat[1]), c(theta[2],t.hat[2]), lty=2) lines(c(theta[1],y[1]), c(theta[2],y[2]), lty=2) } if (TYP==2) { ellip(t.hat,sigma) } } graphics.off()