######################################################################### ## Particle filtering ### ######################################################################### ## This file contains some simple R-functions and commands for computing and ## illustrating the particle filter. No attempt has been made to optimize ## the code ##### Example of Andrade-Netto ####################### ##### (nonlinear AR(1), observation =squared states + noise) ##### Functions to be loaded ################################### f.andrade <- function(x,u,i) { ## Purpose: simulates one state transition ## x = vector of starting values ## u = vector of innovations ## i = time index 0.5*x + 25*x/(1+x^2) + 8*cos(1.2*i) + sqrt(10)*u } h.square <- function(x,y) {dnorm(y-0.05*x^2)} ## computes the likelihood for observation y and vector x of state values bal.sample <- function(size=length(prob),prob) { ## Purpose: generating a balanced sample from (1,2,..,length(prob)) ## ---------------------------------------------------------------------- ## Arguments: size=sample size, prob=vector of probabilites ## Output: Vector of multiplicites and indices of resampled particles ## ---------------------------------------------------------------------- ## Author: Hans-Ruedi Kuensch, Date: 9 Jun 2005, 14:22 w <- floor(size*cumsum(prob) + runif(1)) w <- w-c(0,w[-length(prob)]) list(mult=w,index=rep(1:length(prob),w)) } wtd.quantile <- function(x,wts) { ## Purpose: Computing the 10%, 50% and 90% quantile of a weighted ensemble ## ---------------------------------------------------------------------- ## Arguments: x=particles, wts=weights ## ---------------------------------------------------------------------- ## Author: Hans-Rudolf Kuensch, Date: 7 Apr 2009, 18:45 n <- length(x) index <- sort.list(x) x <- x[index] wts <- cumsum(wts[index]) i1 <- min((1:n)[wts>0.1*n]) i2 <- min((1:n)[wts>0.5*n]) i3 <- min((1:n)[wts>0.9*n]) c(x[i1],x[i2],x[i3]) } f.particle <- function(x0,y, f.state, h) { ## Purpose: Computing the basic particle filter with resampling only ## before propagation ## ---------------------------------------------------------------------- ## Arguments: x0=sample from the starting density ## y=observed values ## f.state=function describing the evolution of states ## x_t = f(x_{t-1},u_t,t) with u_t standard normal ## h=observation density (time homogeneous) ## Output: x[,t] = filter ensemble at time t-1, 1<=t<=length(y)+1 ## lambda[t,] = weights at time t-1 (to be used with x[t,]) ## ---------------------------------------------------------------------- ## Author: H. R. Kuensch, Date: 13 Mar 2005, 22:11 n.rep <- length(x0) # number of particles n <- length(y) # number of time points - 1 x <- matrix(x0,nrow=n.rep,ncol=n+1) # matrix to store all filter samples lambda <- matrix(1/n.rep,nrow=n.rep,ncol=n+1) #matrix to store all weights for (i in (1:n)){ # iteration over time index <- bal.sample(n.rep,lambda[,i])$index #resampling z <- x[index,i] x[,i+1] <- f.state(z,rnorm(n.rep),i) #propagation w <- h(x[,i+1],y[i]) #reweighting lambda[,i+1] <- w/sum(w) } list(x=x,lambda=lambda) } ########### Simulation of the model and plotting ################## n.t <- 100 #number of time points x <- rnorm(n.t+1) x[1] <- 5*x[1] for (i in (1:n.t)) x[i+1] <- f.andrade(x[i],x[i+1],i) y <- 0.05*x[-1]^2 + rnorm(n.t) plot(x[-1],type="l",xlab="time",ylab="state/observations",main= "Simulated states and observations of Andrade-Netto") lines(y,col=2) ############# Particle filter, step by step: #################### ## Each plot shows filter particles at current time, prediction ## particles at next time, likelihood and weights from next ## observation and filter particles at next time ## initialize n.rep <- 20 #number of particles (small for plotting reasons) x.filt <- 5*rnorm(n.rep) #starting values i <- 0 # current time i.filt <- rep(0,n.rep) # For plotting purposes only i.pred <- i.filt+1 ## propagate x.pred <- f.andrade(x.filt,rnorm(n.rep),i) x.mode <- sqrt(20*max(0,y[i+1])) # arg max h(y|x) x.range <- range(c(x.filt,x.pred,x.mode,-x.mode)) # plotting range plot(x.filt,i.filt,ylim=c(0,3),xlim=x.range,xlab="state",ylab="") #filter particles at time i arrows(x.filt,i.filt,x.pred,i.pred,length=0.1) # movement of filter particles from time i to time i+1 points(x.pred,i.pred) # prediction particles at time i+1 ## weight w <- h.square(x.pred,y[i+1]) w <- w/sum(w) # weights xx <- seq(x.range[1],x.range[2],by=0.1) lines(xx,1+0.8*h.square(xx,y[i+1])) # observation likelihood segments(x.pred,i.pred,x.pred,i.pred+w) # weights of filter particles at time i+1 ## resample resample <- bal.sample(n.rep,w) x.filt <- x.pred[resample$index] #equally weighted filter particles at time i+1 for (j in (1:max(resample$mult))) { draw <- resample$mult>=j points(x.pred[draw],rep(2+(j-1)*0.05,sum(draw))) } # filter particles stacked to show identical particles points(x[i+2],2,col=2) #true value of the state ## iterate i <- i+1 # then go back to the propagation step ##### Particle filter for all observation times ###### ## The particle approximation of the filtering distribution is computed and ## analyzed for all time points n.t <- length(y) n.rep <- 200 #number of particles x0 <- 5*rnorm(n.rep) #starting values x.filt <- f.particle(x0,y,f.state=f.andrade,h=h.square) x.f <- x.filt$x lambda.f <- x.filt$lambda ## Filter distributions for the first 20 time points: ## Weighted filter distribution and posterior mode based on observation only ## Observe varying shape of filter distribution and balance of weights par(mfrow=c(4,5)) for (i in (2:21)) { plot(x.f[,i],lambda.f[,i],type="h") abline(v=sqrt(max(0,20*y[i-1]))*c(1,-1),col=2) } ## Plotting the 10%, 50%, 90% quantiles x.filtq <- matrix(0,nrow=3,ncol=n.t+1) for (i in (1:(n.t+1))) x.filtq[,i] <- wtd.quantile(x.f[,i],n.rep*lambda.f[,i]) par(mfrow=c(1,1)) plot(0:n.t,x,type="l",xlab="time",ylab="state",lwd=3, main="True state, filter median and 80% filter interval") polygon(c(0:n.t,n.t:0),c(x.filtq[1,],rev(x.filtq[3,])),col="pink",border=NA) lines(0:n.t,x,type="l",lwd=3) lines(0:n.t,x.filtq[2,],col=2,lwd=2) ## Checking the filter quantiles sum(x