Skip to content

Commit

Permalink
more coverage improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
kingaa committed Jun 20, 2018
1 parent d09822d commit 630a8b9
Show file tree
Hide file tree
Showing 12 changed files with 168 additions and 74 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
Version: 1.16.3.1
Date: 2018-06-19
Version: 1.16.3.2
Date: 2018-06-20
Authors@R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="[email protected]"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
Expand Down
68 changes: 24 additions & 44 deletions R/mif.R
Original file line number Diff line number Diff line change
Expand Up @@ -323,44 +323,33 @@ mif.internal <- function (object, Nmif,

transform <- as.logical(transform)

if (length(start)==0)
stop(
ep,sQuote("start"),
" must be specified if ",
sQuote("coef(object)")," is NULL",
call.=FALSE
)
if (length(start)==0) stop(ep,sQuote("start")," must be specified if ",
sQuote("coef(object)")," is NULL",call.=FALSE)

if (transform)
start <- partrans(object,start,dir="toEstimationScale")

start.names <- names(start)
if (is.null(start.names))
stop(ep,sQuote("start"),
" must be a named vector",call.=FALSE)
if (is.null(start.names) || !is.numeric(start))
stop(ep,sQuote("start")," must be a named numeric vector",call.=FALSE)

rw.names <- names(rw.sd)
if (is.null(rw.names) || any(rw.sd<0))
stop(ep,sQuote("rw.sd"),
" must be a named non-negative numerical vector",call.=FALSE)
stop(ep,sQuote("rw.sd")," must be a named non-negative numerical vector",call.=FALSE)
if (!all(rw.names%in%start.names))
stop(ep,"all the names of ",sQuote("rw.sd"),
" must be names of ",sQuote("start"),call.=FALSE)
stop(ep,"all the names of ",sQuote("rw.sd")," must be names of ",sQuote("start"),call.=FALSE)
rw.names <- rw.names[rw.sd>0]
rw.sd <- rw.sd[rw.sd>0]
if (length(rw.names) == 0)
stop(ep,sQuote("rw.sd"),
" must have one positive entry for each parameter to be estimated",call.=FALSE)
stop(ep,sQuote("rw.sd")," must have one positive entry for each parameter to be estimated",call.=FALSE)

pars <- rw.names[!(rw.names%in%ivps)]

if (!is.character(ivps) || !all(ivps%in%start.names))
stop(ep,sQuote("ivps"),
" must name model parameters",call.=FALSE)
stop(ep,sQuote("ivps")," must name model parameters",call.=FALSE)

ntimes <- length(time(object))
if (is.null(Np)) stop(ep,sQuote("Np"),
" must be specified",call.=FALSE)
if (is.null(Np)) stop(ep,sQuote("Np")," must be specified",call.=FALSE)
if (is.function(Np)) {
Np <- tryCatch(
vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
Expand All @@ -377,15 +366,12 @@ mif.internal <- function (object, Nmif,
if (any(Np<=0))
stop(ep,"number of particles, ",sQuote("Np"),", must always be positive",call.=FALSE)
if (!is.numeric(Np))
stop(ep,sQuote("Np"),
" must be a number, a vector of numbers, or a function",
call.=FALSE)
stop(ep,sQuote("Np")," must be a number, a vector of numbers, or a function",call.=FALSE)
Np <- as.integer(Np)

ic.lag <- as.integer(ic.lag)
if ((length(ic.lag)!=1)||(ic.lag<1))
stop(ep,sQuote("ic.lag"),
" must be a positive integer",call.=FALSE)
if ((length(ic.lag)!=1)||(!is.finite(ic.lag))||(ic.lag<1))
stop(ep,sQuote("ic.lag")," must be a positive integer",call.=FALSE)
if (ic.lag>ntimes) {
warning(
ep,sQuote("ic.lag")," = ",ic.lag," > ",ntimes,
Expand All @@ -406,25 +392,24 @@ mif.internal <- function (object, Nmif,
}

if (missing(cooling.fraction.50))
stop(ep,sQuote("cooling.fraction.50"),
" must be specified",call.=FALSE)
stop(ep,sQuote("cooling.fraction.50")," must be specified",call.=FALSE)
cooling.fraction.50 <- as.numeric(cooling.fraction.50)
if ((length(cooling.fraction.50)!=1)||(cooling.fraction.50<0)||(cooling.fraction.50>1))
stop(ep,sQuote("cooling.fraction.50"),
" must be a number between 0 and 1",call.=FALSE)
if ((length(cooling.fraction.50)!=1)||(!is.finite(cooling.fraction.50))||
(cooling.fraction.50<0)||(cooling.fraction.50>1))
stop(ep,sQuote("cooling.fraction.50")," must be a number between 0 and 1",call.=FALSE)

cooling <- mif.cooling(
type=cooling.type,
fraction=cooling.fraction.50,
ntimes=ntimes
)

if ((length(var.factor)!=1)||(var.factor < 0))
var.factor <- as.numeric(var.factor)
if ((length(var.factor)!=1)||(!is.finite(var.factor))||(var.factor < 0))
stop(ep,sQuote("var.factor")," must be a positive number",call.=FALSE)

Nmif <- as.integer(Nmif)
if (Nmif<0)
stop(ep,sQuote("Nmif")," must be a positive integer",call.=FALSE)
if (Nmif<0) stop(ep,sQuote("Nmif")," must be a positive integer",call.=FALSE)

theta <- start

Expand Down Expand Up @@ -509,7 +494,7 @@ mif.internal <- function (object, Nmif,
fp={ # fixed-point iteration
theta[pars] <- pfp@filter.mean[pars,ntimes,drop=FALSE]
},
stop(ep,"unrecognized method ",sQuote(method),call.=FALSE)
stop(ep,"unrecognized method ",sQuote(method),call.=FALSE) # nocov
)
theta[ivps] <- pfp@filter.mean[ivps,ic.lag]
conv.rec[n+1,-c(1,2)] <- theta
Expand Down Expand Up @@ -564,22 +549,17 @@ setMethod(
method <- match.arg(method)

if (missing(start)) start <- coef(object)
if (missing(rw.sd))
stop(ep,sQuote("rw.sd")," must be specified.",call.=FALSE)
if (!is.numeric(rw.sd))
stop(ep,sQuote("rw.sd")," must be a named numeric vector.",call.=FALSE)
if (missing(rw.sd)) stop(ep,sQuote("rw.sd")," must be specified.",call.=FALSE)
if (!is.numeric(rw.sd)) stop(ep,sQuote("rw.sd")," must be a named numeric vector.",call.=FALSE)
if (missing(ic.lag)) {
if (length(ivps)>0) {
stop(ep,sQuote("ic.lag"),
" must be specified if ",sQuote("ivps"),
" are.",call.=FALSE)
stop(ep,sQuote("ic.lag")," must be specified if ",sQuote("ivps")," are.",call.=FALSE)
} else {
ic.lag <- length(time(object))
}
}

if (missing(Np))
stop(ep,sQuote("Np")," must be specified",call.=FALSE)
if (missing(Np)) stop(ep,sQuote("Np")," must be specified",call.=FALSE)

cooling.type <- match.arg(cooling.type)

Expand Down
3 changes: 1 addition & 2 deletions R/mif2.R
Original file line number Diff line number Diff line change
Expand Up @@ -261,8 +261,7 @@ mif2.internal <- function (object, Nmif, start, Np, rw.sd, transform = FALSE,
gnsi <- as.logical(.getnativesymbolinfo)
Np <- c(Np,Np[1L])

if (Nmif <= 0)
stop(ep,sQuote("Nmif")," must be a positive integer",call.=FALSE)
if (Nmif <= 0) stop(ep,sQuote("Nmif")," must be a positive integer",call.=FALSE)

cooling.fn <- mif2.cooling(
type=cooling.type,
Expand Down
Binary file modified tests/mif-01.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified tests/mif-02.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified tests/mif-03.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified tests/mif-04.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
34 changes: 34 additions & 0 deletions tests/mif.R
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,40 @@ guess2 <- guess1 <- p.truth
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')]

po <- ou2
coef(po) <- numeric(0)
try(mif(po))
try(mif(po,rw.sd=c(3,2)))
try(mif(po,Np=10,rw.sd=c(3,2)))
try(mif(po,Np=10,start=numeric(0),rw.sd=c(3,2)))
try(mif(po,Np=10,start="yes",rw.sd=c(3,2)))
try(mif(po,Np=10,start=c(a="yes"),rw.sd=c(3,2)))
try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3,b=2)))
try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3)))
try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3),cooling.fraction.50=3))
try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3),cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=-3),cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(1),cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1,b.0=2),ivps="c",cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1,b.0=2),ivps="c",ic.lag=10,cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1,b.0=2),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=0),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=1),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
try(mif(po,Np="10",start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
try(mif(po,Np=NULL,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=-10,cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=NA,cooling.fraction.50=1))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=NA))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=NULL))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=1,
var.factor=-100))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=1,
var.factor=NA))
try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=1,
var.factor=NULL))
try(mif(ou2,Np=c(10,20,30),rw.sd=c(alpha.1=1),ivps="x1.0",ic.lag=5,cooling.fraction.50=1))

