[Pomp-commits] r250 - in pkg: R data inst inst/data-R inst/doc inst/examples inst/include man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 20 22:58:12 CEST 2010
Author: kingaa
Date: 2010-05-20 22:58:12 +0200 (Thu, 20 May 2010)
New Revision: 250
Modified:
pkg/R/aaa.R
pkg/R/pomp.R
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/data-R/euler.sir.R
pkg/inst/data-R/gillespie.sir.R
pkg/inst/data-R/ou2.R
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/ou2-first-mif.rda
pkg/inst/doc/ou2-trajmatch.rda
pkg/inst/examples/sir.R
pkg/inst/examples/sir.c
pkg/inst/include/pomp.h
pkg/man/pomp-class.Rd
pkg/man/pomp.Rd
pkg/src/dmeasure.c
pkg/src/ou2.c
pkg/src/pomp.h
pkg/src/rmeasure.c
pkg/src/sir.c
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- add 'obsnames' slot to the 'pomp' class
- change 'dmeasure' and 'rmeasure' functions to make use of 'obsnames' in the case of user-supplied native routines
- modify prototypes of user-supplied native routines for measurement models in keeping with this
- modify source code for SIR and OU2 models
Modified: pkg/R/aaa.R
===================================================================
--- pkg/R/aaa.R 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/R/aaa.R 2010-05-20 20:58:12 UTC (rev 250)
@@ -30,6 +30,7 @@
params = 'numeric',
covar = 'matrix',
tcovar = 'numeric',
+ obsnames = 'character',
statenames = 'character',
paramnames = 'character',
covarnames = 'character',
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/R/pomp.R 2010-05-20 20:58:12 UTC (rev 250)
@@ -11,7 +11,7 @@
pomp <- function (data, times, t0, ..., rprocess, dprocess,
rmeasure, dmeasure, measurement.model,
skeleton.map, skeleton.vectorfield, initializer, covar, tcovar,
- statenames, paramnames, covarnames,
+ obsnames, statenames, paramnames, covarnames,
PACKAGE) {
## save the call
this.call <- match.call()
@@ -125,6 +125,7 @@
call.=TRUE
)
+ if (missing(obsnames)) obsnames <- character(0)
if (missing(statenames)) statenames <- character(0)
if (missing(paramnames)) paramnames <- character(0)
if (missing(covarnames)) covarnames <- character(0)
@@ -200,6 +201,7 @@
initializer = initializer,
covar = covar,
tcovar = tcovar,
+ obsnames = obsnames,
statenames = statenames,
paramnames = paramnames,
covarnames = covarnames,
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/ChangeLog 2010-05-20 20:58:12 UTC (rev 250)
@@ -1,5 +1,10 @@
+2010-05-20 kingaa
+
+ * [r249] src/pomp_internal.h: - initialize arrays with NA
+
2010-05-19 kingaa
+ * [r248] inst/ChangeLog:
* [r247] DESCRIPTION, inst/doc/advanced_topics_in_pomp.Rnw,
inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.Rnw,
inst/doc/intro_to_pomp.pdf, inst/examples/euler_sir.R,
Modified: pkg/inst/data-R/euler.sir.R
===================================================================
--- pkg/inst/data-R/euler.sir.R 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/data-R/euler.sir.R 2010-05-20 20:58:12 UTC (rev 250)
@@ -2,9 +2,22 @@
simulate(
pomp(
- times=seq(1/52,4,by=1/52),
- data=rbind(reports=numeric(52*4)),
+ data=data.frame(
+ time=seq(from=1/52,to=4,by=1/52),
+ reports=NA
+ ),
+ times="time",
t0=0,
+ rprocess=euler.sim(
+ step.fun="sir_euler_simulator",
+ delta.t=1/52/20,
+ PACKAGE="pomp"
+ ),
+ skeleton.vectorfield="sir_ODE",
+ rmeasure="sir_binom_rmeasure",
+ dmeasure="sir_binom_dmeasure",
+ PACKAGE="pomp",
+ obsnames = c("reports"),
statenames=c("S","I","R","cases","W"),
paramnames=c(
"gamma","mu","iota",
@@ -13,15 +26,6 @@
),
zeronames=c("cases"),
comp.names=c("S","I","R"),
- rprocess=euler.sim(
- step.fun="sir_euler_simulator",
- delta.t=1/52/20,
- PACKAGE="pomp"
- ),
- skeleton.vectorfield="sir_ODE",
- rmeasure="binom_rmeasure",
- dmeasure="binom_dmeasure",
- PACKAGE="pomp",
initializer=function(params, t0, comp.names, ...){
p <- exp(params)
snames <- c("S","I","R","cases","W")
Modified: pkg/inst/data-R/gillespie.sir.R
===================================================================
--- pkg/inst/data-R/gillespie.sir.R 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/data-R/gillespie.sir.R 2010-05-20 20:58:12 UTC (rev 250)
@@ -8,6 +8,7 @@
beta3=log(490),
gamma=log(24),
iota=log(0.1),
+ rho=log(0.1),
S.0=log(0.05),
I.0=log(1e-4),
R.0=log(0.95),
@@ -46,8 +47,13 @@
rdeath=c(0,0,1,0,0)
)
),
+ obsnames="reports",
statenames=c("S","I","R","N","cases"),
- paramnames=c("gamma","mu","iota","beta1","nu","nbasis","degree","period"),
+ paramnames=c(
+ "gamma","mu","iota",
+ "beta1","nu",
+ "nbasis","degree","period"
+ ),
zeronames=c("cases"),
measurement.model=reports~binom(size=cases,prob=0.1),
initializer=function(params, t0, ...){
Modified: pkg/inst/data-R/ou2.R
===================================================================
--- pkg/inst/data-R/ou2.R 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/data-R/ou2.R 2010-05-20 20:58:12 UTC (rev 250)
@@ -54,8 +54,8 @@
dim=c(nrep,ntimes-1)
)
},
- dmeasure = "normal_dmeasure",
- rmeasure = "normal_rmeasure",
+ dmeasure = "ou2_normal_dmeasure",
+ rmeasure = "ou2_normal_rmeasure",
skeleton.map = function (x, t, params, ...) {
with(
as.list(c(x,params)),
@@ -72,7 +72,8 @@
"sigma.1","sigma.2","sigma.3",
"tau"
),
- statenames = c("x1","x2")
+ statenames = c("x1","x2"),
+ obsnames = c("y1","y2")
),
params=c(
alpha.1=0.8, alpha.2=-0.5, alpha.3=0.3, alpha.4=0.9,
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-05-20 20:58:12 UTC (rev 250)
@@ -64,7 +64,7 @@
This document serves to give some examples of the use of native (C or FORTRAN) codes in \code{pomp}
and to introduce the low-level interface to \code{pomp} objects.
-\section{Acceleration using native codes.}
+\section{Acceleration using native codes: writing \code{rprocess} and \code{dprocess} from scratch.}
Since many of the methods we will use require us to simulate the process and/or measurement models many times, it is a good idea to use native (compiled) codes for the computational heavy lifting.
This can result in many-fold speedup.
@@ -142,26 +142,37 @@
t0=0,
rprocess = ou2.rprocess,
dprocess = ou2.dprocess,
- dmeasure = "normal_dmeasure",
- rmeasure = "normal_rmeasure",
+ dmeasure = "ou2_normal_dmeasure",
+ rmeasure = "ou2_normal_rmeasure",
paramnames=c(
"alpha.1","alpha.2","alpha.3","alpha.4",
"sigma.1","sigma.2","sigma.3",
"tau"
),
statenames = c("x1","x2"),
+ obsnames = c("y1","y2"),
PACKAGE="pomp"
)
@
Notice that the process model is implemented using using \verb+.C+, while the measurement model is specified by giving the names of native C routines.
Read the source to see the definitions of these functions.
-For convenience, the source codes are provided with the package in the \code{examples} directory.\
+For convenience, the source codes are provided with the package in the \code{examples} directory.
Do
<<view-ou2-source,eval=F>>=
edit(file=system.file("examples/ou2.c",package="pomp"))
@
to view the source code.
+There is an important issue that arises when using native codes.
+This has to do with the order in which parameters, states, and observables are passed into the native codes.
+\pomp\ relies on the names (also row-names and column-names) attributes to identify variables in vectors and arrays.
+When you write a C or FORTRAN version of \code{rprocess} or \code{dmeasure} for example, you write a routine that takes parameters, state variables, and/or observables in the form of a vector.
+However, you have no control over the order in which these are given to you.
+Without some means of knowing which element of each vector corresponds to which variable, you cannot write the codes correctly.
+This is where the \code{paramnames}, \code{statenames}, and \code{obsnames} arguments to \code{pomp} come in.
+When you specifying the names of parameters, state variables, and observables (data variables) here, \code{pomp} matches these names against the corresponding names attributes of vectors and passes to your native routine integer vectors which you can use to identify the correct variables.
+See the source code to see how this is done.
+
We'll specify some parameters:
<<>>=
theta <- c(
@@ -179,12 +190,12 @@
print(toc-tic)
@
-In this example, we've written our simulators and density functions ``from scratch''.
+\section{Acceleration using native codes: using plug-ins with native code}
+
+In the preceding example, we've written our simulators and density functions ``from scratch''.
\code{pomp} provides ``plug-in'' facilities to make it easier to define certain kinds of models.
These plug-ins can be used with native codes as well, as we'll see in the next example.
-\section{The seasonal SIR model implemented using compiled native routines}
-
In the ``intro\_to\_pomp'' vignette, we looked at the SIR model, which we implemented using an Euler-multinomial approximation to the continuous-time Markov process.
Here is the same model implemented using native C codes:
<<sir-def,keep.source=T>>=
@@ -203,10 +214,12 @@
## native routine for the skeleton:
skeleton.vectorfield="sir_ODE",
## binomial measurement model:
- rmeasure="binom_rmeasure",
- ## binomial measurement model:
- dmeasure="binom_dmeasure",
- PACKAGE="pomp", ## name of the shared-object library
+ rmeasure="sir_binom_rmeasure",
+ dmeasure="sir_binom_dmeasure",
+ ## name of the shared-object library containing the native routines:
+ PACKAGE="pomp",
+ ## the order of the observable assumed in the native routines:
+ obsnames=c("reports"),
## the order of the state variables assumed in the native routines:
statenames=c("S","I","R","cases","W"),
## the order of the parameters assumed in the native routines:
@@ -230,12 +243,16 @@
}
) -> sir
@
-The source code for the native routines \verb+sir_euler_simulator+, \verb+sir_ODE+, \verb+binom_rmeasure+, and \verb+binom_dmeasure+ is provided with the package (in the \code{examples} directory).
+The source code for the native routines \verb+sir_euler_simulator+, \verb+sir_ODE+, \verb+sir_binom_rmeasure+, and \verb+sir_binom_dmeasure+ is provided with the package (in the \code{examples} directory).
To see the source code, do
<<view-sir-source,eval=F>>=
edit(file=system.file("examples/sir.c",package="pomp"))
@
Also in the \code{examples} directory is an \R\ script that shows how to compile \verb+sir.c+ into a shared-object library and link it with \R.
+Note that the native routines for this model are included in the package, which is why we give the \verb+PACKAGE="pomp"+ argument to \pomp.
+When you write your own model using native routines, you'll compile them into a dynamically-loadable library.
+In this case, you'll want to specify the name of that library using the \code{PACKAGE} argument.
+Again, refer to the SIR example included in the \code{examples} directory to see how this is done.
Let's specify some parameters and simulate:
<<sir-sim>>=
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ou2-first-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ou2-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/examples/sir.R 2010-05-20 20:58:12 UTC (rev 250)
@@ -43,14 +43,16 @@
delta.t=1/52/20
),
skeleton.vectorfield="sir_ODE", # native routine for the skeleton
- rmeasure="binom_rmeasure", # binomial measurement model
- dmeasure="binom_dmeasure", # binomial measurement model
+ rmeasure="sir_binom_rmeasure", # binomial measurement model
+ dmeasure="sir_binom_dmeasure", # binomial measurement model
PACKAGE="sir", ## name of the shared-object library
+ ## the order of the observables assumed in the native routines:
+ obsnames="reports",
## the order of the state variables assumed in the native routines:
statenames=c("S","I","R","cases","W"),
## the order of the parameters assumed in the native routines:
paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho","nbasis","degree","period"),
- ## reset cases to zero at each new observation:
+ ## reset cases to zero after each new observation:
zeronames=c("cases"),
initializer=function(params,t0,...){
p <- exp(params)
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/examples/sir.c 2010-05-20 20:58:12 UTC (rev 250)
@@ -41,17 +41,17 @@
#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
#define W (x[stateindex[4]]) // integrated white noise
-#define REPORTS (y[0])
+#define REPORTS (y[obsindex[0]])
-void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
+void sir_binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
*lik = dbinom(REPORTS,CASE,exp(LOGRHO),give_log);
}
-void binom_rmeasure (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
+void sir_binom_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
REPORTS = rbinom(CASE,exp(LOGRHO));
}
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/inst/include/pomp.h 2010-05-20 20:58:12 UTC (rev 250)
@@ -115,12 +115,14 @@
// Prototype for measurement model simulation
typedef void pomp_measure_model_simulator (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t);
// Description:
// on input:
// x = pointer to state vector at time t
// p = pointer to parameter vector
+// obsindex = pointer to vector of integers indexing the variables in 'y' in the order specified by
+// the 'obsnames' slot
// stateindex = pointer to vector of integers indexing the states in 'x' in the order specified by
// the 'statenames' slot
// parindex = pointer to vector of integers indexing the parameters in 'p' in the order specified by
@@ -141,7 +143,7 @@
// Prototype for measurement model density evaluator
typedef void pomp_measure_model_density (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t);
// Description:
// on input:
@@ -149,6 +151,8 @@
// x = pointer to state vector at time t
// p = pointer to parameter vector
// give_log = should the log likelihood be returned?
+// obsindex = pointer to vector of integers indexing the variables in 'y' in the order specified by
+// the 'obsnames' slot
// stateindex = pointer to vector of integers indexing the states in 'x' in the order specified by
// the 'statenames' slot
// parindex = pointer to vector of integers indexing the parameters in 'p' in the order specified by
Modified: pkg/man/pomp-class.Rd
===================================================================
--- pkg/man/pomp-class.Rd 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/man/pomp-class.Rd 2010-05-20 20:58:12 UTC (rev 250)
@@ -52,8 +52,8 @@
\item{covar, tcovar}{
An optional table of covariates.
}
- \item{statenames, paramnames, covarnames}{
- Names state variables, parameters, and covariates, respectively.
+ \item{obsnames, statenames, paramnames, covarnames}{
+ Names of observables, state variables, parameters, and covariates, respectively.
These are used to make the interface with native routines more robust.
}
\item{PACKAGE}{
Modified: pkg/man/pomp.Rd
===================================================================
--- pkg/man/pomp.Rd 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/man/pomp.Rd 2010-05-20 20:58:12 UTC (rev 250)
@@ -7,7 +7,8 @@
\usage{
pomp(data, times, t0, \dots, rprocess, dprocess, rmeasure, dmeasure,
measurement.model, skeleton.map, skeleton.vectorfield,
- initializer, covar, tcovar, statenames, paramnames, covarnames,
+ initializer, covar, tcovar,
+ obsnames, statenames, paramnames, covarnames,
PACKAGE)
}
\arguments{
@@ -86,11 +87,12 @@
If a covariate table is supplied, then the value of each of the covariates is interpolated as needed, i.e., whenever \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, or \code{init.state} is evaluated.
The resulting interpolated values are passed to the corresponding functions as a numeric vector named \code{covars}.
}
- \item{statenames, paramnames, covarnames}{
- Optional character vectors specifying the names of state variables, parameters, or covariates, respectively.
+ \item{obsnames, statenames, paramnames, covarnames}{
+ Optional character vectors specifying the names of observables, state variables, parameters, or covariates, respectively.
These are only used in the event that one or more of the basic functions (\code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, \code{skeleton}) are defined using native routines.
In that case, these name vectors are matched against the corresponding names and the indices of the names are passed to the native routines.
Using this facility allows one to write one or more of \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, \code{skeleton} in native code in a way that does not depend on the order of states, parameters, and covariates at run time.
+ See the \dQuote{advanced_topic_in_pomp} vignette for more on this topic and examples.
}
\item{PACKAGE}{
An optional string giving the name of the dynamically loaded library in which any native routines are to be found.
Modified: pkg/src/dmeasure.c
===================================================================
--- pkg/src/dmeasure.c 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/src/dmeasure.c 2010-05-20 20:58:12 UTC (rev 250)
@@ -11,7 +11,7 @@
double *lik, double *y,
double *x, double *times, double *params,
int *give_log, int *ndim,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
double *time_table, double *covar_table)
{
double t, *lp, *xp, *pp, *yp;
@@ -45,7 +45,7 @@
xp = &x[nvar*(p+nrep*k)];
pp = ¶ms[npar*p];
- (*f)(lp,yp,xp,pp,*give_log,stateindex,parindex,covindex,covdim,covar_fn,t);
+ (*f)(lp,yp,xp,pp,*give_log,obsindex,stateindex,parindex,covindex,covdim,covar_fn,t);
}
}
@@ -77,9 +77,9 @@
// this is the measurement pdf evaluated when the user supplies an R function
// (and not a native routine)
-// Note that stateindex, parindex, covindex are ignored.
+// Note that obsindex, stateindex, parindex, covindex are ignored.
static void default_meas_dens (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
int covdim, double *covar, double t)
{
int nprotect = 0;
@@ -112,12 +112,12 @@
SEXP F, fn;
SEXP dimP, dimX, dimY;
SEXP tcovar, covar;
- SEXP statenames, paramnames, covarnames;
- SEXP sindex, pindex, cindex;
- int *sidx, *pidx, *cidx;
+ SEXP statenames, paramnames, covarnames, obsnames;
+ SEXP sindex, pindex, cindex, oindex;
+ int *sidx, *pidx, *cidx, *oidx;
SEXP Xnames, Ynames, Pnames, Cnames;
int use_native;
- int nstates, nparams, ncovars;
+ int nstates, nparams, ncovars, nobsers;
pomp_measure_model_density *ff = NULL;
PROTECT(times = AS_NUMERIC(times)); nprotect++;
@@ -149,9 +149,11 @@
PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
+ PROTECT(obsnames = GET_SLOT(object,install("obsnames"))); nprotect++;
PROTECT(statenames = GET_SLOT(object,install("statenames"))); nprotect++;
PROTECT(paramnames = GET_SLOT(object,install("paramnames"))); nprotect++;
PROTECT(covarnames = GET_SLOT(object,install("covarnames"))); nprotect++;
+ nobsers = LENGTH(obsnames);
nstates = LENGTH(statenames);
nparams = LENGTH(paramnames);
ncovars = LENGTH(covarnames);
@@ -206,6 +208,13 @@
ndim[0] = nreps; ndim[1] = ntimes;
PROTECT(F = makearray(2,ndim)); nprotect++;
+ if (nobsers > 0) {
+ PROTECT(oindex = matchnames(Ynames,obsnames)); nprotect++;
+ oidx = INTEGER(oindex);
+ } else {
+ oidx = 0;
+ }
+
if (nstates > 0) {
PROTECT(sindex = matchnames(Xnames,statenames)); nprotect++;
sidx = INTEGER(sindex);
@@ -231,7 +240,7 @@
ndim[4] = covlen; ndim[5] = covdim; ndim[6] = nobs;
dens_meas(ff,REAL(F),REAL(y),REAL(x),REAL(times),REAL(params),INTEGER(log),ndim,
- sidx,pidx,cidx,
+ oidx,sidx,pidx,cidx,
REAL(tcovar),REAL(covar));
UNPROTECT(nprotect);
Modified: pkg/src/ou2.c
===================================================================
--- pkg/src/ou2.c 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/src/ou2.c 2010-05-20 20:58:12 UTC (rev 250)
@@ -8,13 +8,12 @@
// prototypes
-void normal_rmeasure (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int ncovar, double *covar,
- double t);
-void normal_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
- int covdim, double *covar, double t);
+void ou2_normal_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovar, double *covar, double t);
+void ou2_normal_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int covdim, double *covar, double t);
void ou2_adv (double *x, double *xstart, double *par, double *times, int *n, int *parindex);
void ou2_pdf (double *d, double *X, double *par, double *times, int *n, int *parindex, int *give_log);
static void sim_ou2 (double *x,
@@ -85,13 +84,13 @@
#define X1 (x[stateindex[0]])
#define X2 (x[stateindex[1]])
#define TAU (p[parindex[7]])
-#define Y1 (y[0])
-#define Y2 (y[1])
+#define Y1 (y[obsindex[0]])
+#define Y2 (y[obsindex[1]])
// bivariate normal measurement error density
-void normal_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
- int covdim, double *covar, double t)
+void ou2_normal_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int covdim, double *covar, double t)
{
double sd = fabs(TAU);
double f = 0.0;
@@ -101,10 +100,10 @@
}
// bivariate normal measurement error simulator
-void normal_rmeasure (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int ncovar, double *covar,
- double t)
+void ou2_normal_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovar, double *covar,
+ double t)
{
double sd = fabs(TAU);
Y1 = rnorm(X1,sd);
Modified: pkg/src/pomp.h
===================================================================
--- pkg/src/pomp.h 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/src/pomp.h 2010-05-20 20:58:12 UTC (rev 250)
@@ -115,12 +115,14 @@
// Prototype for measurement model simulation
typedef void pomp_measure_model_simulator (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t);
// Description:
// on input:
// x = pointer to state vector at time t
// p = pointer to parameter vector
+// obsindex = pointer to vector of integers indexing the variables in 'y' in the order specified by
+// the 'obsnames' slot
// stateindex = pointer to vector of integers indexing the states in 'x' in the order specified by
// the 'statenames' slot
// parindex = pointer to vector of integers indexing the parameters in 'p' in the order specified by
@@ -141,7 +143,7 @@
// Prototype for measurement model density evaluator
typedef void pomp_measure_model_density (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t);
// Description:
// on input:
@@ -149,6 +151,8 @@
// x = pointer to state vector at time t
// p = pointer to parameter vector
// give_log = should the log likelihood be returned?
+// obsindex = pointer to vector of integers indexing the variables in 'y' in the order specified by
+// the 'obsnames' slot
// stateindex = pointer to vector of integers indexing the states in 'x' in the order specified by
// the 'statenames' slot
// parindex = pointer to vector of integers indexing the parameters in 'p' in the order specified by
Modified: pkg/src/rmeasure.c
===================================================================
--- pkg/src/rmeasure.c 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/src/rmeasure.c 2010-05-20 20:58:12 UTC (rev 250)
@@ -11,7 +11,7 @@
double *y,
double *x, double *times, double *params,
int *ndim,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
double *time_table, double *covar_table)
{
double t, *xp, *pp, *yp;
@@ -44,7 +44,7 @@
xp = &x[nvar*(p+nrep*k)];
pp = ¶ms[npar*p];
- (*f)(yp,xp,pp,stateindex,parindex,covindex,covdim,covar_fn,t);
+ (*f)(yp,xp,pp,obsindex,stateindex,parindex,covindex,covdim,covar_fn,t);
}
}
@@ -80,9 +80,9 @@
// this is the measurement simulator evaluated when the user supplies an R function
// (and not a native routine)
-// Note that stateindex, parindex, covindex are ignored.
+// Note that obsindex, stateindex, parindex, covindex are ignored.
static void default_meas_sim (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
int covdim, double *covar, double t)
{
int nprotect = 0;
@@ -127,12 +127,12 @@
SEXP Y, fn;
SEXP dimP, dimX, dimD;
SEXP tcovar, covar;
- SEXP statenames, paramnames, covarnames;
- SEXP sindex, pindex, cindex;
- int *sidx, *pidx, *cidx;
+ SEXP statenames, paramnames, covarnames, obsnames;
+ SEXP sindex, pindex, cindex, oindex;
+ int *sidx, *pidx, *cidx, *oidx;
SEXP Xnames, Pnames, Cnames;
int use_native;
- int nstates, nparams, ncovars;
+ int nstates, nparams, ncovars, nobsers;
pomp_measure_model_simulator *ff = NULL;
PROTECT(times = AS_NUMERIC(times)); nprotect++;
@@ -160,9 +160,11 @@
PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
- PROTECT(statenames = GET_SLOT(object,install("statenames"))); nprotect++;
- PROTECT(paramnames = GET_SLOT(object,install("paramnames"))); nprotect++;
- PROTECT(covarnames = GET_SLOT(object,install("covarnames"))); nprotect++;
+ PROTECT(obsnames = GET_SLOT(object,install("obsnames"))); nprotect++;
+ PROTECT(statenames = GET_SLOT(object,install("statenames"))); nprotect++;
+ PROTECT(paramnames = GET_SLOT(object,install("paramnames"))); nprotect++;
+ PROTECT(covarnames = GET_SLOT(object,install("covarnames"))); nprotect++;
+ nobsers = LENGTH(obsnames);
nstates = LENGTH(statenames);
nparams = LENGTH(paramnames);
ncovars = LENGTH(covarnames);
@@ -215,6 +217,13 @@
PROTECT(Y = makearray(3,ndim)); nprotect++;
setrownames(Y,OBNM,3);
+ if (nobsers > 0) {
+ PROTECT(oindex = matchnames(OBNM,obsnames)); nprotect++;
+ oidx = INTEGER(oindex);
+ } else {
+ oidx = 0;
+ }
+
if (nstates > 0) {
PROTECT(sindex = matchnames(Xnames,statenames)); nprotect++;
sidx = INTEGER(sindex);
@@ -242,13 +251,14 @@
if (use_native) GetRNGstate();
simul_meas(ff,REAL(Y),REAL(x),REAL(times),REAL(params),ndim,
- sidx,pidx,cidx,
+ oidx,sidx,pidx,cidx,
REAL(tcovar),REAL(covar));
if (use_native) PutRNGstate();
if (OIDX != 0) Free(OIDX);
OIDX = 0;
+ FIRST = 0;
UNPROTECT(nprotect);
return Y;
Modified: pkg/src/sir.c
===================================================================
--- pkg/src/sir.c 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/src/sir.c 2010-05-20 20:58:12 UTC (rev 250)
@@ -41,17 +41,17 @@
#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
#define W (x[stateindex[4]]) // integrated white noise
-#define REPORTS (y[0])
+#define REPORTS (y[obsindex[0]])
-void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
+void sir_binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
*lik = dbinom(REPORTS,CASE,exp(LOGRHO),give_log);
}
-void binom_rmeasure (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
+void sir_binom_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
REPORTS = rbinom(CASE,exp(LOGRHO));
}
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2010-05-20 11:48:58 UTC (rev 249)
+++ pkg/tests/sir.R 2010-05-20 20:58:12 UTC (rev 250)
@@ -18,7 +18,7 @@
## the C codes "sir_euler_simulator" and "sir_euler_density" are included in the "examples" directory (file "sir.c")
po <- pomp(
times=1/52*seq.int(length=4*52),
- data=rbind(measles=numeric(52*4)),
+ data=rbind(reports=numeric(52*4)),
t0=0,
tcovar=tbasis,
covar=basis,
@@ -107,7 +107,7 @@
}
)
},
- measurement.model=measles~binom(size=cases,prob=exp(rho)),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 250
More information about the pomp-commits
mailing list