tests/process.R

In kidusasfaw/pomp: Statistical Inference for Partially Observed Markov Processes

options(digits=3)

library(pomp2)
library(magrittr)

ou2() %>% window(end=10) -> po

set.seed(293982095)

po %>%
  simulate(format="arrays",nsim=3,seed=3434388L) %>%
  extract2("states") -> x
po %>% time() -> t
coef(po) %>% parmat(7) -> p
p["sigma_1",] <- seq(from=1,to=7,by=1)

try(dprocess("ou2",x=x,times=t,params=p))
try(dprocess(x=x,times=t,params=p))
try(po %>% dprocess(x=x,times=t,params=p,log=TRUE))
try(po %>% dprocess(x=x,times=t,params=p[,1:2],log=TRUE))
try(po %>% dprocess(x=x,times=t[1:5],params=p,log=TRUE))
po %>% dprocess(x=x,times=t,params=p[,1:3],log=TRUE) -> d1
po %>% dprocess(x=x,times=t,params=p[,2],log=TRUE) -> d2
stopifnot(d1[2,]==d2[2,])
try(po %>% dprocess(x=x[,,2],times=t[2],params=p[,1:3],log=FALSE))
try(po %>% dprocess(x=x[,,2:5],times=t[2:5],params=p[,1:2],log=FALSE))
po %>% dprocess(x=x[,,2:5],times=t[2:5],params=p[,1:3],log=TRUE) %>%
  apply(1,sum)
try(po %>% dprocess(x=x[,1:2,2:5],times=t[2:5],params=p[,1:3],log=TRUE) %>%
    apply(1,sum))
po %>% dprocess(x=x[,1,2:5],times=t[2:5],params=p[,1:3],log=TRUE) %>%
  apply(1,sum)
po %>% pomp(dprocess=NULL) %>%
  dprocess(x=x[,1,2:5],times=t[2:5],params=p[,1:3],log=TRUE) %>% is.na() %>%
  stopifnot()

po %>% rinit(params=coef(po)) -> x0
freeze(po %>%
    rprocess(params=coef(po),xstart=parmat(x0,3),times=time(po,t0=TRUE),
      offset=1),
  seed=3434388L) -> x1

stopifnot(max(abs(x-x1))==0)

po %>% rinit(nsim=6) -> x0
try(rprocess("ou2",xstart=x0,times=t,params=p))
try(rprocess(xstart=x0,times=t,params=p))
try(po %>% rprocess(times=t,params=p))
try(po %>% rprocess(xstart=x0,params=p))
try(po %>% rprocess(xstart=x0,times=t))
try(po %>% rprocess(xstart=x0,times=t,params=p))
try(po %>% rprocess(xstart=x0,times=t,params=p[,-1],offset=-3))
try(po %>% rprocess(xstart=x0,times=t,params=p[,-1],offset=500))
try(po %>% rprocess(xstart=x0,times=t,params=p[,-1],offset=NA))
try(po %>% rprocess(xstart=x0,times=t,params=p[,-1],offset=Inf))
po %>% rprocess(xstart=x0,times=t,params=p[,1:3]) -> x
stopifnot(
  dim(x)==c(2,6,10),
  names(dimnames(x))==c("variable","rep","time")
)
po %>% rprocess(xstart=x0[,2],times=t,params=p[,1:3]) -> x
stopifnot(
  dim(x)==c(2,3,10),
  names(dimnames(x))==c("variable","rep","time")
)

try(po %>% rprocess(xstart=x0,times=t[2],params=p))
try(po %>% rprocess(xstart=x0[,2],times=t[2],params=p[,1:3]))
try(po %>% rprocess(xstart=x0[,2:4],times=t[2:5],params=p[,1:2]))
po %>% rprocess(xstart=x0[,2:4],times=t[2:5],params=p[,1:3]) %>%
  apply(1,sum)

simulate(
  times=seq(0,10), t0=0,
  params=c(s=3,x_0=0,tau=1),
  covar=covariate_table(z=c(1,1),times=c(0,10)),
  rprocess = onestep(
    function (t, x, s, delta.t, ...) {
      c(x=rnorm(n=1,mean=x,sd=s*sqrt(delta.t)))
    }
  ),
  dprocess = function (x_1, x_2, s, t_1, t_2, ...) {
    delta.t <- t_2-t_1
    dnorm(x=x_2,mean=x_1,sd=s*sqrt(delta.t),log=TRUE)
  },
  seed=3434388L
) -> rw

dprocess(rw,x=states(rw),params=coef(rw),times=time(rw),log=TRUE) -> d
stopifnot(
  round(sum(d),1)==-23.2,
  dim(d)==c(1,10),
  names(dimnames(d))==c("rep","time")
)