tests/ou2-mif2.R

In pomp: Statistical Inference for Partially Observed Markov Processes

library(pomp)

pompExample(ou2)

pdf(file="ou2-mif2.pdf")

p.truth <- coef(ou2)
guess2 <- guess1 <- p.truth
guess1[c('x1.0','x2.0','alpha.2','alpha.3')] <- 0.9*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
guess2[c('x1.0','x2.0','alpha.2','alpha.3')] <- 1.2*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]

set.seed(64857673L)
mif1a <- mif(ou2,Nmif=100,start=guess1,
             pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
             rw.sd=c(
               x1.0=.5,x2.0=.5,
               alpha.2=0.1,alpha.3=0.1),
             transform=F,
             Np=1000,
             var.factor=1,
             ic.lag=10,
             cooling.type="hyperbolic",
             cooling.fraction=0.05,
             method="mif2",
             tol=1e-8
             )

mif2a <- mif(ou2,Nmif=100,start=guess1,
             pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
             rw.sd=c(
               x1.0=0.5,x2.0=.5,
               alpha.2=0.1,alpha.3=0.1),
             transform=F,
             Np=1000,
             var.factor=1,
             ic.lag=10,
             cooling.type="geometric",
             cooling.fraction=0.95^50,
             max.fail=100,
             method="mif",
             tol=1e-8
             )

plot(c(mif1a,mif2a))

set.seed(64857673L)
mif1b <- mif(ou2,Nmif=50,start=guess1,
             pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
             rw.sd=c(
               x1.0=.5,x2.0=.5,
               alpha.2=0.1,alpha.3=0.1),
             transform=F,
             Np=1000,
             var.factor=1,
             ic.lag=10,
             cooling.type="hyperbolic",
             cooling.fraction=0.05,
             method="mif2"
             )
mif1b <- continue(mif1b,Nmif=50)

mif2b <- mif(ou2,Nmif=50,start=guess1,
             pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
             rw.sd=c(
               x1.0=0.5,x2.0=.5,
               alpha.2=0.1,alpha.3=0.1),
             transform=F,
             Np=1000,
             var.factor=1,
             ic.lag=10,
             cooling.whatsit=200,
             cooling.type="geometric",
             cooling.fraction=0.95^50,
             max.fail=100,
             method="mif"
             )
mif2b <- continue(mif2b,Nmif=50)

mif2c <- mif(ou2,Nmif=50,start=guess1,
             pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
             rw.sd=c(
               x1.0=0.5,x2.0=.5,
               alpha.2=0.1,alpha.3=0.1),
             transform=F,
             Np=1000,
             var.factor=1,
             cooling.type="hyperbolic",
             cooling.fraction=0.05,
             max.fail=100,
             method="mif2"
             )
mif2c <- continue(mif2c,Nmif=50)

plot(c(mif1b,mif2b))

plot(c(mif1a,mif1b))
plot(c(mif2a,mif2b))

plot(mfl1 <- c(mif1b,mif2c))

mfl2 <- c(mfl1,mif2c)
mfl3 <- c(mif2a,mfl1)

try(c(mif2a,continue(mif2b,Nmif=1)))

guess2 <- guess1 <- coef(ou2)
guess1[c('x1.0','x2.0','alpha.2','alpha.3')] <- 0.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
guess2[c('x1.0','x2.0','alpha.2','alpha.3')] <- 1.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]

m1 <- mif2(ou2,Nmif=100,start=guess1,Np=1000,
           cooling.type="hyperbolic",cooling.fraction.50=0.05,
           rw.sd=rw.sd(
             x1.0=ivp(0.5),x2.0=ivp(0.5),
             alpha.2=0.1,alpha.3=0.1))

m2 <- mif2(ou2,Nmif=50,start=guess2,Np=1000,
           cooling.type="hyperbolic",cooling.fraction.50=0.05,
           rw.sd=rw.sd(
             x1.0=ivp(0.5),x2.0=ivp(0.5),
             alpha.2=0.1,alpha.3=0.1))
m2 <- continue(m2,Nmif=50)

plot(c(m1,m2))

rbind(mle1=c(coef(m1),loglik=logLik(pfilter(m1,Np=1000))),
      mle2=c(coef(m2),loglik=logLik(pfilter(m1,Np=1000))),
      truth=c(coef(ou2),loglik=logLik(pfilter(m1,Np=1000))))

m3 <- mif2(ou2,Nmif=3,start=guess1,Np=200,
           cooling.fraction=0.2,
           rw.sd=rw.sd(
             x1.0=c(0.5,rep(0.2,99)),
             x2.0=ivp(0.5),
             alpha.2=ifelse(time==1,0.2,0.1),
             alpha.3=0.2*(time<10)))

m4 <- mif2(ou2,Nmif=3,start=guess1,
           Np=function(k)if(k<20) 200 else 100,
           cooling.fraction=0.2,
           rw.sd=rw.sd(
             x1.0=c(0.5,rep(0.2,99)),
             x2.0=ivp(0.5),
             alpha.2=ifelse(time==1,0.2,0.1),
             alpha.3=0.2*(time<10)))

m4 <- mif2(m4,Nmif=2,Np=c(rep(200,20),rep(100,80),200))
m4 <- continue(m4,Nmif=2,cooling.fraction=0.1)
pf <- pfilter(m4)
m4 <- mif2(pf,Nmif=2,
           cooling.fraction=0.2,
           rw.sd=rw.sd(
             x1.0=c(0.5,rep(0.2,99)),
             x2.0=ivp(0.5),
             alpha.2=ifelse(time==1,0.2,0.1),
             alpha.3=0.2*(time<10)))

dev.off()