gfortran
when I try to install pomp from source on my Mac?pomp
for multivariate data?dmeasure
returning illegal values. What does it mean?pomp
?c
operator to concatenate certain kinds of pomp structures. How do I fix this?You are always welcome to ask for help, but maintainer time is limited. Therefore, before submitting a request, first see if you can solve your own problem. The following steps will help you to do so.
If you find that your question has not been answered, or you believe you have found a bug, then submit a request for help via the Issues page if possible or to Aaron King (kingaa at umich dot edu) if necessary. Include in your request—at a minimum—the code that you ran that produced the error, together with a transcript of the R session that gave the errors. Better still, construct a minimal example that will reproduce the error: this allows for the most efficient solution of problems.
Make sure to include the output of following, so that I can see what version of pomp you are using, what version of R, what kind of machine you have, any customizations you have that might be affecting performance, etc.
source("https://kingaa.github.io/scripts/diagnostics.R")
To cite the ‘pomp’ package in publications, use:
A. A. King, D. Nguyen, E. L. Ionides (2016). Statistical Inference for Partially Observed Markov Processes via the R Package pomp. Journal of Statistical Software 69(12): 1–43. DOI:10.18637/jss.v069.i12.
Additionally, consider citing the package itself:
A. A. King, E. L. Ionides, C. M. Breto, S. P. Ellner, M. J. Ferrari, S. Funk, S. G. Johnson, B. E. Kendall, M. Lavine, D. Nguyen, E. B. O’Dea, D. C. Reuman, H. Wearing, and S. N. Wood (2024). pomp: Statistical Inference for Partially Observed Markov Processes. (R package, version 5.11.0.0). https://kingaa.github.io/pomp/
You can get this information by executing citation("pomp")
in an R session. LaTeX users can get a BibTeX entry for these publications by executing toBibtex(citation("pomp"))
.
If you use pomp in a publication, please take a moment to let us know! Send an email with the citation, and we’ll include it in the bibliography.
There are many resources available that can help. First, check the package manual to see if it describes the situation you have encountered. Second, check the discussions forum. Many conversations about various situations that have arisen are archived there. Third, search the issues (including the closed issues) to see if your situation has been encountered and resolved already. Finally, feel free to submit an issue. In doing so, display the error message you get and document how it arises. Please follow the instructions in FAQ 1.1 above.
The pomp developers have a goal of replacing all unhelpful error messages with helpful ones. By letting us know when messages are confusing or uninformative, you help us to make the package more usable.
Of course, there could be many reasons for this. However, you should know that pomp version 2 made many changes to the pomp interface so that older code will not run under versions ≥ 2. See the pomp version 2 upgrade guide for information on how to fix broken code from earlier versions.
The paper announcing pomp was published in the Journal of Statistical Software (King, Nguyen, and Ionides 2016). An updated version of that paper can be found on the pomp website.
The video and visual aids from a talk entitled “Two perspectives on the pomp project”, which describes the goals and organization of the package from both the user’s and developer’s persectives is also linked from the website.
See the installation instructions.
gfortran
when I try to install pomp from source on my Mac?Compilation of pomp requires a FORTRAN compiler, which has not been installed (or not installed properly) on your system. See the installation instructions.
Firstly, as a general rule, it’s very probably a good idea to simply install the latest version of R, since each new release fixes bugs, adds and improves features, etc.
Usually, this warning arises because you are trying to install binaries from a repository on which binaries built for your version of R are no longer available. The preferred solution in this case is, again, to upgrade R, but if for some reason that is impossible or undesirable, you can always build the package from source using the source code retrieved from the package website or from one of the CRAN mirrors. For more, see the installation instructions.
To be clear: the statement that pomp will work for R(>=Y) is a guarantee that the code is compatible with those versions of R. It is not a promise that binaries will always be available for old versions of R.
A C snippet can make use of any feature of the C language. In particular, we can use pointers to give access to arrays. For example, consider the following, which implement a simple chain of random variables.
Csnippet("
rSim <- double *x = &X1;
int i, n = (int) N;
for (i = 0; i < n-1; i++) x[i] = x[i+1];
x[n-1] = rnorm(0,1);
")
Csnippet("
rInit <- double *x = &X1;
int i, n = (int) N;
for (i=0; i < n; i++) x[i] = 0.0;
")
Csnippet("
rMeas <- Y1 = rnorm(X2,2);
")
The following simulates and plots using the above.
library(ggplot2)
library(tidyr)
simulate(
t0=0,
times=1:20,
rprocess=euler(rSim,delta.t=1),
rmeasure=rMeas,
rinit=rInit,
obsnames="Y1",
statenames=sprintf("X%d",1:5),
paramnames="N",
params=c(N=5),
format="data.frame"
|>
) pivot_longer(-c(.id,time)) |>
ggplot(aes(x=time,y=value,group=name,color=name))+
geom_line()+
theme_bw()
pomp
for multivariate data?Yes, this is no problem. The data
you supply to pomp
can contain multiple variables. You simply refer to these variables by name in writing the rmeasure
and dmeasure
C snippets. The following example illustrates this for a multivariate Ornstein-Uhlenbeck process (as in pompExample("ou2")
) with a slightly peculiar measurement model.
simulate(
times=1:100, t0=0,
rmeasure=Csnippet("
y1 = rnorm(x1+x2,sigma1);
y2 = rnorm(x2-x1,sigma2);"
),dmeasure=Csnippet("
lik = dnorm(y1,x1+x2,sigma1,1)+dnorm(y2,x2-x1,sigma2,1);
lik = (give_log) ? lik : exp(lik);"
),rprocess=discrete_time(
step.fun=Csnippet("
double tx1, tx2;
tx1 = rnorm(a11*x1 + a12*x2,nu1);
tx2 = rnorm(a21*x1 + a22*x2,nu2);
x1 = tx1; x2 = tx2;"
),delta.t=1),
rinit=Csnippet("
x1 = 0;
x2 = 0;
"),
obsnames=c("y1","y2"),
statenames=c("x1","x2"),
paramnames=c("a11","a12","a21","a22","sigma1","sigma2","nu1","nu2"),
params=c(a11=0.5,a12=-0.1,a21=0.2,a22=-1,nu1=0.3,nu2=0.1,sigma1=0.1,sigma2=0.3)
obj ) ->
Missing data are a common complication. pomp can handle NAs in the data, but if it is needed, the measurement model probability density function, dmeasure
must be written so as to deal with NAs appropriately. For example, look at the following variant of the SIR model describing the influenza outbreak in a boarding school:
library(ggplot2)
read.csv(text="
B,day
NA,0
1,1
6,2
26,3
73,4
NA,5
293,6
258,7
236,8
191,9
124,10
69,11
26,12
11,13
4,14
") -> dat
|>
dat ggplot(aes(x=day, y=B)) +
geom_line() + geom_point()
Data are missing at days 0 and 5.
We create a pomp
object implementing a simple SIR process model and a binomial measurement model, as in the original example. The only difference is in the dmeasure
:
library(pomp)
Csnippet("
sir_step <- double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
double dN_IR = rbinom(I,1-exp(-gamma*dt));
S -= dN_SI;
I += dN_SI - dN_IR;
R += dN_IR;
H += dN_IR;
")
Csnippet("
sir_init <- S = N-1;
I = 1;
R = 0;
H = 0;
")
Csnippet("B = rbinom(H,rho);")
rmeas <-
Csnippet("
dmeas <- if (ISNA(B)) {
lik = (give_log) ? 0 : 1;
} else {
lik = dbinom(B,H,rho,give_log);
}
")
|>
dat pomp(
time="day",t0=0,
rprocess=euler(sir_step,delta.t=1/12),
rinit=sir_init,
rmeasure=rmeas,
dmeasure=dmeas,
accumvars="H",
paramnames=c("Beta","gamma","N", "rho"),
statenames=c("S","I","R","H")
sir_na ) ->
In the above, to check for missing data, we use the ISNA
macro, which is part of R’s C API. Note that the dmeasure
returns a likelihood of 1 (log likelihood 0) for the missing data. [What’s the probability of seeing something if you don’t look?] When there is an observation, it returns a binomial likelihood as usual.
Our simulations now fill in simulated values for the missing data,
|>
sir_na simulate(
params=c(Beta=3,gamma=2,rho=0.9,N=2600),
nsim=10,
format="data.frame",
include.data=TRUE
|>
) ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
geom_line()+
guides(color="none")+
theme_bw()
and the particle filter handles the missing data correctly:
|>
sir_na pfilter(Np=1000,params=c(Beta=3,gamma=2,rho=0.9,N=2600)) |>
as.data.frame() |>
pivot_longer(-day) |>
ggplot(aes(x=day,y=value,color=name))+
guides(color="none")+
geom_line()+
facet_wrap(~name,ncol=1,scales='free_y')+
theme_bw()
In the above particle filter computation, notice that the effective sample size (ESS) is full, as it should be, when the missing data contribute no information.
As described in the documentation (?pomp
), pomp allows you to define “accumulator variables” that collect events occurring between observations. The example above illustrates this.
In that example, we took t0
equal to the first observation time. Sometimes we want to take t0
to be much earlier. For example, if we wish to posit that the initial state of the unobserved Markov process is distributed at t0
according to its stationary distribution, one way to achieve this is to put t0
so early that simulations will equilibrate before the first observation is made.
This is straightforward, but the presence of accumulator variables leads to a small twist. Let’s look at the boarding-school flu example to illustrate.
library(pomp)
Csnippet("
sir_step <- double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
double dN_IR = rbinom(I,1-exp(-gamma*dt));
S -= dN_SI;
I += dN_SI - dN_IR;
R += dN_IR;
H += dN_IR;
")
Csnippet("
sir_init <- S = N-1;
I = 1;
R = 0;
H = 0;
")
Csnippet("B = rbinom(H,rho);")
rmeas <-
Csnippet("
dmeas <- if (ISNA(B)) {
lik = (give_log) ? 0 : 1;
} else {
lik = dbinom(B,H,rho,give_log);
}")
read.csv(text="
day,B
1,3
2,8
3,28
4,76
5,222
6,293
7,257
8,237
9,192
10,126
11,70
12,28
13,12
14,5
") -> dat
|>
dat pomp(
time="day",t0=-8,
rprocess=euler(sir_step,delta.t=1/12),
rinit=sir_init,
rmeasure=rmeas,
dmeasure=dmeas,
accumvars="H",
paramnames=c("Beta","gamma","N", "rho"),
statenames=c("S","I","R","H")
bsflu ) ->
Note that, as above, we’ve allowed for the possibility of NAs in the data.
Now let’s look at the data and some simulations from the model.
library(ggplot2)
|>
bsflu as.data.frame() |>
ggplot(aes(x=day,y=B))+
geom_point()+geom_line()+
theme_bw()
|>
bsflu simulate(
params=c(Beta=1.5,gamma=1,rho=0.9,N=2600),
nsim=10,format="d",include.data=TRUE
|>
) ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
geom_line()+
guides(color="none")+
theme_bw()
The data plot looks fine, but what’s going on with the simulations? It’s easy to understand: we are assuming that there is one infectious host at t = t0 = -8
days. In most of the simulations, this leads to a number of new infections, which the H
variable accumulates from t=-8
to t=1
, the time of our first observation. But our first observation is not the number of new cases to have occurred in the last 9 days, but the number that occurred (and were reported) in the last 1 day.
We can fix this by introducing a fictitious “observation” at t=0
, with missing data, i.e., by assuming we observed nothing at t=0
:
library(dplyr)
|>
dat bind_rows(c(day=0,B=NA)) |>
arrange(day) |>
pomp(
time="day",t0=-8,
rprocess=euler(sir_step,delta.t=1/12),
rinit=sir_init,
rmeasure=rmeas,
dmeasure=dmeas,
accumvars="H",
paramnames=c("Beta","gamma","N", "rho"),
statenames=c("S","I","R","H")
bsflu
) ->
|>
bsflu simulate(
params=c(Beta=1.5,gamma=1,rho=0.9,N=2600),
nsim=10,format="d",include.data=TRUE
|>
) filter(day>0) |>
ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
geom_line()+
guides(color="none")+
theme_bw()
The answer to this question depends on the nature of the data. The measurement model is as much a part of the model as the latent process model. Accordingly, the choice of the measurement model is governed by the same considerations: Which model fits the data best? Which model has the greater predictive power? Which model makes the most sense?
Nor is the choice of the measurement model without consequences for the choice of the latent process model. Despite the fact that one typically is much more interested in the latent process than in the measurement, if one does not pay careful attention to both, one risks making incorrect or sub-optimal inference.
The Euler multinomial approximation of a continuous-time, stochastic compartmental model is as follows. Suppose a compartment has occupancy Nt at time t and that there are K ways of exiting the compartment, with per capita rates (hazards) μ1,…,μK, respectively.
To make the Euler multinomial approximation, we approximate the total exit rates as constant over a small interval [t,t+Δt). Let the random variable Δnk, k=1,…,K, be the number that exit by path k in this time interval and Δn0 be the number that remain. Under this assumption, the vector of numbers of exits, (Δn0,Δn1,…,ΔnK) is multinomially distributed with size Nt and probabilities (pk)Kk=0, where p0=exp(−∑μiΔt), and pk=μk∑μi(1−p0),k=1,…,K. By way of shorthand, we say that Δn=(Δnk)Kk=1 is Euler-multinomially distributed with size Nt, rates μ=(μk)Kk=1, and time-step Δt and we write Δn∼Eulermultinom(Nt,μ,Δt).
When constructing a POMP for a compartmental model using the Euler multinomial approximation, to write the rprocess basic component, one uses the euler
plugin. One chooses the time step using the delta.t
argument. In writing the step.fun
C snippet, one needs exactly one call to reulermultinom
for each compartment of the model. Using the notation above, one has to pack the K rates μ1,…,μK into contiguous memory locations and retrieve the results in (a different set of) contiguous memory locations. For example, if rate
is a pointer to K contiguous memory locations holding the rates and dn
is a pointer to K contiguous memory locations ready to hold the results, then
reulermultinom(K,N,rate,dt,dn);
will result in a random sample from the Euler multinomial distribution (with timestep dt) being stored in dn[0]
, …, dn[K-1]
. In the foregoing, we’ve assumed that the quantities Nt and K are stored in the memory locations denoted N
and K
, respectively. Note that the rules for step.fun
C snippets guarantee that the operative time step (Δt) is stored in dt
and that this will not exceed the user-specified delta.t
. See ?rprocess_spec
and ?Csnippet
for more information.
When the step function is provided as an R function, the corresponding call is
reulermultinom(n=1,size=N,rate=mu,dt=delta.t) dn <-
where mu
is a vector containing the rates.
The manual page contains more information on the Euler-multinomial distribution.
Yes, this is straightforward. One can make use of the globals
and shlib.args
arguments to do so.
Any valid C code passed to globals
is included at the top level of the C file containing C snippets. For example, any variables or functions defined there will be visible to every C snippet in that file. If one wishes to call functions defined in a precompiled C library, one can include the library’s header file via the usual #include
statement, passed via the globals
argument.
Any text passed to the shlib.args
argument is included in the command-line call to R CMD SHLIB
. To link to a precompiled library, then, one can do something like
"<library>" shlib.args =
where <library>
is the path of the library.
Consider the following simple example. The function truncated_normal_a_pdf
is defined in the library libtn.so
, located in directory libs
. It has a header file, located in include/tn.h
.
simulate(
t0 = 0,
times = 0:100,
obsnames = "y",
rprocess=discrete_time(
step.fun=Csnippet("
double e = rnorm(0,sigma);
N = exp(log(r)+log(N)-N+e);"),
delta.t=1
),rmeasure = Csnippet("
y = rnorm(N,4);
y = (y > 0) ? y : 0;"
),dmeasure = Csnippet("
lik = truncated_normal_a_pdf(y,N,4,0);
if (give_log) lik = log(lik);"
),statenames = "N",
paramnames = c("r", "sigma"),
params = c(N_0 = 7, r = 45, sigma = 0.3),
globals = r"{#include "include/tn.h"}",
shlib.args = "libs/libtn.so"
rick
) ->
|>
rick pfilter(Np=1000) |>
logLik()
In Discussion #156 there is another simple example.
dmeasure
returning illegal values. What does it mean?A common error message is
: in ‘pfilter’: ‘dmeasure’ returns illegal value Error
This error arises because your dmeasure
function is giving a negative, infinite, NA
, or otherwise illegal value of the likelihood. Common causes are:
dmeasure
function or C snippet returning the log likelihood when the likelihood is requested, or vice versa.lik
has not been set properly. Note that, when the snippet is executed, lik
will be undefined. In particular, you cannot assume that it has value 0.See the help (?Csnippet
) for information on how to properly write dmeasure
C snippets. As of version 1.5.6.3, when this error occurs, the values of the likelihood, parameters, and states are printed. If this alone doesn’t make the underlying cause clear, explore the model behavior further. A useful way to do this is to put Rprintf
statements in your dmeasure
and rprocess
C snippets to view the values of the data, state variables, and likelihoods. Issue #13 has an example.
Remember, dmeasure
has to be able to return the likelihood—or the log likelihood—upon request!
If your system does not support R CMD SHLIB
, you won’t be able to use C snippets. R compiles native codes (C, C++, and FORTRAN) into shared-object libraries using SHLIB
. To check whether your R installation is configured to use SHLIB
, run the following code snippet in an R session:
cat(
"#include <R.h>
void hello (void) {
Rprintf(\"hello world!\\n\");
}",
file="hello.c")
system("R CMD SHLIB hello.c")
dyn.load(paste0("hello",.Platform$dynlib.ext))
.C("hello",PACKAGE="hello")
If it works, you should see
!
hello worldlist()
If it doesn’t work, consult the installation instructions.
pomp
?See the preceding FAQ.
Also, some Windows users have reported trouble with certain parallel codes. Specifically, pomp codes that require compilation of C snippets trigger errors when executed in parallel on certain Windows machines.
The workaround is to execute serially all codes that require compilation; that is, to execute them outside of all parallel code blocks.
Calls that induce compilation include:
rmeasure
, dmeasure
, rprocess
, dprocess
, skeleton
, rprior
, dprior
, partrans
, and rinit
. Whenever one of these arguments is furnished a Csnippet
, the snippet is compiled.Csnippet
in such cases, each such call results in a Csnippet
being written internally, and then compiled.In addition, the use of temporary directories to hold the C files and dynamically loadable libraries created from them that underlie the C snippet facility causes problems on certain Windows machines. We hypothesize that this is due to Windows securities features.
Again, there is a workaround. One can specify to pomp that all compilation should be done in the current working directory.
An illustration is given in the Simulation-based Inference for Epidemiological Dynamics short course.
c
operator to concatenate certain kinds of pomp structures. How do I fix this?This appears to arise from incompatibility between certain R installations and binary versions of the pomp package from CRAN. Try installing the latest version from the github site:
install.packages("pomp",repos="https://kingaa.github.io/")
If you continue to see the error, try installing from source. The easiest way to do this is to use devtools:
::install_github("kingaa/pomp") devtools
For more help with installation, see the installation instructions.
If these steps do not solve the problem, please open an issue. In doing so, please follow the instructions in FAQ 1.1 above.
pomp incorporates time-varying covariates via a lookup table. Whenever your model needs the value of a covariate at a particular time, it uses the lookup table provided (via the covar
and tcovar
options of pomp
), linearly interpolating if necessary. It can happen that the required computation depends on values of the covariates outside the time interval covered by the lookup table. When this happens, pomp does linear extrapolation, with a warning.
This can be a problem if the linear extrapolation is not accurate. To avoid this, be sure to give values for the covariates over a time interval that covers the full range of relevant times, i.e., t0
to the time of the last observation.
Note that the ODE integration algorithms underlying trajectory
can sometimes evaluate the covariates beyond the last observation time. To avoid warnings (and the potential inaccuracy they signal) in this case, it may be necessary to supply values for the covariates somewhat beyond the final observation time.
See Issue 115 for an example. This can happen if you use file or directory names that contain whitespace. For the tools that compile the C snippets—as with (all?) other command-line tools—whitespace is syntactic. By far the best solution to this problem is to avoid using whitespace altogether in file and directory names! The standard choice is to substitute the underscore "_" for whitespace.
Error messages containing the following elements often arise.
...: error in building shared-object library from C snippets:
Errorin ‘Cbuilder’: compilation error: cannot compile shared-object library
...:
compiler messages
...: "beta" redefined
warning ...
This is due to a name conflict between a quantity beta
defined as a model parameter, covariate, data variable, or in some other way defined in a C snippet, and the Euler beta function, which is defined as beta
in Rmath.h
and documented in the R extensions manual. To work around this, simply rename the parameter. For example, it could be called Beta
instead.
The error message displays the “prototype” of the function pomp expects. For example, in the message
: in 'sannbox': 'candidate.dist' must be a function of the form 'candidate.dist(par, temp, scale, ...)'. Error
the prototype expected is candidate.dist(par, temp, scale, ...)
. This means that a function supplied to the candidate.dist
argument must have at least the arguments par
, temp
, scale
, and ...
. Note that the ellipses (...
) are necessary. Similarly, the error message
: in ‘rmeasure’: ‘rmeasure’ must be a function of the form ‘rmeasure(...)’ Error
tells us that we have left the ...
out of the argument list in the function we supplied for rmeasure
.
Error messages like the following
: in 'pmcmc': in 'dprior': variable 'r' not found among the parameters. Error
or
: in 'dmeasure': variable 'y2' not found among the observables. Error
or similar messages, wherein a named variable is not found among the parameters, state variables, covariates, or observables, arise when a C snippet-based workhorse function, i.e., one that executes a basic model component, cannot find the named variable among the supplied parameters, state variables, covariates, or observables. Typically this means that either (a) the full vector of model parameters has been supplied neither via the params
argument nor in the pomp object’s coef
slot; (b) a basic model component that computes state variables (e.g., rinit
or rprocess
) is not returning the full vector of state variables; (c) the covariate table does not contain the full set of covariates; or (d) some data variable is missing.
It may seem puzzling when, for example, the rinit
workhorse complains about the lack of a parameter upon which the rinit basic component does not in fact depend. To understand this, it is necessary to know that, whenever one or more C snippets are passed to a pomp function, the names of the parameters, state variables, covariates, and observables that are needed by any of them will be expected by all of them. Thus, when rinit
complains that it cannot find a certain parameter that it actually doesn’t need to know, it is effectively doing so on behalf of another basic model component which will depend on that parameter.
Finally, note that, while one informs pomp functions about the names of parameters using the paramnames
argument and about the names of state variables using statenames
, the names of covariates are by default drawn from the covariate table itself, and the names of the observables are taken from the data. This behavior can be over-ridden using the optional covarnames
and obsnames
arguments.
mif2
gives me the same as the likelihood I get from pfilter
?No. As its name suggests, the iterated filtering (IF) algorithm (that implemented in pomp as mif2
) work by iteratively performing particle filtering on a model, but that model is slightly different from the focal model. In particular, IF pertubs the parameters of the model by adding random noise to them. Actually, if one were to simulate from this model, the parameters would be taking random walks. This effectively smooths the likelihood surface, combats particle depletion, and—critically for IF—allows local exploration of the parameter space. The likelihood that mif2
reports is the estimated likelihood of this perturbed model. It can be either higher or lower than that of the focal model. Iterated filtering works by gradually reducing the magnitude of the perturbations with successive iterations (“cooling”), so that the perturbed model gradually becomes more similar to the focal model.
To compute the likelihood at a point returned by mif2
, it is necessary to re-run the particle filter. By doing the latter several times, one can obtain a more accurate estimate of the likelihood along with a measure of the associated uncertainty.
mif2
iterations?In Issue 27, @theresasophia asked this question. She provided codes for an estimation procedure based on mif2
. The plot below illustrates what she found.
The plots above show that, after a few iterations during which it rapidly increases, the mif2
likelihood displays a slow downward trend. However, the mif2
likelihood is not identical to the true likelihood (see FAQ 5.1). This is because mif2
adds extra variability to the model in the form of random perturbations to the parameters. In effect, the model mif2
sees is different from the model of interest. Indeed, mif2
works by gradually reducing the intensity of this extra variability, allowing one model to converge to the other. We refer to this gradual reduction as “cooling”. In the following, “MOI” refers to the model of interest and “MM” to the corresponding model that mif2
sees. You can estimate the true likelihood of the MOI using pfilter
. For example, it is informative to compare the MOI likelihood at the end of the 200 mif2
iterations above to that after, say, iteration 20.
Doing so will reveal that either (1) the MOI likelihood is stable or growing or (2) the MOI likelihood, like the MM likelihood, is declining. These two cases have different interpretations.
mif2
is doing its job, which is to maximize the likelihood of the MOI.mif2
cooling is happening too quickly.See also Issue 27.
partrans
, parameter_trans
) in pomp and when should I use them?It is frequently the case that the parameters of a model make sense only for some restricted set of values. For example, rates are typically positive numbers, and probabilities must be between 0 and 1. Put another way, the parameter space that is to be searched is constrained. Many powerful estimation algorithms, however, are unconstrained optimization algorithms, i.e., they assume that the parameter space to be searched is Euclidean space Rn. These algorithms will therefore potentially explore regions of parameter space for which the model constraints are violated (e.g., negative rates and probabilities >1.) One way to enforce the constraints is to transform the model parameter space so as to remove the constraints. Thus for example, the logarithm transforms the space of positive rates onto the real line and the logit function (the log-odds-ratio) transforms the unit interval (0,1) onto the real line.
To accommodate this approach, pomp provides a facility for the user to supply such parameter transformations. These transformations take the parameters on the “natural scale” to another “estimation scale”. One specifies parameter transformations using the parameter_trans
function, passing its output to the partrans
argument to any pomp estimation algorithm (or to the basic constructor, pomp
). For convenience, shortcuts for the commonly used log, logit, and log-barycentric transformations are provided. [Note that the latter use C snippets internally, so that the paramnames
argument must be specified when the log
, logit
, or barycentric
arguments of parameter_trans
are used.] See the parameter_trans
manual page for more detail.
mif2
)rw.sd
)?First, some background on how the iterated filtering (IF2) algorithm works. As its name suggests, IF2 (implemented as mif2
in pomp) essentially iterates a particle filter. It does this Nmif
times. Crucially, however, it does this not using the exact model you’ve specified in your pomp object, which has parameters that are constant in time. Instead, it filters a model in which the time-constant parameters have been replaced by parameters that take a random walk with intensities given by the rw.sd
argument. In other words, at each time-step, a random value (normally distributed in the current implementation) with mean zero and s.d. set by rw.sd
is added to the parameter’s value. As the iterations proceed, the intensities are gradually drawn down (according to the “cooling schedule” specified by cooling.type
and cooling.factor.50
). The theory tells us that the swarm of particles will then concentrate around the maximum likelihood estimate (MLE) eventually, i.e., with sufficiently many iterations if the cooling is not too quick. All this is explained in the IF2 paper (Ionides, Nguyen, Atchadé, Stoev, and King 2015).
Some insight into the algorithm’s effectiveness is afforded by considering that the effect of the added perturbations is to smooth the likelihood surface. This smoothing out of the local structure is a double-edged sword. On the one hand, it makes the large-scale structure of the surface plainer, which is useful when one is far from the MLE. On the other hand, it makes it impossible to find the exact MLE when it is close by.
The key thing to understand is that there are two parameter-space scales at work. First, there is the scale dictated by the distance between your starting guess and the MLE. This can be hard to know a priori, obviously. (I.e., if you knew where the MLE was, you wouldn’t be searching for it.) Second, there is the scale over which the log likelihood function changes appreciably. (Recall that a unit change in the log likelihood is considered “appreciable”.) Again, this can be hard to know a priori and it can be very different in different regions of parameter space. In particular, it can be quite different in the lowlands far from the MLE than it is in the high country near the maximum of the likelihood. Moreover, there is a tension between these scales: If you choose the random-walk intensities too small, it will take more IF2 iterations to reach the MLE; If you choose it too large, the random-walk noise will obscure the fine-structure of the likelihood surface.
So, back to the question: What to do? The following gives some heuristics. When model parameters are rates (or probabilities), it is plausible to imagine that multiplicative perturbations on the order of a few percent to the rates (or odds ratios) will lead to relatively small effects on the model behavior. Of course, this is famously not the case in general! Nevertheless, it suggests perturbing the rates (or odds ratios) additively on the log scale with random increments on the order of a few percent. The idea here is to err on the side of small perturbations, counting on cheap computing power to perform the IF2 iterations needed to achieve the MLE. At any rate, this is the reasoning behind the common choice of setting rw.sd
to a few percent (say 0.01–0.05). [When parameter transformations are in use, the perturbations will be applied on the estimation scale.] It frequently happens, then, that IF2 rapidly climbs the likelihood surface, despite the choice of a “small” rw.sd
. In subsequent iterations, it is often useful to reduce the random-walk intensities further, as attention shifts to the fine structure in the highlands of the likelihood surface.
Special considerations apply in the case of so-called initial value parameters (IVPs). Commonly, a model has some parameters that affect only the initial values (i.e., the values at t=t0) of the latent Markovian state process. Clearly, it is useless to perturb these parameters when t>t0. Indeed, it’s worse than useless since the perturbations will impede the progress of such parameters toward their maximum-likelihood values. Declaring such parameters to be initial value parameters using ivp
results in perturbations that are only applied at t=t0. More generally, there may be parameters for which perturbations should only be applied at certain times. The rw.sd
function allows one to specify (via a simple R expression in terms of time t
) which parameters fall into this category and when and how much they are to be perturbed. See the mif2
manual page for details.
To see what is happening during the execution of a C snippet, one can embed calls to Rprintf
. Such calls cause information to be printed to the R console. In the following example, we have a dmeasure snippet that implements a measurement model wherein D1 is negative-binomially distributed with mean ρ1N1 and variance ρ1N1(1+kρ1/N1) and D2 is binomially distributed with size N2 and probability ρ2. Whenever a non-finite log likelihood results (i.e., -Inf
, Inf
, NA
, or NaN
), the values of some relevant quantities are printed.
Csnippet("
dmeas <- lik = dnbinom_mu(D1,1/k,rho1*N1,1) + dbinom(D2,N2,rho2,1);
if (!R_FINITE(lik)) {
Rprintf(\"t = %lg, D1 = %lg, loglik = %lg\\n\",t,D1,lik);
}
lik = (log) ? lik : exp(lik);
")
More generally, one can make use of any aspect of R’s C API in writing C snippets.
Note that a call to Rprintf
involves passing a quoted string. To avoid having to escape those quotes, and the linefeed character "\n"
, as in the above code, one can use R’s raw string literals, as in the following.
Csnippet(r"{
dmeas <- lik = dnbinom_mu(D1,1/k,rho1*N1,1) + dbinom(D2,N2,rho2,1);
if (!R_FINITE(lik)) {
Rprintf("t = %lg, D1 = %lg, loglik = %lg\n",t,D1,lik);
}
lik = (log) ? lik : exp(lik);
}")
Ionides EL, Nguyen D, Atchadé Y, Stoev S, King AA (2015). “Inference for Dynamic and Latent Variable Models via Iterated, Perturbed Bayes Maps.” Proc Natl Acad Sci, 112(3), 719–724. https://doi.org/10.1073/pnas.1410597112.
King AA, Nguyen D, Ionides EL (2016). “Statistical Inference for Partially Observed Markov Processes via the R Package Pomp.” J Stat Softw, 69(12), 1–43. https://doi.org/10.18637/jss.v069.i12.