tests/ou2-bsmc2.R

In pomp: Statistical Inference for Partially Observed Markov Processes

library(pomp)

set.seed(398585L)
pompExample(ou2)

time(ou2) <- 1:10

Np <- 50000

prior.bounds <- rbind(
                      alpha.2=c(-0.55,-0.45),
                      alpha.3=c(0.25,0.35)
                      )
colnames(prior.bounds) <- c("lower","upper")

estnames <- rownames(prior.bounds)

prior <- matrix(data=coef(ou2),nrow=length(coef(ou2)),ncol=Np)
rownames(prior) <- names(coef(ou2))
for (n in estnames) {
  prior[n,] <- runif(n=Np,min=prior.bounds[n,1],max=prior.bounds[n,2])
}

##Run Liu & West particle filter
tic <- Sys.time()
smc <- bsmc2(
             ou2,
             params=prior,
             est=estnames,
             smooth=0.02
             )
toc <- Sys.time()

prior <- smc$prior
post <- smc$post

print(etime <- toc-tic)

print(
      cbind(
            prior.mean=apply(prior,1,mean),
            posterior.mean=apply(post,1,mean),
            truth=coef(ou2),
            t(apply(post,1,quantile,c(0.025,0.5,0.975)))
            )
      )

print(min(smc$eff.sample.size))
print(smc$log.evidence)

ou2 <- pomp(ou2,
            rprior=function(params,...){
              params
            }
            )

smc <- bsmc2(ou2,Np=25000,smooth=0.1,est=estnames,seed=648651945L)
print(smc$eff.sample.size)
print(smc$log.evidence)