seqMC: seqMC creates a seqMC object

In chuanwen/seqMC: Sequential Monte Carlo

Description

seqMC creates a seqMC object

Usage

seqMC(f, logprob_y_given_x, x0, y0, sample_method = c("systematic",
  "residual", "bootstrap"))

Arguments

f	function, when called with parameter t (time point) and x_t (state vector at time t), it would return x_(t+1)
logprob_y_given_x	function, when called with parameter t (time point), y_t (observation vector at time t) and x_t (state vector at time t), it would return the conditional log_probability: log(Prob(y_t | x_t ))
x0	matrix, sample of state vector at time 0, each col is a sample of state at time 0.
y0	observation at time 0 (can be missing).
sample_method	character, specify sample method in the resample stage. Default systematic, means "systematic resampling".

Value

a seqMC object, which can be updated at each time point. e.g. obj = update(ojb, y)

Examples

f <- function(t, x) {
 0.5 * x + 25 * x / (1 + x * x) + 8.0 * cos(1.2 * (t-1)) + rnorm(length(x), sd=sqrt(10.0))
}

logprob_y_given_x <- function(t, y, x) {
	 as.numeric(-(y - x * x / 20.0)**2/2.0)
}

x0 = matrix(rnorm(4000, sd=2), nrow=1, ncol=4000)

mod = seqMC(f, logprob_y_given_x, x0)

### simulate true path ###
x0 = 0.1
T = 50
x = rep(0.0, T)
x[1] = f(0, x0)
for (t in 1:(T-1)) {
	 x[t+1] = f(t, x[t])
}
y = x * x / 20 + rnorm(length(x))

### estimate the posterior of state vector given y[t] ####
xhat = sapply(1:T, function(t) {
   mod <<- update(mod, y[t])
   estimate.seqMC(mod)
})

plot(x, ylim=c(-40, 40), pch='*')
lines(xhat[1,])
lines(xhat[2,], lty='dotted')
lines(xhat[3,], lty='dotted')

chuanwen/seqMC documentation built on May 4, 2019, 3:19 p.m.