mif1 <- mif(ou2,Nmif=25,start=guess1,
ivps=c('x1.0','x2.0'),
rw.sd=c(
Expand Down
64 changes: 63 additions & 1 deletion tests/mif.Rout.save
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,68 @@ particle filter log likelihood at truth
> 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')]
>
> po <- ou2
> coef(po) <- numeric(0)
> try(mif(po))
Error : in 'mif': 'rw.sd' must be specified.
> try(mif(po,rw.sd=c(3,2)))
Error : in 'mif': 'Np' must be specified
> try(mif(po,Np=10,rw.sd=c(3,2)))
Error : in 'mif': 'start' must be specified if 'coef(object)' is NULL
> try(mif(po,Np=10,start=numeric(0),rw.sd=c(3,2)))
Error : in 'mif': 'start' must be specified if 'coef(object)' is NULL
> try(mif(po,Np=10,start="yes",rw.sd=c(3,2)))
Error : in 'mif': 'start' must be a named numeric vector
> try(mif(po,Np=10,start=c(a="yes"),rw.sd=c(3,2)))
Error : in 'mif': 'start' must be a named numeric vector
> try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3,b=2)))
Error : in 'mif': all the names of 'rw.sd' must be names of 'start'
> try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3)))
Error : in 'mif': 'cooling.fraction.50' must be specified
> try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3),cooling.fraction.50=3))
Error : in 'mif': 'cooling.fraction.50' must be a number between 0 and 1
> try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=3),cooling.fraction.50=1))
Error : in 'mif': in default 'initializer': there are no parameters with suffix '.0'. See '?pomp'.
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=-3),cooling.fraction.50=1))
Error : in 'mif': 'rw.sd' must be a named non-negative numerical vector
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),cooling.fraction.50=1))
Error : in 'mif': in 'mif.pfilter': process simulation error: in 'discrete.time.sim' plugin: variable 'x1' not found among the state variables
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(1),cooling.fraction.50=1))
Error : in 'mif': 'rw.sd' must be a named non-negative numerical vector
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1,b.0=2),ivps="c",cooling.fraction.50=1))
Error : in 'mif': 'ic.lag' must be specified if 'ivps' are.
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1,b.0=2),ivps="c",ic.lag=10,cooling.fraction.50=1))
Error : in 'mif': 'ivps' must name model parameters
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1,b.0=2),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
Error : in 'mif': in 'mif.pfilter': process simulation error: in 'discrete.time.sim' plugin: variable 'x1' not found among the state variables
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=0),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
Error : in 'mif': 'rw.sd' must have one positive entry for each parameter to be estimated
> try(mif(po,Np=10,start=c(a=3),rw.sd=c(a=1),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
Error : in 'mif': 'ivps' must name model parameters
> try(mif(po,Np="10",start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
Error : in 'mif': 'Np' must be a number, a vector of numbers, or a function
> try(mif(po,Np=NULL,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=10,cooling.fraction.50=1))
Error : in 'mif': 'Np' must be specified
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=-10,cooling.fraction.50=1))
Error : in 'mif': 'ic.lag' must be a positive integer
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=NA,cooling.fraction.50=1))
Error : in 'mif': 'ic.lag' must be a positive integer
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=NA))
Error : in 'mif': 'cooling.fraction.50' must be a number between 0 and 1
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=NULL))
Error : in 'mif': 'cooling.fraction.50' must be a number between 0 and 1
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=1,
+ var.factor=-100))
Error : in 'mif': 'var.factor' must be a positive number
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=1,
+ var.factor=NA))
Error : in 'mif': 'var.factor' must be a positive number
> try(mif(po,Np=10,start=c(a=3,b.0=1),rw.sd=c(a=1),ivps="b.0",ic.lag=5,cooling.fraction.50=1,
+ var.factor=NULL))
Error : in 'mif': 'var.factor' must be a positive number
> try(mif(ou2,Np=c(10,20,30),rw.sd=c(alpha.1=1),ivps="x1.0",ic.lag=5,cooling.fraction.50=1))
Error : in 'mif': 'Np' must have length 1 or length 101
>
> mif1 <- mif(ou2,Nmif=25,start=guess1,
+ ivps=c('x1.0','x2.0'),
+ rw.sd=c(
Expand Down Expand Up @@ -101,7 +163,7 @@ in 'plot-mif': 'y' is ignored
> mif1a <- c(mif1)
> coef(mif2)
alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau x1.0 x2.0
0.800 -0.481 0.343 0.900 3.000 -0.500 2.000 1.000 -3.294 2.481
0.800 -0.508 0.344 0.900 3.000 -0.500 2.000 1.000 -3.462 2.563
> dim(coef(mif12))
[1] 2 10
> dim(coef(c(mif12,mif2)))
Expand Down
10 changes: 5 additions & 5 deletions tests/sir.Rout.save
Original file line number Diff line number Diff line change
Expand Up @@ -193,8 +193,8 @@ process model simulator, rprocess = function (xstart, times, params, ..., zerona
stop(ep, conditionMessage(e), call. = FALSE)
})
}
<bytecode: 0x38b87e0>
<environment: 0x2ac9c20>
<bytecode: 0x3449590>
<environment: 0x2659230>
process model density, dprocess = function (x, times, params, ..., tcovar, covar, log = FALSE,
.getnativesymbolinfo = TRUE)
{
Expand All @@ -204,8 +204,8 @@ process model density, dprocess = function (x, times, params, ..., tcovar, covar
stop(ep, conditionMessage(e), call. = FALSE)
})
}
<bytecode: 0x35b3540>
<environment: 0x276b100>
<bytecode: 0x31442f0>
<environment: 0x2305858>
measurement model simulator, rmeasure = native function '__pomp_rmeasure', defined by a Csnippet
measurement model density, dmeasure = function (y, x, t, params, log, covars, ...)
{
Expand Down Expand Up @@ -236,7 +236,7 @@ skeleton (vectorfield) = function (x, t, params, covars, ...)
xdot
})
}
<environment: 0x3d320b8>
<environment: 0x38c4590>
initializer = function (params, t0, ...)
{
with(as.list(params), {
Expand Down
7 changes: 7 additions & 0 deletions tests/verbosity.R
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,14 @@ capture.output(simulate(pomp(ricker,rmeasure=Csnippet("y=rpois(N);"),statenames=
cfile="bob",verbose=TRUE),verbose=TRUE) -> po) -> out
gsub("(\\w+)\\s.+","\\1",out,perl=TRUE)

set.seed(3949586L)
capture.output(invisible(mif2(window(ricker,end=10),Nmif=1,Np=1,rw.sd=rw.sd(r=1),
transform=TRUE,cooling.fraction.50=1,verbose=TRUE)),
type="message") -> out
stopifnot(sum(grepl("filtering failure at time",out))==5)

set.seed(3949586L)
capture.output(invisible(mif(window(ricker,end=10),Nmif=1,Np=1,rw.sd=c(r=1),
transform=TRUE,cooling.fraction.50=1,verbose=TRUE)),
type="message") -> out
stopifnot(sum(grepl("filtering failure at time",out))==5)
Loading

0 comments on commit 630a8b9

Please sign in to comment.